random-utility-models.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# <center>Block 3b: Random utility models</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": [
    "### Learning objectives\n",
    "\n",
    "* Emax operator and generalized entropy of choice\n",
    "\n",
    "* The Daly-Zachary-Williams theorem\n",
    "\n",
    "* The MEV class"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### References\n",
    "\n",
    "* Galichon (2016). *Optimal Transport Methods in Economics*, App. E. Princeton University Press.\n",
    "\n",
    "* McFadden (1981). \"Econometric Models of Probabilistic Choice\", in C.F. Manski and D. McFadden (eds.), *Structural Analysis of Discrete Data with Econometric Applications*, MIT Press.\n",
    "\n",
    "* McFadden (1989). \"A Method of Simulated Moments for Estimation of Discrete Response Models Without Numerical\n",
    "Integration\". *Econometrica*.\n",
    "\n",
    "* Berry, Levinsohn, and Pakes (1995). \"Automobile Prices in Market Equilibrium,\" *Econometrica*.\n",
    "\n",
    "* Train. (2009). *Discrete Choice Methods with Simulation*. 2nd Edition. Cambridge University Press.\n",
    "\n",
    "* Galichon and Salanie (2020). \"Cupid's Invisible Hands\". Preprint (first version 2011).\n",
    "\n",
    "* Chiong, Galichon and Shum (2016), \"Duality in Discrete Choice Models\". *Quantitative Economics*.\n",
    "\n",
    "* Galichon. \"On the representation of the nested logit model\". *Econometric Theory* (forthcoming). \n",
    "\n",
    "* Greene and Hensher (1997), *Multinomial Logit and Discrete Choice Models*."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## A look at our data\n",
    "The data is taken from Greene and Hensher (1997). 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",
    "---\n",
    "\n",
    "First, let's load the libraries we will need and let us get the data files."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "\n",
    "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['choice'] = np.where(travelmode['choice'] =='yes' , 1, 0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our data looks as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "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>0</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>0</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>0</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>1</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>0</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       0    69     59     100     70      35     1\n",
       "1           1  train       0    34     31     372     71      35     1\n",
       "2           1    bus       0    35     25     417     70      35     1\n",
       "3           1    car       1     0     10     180     30      35     1\n",
       "4           2    air       0    64     58      68     68      30     2"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "travelmode.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Each individual is surveyed about their perceiption of various costs associated with each of the four alternatives. Hence each individual appears four times in the data set, and there are 210 observations in total:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "210"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "nobs,ncols = travelmode.shape\n",
    "nbchoices = 4\n",
    "ninds = int(nobs/nbchoices)\n",
    "ninds"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We consutruct a matrix `muhat_i_y` of dummies for various choices (1 if individual `i` chooses alternative `y`, 0 else)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "muhat_i_y = travelmode['choice'].values.reshape(ninds,4).T\n",
    "muhat_iy = muhat_i_y.flatten()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We compute the inconditional market shares associated with each of the four travel modes:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "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>mode</th>\n",
       "      <th>air</th>\n",
       "      <th>bus</th>\n",
       "      <th>car</th>\n",
       "      <th>train</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>choice</th>\n",
       "      <td>0.27619</td>\n",
       "      <td>0.142857</td>\n",
       "      <td>0.280952</td>\n",
       "      <td>0.3</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "mode        air       bus       car  train\n",
       "choice  0.27619  0.142857  0.280952    0.3"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "q_y = travelmode.groupby(['mode']).mean()['choice'].to_frame().sort_index()\n",
    "q_y.transpose()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Emax operators and demand maps"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Discrete choice models\n",
    "\n",
    "\n",
    "Assume a consumer is facing a number of options $y\\in\\mathcal{Y}_{0}=\\mathcal{Y}\\cup\\left\\{  0\\right\\}$, where $y=0$ is a default option. The consumer is drawing a utility shock which is a vector $\\varepsilon =(\\varepsilon_{0},\\ldots,\\varepsilon_{\\left\\vert \\mathcal{Y}\\right\\vert })\\sim\\mathbf{P}$ such that the utility of option $y$ is $U_{y}+\\varepsilon _{y}$, while the outside option yields utility $\\varepsilon_{0}$.\n",
    "\n",
    "$U$ is called vector of *systematic utilities*; $\\varepsilon$ is called vector of *utility shocks*.\n",
    "\n",
    "We assume throughout that $\\mathbf{P}$ has a density with respect to the Lebesgue measure, and has full support.\n",
    "\n",
    "The preferred option is the one which attains the maximum in\n",
    "\n",
    "\\begin{align*}\n",
    "\\max_{y\\in\\mathcal{Y}}\\left\\{  U_{y}+\\varepsilon_{y},\\varepsilon_{0}\\right\\}.\n",
    "\\end{align*}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Demand map\n",
    "\n",
    "Let $q_y=\\sigma_{y}\\left(  U\\right)  $ be the probability of choosing option $y$, where $\\sigma$ is given by \n",
    "\n",
    "\\begin{align*}\n",
    "\\sigma_{y}\\left(  U\\right)  =\\Pr(U_{y}+\\varepsilon_{y}\\geq U_{z}+\\varepsilon_{z}\\mbox{ for all }z\\in\\mathcal{Y}_{0}).\n",
    "\\end{align*}\n",
    "\n",
    "The map $\\sigma$ is called *demand map*, and the vector $s$ is called vector of market shares, or vector of choice probabilities.\n",
    "\n",
    "Note that if $s=\\sigma(U)$, then $q_y>0$ for all $y\\in\\mathcal{Y}_{0}$ and $\\sum_{y\\in\\mathcal{Y}_{0}}q_y=1$.\n",
    "\n",
    "Note that because the distribution $\\mathbf{P}$ of $\\varepsilon$ is continuous, the probability of being indifferent between two options is zero, and hence we could have indifferently replaced weak preference $\\geq$ by strict preference $>$. Without this, choice probabilities may not have been well defined."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Properties\n",
    "\n",
    "* $\\sigma_{y}\\left(  U\\right)  $ is increasing in $U_{y}$.\n",
    "\n",
    "* $\\sigma_{y}\\left(  U\\right)  $ is weakly decreasing in $U_{y^{\\prime}}$ for $y^{\\prime}\\neq y$.\n",
    "\n",
    "* If one replaces $\\left(  U_{y}\\right)  $ by $\\left(  U_{y}+c\\right)  $, for a constant $c$, one has $\\sigma\\left(  U+c\\right)  =\\sigma\\left( U\\right)  .$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Normalization\n",
    "\n",
    "Because of the last property, we can normalize the utility of one of the alternatives. We will normalize the utility of the utility associated to $y=0$, and hence take\n",
    "\n",
    "\\begin{align*}\n",
    "U_{0}=0.\n",
    "\\end{align*}\n",
    "\n",
    "Thus in the sequel, $\\sigma$ will be seen as a mapping from $\\mathbb{R}^{\\mathcal{Y}}$ to the set of $\\left(  q_y\\right)_{y\\in \\mathcal{Y}}$ such that $q_y>0$ and $\\sum_{y\\in\\mathcal{Y}}q_y<1$, and the choice probability of alternative $y=0$ is recovered by\n",
    "\n",
    "\\begin{align*}\n",
    "s_{0}=1-\\sum_{y\\in\\mathcal{Y}}q_y.\n",
    "\\end{align*}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The Daly-Williams-Zachary theorem\n",
    "\n",
    "Define the expected indirect utility of consumers by\n",
    "\n",
    "\\begin{align*}\n",
    "G(U)=\\mathbb{E}\\left[  \\max_{y\\in\\mathcal{Y}}(U_{y}+\\varepsilon_{y}%\n",
    ",\\varepsilon_{0})\\right]\n",
    "\\end{align*}\n",
    "\n",
    "This is called *Emax operator*, a.k.a. *McFadden's surplus function*.\n",
    "\n",
    "As the expectation of the maximum of terms which are linear in $U$, $G$ is convex function in $U$ (strictly convex in fact), and\n",
    "\n",
    "\\begin{align*}\n",
    "\\frac{\\partial G}{\\partial U_{y}}(U)=\\Pr(U_{y}+\\varepsilon_{y}\\geq\n",
    "U_{z}+\\varepsilon_{z}\\mbox{ for all }z\\in\\mathcal{Y}_{0}).\n",
    "\\end{align*}\n",
    "\n",
    "But the right-hand side is simply the probability $q_y$ of choosing option $y$; therefore, we get:\n",
    "\n",
    "---\n",
    "\n",
    "**Theorem (Daly-Zachary-Williams)**\n",
    "\n",
    "The demand map $\\sigma$ is the gradient of the Emax operator $G$, that is\n",
    "\n",
    "<a name='eq:derGU'></a>\n",
    "\\begin{align*}\n",
    "\\sigma\\left(  U\\right)  =\\nabla G(U).\n",
    "\\end{align*}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Examples"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Model 1: the logit model\n",
    "\n",
    "Assume that $\\mathbf{P}$ is the distribution of i.i.d. *centered type $1$ extreme value* a.k.a. *centered Gumbel* terms, which has c.d.f.\n",
    "\n",
    "\\begin{align*}\n",
    "F\\left(  z\\right)  =\\exp\\left(  -\\exp\\left(  -x+\\gamma\\right)  \\right)\n",
    "\\end{align*}\n",
    "\n",
    "where $\\gamma=0.5772...$ (Euler's constant). The mean of this distribution is zero.\n",
    "\n",
    "Basic fact from extreme value theory: if $(\\varepsilon_{y})$ have i.i.d. Gumbel distributions, then $\\max\\left\\{U_{y}+\\varepsilon_{y}\\right\\}  $ has the same distribution as $\\log\\left(\\sum_{y=1}^{n}\\exp U_{y}\\right)  +\\epsilon$, where $\\epsilon$ is also a Gumbel. (Proof of this fact later).\n",
    "\n",
    "Notes:\n",
    "\n",
    "* This distribution is sometimes called the \"Gumbel max\" distribution, to contrast it with the distribution of its opposite, which is then called *Gumbel min*.\n",
    "\n",
    "* The literature usually calls *standard Gumbel* the distribution with c.d.f. $\\exp\\left(-\\exp\\left(  -x\\right)  \\right)  $; but that distribution has mean $\\gamma$, which is why we slightly depart from the convention."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Emax operator associated with the logit model can be given in closed form as\n",
    "\n",
    "\\begin{align*}\n",
    "G(U)=\\log\\left(  1+\\sum_{y\\in\\mathcal{Y}}\\exp(U_{y})\\right)\n",
    "\\end{align*}\n",
    "\n",
    "where $s_{0}=1-\\sum_{y\\in\\mathcal{Y}}q_y$. This is called a *log-partition function*\n",
    "\n",
    "As a result, the choice probability of alternative $y$ is proportional to the exponential of the systematic utility associated with $U$, that is \n",
    "\n",
    "\\begin{align*}\n",
    "\\sigma_{y}\\left(  U\\right)  =\\frac{\\exp U_{y}}{1+\\sum_{y^{\\prime}%\n",
    "\\in\\mathcal{Y}}\\exp(U_{y^{\\prime}})}\n",
    "\\end{align*}\n",
    "\n",
    "which is called a *Gibbs distribution*.\n",
    "\n",
    "* Assume that the random utility shock is scaled by a factor $T$. Then\n",
    "\n",
    "\\begin{align*}\n",
    "\\sigma_{y}\\left(  U\\right)  =\\frac{\\exp\\left(  U_{y}/T\\right)  }{1+\\sum_{y^{\\prime}\\in\\mathcal{Y}}\\exp(U_{y^{\\prime}}/T)}\n",
    "\\end{align*}\n",
    "\n",
    "which is sometimes called the *soft-max operator* and converges as $T\\rightarrow0$ toward\n",
    "\n",
    "\\begin{align*}\n",
    "\\max_{y\\in\\mathcal{Y}}\\left\\{  U_{y},0\\right\\}  .\n",
    "\\end{align*}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Generalizing the logit model: the multivariate extreme value (MEV) class\n",
    "\n",
    "**See class notes for an introduction using Pickand's representation.**\n",
    "\n",
    "Let $\\mathbf{F}$ be a cumulative distribution such that function $g$ defined by\n",
    "\n",
    "\\begin{align*}\n",
    "g\\left(  x_{1},...,x_{n}\\right)  =-\\log\\mathbf{F}\\left(  -\\log x_{1},...,-\\log\n",
    "x_{n}\\right)  \\label{cdfMEV}\n",
    "\\end{align*}\n",
    "\n",
    "is positive homogeneous of degree $1$. (This inverts into $\\mathbf{F}\\left(u_{1},...,u_{n}\\right)  =\\exp\\left(  -g\\left(  e^{-u_{1}},...,e^{-u_{n} }\\right)  \\right)  $). We have by a theorem of McFadden (1978):\n",
    "\n",
    "---\n",
    "**Theorem** \n",
    "\n",
    "Let $\\left(  \\varepsilon_{y}\\right)  _{1\\leq y\\leq n}$ be a random vector with c.d.f. $\\mathbf{F}$, and define\n",
    "\n",
    "\\begin{align*}\n",
    "Z=\\max_{y=1,...,n}\\left\\{  U_{y}+\\varepsilon_{y}\\right\\}  .\n",
    "\\end{align*}\n",
    "\n",
    "Then $Z$ has the same distribution as $\\log g\\left(  e^{U_{1}},...,e^{U_{n}}\\right)  +\\gamma+\\epsilon$, where $\\epsilon$ is a standard Gumbel. In particular,\n",
    "\\begin{align*}\n",
    "\\mathbb{E}\\left[  \\max_{y=1,...,n}\\left\\{  U_{y}+\\varepsilon_{y}\\right\\}\n",
    "\\right]  =\\log g\\left(  e^{U_{1}},...,e^{U_{n}}\\right)  +\\gamma\n",
    "\\end{align*}\n",
    "where $\\gamma$ is the Euler constant $\\gamma\\simeq0.5772$.\n",
    "\n",
    "---\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Proof**\n",
    "\n",
    "Let $F_{Z}$ be the c.d.f. of $Z=\\max_{y=1,\\ldots,n}\\left\\{  U_{y}+\\varepsilon_{y}\\right\\}  $. One has \n",
    "\n",
    "\\begin{align*}\n",
    "F_{Z}\\left(  z\\right)   &  =\\Pr\\left(  \\max_{y=1,\\ldots,n}\\left\\{\n",
    "U_{y}+\\varepsilon_{y}\\right\\}  \\leq z\\right)  =\\Pr\\left(  \\forall\n",
    "y:~\\varepsilon_{y}\\leq z-U_{y}\\right) \\\\\n",
    "&  =\\mathbf{F}\\left(  z-U_{1},...,z-U_{n}\\right)  =\\exp\\left(  -g\\left(\n",
    "e^{U_{1}-z},...,e^{U_{n}-z}\\right)  \\right) \\\\\n",
    "&  =\\exp\\left(  -e^{-z}g\\left(  e^{U_{1}},...,e^{U_{n}}\\right)  \\right)\n",
    "=\\varphi\\left(  z-\\log g\\left(  e^{U_{1}},...,e^{U_{n}}\\right)  -\\gamma\n",
    "\\right)\n",
    "\\end{align*}\n",
    "\n",
    "where $\\varphi\\left(  z\\right)  :=\\exp\\left(  -e^{-\\left(  z-\\gamma\\right) }\\right)  $ is the cdf of the standard Gumbel distribution. Hence $Z$ has the distribution of $\\log g\\left(  e^{U_{1}},...,e^{U_{n}}\\right)  +\\gamma +\\epsilon$, where $\\epsilon$ is a standard Gumbel."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As a result, the choice probability of alternative $y$ is\n",
    "\n",
    "\\begin{align*}\n",
    "\\sigma_{y}\\left(  U\\right)  =\\frac{\\frac{\\partial g}{\\partial x_{y}}\\left(e^{U_{1}},...,e^{U_{n}}\\right)  }{g\\left(  e^{U_{1}},...,e^{U_{n}}\\right)\n",
    "}e^{U_{y}}.\n",
    "\\end{align*}\n",
    "\n",
    "The MEV\\ framework has several commonly used examples: logit, nested logit, mixture of logit...\n",
    "\n",
    "We just saw the logit model, in which $g\\left(  x_{1},...,x_{n}\\right) =e^{-\\gamma}\\sum_{y=1}^{n}x_{y}$. In this case, the distribution of\n",
    "\n",
    "\\begin{align*}\n",
    "Z=\\max_{y=1,\\ldots,n}\\left\\{  U_{y}+\\varepsilon_{y}\\right\\}\n",
    "\\end{align*}\n",
    "\n",
    "is $\\log\\sum_{y=1}^{n}e^{U_{y}}+\\epsilon$, where $\\epsilon$ is a standard Gumbel."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Model 2: the nested logit model\n",
    "\n",
    "The nested logit model is an instance of MEV model where alternatives can be grouped in nests. Eg, people choose their means of transportation (nest), and within this nest, a particular operator.\n",
    "\n",
    "**See class notes for an introduction using positive-stable distributions representation.**\n",
    "\n",
    "\n",
    "Let $\\mathcal{X}$ be the set of nests and assume that for each nest $x$, there is a set $\\mathcal{Y}_{x}$ alternatives. Let $U_{xy}$ be utility from alternative $y$ in nest $x$, and $\\lambda_{x}\\in\\left[  0,1\\right]  $ and\n",
    "\n",
    "\\begin{align*}\n",
    "g\\left(  U_{xy}\\right)  =e^{-\\gamma}\\sum_{x\\in\\mathcal{X}}\\left(  \\sum\n",
    "_{y\\in\\mathcal{Y}_{x}}U_{xy}^{1/\\lambda_{x}}\\right)  ^{\\lambda_{x}}.\n",
    "\\end{align*}\n",
    "\n",
    "In this case,\n",
    "\n",
    "\\begin{align*}\n",
    "G\\left(  U\\right)   &  =\\mathbb{E}\\left[  \\max_{x\\in\\mathcal{X}}\\max\n",
    "_{y\\in\\mathcal{Y}_{x}}\\left\\{  U_{xy}+\\varepsilon_{xy}\\right\\}  \\right] =\\log\\sum_{x\\in\\mathcal{X}}\\left(  \\sum_{y\\in\\mathcal{Y}_{x}}e^{U_{xy}/\\lambda_{x}}\\right)  ^{\\lambda_{x}}\\\\\n",
    "\\sigma_{xy}\\left(  U\\right)   &  =\\frac{\\left(  \\sum_{y\\in\\mathcal{Y}_{x}}e^{U_{xy}/\\lambda_{x}}\\right)  ^{\\lambda_{x}}}{\\sum_{x\\in\\mathcal{X}}\\left(\\sum_{y\\in\\mathcal{Y}_{x}}e^{U_{xy}/\\lambda_{x}}\\right)  ^{\\lambda_{x}}}\\frac{e^{U_{xy}/\\lambda_{x}}}{\\left(  \\sum_{y\\in\\mathcal{Y}_{x}}e^{U_{xy}/\\lambda_{x}}\\right)}\n",
    "\\end{align*}\n",
    "\n",
    "so the demand map has an interesting interpretation as \"choice of nest then choice of alternative\"."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Assume that $\\left(  \\varepsilon_{1},\\varepsilon_{2}\\right)  $ have a nested logit distribution with two nests, that is, their cdf is given by \n",
    "\n",
    "\\begin{align*}\n",
    "\\mathbf{F}\\left(  u_{1},u_{2}\\right)  =\\exp\\left(  -e^{-\\gamma}\\left(e^{-u_{1}/\\lambda}+e^{-u_{2}/\\lambda}\\right)  ^{\\lambda}\\right)  .\n",
    "\\end{align*}\n",
    "\n",
    "Particular cases:\n",
    "\n",
    "* When $\\lambda=1$, $\\varepsilon_{1}$ and $\\varepsilon_{2}$ are independent and one recovers the logit model.\n",
    "\n",
    "* When $\\lambda\\rightarrow0$, $\\mathbf{F}\\left(  u_{1},u_{2}\\right) =\\exp\\left(  -e^{-\\gamma}e^{\\max\\left\\{  -u_{1},-u_{2}\\right\\}  }\\right)=\\min\\left\\{  \\mathbf{F}\\left(  u_{1}\\right)  ,\\mathbf{F}\\left(  u_{2}\\right) \\right\\}  $ and therefore $\\varepsilon_{1}$ and $\\varepsilon_{2}$ are perfectly correlated.\n",
    "\n",
    "In general one can show that\n",
    "\n",
    "\\begin{align*}\n",
    "\\lambda=\\sqrt{1-cor\\left(  \\varepsilon_{1},\\varepsilon_{2}\\right)  }%\n",
    "\\end{align*}\n",
    "\n",
    "This formula, due to Tiago de Oliviera, is now easy to prove using the positive-stable representation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Other popular models (beyond MEV)\n",
    "\n",
    "* Probit model (later)\n",
    "\n",
    "* Berry-Pakes' pure characteristics model (later)\n",
    "\n",
    "* Berry-Levinsohn-Pakes' mixed logit coefficient model (later)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Demand inversion"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In many settings, the econometrician observes the market shares $q_y$ and wants to deduce the corresponding vector of systematic utilities. That is, we would like to solve:\n",
    "\n",
    "---\n",
    "**Problem**. \n",
    "\n",
    "Given a vector $s$ with positive entries satisfying $\\sum_{y\\in\\mathcal{Y}}q_y<1$, characterize and compute the set\n",
    "\n",
    "\\begin{align*}\n",
    "\\sigma^{-1}\\left(  s\\right)  =\\left\\{  U\\in\\mathbb{R}^{\\mathcal{Y}}%\n",
    ":\\sigma\\left(  U\\right)  =s\\right\\}  .\n",
    "\\end{align*}\n",
    "\n",
    "---\n",
    "\n",
    "This problem is called \"demand inversion\", or \"conditional choice\n",
    "probability inversion\", or \"identification problem\". It is a central issue in econometrics/industrial organization and will be a key building block for matching models."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Demand inversion via convex analysis\n",
    "\n",
    "Galichon and Salanié (2011) suggested to use convex analysis for the inversion of discrete choice problems. Indeed, if $G$ is strictly convex and $C^{1}$, then\n",
    "\n",
    "\\begin{align*}\n",
    "\\sigma^{-1}\\left(  s\\right)  =\\nabla G^{-1}(s)=\\nabla G^{\\ast}\\left(\n",
    "s\\right),\n",
    "\\end{align*}\n",
    "\n",
    "which leads Galichon and Salanié to define the *entropy of choice*  by\n",
    "\n",
    "\\begin{align*}\n",
    "G^{\\ast}(s)=\\max_{U}\\left\\{  \\sum_{y\\in\\mathcal{Y}}q_yU_{y}-G(U)\\right\\}.\n",
    "\\label{eq:constrG}%\n",
    "\\end{align*}\n",
    "\n",
    "In terms of convex analysis, $G^{\\ast}$ is the Legendre-Fenchel transform of $G$. It follows from the envelope theorem that $\\sigma^{-1}\\left(  s\\right)  $ is the vector $U$ such that\n",
    "\n",
    "\\begin{align*}\n",
    "U\\in\\arg\\max_{U}\\left\\{  \\sum_{y\\in\\mathcal{Y}}q_yU_{y}-G(U)\\right\\} .\n",
    "\\end{align*}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Entropy of choice\n",
    "\n",
    "Convex duality implies that if $s$ and $U$ are related by $s\\in\\partial G\\left(  U\\right)$, then\n",
    "\n",
    "\\begin{align*}\n",
    "G(U)=\\sum_{y\\in\\mathcal{Y}}q_yU_{y}-G_{x}^{\\ast}(s). \\label{fenchel}%\n",
    "\\end{align*}\n",
    "\n",
    "But letting $Y=\\arg\\max_{y}\\left\\{  U_{y}+\\varepsilon_{y}\\right\\}$, $G\\left(  U\\right)  =\\mathbb{E}\\left[  U_{Y}+\\varepsilon_{Y}\\right]$ implies\n",
    "\n",
    "\\begin{align*}\n",
    "G(U)=\\sum_{y\\in\\mathcal{Y}}q_yU_{y}+\\mathbb{E}\\left[  \\varepsilon\n",
    "_{Y}\\right]  ,\n",
    "\\end{align*}\n",
    "\n",
    "thus one has\n",
    "\n",
    "\\begin{align*}\n",
    "G^{\\ast}(s)=-\\mathbb{E}\\left[  \\varepsilon_{Y}\\right]  . \\label{Fenchel}%\n",
    "\\end{align*}\n",
    "\n",
    "Hence, the entropy of choice $G^{\\ast}\\left(  s\\right)  $ is interpreted as minus\\ the expected amount of heterogeneity needed to rationalize the choice probabilities $s$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Examples of entropy of choice and identification\n",
    "### Logit model\n",
    "\n",
    "Then\n",
    "\n",
    "\\begin{align*}\n",
    "G^{\\ast}\\left(  s\\right)  =s_{0}\\log(s_{0})+\\sum_{y\\in\\mathcal{Y}}q_y\\log q_y\n",
    "\\end{align*}\n",
    "\n",
    "where $s_{0}=1-\\sum_{y\\in\\mathcal{Y}}q_y$. Hence, $G^{\\ast}$ is a bona fide entropy function when $\\mathbf{P}$ is Gumbel--hence the name of *entropy  of choice* in general.\n",
    "\n",
    "As a result,\n",
    "\n",
    "\\begin{align*}\n",
    "\\sigma_{y}^{-1}\\left(  s\\right)  =\\log\\frac{q_y}{s_{0}}\n",
    "\\end{align*}\n",
    "\n",
    "which is the celebrated *log-odds ratio formula*: the log of the odds of alternatives $y$ and $0$ identify the difference between the systematic utilities of these alternatives."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Back to our dataset**. Ignoring any of the information we have on individual, and thus assuming there is a representative consumer with logit random utility, we shall retrieve the systematic utilities from the market shares. \n",
    "\n",
    "Define \"car\" as the default alternative. The utilities in the logit model are obtained by the log-odds ratio formula:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "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>mode</th>\n",
       "      <th>air</th>\n",
       "      <th>bus</th>\n",
       "      <th>car</th>\n",
       "      <th>train</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>choice</th>\n",
       "      <td>-0.017094</td>\n",
       "      <td>-0.67634</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.065597</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "mode         air      bus  car     train\n",
       "choice -0.017094 -0.67634  0.0  0.065597"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Ulogit = np.log(q_y/q_y.loc[['car']].values[0])\n",
    "Ulogit.transpose()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Nested logit model\n",
    "The entropy of choice $G^{\\ast}$ in the nested logit model is given by\n",
    "\n",
    "\\begin{align*}\n",
    "G^{\\ast}\\left(  s\\right)  =\\sum_{x\\in\\mathcal{X}}\\sum_{y\\in\\mathcal{Y}_{x}}\\lambda_{x}s_{xy}\\ln s_{xy}+\\sum_{x\\in\\mathcal{X}}\\left(  1-\\lambda_{x}\\right)  \\left(  \\sum_{z\\in\\mathcal{Y}_{x}}s_{xz}\\right)  \\ln\\left(\\sum_{z\\in\\mathcal{Y}_{x}}s_{xz}\\right)  \\label{Gstarnestedlogit}\n",
    "\\end{align*}\n",
    "\n",
    "if $s_{xy}\\geq0$ and $\\sum_{x\\in\\mathcal{X}}\\sum_{y\\in\\mathcal{Y}_{x}}s_{xy}=1$, $G^{\\ast}\\left(  s\\right)  =+\\infty$ otherwise.\n",
    "\n",
    "Identification in the nested logit model: with normalization $\\sum _{x\\in\\mathcal{X}}\\left(  \\sum_{y\\in\\mathcal{Y}_{x}}e^{U_{xy}/\\lambda_{x}}\\right)  ^{\\lambda_{x}}=1$, one has $s_{xy}=\\left(  \\sum_{y\\in\\mathcal{Y}_{x}}e^{U_{xy}/\\lambda_{x}}\\right)  ^{\\lambda_{x}-1}e^{U_{xy}/\\lambda_{x}}$, thus $\\sum_{y\\in\\mathcal{Y}_{x}}e^{U_{xy}/\\lambda_{x}}=\\left(  \\sum _{y\\in\\mathcal{Y}_{x}}s_{xy}\\right)  ^{1/\\lambda_{x}}$, therefore\n",
    "\n",
    "\\begin{align*}\n",
    "U_{xy}=\\lambda_{x}\\log s_{xy}-\\left(  \\lambda_{x}-1\\right)  \\log\\sum\n",
    "_{y\\in\\mathcal{Y}_{x}}s_{xy}.\n",
    "\\end{align*}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Back to our dataset**. Still ignoring any observation about individual agents, let us retrieve the systematic utilities from the market shares in the nested logit model with two nests, \"noncar\" and \"car\", and taking $\\lambda=0.5$ in both nests. Do:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "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>mode</th>\n",
       "      <th>air</th>\n",
       "      <th>bus</th>\n",
       "      <th>train</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>choice</th>\n",
       "      <td>-0.808246</td>\n",
       "      <td>-1.137869</td>\n",
       "      <td>-0.7669</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "mode         air       bus   train\n",
       "choice -0.808246 -1.137869 -0.7669"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "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>mode</th>\n",
       "      <th>car</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>choice</th>\n",
       "      <td>-1.26957</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "mode        car\n",
       "choice -1.26957"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "l = [0.5, 0.5]\n",
    "nest = Ulogit.index.isin(['car'])\n",
    "Unocar = q_y[~nest].apply(np.log).multiply(l[0]).add((1-l[0])*np.sum(q_y[~nest]).apply(np.log))\n",
    "display(Unocar.transpose())\n",
    "Ucar = q_y[nest].apply(np.log).multiply(l[1]).add((1-l[1])*np.sum(q_y[nest]).apply(np.log))\n",
    "display(Ucar.transpose())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This yields the systematic utilities (nested logit):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "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>mode</th>\n",
       "      <th>air</th>\n",
       "      <th>bus</th>\n",
       "      <th>car</th>\n",
       "      <th>train</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>choice</th>\n",
       "      <td>0.461324</td>\n",
       "      <td>0.131701</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.50267</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "mode         air       bus  car    train\n",
       "choice  0.461324  0.131701  0.0  0.50267"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Unested = pd.concat([Unocar,Ucar]).add(-Ucar.values[0]).sort_index()\n",
    "Unested.transpose()"
   ]
  }
 ],
 "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
}