characteristics-models.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# <center>Block 7b: Characteristics-based models of demand</center>\n",
    "### <center>Alfred Galichon (NYU & Sciences Po)</center>\n",
    "## <center>'math+econ+code' masterclass on optimal transport and economic applications</center>\n",
    "<center>© 2018-2022 by Alfred Galichon. Past and present support from NSF grant DMS-1716489, ERC grant CoG-866274, and contributions by Jules Baudet, Pauline Corblet, Gregory Dannay, and James Nesbit are acknowledged.</center>\n",
    "\n",
    "#### <center>With python code examples</center>\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 2021. https://github.com/math-econ-code/mec_optim"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Learning objectives\n",
    "\n",
    "* Beyond GEV: the pure characteristics models, the random coefficient logit model, the probit model\n",
    "\n",
    "* Simulation methods: Accept-Reject and SARS\n",
    "\n",
    "* Demand inversion: The inversion theorem"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### References\n",
    "\n",
    "* Galichon (2016). *Optimal Transport Methods in Economics*, Chapter 9.2, 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",
    "* 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",
    "* Bonnet, Galichon, O'Hara and Shum (2017). \"Yogurts choose consumers? Identification of Random Utility Models via Two-Sided Matching\". Working paper."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Libraries\n",
    "\n",
    "Let's start by loading the libraries we shall need for this course."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy.sparse as spr\n",
    "# !python -m pip install -i https://pypi.gurobi.com gurobipy ## only if Gurobi not here\n",
    "import gurobipy as grb"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Datasets\n",
    "\n",
    "Now we load the dataset we shall use. \n",
    "\n",
    "XXX Load BLP's dataset here XXX."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Choice models beyond GEV"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The need for further models\n",
    "\n",
    "The GEV models are convenient analytically, but not very flexible.\n",
    "\n",
    "The logit model imposes zero correlation across alternatives\n",
    "\n",
    "The nested logit allows for nonzero correlation, but in a very rigid way (needs to define nests).\n",
    "\n",
    "A good example is the probit model, where $\\varepsilon$ is a Gaussian vector. For this model, there is no close-form solution neither for $G$ nor for $G^*$.\n",
    "\n",
    "More recently, a number of modern models don't have closed-form either. These models require simulation methods in order to approximate them by discrete models."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The pure characteristics model\n",
    "\n",
    "#### Motivation\n",
    "\n",
    "The pure characteristics model (Berry and Pakes, 2007) can be motivated as follows. Assume $y$ stands for the number of bedrooms. The logit model would assume that the random utility associated with a 2-BR is uncorrelated with a 3-BR, which is not realistic.\n",
    "\n",
    "Let $\\xi_{y}$ is the typical size of a bedroom of size $y$, one may introduce $\\epsilon$ as the valuation of size; in which case the utility shock associated with $y$ should be $\\varepsilon_{y}=\\epsilon\\xi_{y}$. More generally, the characteristics $\\xi_{y}$ is a $d$-dimensional (deterministic) vector, and $\\epsilon\\sim\\mathbf{P}_{\\epsilon}$ is a (random) vector of the same size standing for the valuations of the respective dimensions, so that\n",
    "\n",
    "\\begin{align*}\n",
    "\\varepsilon_{y}=\\epsilon^{\\intercal}\\xi_{y}.\n",
    "\\end{align*}\n",
    "\n",
    "For example, if each alternative $y$ stands for a model of car, the first component of $\\xi_{y}$ may be the price of car $y$; the other components may be other characteristics such as number of seats, fuel efficiency, size, etc. In that case, for a given dimension $y\\in\\mathcal{Y}_{0}$, $\\epsilon_{y}$ is the (random) valuation of this dimension by the consumer with taste vector $\\epsilon$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Definition\n",
    "\n",
    "Assume without loss of generality that $\\varepsilon_{y}=0$, that is $\\xi_{0}=0$ as we can always reduce the setting to this case by replacing $\\xi_{y}$ by $\\xi_{y}-\\xi_{0}$.\n",
    "\n",
    "Letting $Z$ be the $\\left\\vert \\mathcal{Y}\\right\\vert \\times d\\,$\\ matrix of $\\left(  y,k\\right)  $-term $\\xi_{y}^{k}$, this rewrites as\n",
    "\n",
    "\\begin{align*}\n",
    "\\varepsilon = Z\\epsilon.\n",
    "\\end{align*}\n",
    "\n",
    "Hence, we have\n",
    "\n",
    "\\begin{align*}\n",
    "G\\left(  U\\right)  =\\mathbb{E}\\left[  \\max\\left\\{  U+Z\\epsilon,0\\right\\}\\right].\n",
    "\\end{align*}\n",
    "\n",
    "and\n",
    "\n",
    "\\begin{align*}\n",
    "\\sigma_{y}\\left(  U\\right)  =\\Pr\\left(  U_{y}-U_{z}\\geq\\left(  Z\\epsilon\\right)_{y}-\\left(  Z\\epsilon\\right)_{z}, \\quad forall z\\in\\mathcal{Y}_{0}\\backslash\\left\\{  y\\right\\}  \\right).\n",
    "\\end{align*}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### In dimension 1\n",
    "\n",
    "When $d=1$ (scalar characteristics), one has $\\sigma_{y}\\left(U\\right)  =\\Pr\\left(  U_{y}-U_{z}\\geq\\left(  \\xi_{y}-\\xi_{z}\\right)\\epsilon~\\forall z\\in\\mathcal{Y}_{0}\\backslash\\left\\{  y\\right\\}  \\right)  $, and thus\n",
    "\n",
    "\\begin{align*}\n",
    "\\sigma_{y}\\left(  U\\right)  =\\Pr\\left(  \\max_{z:\\xi_{y}>\\xi_{z}}\\left\\{\\frac{U_{y}-U_{z}}{\\xi_{y}-\\xi_{z}}\\right\\}  \\leq\\epsilon\\leq\\min_{z:\\xi_{y}<\\xi_{z}}\\left\\{  \\frac{U_{y}-U_{z}}{\\xi_{y}-\\xi_{z}}\\right\\}  \\right)\n",
    "\\end{align*}\n",
    "\n",
    "with the understanding that $\\max_{z\\in\\emptyset}f_{z}=-\\infty$ and \n",
    "$\\min_{z\\in\\emptyset}f_{z}=+\\infty$.\n",
    "\n",
    "Therefore, letting $\\mathbf{F}_{\\epsilon}$ be the c.d.f. associated with the distribution of $\\epsilon$, one has a closed-form expression for $\\sigma_{y}$:\n",
    "\n",
    "\\begin{align*}\n",
    "\\sigma_{y}\\left(  U\\right)  =\\mathbf{F}_{\\epsilon}\\left(  \\left[  \\max_{z:\\xi_{y}>\\xi_{z}}\\left\\{  \\frac{U_{y}-U_{z}}{\\xi_{y}-\\xi_{z}}\\right\\},\\min_{z:\\xi_{y}<\\xi_{z}}\\left\\{  \\frac{U_{y}-U_{z}}{\\xi_{y}-\\xi_{z}}\\right\\}\\right]  \\right)\n",
    "\\end{align*}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The probit model\n",
    "\n",
    "When $\\mathbf{P}_{\\epsilon}$ is the $\\mathcal{N}\\left(  0,S\\right)  $\n",
    "distribution, then the pure characteristics model is called a Probit model; in this case,\n",
    "\n",
    "\\begin{align*}\n",
    "\\varepsilon\\sim\\mathcal{N}\\left(  0,\\Sigma\\right)  \\text{ where }%\n",
    "\\Sigma=ZSZ^{\\intercal}.\n",
    "\\end{align*}\n",
    "\n",
    "Note the distribution $\\varepsilon\\,$will not have full support unless $d\\geq\\left\\vert \\mathcal{Y}\\right\\vert $ and $Z$ is of full rank.\n",
    "\n",
    "Computing $\\sigma$ in the Probit model thus implies computing the mass assigned by the Gaussian distribution to rectangles of the type \n",
    "\n",
    "\\begin{align*}\n",
    "\\left[  l_{y},u_{y}\\right]  .\n",
    "\\end{align*}\n",
    "\n",
    "When $\\Sigma$ is diagonal (random utility terms are i.i.d. across alternatives), this is numerically easy. However, this is computationally difficult in general (more on this later)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The random coefficient logit model\n",
    "\n",
    "The random coefficient logit model (Berry, Levinsohn and Pakes, 1995) may be viewed as an interpolant between the random characteristics model and the logit model. In this case,\n",
    "\n",
    "\\begin{align*}\n",
    "\\varepsilon=\\left(  1-\\lambda\\right)  Z\\epsilon+\\lambda\\eta\n",
    "\\end{align*}\n",
    "\n",
    "where $\\epsilon\\sim\\mathbf{P}_{\\epsilon}$, $\\eta$ is an EV1 distribution independent from the previous term, and $\\lambda$ is a interpolation parameter ($\\lambda=1$ is the logit model, and $\\lambda=0$ is the pure characteristics model).\n",
    "\n",
    "In this case, one may compute the Emax operator as \n",
    "\n",
    "\\begin{align*}\n",
    "G\\left(  U\\right)   &  =\\mathbb{E}\\left[  \\max_{y\\in\\mathcal{Y}_{0}}\\left\\{U_{y}+\\left(  1-\\lambda\\right)  \\left(  Z\\epsilon\\right)  _{y}+\\lambda\\eta_{y}\\right\\}  \\right] \\\\\n",
    "& =\\mathbb{E}\\left[  \\mathbb{E}\\left[  \\max_{y\\in\\mathcal{Y}_{0}}\\left\\{U_{y}+\\left(  1-\\lambda\\right)  \\left(  Z\\epsilon\\right)  _{y}+\\lambda\\eta_{y}\\right\\}  |\\epsilon\\right]  \\right] \\\\\n",
    "&  =\\mathbb{E}\\left[  \\lambda\\log\\sum_{y\\in\\mathcal{Y}_{0}}\\exp\\left(\n",
    "\\frac{U_{y}+\\left(  1-\\lambda\\right)  \\left(  Z\\epsilon\\right)  _{y}}{\\lambda}\\right)  \\right]\n",
    "\\end{align*}\n",
    "\n",
    "Recall\n",
    "\n",
    "\\begin{align*}\n",
    "G\\left(  U\\right)  =\\mathbb{E}\\left[  \\lambda\\log\\sum_{y\\in\\mathcal{Y}_{0}}\\exp\\left(  \\frac{U_{y}+\\left(  1-\\lambda\\right)  \\left(  Z\\epsilon\\right){y}}{\\lambda}\\right)  \\right]  .\n",
    "\\end{align*}\n",
    "\n",
    "The demand map in the random coefficients logit model is obtained by derivation of the expression of the Emax, i.e.\n",
    "\n",
    "\\begin{align*}\n",
    "\\sigma_{y}\\left(  U\\right)  =\\mathbb{E}\\left[  \\frac{\\exp\\left(  \\frac{U_{y}+\\left(  1-\\lambda\\right)  \\left(  Z\\epsilon\\right)  _{y}}{\\lambda}\\right)  }{\\sum_{y^{\\prime}\\in\\mathcal{Y}_{0}}\\exp\\left(  \\frac{U_{y^{\\prime}}+\\left(  1-\\lambda\\right)  \\left(  Z\\epsilon\\right)  _{y^{\\prime}}}{\\lambda}\\right)  }\\right]  .\n",
    "\\end{align*}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Simulation methods\n",
    "\n",
    "In a number of cases, one cannot compute the choice probabilities $\\sigma\\left(  U\\right)$ using a closed-form expression. In this case, we need to resort to simulation to compute $G$, $G^{\\ast}$, $\\sigma$ and $\\sigma^{-1}$.\n",
    "\n",
    "The idea is that:\n",
    "\n",
    "* One is able to compute $G$ and $G^{\\ast}$ for discrete distributions (more on this later)\n",
    "\n",
    "* The sampled versions of $G$, $G^{\\ast}$, $\\sigma$ and $\\sigma^{-1}$ converge to the populations objects when the sample size is large."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Accept-reject simulator\n",
    "\n",
    "One simulates $N$ points $\\varepsilon^{i}\\sim P$. The Emax operator associated with the empirical sample distribution $P_{N}$ is\n",
    "\n",
    "\\begin{align*}\n",
    "G_{N}=N^{-1}\\sum_{i=1}^{N}\\max_{y\\in\\mathcal{Y}}\\left\\{  U_{y} + \\varepsilon_{y}^{i}\\right\\}\n",
    "\\end{align*}\n",
    "\n",
    "and the demand map is given by\n",
    "\n",
    "\\begin{align*}\n",
    "\\sigma_{N,y}\\left(  U\\right)  =N^{-1}\\sum_{i=1}^{N}1\\left\\{  U_{y} + \\varepsilon_{y}^{i}\\geq U_{z}+\\varepsilon_{z}^{i}, \\quad \\forall z\\in\\mathcal{Y}\n",
    "_{0}\\right\\}\n",
    "\\end{align*}\n",
    "\n",
    "In the literature, $\\sigma_{N}$ is called the *accept-reject simulator*."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Example**. We shall code the AR simulator for the probit model and generate choice probabilities. Take a vector of systematic utilities:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "seed = 777\n",
    "nbDraws = 1000\n",
    "U_y = np.array([0.4, 0.5, 0.2, 0.3, 0.1, 0])\n",
    "nbY = len(U_y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We shall specify a probit model where alternatives have correlation $\\rho = .5$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "rho = 0.5\n",
    "Covar = rho * np.ones((nbY, nbY)) + (1 - rho) * np.eye(nbY)\n",
    "\n",
    "E = np.linalg.eigh(Covar)\n",
    "V = E[0]\n",
    "Q = E[1]\n",
    "SqrtCovar = Q.dot(np.diag(np.sqrt(V))).dot(Q.T)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now generate the $\\varepsilon_{iy}$'s:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "epsilon_iy = np.random.normal(0,1,nbDraws*nbY).reshape(nbDraws,nbY).dot(SqrtCovar)\n",
    "u_iy = epsilon_iy + U_y\n",
    "\n",
    "ui = np.max(u_iy, axis=1)\n",
    "s_y = np.sum((u_iy.T - ui).T == 0, axis=0) / nbDraws"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### McFadden's SARS\n",
    "\n",
    "McFadden's smoothed accept-reject simulator (SARS) consists in sampling $\\varepsilon\\sim P$: $\\varepsilon^{1},...,\\varepsilon^{N}$, and replacing the max by the smooth-max\n",
    "\n",
    "\\begin{align*}\n",
    "\\sigma_{N,T,y}\\left(  U\\right)  =\\sum_{i=1}^{N}\\frac{1}{N}\\frac{\\exp\\left((U_{y}+\\varepsilon_{y}^{i})/T\\right)  }{\\sum_{z}\\exp\\left(  (U_{z}+\\varepsilon_{z}^{i})/T\\right)}\n",
    "\\end{align*}\n",
    "\n",
    "One seeks $U$ so that the induced choice probabilities are $s$, that is\n",
    "\n",
    "\\begin{align*}\n",
    "s_{y}=\\sum_{i=1}^{N}\\frac{1}{N}\\frac{\\exp\\left(  (U_{y}+\\varepsilon_{y}^{i})/T\\right)  }{\\sum_{z}\\exp\\left(  (U_{z}+\\varepsilon_{z}^{i})/T\\right)}.\n",
    "\\end{align*}\n",
    "\n",
    "The associated Emax operator is\n",
    "\n",
    "\\begin{align*}\n",
    "G_{N,T}\\left(  U\\right)  =\\mathbb{E}_{\\mathbf{P}_{N}}\\left[G_{\\operatorname{logit}}\\left(  U+\\varepsilon^{i}\\right)  \\right]\n",
    "\\end{align*}\n",
    "\n",
    "so the underlying random utility structure is a random coefficient logit."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# The inversion theorem"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following theorem, first shown in Galichon and Salanié (2011), shows that the inversion of discrete choice models is an optimal transport problem.\n",
    "\n",
    "---\n",
    "**Theorem** [Galichon and Salanie]\n",
    "\n",
    "Consider a solution $\\left(u\\left(  \\varepsilon\\right),v_{y}\\right)  $ to the dual Monge-Kantorovich problem with cost $\\Phi\\left(\\varepsilon,y\\right)  =\\varepsilon_{y}$, that is:\n",
    "\n",
    "<a name='MKDualDC'></a>\n",
    "\\begin{align*}\n",
    "\\min_{u,v}  &  \\int u\\left(  \\varepsilon\\right)  d\\mathbf{P}\\left(\n",
    "\\varepsilon\\right)  +\\sum_{y\\in\\mathcal{Y}_{0}}v_{y}s_{y}\\\\\n",
    "s.t.~  &  u\\left(  \\varepsilon\\right)  +v_{y}\\geq\\Phi\\left(  \\varepsilon\n",
    ",y\\right)\n",
    "\\end{align*}\n",
    "\n",
    "Then:\n",
    "\n",
    "* $U=\\sigma^{-1}\\left(  s\\right)  $ is given by $U_{y}=v_{0}-v_{y}$.\n",
    "\n",
    "* The value of the [MK dual](#MKDualDC) is $-G^{\\ast}\\left(  s\\right)  $.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Proof**\n",
    "\n",
    "$\\sigma^{-1}\\left(  s\\right)  =\\arg\\max_{U:U_{0}=0}\\left\\{  \\sum_{y\\in\\mathcal{Y}}s_{y}U_{y}-G(U)\\right\\}  $, thus, letting $v=-U$, $v$ is the solution to\n",
    "\n",
    "\\begin{align*}\n",
    "\\min_{v:v_{0}=0}\\left\\{  \\sum_{y\\in\\mathcal{Y}_{0}}s_{y}v_{y}+G(-v)\\right\\}\n",
    "\\end{align*}\n",
    "which is exactly the [MK dual](#MKDualDC).\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Inversion of the pure characteristics model\n",
    "\n",
    "It follows from the inversion theorem that the problem of demand inversion in the pure characteristics model is a semi-discrete transport problem, a point made in Bonnet, Galichon, O'Hara, and Shum (2017).\n",
    "\n",
    "Indeed, the correspondence is:\n",
    "\n",
    "* an alternative $y$ is a fountain\n",
    "\n",
    "* the characteristics of an alternative is a fountain location\n",
    "\n",
    "* the systematic utility associated with alternative $y$ is minus the price of fountain $y$\n",
    "\n",
    "* the market share of altenative $y$ coindides with the capacity of fountain $y$\n",
    "\n",
    "* the random vector $\\epsilon$ is the location of an inhabitant"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Example**. As an example, we invert the probit model above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Using license file /opt/gurobi/gurobi.lic\n",
      "Set parameter GURO_PAR_WLSACCESSID\n",
      "Set parameter GURO_PAR_WLSSECRET\n",
      "Set parameter GURO_PAR_WLSTOKEN\n",
      "Set parameter GURO_PAR_LICENSEID to value 557591\n",
      "Gurobi Optimizer version 9.1.1 build v9.1.1rc0 (linux64)\n",
      "Thread count: 2 physical cores, 4 logical processors, using up to 4 threads\n",
      "Optimize a model with 1006 rows, 6000 columns and 12000 nonzeros\n",
      "Model fingerprint: 0x877b70a7\n",
      "Coefficient statistics:\n",
      "  Matrix range     [1e+00, 1e+00]\n",
      "  Objective range  [4e-05, 4e+00]\n",
      "  Bounds range     [0e+00, 0e+00]\n",
      "  RHS range        [1e-03, 3e-01]\n",
      "Presolve time: 0.02s\n",
      "Presolved: 1006 rows, 6000 columns, 12000 nonzeros\n",
      "\n",
      "Iteration    Objective       Primal Inf.    Dual Inf.      Time\n",
      "       0    2.4093616e+33   6.004000e+33   1.927489e+04      0s\n",
      "    1251    8.6881256e-01   0.000000e+00   0.000000e+00      0s\n",
      "\n",
      "Solved in 1251 iterations and 0.06 seconds\n",
      "Optimal objective  8.688125603e-01\n",
      "U_y (true and recovered)\n",
      "[0.4 0.5 0.2 0.3 0.1 0. ]\n",
      "[ 0.39468612  0.49545432  0.17618765  0.29614658  0.08972953 -0.        ]\n"
     ]
    }
   ],
   "source": [
    "A1 = spr.kron(np.ones((1, nbY)), spr.identity(nbDraws))\n",
    "A2 = spr.kron(spr.identity(nbY), np.ones((1, nbDraws)))\n",
    "A = spr.vstack([A1, A2])\n",
    "obj = epsilon_iy.flatten(order='F') # change this\n",
    "rhs = np.ones(nbDraws)/nbDraws\n",
    "rhs = np.append(rhs, s_y)\n",
    "m = grb.Model('optimal')\n",
    "x = m.addMVar(len(obj), name='couple')\n",
    "m.setObjective(obj @ x, grb.GRB.MAXIMIZE)\n",
    "m.addConstr(A @ x == rhs, name=\"Constr\")\n",
    "\n",
    "m.optimize()\n",
    "if m.Status == grb.GRB.Status.OPTIMAL:\n",
    "    pi = m.getAttr('pi')\n",
    "    Uhat_y = -np.subtract(pi[nbDraws:nbY+nbDraws], pi[nbY + nbDraws - 1])\n",
    "    print('U_y (true and recovered)')\n",
    "    print(U_y)\n",
    "    print(Uhat_y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### McFadden's SARS and regularized Optimal Transport\n",
    "\n",
    "Cf. Bonnet, Galichon, O'Hara and Shum (2017). Let $u_{i}=T\\log\\sum_{z}\\exp\\left((U_{z}+\\varepsilon_{z}^{i})/T\\right)  $. One has\n",
    "\n",
    "\\begin{align*}\n",
    "\\left\\{\n",
    "\\begin{array}\n",
    "[c]{l}%\n",
    "s_{y}=\\sum_{i=1}^{N}\\frac{1}{N}\\exp\\left(  (U_{y}-u_{i}+\\varepsilon_{y}^{i})/T\\right) \\\\\n",
    "\\frac{1}{N}=\\sum_{y}\\frac{1}{N}\\exp\\left(  (U_{y}-u_{i}+\\varepsilon_{y}^{i})/T\\right)\n",
    "\\end{array}\n",
    "\\right. .\n",
    "\\end{align*}\n",
    "\n",
    "As a result, $\\left(  u_{i},U_{y}\\right)  $ are the solution of the regularized OT problem\n",
    "\n",
    "\\begin{align*}\n",
    "\\min_{u,U}\\sum_{i=1}^{N}\\frac{1}{N}u_{i}-\\sum s_{y}U_{y}+\\sum_{i,y}\\frac{1}{N}\\exp\\left(  (U_{y}-u_{i}+\\varepsilon_{y}^{i})/T\\right)  .\n",
    "\\end{align*}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### BLP's contraction mapping\n",
    "\n",
    "Consider the IPFP algorithm for solving the latter problem:\n",
    "\n",
    "\\begin{align*}\n",
    "\\left\\{\n",
    "\\begin{array}\n",
    "[c]{l}\n",
    "\\exp\\left(  u_{i}^{k+1}/T\\right)  =\\sum_{z}\\exp\\left(  (U_{z}^{k}%\n",
    "+\\varepsilon_{z}^{i})/T\\right) \\\\\n",
    "\\exp U_{y}^{k+1}/T=\\frac{Ns_{y}}{\\sum_{i=1}^{N}\\exp\\left(  (-u_{i}%\n",
    "^{k+1}+\\varepsilon_{y}^{i})/T\\right)  }\n",
    "\\end{array}\n",
    "\\right.\n",
    "\\end{align*}\n",
    "\n",
    "This rewrites as\n",
    "\n",
    "\\begin{align*}\n",
    "\\exp U_{y}^{k+1}/T  &  =\\frac{Ns_{y}}{\\sum_{i=1}^{N}\\frac{\\exp\\left(\n",
    "\\varepsilon_{y}^{i}/T\\right)  }{\\sum_{z}\\exp\\left(  (U_{z}^{k}+\\varepsilon\n",
    "_{z}^{i})/T\\right)  }},\\text{ i.e.}\\\\\n",
    "U_{y}^{k+1}  &  =T\\log s_{y}-T\\log\\sum_{i=1}^{N}\\frac{1}{N}\\frac{\\exp\\left(\n",
    "\\varepsilon_{y}^{i}/T\\right)  }{\\sum_{z}\\exp\\left(  (U_{z}^{k}+\\varepsilon\n",
    "_{z}^{i})/T\\right)  }\n",
    "\\end{align*}\n",
    "\n",
    "which is exactly the contraction mapping algorithm of Berry, Levinsohn and Pakes (1995, appendix 1)."
   ]
  }
 ],
 "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
}