discrete-choice-estimation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# <center>Block 4a: Discrete choice estimation</center>\n",
    "### <center>Alfred Galichon (NYU & Sciences Po)</center>\n",
    "## <center>'math+econ+code' masterclass on optimal transport and economic applications</center>\n",
    "#### <center>With python code examples</center>\n",
    "© 2018-2022 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/theteam).\n",
    "\n",
    "**If you reuse material from this masterclass, please cite as:**<br>\n",
    "Alfred Galichon, 'math+econ+code' masterclass on optimal transport and economic applications, January 2022. https://github.com/math-econ-code/mec_optim"
   ]
  },
  {
   "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": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# !python -m pip install -i https://pypi.gurobi.com gurobipy ## only if Gurobi not here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "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": null,
   "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. 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": null,
   "metadata": {},
   "outputs": [],
   "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": null,
   "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": null,
   "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": null,
   "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": null,
   "metadata": {},
   "outputs": [],
   "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": null,
   "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": null,
   "metadata": {},
   "outputs": [],
   "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": null,
   "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": null,
   "metadata": {},
   "outputs": [],
   "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": null,
   "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": null,
   "metadata": {},
   "outputs": [],
   "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": null,
   "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": null,
   "metadata": {},
   "outputs": [],
   "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": null,
   "metadata": {},
   "outputs": [],
   "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",
   "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.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}