dc02_discrete-choice-estimation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# <center>Discrete choice estimation</center>\n",
    "### <center>Alfred Galichon (NYU & Sciences Po)</center>\n",
    "## <center>'math+econ+code' masterclass series</center>\n",
    "#### <center>With python code examples</center>\n",
    "© 2018–2023 by Alfred Galichon. Past and present support from NSF grant DMS-1716489, ERC grant CoG-866274 are acknowledged, as well as inputs from contributors listed [here](http://www.math-econ-code.org/team).\n",
    "\n",
    "**If you reuse material from this masterclass, please cite as:**<br>\n",
    "Alfred Galichon, 'math+econ+code' masterclass series. https://www.math-econ-code.org/"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## References\n",
    "\n",
    "* Savage, L. (1951). The theory of statistical decision. JASA.\n",
    "* Bonnet, Fougère, Galichon, Poulhès (2021). Minimax estimation of hedonic models. Preprint."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Loading the libraries\n",
    "\n",
    "First, let's load the libraries we shall need."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# !python -m pip install -i https://pypi.gurobi.com gurobipy ## only if Gurobi not here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import scipy.sparse as spr\n",
    "from scipy import optimize, special\n",
    "import gurobipy as grb\n",
    "from sklearn import linear_model\n",
    "from tabulate import tabulate\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will also import objects created in the previous lectures, which are now stored in the `objects_D3`module."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# from objects_D3 import *"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Our data\n",
    "We will go back to the dataset of Greene and Hensher (1997). As a reminder, 210 individuals are surveyed about their choice of travel mode between Sydney, Canberra and Melbourne, and the various costs (time and money) associated with each alternative. More precisely, recall that the variable `wait` is the waiting time of the plane, train or bus (waiting time is zero for car). `vcost` is the in-vehicle cost, without accounting for the value of the time savings. `travel` is the in-vehicle travel time. `gcost` is the generalized cost, which is the sum of ` vcost` and the value associated with time savings. `income` is the income of the household, in thousands of dollars. `size` is the size of the travel group.\n",
    "\n",
    "Therefore there are 840 = 4 x 210 observations, which we can stack into `travelmodedataset` a 3 dimensional array whose dimensions are mode,individual,dummy for choice+covariates.\n",
    "\n",
    "Let's load the dataset and take a first glance at it:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>individual</th>\n",
       "      <th>mode</th>\n",
       "      <th>choice</th>\n",
       "      <th>wait</th>\n",
       "      <th>vcost</th>\n",
       "      <th>travel</th>\n",
       "      <th>gcost</th>\n",
       "      <th>income</th>\n",
       "      <th>size</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1</td>\n",
       "      <td>air</td>\n",
       "      <td>no</td>\n",
       "      <td>69</td>\n",
       "      <td>59</td>\n",
       "      <td>100</td>\n",
       "      <td>70</td>\n",
       "      <td>35</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1</td>\n",
       "      <td>train</td>\n",
       "      <td>no</td>\n",
       "      <td>34</td>\n",
       "      <td>31</td>\n",
       "      <td>372</td>\n",
       "      <td>71</td>\n",
       "      <td>35</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>1</td>\n",
       "      <td>bus</td>\n",
       "      <td>no</td>\n",
       "      <td>35</td>\n",
       "      <td>25</td>\n",
       "      <td>417</td>\n",
       "      <td>70</td>\n",
       "      <td>35</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>1</td>\n",
       "      <td>car</td>\n",
       "      <td>yes</td>\n",
       "      <td>0</td>\n",
       "      <td>10</td>\n",
       "      <td>180</td>\n",
       "      <td>30</td>\n",
       "      <td>35</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>2</td>\n",
       "      <td>air</td>\n",
       "      <td>no</td>\n",
       "      <td>64</td>\n",
       "      <td>58</td>\n",
       "      <td>68</td>\n",
       "      <td>68</td>\n",
       "      <td>30</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   individual   mode choice  wait  vcost  travel  gcost  income  size\n",
       "0           1    air     no    69     59     100     70      35     1\n",
       "1           1  train     no    34     31     372     71      35     1\n",
       "2           1    bus     no    35     25     417     70      35     1\n",
       "3           1    car    yes     0     10     180     30      35     1\n",
       "4           2    air     no    64     58      68     68      30     2"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "thepath = 'https://raw.githubusercontent.com/math-econ-code/mec_optim_2021-01/master/data_mec_optim/demand_travelmode/'\n",
    "travelmode =  pd.read_csv(thepath+'travelmodedata.csv')\n",
    "travelmode.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# GLM to estimate discrete choice models\n",
    "\n",
    "## Estimation with observed heterogeneity\n",
    "\n",
    "We assume that we observe individual characteristics that are relevant for individual choices, that is $U_{iy}=\\sum_k \\Phi_{iyk} \\beta_k$, or in matrix form<br>\n",
    "$U = \\Phi \\beta,$<br>\n",
    "where $\\beta\\in\\mathbb{R}^{p}$ is a parameter, and $\\Phi$ is a $\\left(\\left\\vert \\mathcal{I}\\left\\vert\\right\\vert\\mathcal{Y}\\right\\vert \\right) \\times p$ matrix.\n",
    "\n",
    "Assume $u_{iy}=U_{iy} + \\varepsilon _{iy}= \\sum_{k}\\Phi _{iyk} \\beta _{k}+\\varepsilon _{iy}$.\n",
    "\n",
    "Let $\\hat{\\mu}_{iy}$ be the indicator that $i$ chooses alternative $y$. \n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We create a `discreteChoicePb` class where we store that data. An arc $a$ is a pair $iy$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "class DiscreteChoicePb():\n",
    "    def __init__(self,Φ_i_y_k, μhat_i_y):\n",
    "        self.nbi,self.nby,self.nbk = Φ_i_y_k.shape\n",
    "        self.nba = self.nbi * self.nby\n",
    "        self.Φ_a_k = Φ_i_y_k.reshape((self.nba,-1))\n",
    "        self.μhat_a = μhat_i_y.flatten()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We build an object `travelEx` based on our travel data and a parametric model where the regressors are `travel`, `-travel*income` and `gcost`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def prepare_travel_data():\n",
    "    μhat_a = np.where(travelmode['choice'] =='yes' , 1, 0)\n",
    "    #nobs,ncols = travelmode.shape\n",
    "    nby = travelmode['mode'].nunique()\n",
    "    nbi = travelmode.shape[0] // nby\n",
    "    covariates = travelmode[['travel', 'income', 'gcost']].values\n",
    "    Φ_a_k = np.column_stack([ covariates[:,0] , - (covariates[:,0] * covariates[:,1] ), - covariates[:,2] ])\n",
    "    _,nbk = Φ_a_k.shape \n",
    "    Φbar_k = Φ_a_k.mean(axis = 0)\n",
    "    Φstdev_k = Φ_a_k.std(axis = 0, ddof = 1)\n",
    "    Φ_i_y_k = ((Φ_a_k - Φbar_k[None,:]) / Φstdev_k[None,:]).reshape((nbi,nby,nbk))\n",
    "    return DiscreteChoicePb(Φ_i_y_k,μhat_a)\n",
    "\n",
    "travelEx = prepare_travel_data()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Maximum likelihood estimation\n",
    "\n",
    "The probability $\\mu_{iy}$ that individual $i$ chooses alternative $y$ is given by $\\partial G_i / \\partial U_y (\\Phi \\beta)$. \n",
    "\n",
    "\n",
    "The log-likelihood function is given by\n",
    "\n",
    "$\n",
    "l\\left(  \\beta\\right)  =\\sum_{y}\\hat{\\mu}_{iy}\\log \\mu_{iy}\\left(\\Phi \\beta\\right) = \\sum_{y}\\hat{\\mu}_{iy}\\log  \n",
    "\\frac {\\partial G_i}  {\\partial U_y} (\\Phi \\beta)\n",
    "$\n",
    "\n",
    "A common estimation method of $\\beta$ is by maximum likelihood: $\\max_{\\beta}l\\left(  \\beta\\right) $; MLE is statistically efficient; the problem is that the problem is not guaranteed to be convex, so there may be computational difficulties (e.g. local optima).\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### MLE, logit case\n",
    "\n",
    "In the logit case, the log-likelihood associated with observation $i\\in\\mathcal{I}$ is\n",
    "\n",
    "$\n",
    "l_{i}\\left( \\beta \\right) =\\sum_{y\\in \\mathcal{Y}}\\hat{\\mu}_{iy}\\left( \\Phi\n",
    "\\beta \\right) _{iy}-\\log \\sum_{y\\in \\mathcal{Y}}\\exp \\left( \\Phi \\beta\n",
    "\\right) _{iy}$\n",
    "\n",
    "and the max-likelihood rewrites as\n",
    "\n",
    "$\\max_{\\beta }\\left\\{ \\sum_{i\\in \\mathcal{I},y\\in \\mathcal{Y}}\\hat{\\mu}%\n",
    "_{iy}\\left( \\Phi \\beta \\right) _{iy}-\\sum_{i\\in \\mathcal{I}}\\log \\sum_{y\\in \n",
    "\\mathcal{Y}}\\exp \\left( \\Phi \\beta \\right) _{iy}\\right\\}$\n",
    "\n",
    "so that the max-likehood boils down to\n",
    "\n",
    "\\begin{align*}\n",
    "\\max_{\\beta}\\left\\{  \\hat{\\mu}^{\\intercal} \\Phi \\beta- \\sum_i G_i\\left( \\Phi \\beta\\right)\\right\\}\n",
    "\\end{align*}\n",
    "\n",
    "whose value is the Legendre-Fenchel transform of $\\beta\\rightarrow \\sum_i G_i\\left( \\Phi \\beta\\right)$ evaluated at $\\Phi ^{^{\\intercal}}\\hat{\\mu}$.\n",
    "\n",
    "Note that the vector $\\Phi^{^{\\intercal}}\\hat{\\mu}$ is the vector of empirical moments, which is a sufficient statistics in the logit model.\n",
    "\n",
    "As a result, in the logit case, the MLE is a convex optimization problem, and it is therefore both statistically efficient and computationally efficient.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Moment estimation\n",
    "\n",
    "The previous remark will inspire an alternative procedure based on the moments statistics $\\Phi^{^{\\intercal}}\\hat{\\mu}$.\n",
    "\n",
    "The social welfare is given in general by $W\\left(  \\beta\\right) =\\sum_i G_i\\left(  \\Phi\\beta\\right)  $. One has <br>$\\partial_{\\beta^{k}}W\\left(\\beta\\right)  =\\sum_{iy} \\frac {\\partial G_i} {\\partial U_y}(\\Phi\\beta) \\Phi_{yik}$, that is \n",
    "\n",
    "\\begin{align*}\n",
    "\\nabla W\\left(  \\beta\\right)  = \\Phi^{\\intercal}\\nabla G_i\\left(  \\Phi\\beta\\right)  ,\n",
    "\\end{align*}\n",
    "\n",
    "which is the vector of predicted moments.\n",
    "\n",
    "Therefore the program\n",
    "\n",
    "$\\max_{\\beta }\\left\\{ \\hat{\\mu}^{\\top }\\Phi \\beta -\\sum_{i\\in \\mathcal{I}%\n",
    "}G\\left( \\left( \\Phi \\beta \\right) _{i.}\\right) \\right\\} ,$\n",
    "\n",
    "(where G is the Emax operator associated with the distribution of the random utility), picks up the parameter $\\beta$ which matches the empirical moments $\\Phi^{^{\\intercal}}\\hat{q}$ with the predicted ones $\\nabla W\\left(\\beta\\right)  $. This procedure is not statistically efficient, but is computationally efficient becauses it arises as a convex optimization problem."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Creating a DiscreteChoicePb class"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def DiscreteChoicePb_mle_diy(self):\n",
    "    def minus_log_likelihood(β_k, σ = 1):\n",
    "        Φβ_i_y = (self.Φ_a_k.dot(β_k)).reshape((-1,self.nby)) \n",
    "        maxΦβ_i = Φβ_i_y.max(axis = 1)\n",
    "        d_i = np.sum(np.exp((Φβ_i_y -maxΦβ_i[:,None])/σ ), axis = 1)\n",
    "        return - ((Φβ_i_y.flatten()*self.μhat_a).sum() / σ  -  (maxΦβ_i / σ + np.log(d_i)).sum())\n",
    "\n",
    "    def grad_minus_log_likelihood(β_k, σ = 1):\n",
    "        Φβ_i_y = (self.Φ_a_k.dot(β_k)).reshape((-1,self.nby)) \n",
    "        maxΦβ_i = Φβ_i_y.max(axis = 1)\n",
    "        d_i = np.sum(np.exp((Φβ_i_y - maxΦβ_i[:,None] )/σ ), axis = 1)\n",
    "        μβ_iy = (np.exp((Φβ_i_y - maxΦβ_i[:,None] )/σ ) / d_i[:,None]).flatten()\n",
    "        return  - ((self.μhat_a - μβ_iy).reshape((1,-1)) @ self.Φ_a_k).flatten()\n",
    "\n",
    "    βtilde0_k = np.zeros(self.nbk)\n",
    "    res = optimize.minimize(minus_log_likelihood,method = 'CG',jac = grad_minus_log_likelihood, x0 = βtilde0_k )\n",
    "    βtilde_k = res['x']\n",
    "    print(-minus_log_likelihood (βtilde_k))\n",
    "    return βtilde_k[:-1] / βtilde_k[-1],  1 / βtilde_k[-1], - res['fun']\n",
    "\n",
    "DiscreteChoicePb.mle_diy = DiscreteChoicePb_mle_diy\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compute using:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-277.7052141445852\n",
      "DIY approach. βmle_k =  [0.33826733 0.85179311] ; Tmle =  1.8162760072465682  ; ll= -277.7052141445852 .\n"
     ]
    }
   ],
   "source": [
    "βmle_k,Tmle,llmle = travelEx.mle_diy()\n",
    "print('DIY approach. βmle_k = ',βmle_k,'; Tmle = ',Tmle, ' ; ll=',llmle,'.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Computation as a Poisson regression using GLM in scikit-learn\n",
    "\n",
    "As a reminder, this can be computed as \n",
    "\\begin{align*}\n",
    "\\min_{\\beta,u} \\left\\{ \\sum_{iy} \\hat{\\mu}_{iy} \\left(   \\Phi\\beta - (I_\\mathcal{I} \\otimes 1_\\mathcal{Y}) u \\right)  _{iy} - \\sum_{iy} \\exp\\left(  \\left(  \\Phi\\beta - (I_\\mathcal{I} \\otimes 1_\\mathcal{Y}) u \\right) _{iy} \\right)  \\right\\}\n",
    "\\end{align*}\n",
    "\n",
    "which leads to the following call to `scikit-learn`:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def DiscreteChoicePb_mle_glm(self, max_iter = 1000, tol=0.0001):\n",
    "    poisson = linear_model.PoissonRegressor(alpha = 0, fit_intercept=False,max_iter= max_iter, tol = tol)\n",
    "    X_a_l = spr.hstack([self.Φ_a_k, -spr.kron(spr.identity(self.nbi), np.ones((self.nby,1)))])\n",
    "    poisson.fit(X_a_l, self.μhat_a)\n",
    "    val =  self.μhat_a @ X_a_l @ poisson.coef_ - (np.exp(X_a_l @ poisson.coef_)).sum() + self.nbi\n",
    "    return poisson.coef_[:self.nbk-1] / poisson.coef_[self.nbk-1],1 / poisson.coef_[self.nbk-1], val\n",
    "\n",
    "DiscreteChoicePb.mle_glm = DiscreteChoicePb_mle_glm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We verify that both approaches are indeed equivalent:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "GLM approach. βmle_k =  [0.33826722 0.85179332] ; Tmle =  1.8162751318843122  ; ll= -277.7052141446431 .\n",
      "DIY approach. βmle_k =  [0.33826733 0.85179311] ; Tmle =  1.8162760072465682  ; ll= -277.7052141445852 .\n"
     ]
    }
   ],
   "source": [
    "βmleglm_k,Tmleglm,llmleglm = travelEx.mle_glm(max_iter = 10000, tol = 1e-9)\n",
    "print('GLM approach. βmle_k = ',βmleglm_k,'; Tmle = ',Tmleglm, ' ; ll=',llmleglm,'.')\n",
    "print('DIY approach. βmle_k = ',βmle_k,'; Tmle = ',Tmle, ' ; ll=',llmle,'.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Fixed temperature MLE\n",
    "\n",
    "Back to the logit case. Recall we have\n",
    "\n",
    "\\begin{align*}\n",
    "l\\left(  \\tilde{\\beta}\\right)  =N\\left\\{  \\hat{\\mu}^{\\intercal}\\Phi\\tilde{\\beta}-\\sum_i\\log\\sum_{y} \\exp\\left(  \\Phi\\tilde{\\beta}\\right)  _{iy}\\right\\}\n",
    "\\end{align*}\n",
    "\n",
    "Assume that we restrict ourselves to $\\tilde{\\beta}[k]>0$. Then we can define $\\beta =  \\tilde{\\beta} / \\tilde{\\beta}[k]$ and $T=1/ \\tilde{\\beta}[k]$ so we have $\\tilde{\\beta}=\\beta/T$  and $\\beta[k]=1$. Letting $B=\\left\\{  \\beta\\in\\mathbb{R}^{p},\\beta[k]=1\\right\\}  $, so that $\\beta\\in B$. We have for $\\beta \\in B$ and $T>0$\n",
    "\n",
    "\\begin{align*}\n",
    "l\\left(  \\beta,T\\right)  =\\frac{N}{T}\\left\\{  \\hat{\\mu}^{\\intercal}\n",
    "\\Phi\\beta-T\\sum_i\\log\\sum_{y}\\exp\\left(  \\frac{\\left(  \\Phi\\beta\\right)  _{iy}}{T}\\right)  \\right\\}\n",
    "\\end{align*}\n",
    "\n",
    "and we define the *fixed temperature maximum likelihood estimator* by\n",
    "\n",
    "\\begin{align*}\n",
    "\\beta\\left(  T\\right)  =\\arg\\max_{\\beta \\in B}l\\left(  \\beta,T\\right)\n",
    "\\end{align*}\n",
    "\n",
    " Note that $\\beta\\left(  T\\right)  =\\arg\\max_{\\beta\\in B}Tl\\left(\\beta,T\\right)$ where\n",
    "\n",
    "\\begin{align*}\n",
    "Tl\\left(  \\beta,T\\right)  =N\\left\\{  \\hat{\\mu}^{\\intercal}\\Phi\\beta-T\\sum_i\\log\\sum _{y}\\exp\\left(  \\frac{\\left(  \\Phi\\beta\\right)  _{iy}}{T}\\right)  \\right\\}.\n",
    "\\end{align*}\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Normalization and implementation\n",
    "\n",
    "We denote $\\phi[a] = \\Phi[k,a]$ and we have \n",
    "\\begin{align*}\n",
    "\\min_{\\beta,u} \\left\\{ \\sum_{iy} \\hat{\\mu}_{iy} \\left(  \\frac{\\left(  \\phi + \\Phi\\beta - (I_\\mathcal{I} \\otimes 1_\\mathcal{Y}) u \\right)  }{T}\\right)_{iy} - \\sum_{iy} \\exp\\left(  \\frac{\\left(\\phi+  \\Phi\\beta - (I_\\mathcal{I} \\otimes 1_\\mathcal{Y}) u \\right)  }{T}\\right)_{iy}  \\right\\}\n",
    "\\end{align*}\n",
    "\n",
    "which amounts to a weighted Poisson regression with weights $\\exp(\\phi/T)$ and dependent variable $\\hat{\\mu} \\exp(-\\phi/T)$.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have\n",
    "\n",
    "\\begin{align*}\n",
    "\\frac{Tl\\left(  \\beta,T\\right)  }{N}=\\hat{\\mu}^{\\intercal}( \\phi+ \\Phi\\beta)-T\\sum_i \\log\\sum_{y}\\exp\\left(  \\frac{ \\phi_{iy} + \\left(   \\Phi\\beta\\right)  _{iy}}{T}\\right)\n",
    "\\end{align*}\n",
    "\n",
    "One has \n",
    "\n",
    "\\begin{align*}\n",
    "\\beta\\left(  T\\right)  \\in\\arg\\max\\left\\{  \\hat{\\mu}^{\\intercal}(\\phi+ \\Phi\\beta)- \\sum_i u_i\\left(  \\beta\\right)  -T\\sum_i \\log\\sum_{y}\\exp\\left(  \\frac{\\left( \\phi+ \\Phi\\beta\\right)  _{iy}- u_i\\left(  \\beta\\right)  }{T}\\right)  \\right\\}\n",
    "\\end{align*}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Implement as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "def DiscreteChoicePb_mle_fixed_temp_diy (self, T=1):\n",
    "    def minus_log_likelihood(β_k, *args):\n",
    "        T = args[0]\n",
    "        Φβ_i_y = (self.Φ_a_k.dot(np.append(β_k ,1))).reshape((-1,self.nby)) \n",
    "        maxΦβ_i = Φβ_i_y.max(axis = 1)\n",
    "        d_i = np.sum(np.exp((Φβ_i_y -maxΦβ_i[:,None])/T ), axis = 1)\n",
    "        return - ((Φβ_i_y.flatten()*self.μhat_a).sum() / T  -  (maxΦβ_i / T + np.log(d_i)).sum())\n",
    "\n",
    "    def grad_minus_log_likelihood(β_k, *args):\n",
    "        T = args[0]\n",
    "        Φβ_i_y = (self.Φ_a_k.dot(np.append(β_k ,1))).reshape((-1,self.nby)) \n",
    "        maxΦβ_i = Φβ_i_y.max(axis = 1)\n",
    "        d_i = np.sum(np.exp((Φβ_i_y - maxΦβ_i[:,None] )/T ), axis = 1)\n",
    "        μβ_iy = (np.exp((Φβ_i_y - maxΦβ_i[:,None] )/T ) / d_i[:,None]).flatten()\n",
    "        return  - ((self.μhat_a - μβ_iy).reshape((1,-1)) @ self.Φ_a_k).flatten()[:-1]\n",
    "\n",
    "    β0_k = np.zeros(self.nbk-1)\n",
    "    res = optimize.minimize(minus_log_likelihood,method = 'CG',jac = grad_minus_log_likelihood, args = (T,), x0 = β0_k )\n",
    "    β_k = res['x']\n",
    "    return β_k,T,-res['fun']\n",
    "\n",
    "DiscreteChoicePb.mle_fixed_temp_diy = DiscreteChoicePb_mle_fixed_temp_diy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "DIY approach. β2_k =  [0.35608773 0.9495755 ]  ; ll= -277.7440071232399 .\n"
     ]
    }
   ],
   "source": [
    "β2_k ,_,ll2 = travelEx.mle_fixed_temp_diy(2)\n",
    "print('DIY approach. β2_k = ',β2_k, ' ; ll=',ll2,'.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As before, we can also compute the problem using `glm` as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "def DiscreteChoicePb_mle_fixed_temp_glm(self, T=1,max_iter = 1000, tol=0.0001):\n",
    "    X_a_l = spr.hstack([self.Φ_a_k[:,:-1], -spr.kron(spr.identity(self.nbi), np.ones((self.nby,1)))])\n",
    "    poisson = linear_model.PoissonRegressor(alpha = 0, fit_intercept=False,max_iter= max_iter, tol = tol)\n",
    "    Φ_a = self.Φ_a_k[:,-1]\n",
    "    poisson.fit(X_a_l, self.μhat_a * np.exp(-Φ_a / T), sample_weight = np.exp(Φ_a / T) )\n",
    "    β_k = T * poisson.coef_[:(self.nbk-1)]\n",
    "    val =  self.μhat_a @ (X_a_l @ poisson.coef_ + Φ_a / T)- (np.exp( (X_a_l @ poisson.coef_+ Φ_a / T ) ) ).sum() + self.nbi\n",
    "    return β_k, T,val\n",
    "\n",
    "DiscreteChoicePb.mle_fixed_temp_glm = DiscreteChoicePb_mle_fixed_temp_glm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compare the two:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "GLM approach. β2_k =  [0.3499578  0.94529618]  ; ll= -277.74802551270767 .\n",
      "DIY approach. β2_k =  [0.35608773 0.9495755 ]  ; ll= -277.7440071232399 .\n"
     ]
    }
   ],
   "source": [
    "β2glm_k,_ ,ll2glm = travelEx.mle_fixed_temp_glm(2) \n",
    "print('GLM approach. β2_k = ',β2glm_k, ' ; ll=',ll2glm,'.')\n",
    "print('DIY approach. β2_k = ',β2_k, ' ; ll=',ll2,'.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Minimax-regret estimation\n",
    "\n",
    "Let $\\beta\\left(  0\\right)  =\\lim_{T\\rightarrow0}\\beta\\left(T\\right)  $. Calling $u_i\\left(  \\beta\\right)  =\\max_{y\\in\\mathcal{Y}}\\left\\{\\phi_{iy} + \\left(   \\Phi\\beta\\right)  _{iy}\\right\\}  $, we have\n",
    "\n",
    "\\begin{align*}\n",
    "\\beta\\left(  0\\right)  \\in\\arg\\max_{\\beta}\\left\\{  \\hat{\\mu}^{\\intercal}(\\phi+ \\Phi\\beta)-\\sum_i u_i\\left(  \\beta\\right)  \\right\\},\n",
    "\\end{align*}\n",
    "\n",
    "or\n",
    "\n",
    "\\begin{align*}\n",
    "\\beta\\left(  0\\right)  \\in\\arg\\min_{\\beta}\\left\\{ \\sum_i u_i\\left(  \\beta\\right)-\\hat{\\mu}^{\\intercal}(\\phi+ \\Phi\\beta)\\right\\},\n",
    "\\end{align*}\n",
    "\n",
    "Note that $Tl\\left(  \\beta,T\\right)  \\rightarrow N\\left\\{  \\hat{\\mu}^{\\intercal}\\Phi\\beta-\\sum_i\\max_{y\\in\\mathcal{Y}}\\left\\{  \\left(  \\Phi\\beta\\right)_{iy}\\right\\}  \\right\\}  $ as $T\\rightarrow0$. As a result,\n",
    "\n",
    "\\begin{align*}\n",
    "\\beta\\left(  0\\right)  \\in\\arg\\max\\left\\{  \\hat{\\mu}^{\\intercal} (\\phi+ \\Phi\\beta)\n",
    "-\\sum_i u_i\\left(  \\beta\\right)  \\right\\}  .\n",
    "\\end{align*}\n",
    "\n",
    "Define $R_{i}\\left(  \\beta,y\\right)  =\\left(  \\phi+ \\Phi\\beta\\right)_{iy}-\\left( \\phi+ \\Phi\\beta\\right)  _{iy_{i}}$ the regret associated with observation $i$ with respect to $y$. This is equal to the difference between the payoff given by $y$ and the payoff obtained under observation $i$, denoting $y_{i}$ the action taken in observation $i$. The max-regret associated with observation $i$ is therefore\n",
    "\n",
    "\\begin{align*}\n",
    "\\max_{y\\in\\mathcal{Y}}R_{i}\\left(  \\beta,y\\right)  =\\max_{y\\in\\mathcal{Y}}\\left\\{  \\left(  \\phi+ \\Phi\\beta\\right)_{iy}-\\left(  \\phi+\\Phi\\beta\\right)_{iy_{i}}  \\right\\}= u_i(\\beta) - \\sum_y \\hat{\\mu}_{iy} \\left( \\phi+ \\Phi\\beta\\right)_{iy} \n",
    "\\end{align*}\n",
    "\n",
    "and the max-regret associated with the sample is $\\frac{1}{N}\\sum_i \\max_{y\\in\\mathcal{Y}}\\left\\{  R_{i}\\left(  \\beta,y\\right)  \\right\\}  $, that is $\\sum_i u_i(\\beta)  - \\sum_{iy} \\hat{\\mu}_{iy} \\left( \\phi+ \\Phi\\beta\\right)_{iy}$.\n",
    "\n",
    "This leads to the minimax regret estimator\n",
    "\n",
    "\\begin{align*}\n",
    "\\hat{\\beta}^{MMR}=\\min_{\\beta}\\left\\{  \\sum_i u_i(\\beta)  -\\hat\n",
    "{\\mu}^{\\intercal}(\\phi+\\Phi\\beta)\\right\\}\n",
    "\\end{align*}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Linear programming formulation\n",
    "\n",
    "The minimax regret estimator\n",
    "\n",
    "\\begin{align*}\n",
    "\\hat{\\beta}^{MMR}=\\min_{\\beta}\\left\\{  \\sum_i u_i(\\beta)  -\\hat\n",
    "{\\mu}^{\\intercal}\\Phi\\beta\\right\\}\n",
    "\\end{align*}\n",
    "\n",
    "has a linear programming fomulation\n",
    "\n",
    "\\begin{align*}\n",
    "&  \\min_{u_i,\\beta}\\sum_i u_i -\\hat{\\mu}^{\\intercal}(\\phi+\\Phi\\beta)\\\\\n",
    "s.t.~ &  u_i -\\left( \\Phi\\beta\\right)  _{iy}\\geq  \\phi_{iy} ~\\forall i \\in \\mathcal{I} ~\\forall y\\in\\mathcal{Y}\n",
    "\\end{align*}\n",
    "\n",
    "dropping the unnecessary term from the objective function, this becomes\n",
    "\n",
    "\\begin{align*}\n",
    "&  \\min_{u_i,\\beta}\\sum_i u_i -\\hat{\\mu}^{\\intercal}\\Phi\\beta\\\\\n",
    "s.t.~ &  u_i -\\left( \\Phi\\beta\\right)  _{iy}\\geq  \\phi_{iy} ~\\forall i \\in \\mathcal{I} ~\\forall y\\in\\mathcal{Y}\n",
    "\\end{align*}\n",
    "\n",
    "that is\n",
    "\n",
    "\\begin{align*}\n",
    " \\min_{u_i,\\beta}~ &  1^{\\intercal}_\\mathcal{I} u- \\hat{\\mu}^{\\intercal}  \\Phi\\beta \\\\\n",
    "s.t.~ &  (I_\\mathcal{I} \\otimes 1_\\mathcal{Y}) u  - \\Phi\\beta \\geq  \\phi \n",
    "\\end{align*}\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "def DiscreteChoicePb_minimax_regret(self,OutputFlag = False):\n",
    "    Φ_a_k = self.Φ_a_k[:,:-1]\n",
    "    φ_a = self.Φ_a_k[:,-1]\n",
    "    nba,nbk = Φ_a_k.shape\n",
    "    nbi,nby = self.nbi, self.nby\n",
    "    μhat_a = self.μhat_a\n",
    "    \n",
    "    m = grb.Model()\n",
    "    m.setParam( 'OutputFlag', OutputFlag )\n",
    "    grb_β_k = m.addMVar(nbk , lb=-grb.GRB.INFINITY)\n",
    "    grb_u_i = m.addMVar(nbi , lb=-grb.GRB.INFINITY)\n",
    "    m.setObjective(np.ones((1,nbi)) @ grb_u_i -  μhat_a.reshape(1,nba) @ Φ_a_k @ grb_β_k  ,  grb.GRB.MINIMIZE)\n",
    "\n",
    "    m.addConstr( spr.kron(spr.identity(nbi),np.ones((nby,1))) @ grb_u_i - Φ_a_k @ grb_β_k >= φ_a)\n",
    "    m.optimize()\n",
    "    if m.status == grb.GRB.Status.OPTIMAL:\n",
    "        β_k = np.array(m.getAttr('x'))[:nbk]\n",
    "        return β_k,0\n",
    "    \n",
    "DiscreteChoicePb.minimax_regret = DiscreteChoicePb_minimax_regret"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Set parameter Username\n",
      "Academic license - for non-commercial use only - expires 2024-12-08\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(array([0.13810992, 0.03058339]), 0)"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "travelEx.minimax_regret()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Set-identification\n",
    "\n",
    "Note that the set of $\\theta$ that enter the solution to the problem above is not unique, but is a convex set. Denoting $V$ the value of program, we can look for bounds of $\\theta^{\\intercal}d$ for a chosen direction $d$ by\n",
    "\n",
    "\\begin{align*}\n",
    "& \\min_{\\beta,u}/\\max_{\\beta,u}   \\beta^{\\intercal}d\\\\\n",
    "s.t.~  &  1^{\\intercal}_\\mathcal{I} u- \\hat{\\mu}^{\\intercal}  \\Phi\\beta =V\\\\\n",
    "&  (I_\\mathcal{I} \\otimes 1_\\mathcal{Y}) u  - \\Phi\\beta \\geq  \\phi \n",
    "\\end{align*}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Plotting the whole path"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>T</th>\n",
       "      <th>β_1</th>\n",
       "      <th>β_2</th>\n",
       "      <th>logLikelihood</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.138110</td>\n",
       "      <td>0.030583</td>\n",
       "      <td>-inf</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0.036326</td>\n",
       "      <td>0.142615</td>\n",
       "      <td>0.027921</td>\n",
       "      <td>-2127.900288</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>0.072651</td>\n",
       "      <td>0.145240</td>\n",
       "      <td>0.034438</td>\n",
       "      <td>-1082.761589</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>0.108977</td>\n",
       "      <td>0.152885</td>\n",
       "      <td>0.047561</td>\n",
       "      <td>-743.861319</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>0.145302</td>\n",
       "      <td>0.162268</td>\n",
       "      <td>0.066417</td>\n",
       "      <td>-581.784325</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>0.181628</td>\n",
       "      <td>0.168505</td>\n",
       "      <td>0.074645</td>\n",
       "      <td>-489.953271</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>0.217953</td>\n",
       "      <td>0.175305</td>\n",
       "      <td>0.087783</td>\n",
       "      <td>-432.647199</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>0.254279</td>\n",
       "      <td>0.181414</td>\n",
       "      <td>0.100840</td>\n",
       "      <td>-394.515489</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "          T       β_1       β_2  logLikelihood\n",
       "0  0.000000  0.138110  0.030583           -inf\n",
       "1  0.036326  0.142615  0.027921   -2127.900288\n",
       "2  0.072651  0.145240  0.034438   -1082.761589\n",
       "3  0.108977  0.152885  0.047561    -743.861319\n",
       "4  0.145302  0.162268  0.066417    -581.784325\n",
       "5  0.181628  0.168505  0.074645    -489.953271\n",
       "6  0.217953  0.175305  0.087783    -432.647199\n",
       "7  0.254279  0.181414  0.100840    -394.515489"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "indMax=100\n",
    "tempMax= 2*Tmle\n",
    "outcomemat = np.zeros((indMax+1,travelEx.nbk+1))\n",
    "\n",
    "outcomemat[0,1:-1],_ = travelEx.minimax_regret()\n",
    "outcomemat[0,-1] = -np.inf\n",
    "iterMax = indMax+1\n",
    "for k in range(2,iterMax+1,1):\n",
    "    thetemp = tempMax * (k-1)/ indMax\n",
    "    outcomeFixedTemp,_,ll = travelEx.mle_fixed_temp_diy(thetemp)\n",
    "    outcomemat[k-1,0] = thetemp\n",
    "    outcomemat[k-1,1:-1] = outcomeFixedTemp \n",
    "    outcomemat[k-1,-1] = ll \n",
    "df = pd.DataFrame(outcomemat, columns=['T', 'β_1', 'β_2' , 'logLikelihood'])\n",
    "df.head(8)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}