{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Block 8: Discrete choice estimation\n",
"### Alfred Galichon (NYU & Sciences Po)\n",
"## 'math+econ+code' masterclass on optimal transport and economic applications\n",
"© 2018-2021 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.\n",
"\n",
"#### With python code examples"
]
},
{
"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": 2,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import os\n",
"import pandas as pd\n",
"import string as str\n",
"import math\n",
"import sys\n",
"import time\n",
"import scipy.sparse as spr\n",
"from scipy import optimize, special\n",
"import gurobipy as grb"
]
},
{
"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 represent it conveniently in a similar fashion as in block 6:"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"thepath = 'https://raw.githubusercontent.com/math-econ-code/mec_optim_2021-01/master/data_mec_optim/demand_travelmode/'\n",
"\n",
"travelmode = pd.read_csv(thepath+'travelmodedata.csv')\n",
"\n",
"travelmode['choice'] = np.where(travelmode['choice'] =='yes' , 1, 0)\n",
"\n",
"nobs = travelmode.shape[0]\n",
"ncols = travelmode.shape[1]\n",
"nbchoices = 4\n",
"ninds = int(nobs/nbchoices)\n",
"\n",
"muhat_i_y = travelmode['choice'].to_numpy().reshape(ninds,nbchoices).T\n",
"muhat_iy = muhat_i_y.flatten()\n",
"\n",
"muhat_i_y = travelmode['choice'].to_numpy().reshape(ninds,4).T\n",
"muhat_iy = muhat_i_y.flatten()\n",
"\n",
"s_y = travelmode.groupby(['mode']).mean()['choice'].to_frame().sort_index()\n",
"\n",
"def two_d(X):\n",
" return np.reshape(X,(X.size, 1))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Estimation with no observable heterogeneity\n",
"\n",
"Start with assuming that there is no observable heterogeneity, so the only observation we have at hand are the aggregate market shares $s_y$. Hence the systematic utility will be the same for every agent. However, we wish to write a parametric model for it, namely assume a knwon parametric form for the dependence of $U_y$ with respect to various observed characteristics associated with $y$.\n",
"\n",
"Assume then that the utilities are parameterized as follows: $U = \\Phi \\beta$ where $\\beta\\in\\mathbb{R}^{p}$ is a parameter, and $\\Phi$ is a $\\left\\vert \\mathcal{Y}\\right\\vert \\times p$ matrix.\n",
"\n",
"The log-likelihood function is given by\n",
"\n",
"\\begin{align*}\n",
"l\\left( \\beta\\right) =N\\sum_{y}\\hat{s}_{y}\\log\\sigma_{y}\\left(\\Phi \\beta\\right)\n",
"\\end{align*}\n",
"\n",
"A common estimation method of $\\beta$ is by maximum likelihood%\n",
"\n",
"\\begin{align*}\n",
"\\max_{\\beta}l\\left( \\beta\\right) .\n",
"\\end{align*}\n",
"\n",
"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)."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### MLE, logit case\n",
"\n",
"In the logit case,\n",
"\n",
"\\begin{align*}\n",
"l\\left( \\beta\\right) =N\\left\\{ \\hat{s}^{\\intercal}\\Phi\\beta-\\log\\sum_{y}\\exp\\left( \\Phi\\beta\\right) _{y}\\right\\}\n",
"\\end{align*}\n",
"\n",
"so that the max-likehood amounts to\n",
"\n",
"\\begin{align*}\n",
"\\max_{\\beta}\\left\\{ \\hat{s}^{\\intercal} \\Phi \\beta-G\\left( \\Phi \\beta\\right)\n",
"_{y}\\right\\}\n",
"\\end{align*}\n",
"\n",
"whose value is the Legendre-Fenchel transform of $\\beta\\rightarrow G\\left( \\Phi \\beta\\right)$ evaluated at $\\Phi ^{^{\\intercal}}\\hat{s}$.\n",
"\n",
"Note that the vector $\\Phi^{^{\\intercal}}\\hat{s}$ 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."
]
},
{
"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{s}$.\n",
"\n",
"The social welfare is given in general by $W\\left( \\beta\\right) =G\\left( \\Phi\\beta\\right) $. One has $\\partial_{\\beta^{i}}W\\left(\\beta\\right) =\\sum_{y}\\sigma_{y}\\left( \\Phi\\beta\\right) \\Phi_{yi}$, that is \n",
"\n",
"\\begin{align*}\n",
"\\nabla W\\left( \\beta\\right) = \\Phi^{\\intercal}\\sigma\\left( \\Phi\\beta\\right) ,\n",
"\\end{align*}\n",
"\n",
"which is the vector of predicted moments.\n",
"\n",
"Therefore the program\n",
"\n",
"\\begin{align*}\n",
"\\max_{\\beta}\\left\\{ \\hat{s}^{\\intercal}\\Phi\\beta-G\\left( \\Phi\\beta\\right)\n",
"_{y}\\right\\}\n",
"\\end{align*}\n",
"\n",
"picks up the parameter $\\beta$ which matches the empirical moments $X^{^{\\intercal}}\\hat{s}$ with the predicted ones $\\nabla W\\left(\\beta\\right) $. This procedure is not statistically efficient, but is computationally efficient becauses it arises from a convex optimization problem."
]
},
{
"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( \\beta\\right) =N\\left\\{ \\hat{s}^{\\intercal}\\Phi\\beta-\\log\\sum_{y} \\exp\\left( \\Phi\\beta\\right) _{y}\\right\\}\n",
"\\end{align*}\n",
"\n",
"Assume that we restrict ourselves to $\\beta^{\\top}z>0$. Then we can write $\\beta=\\theta/T$ where $T=1/\\beta^{\\top}z$ and $\\theta=\\beta T$. Call $\\Theta=\\left\\{ \\theta\\in\\mathbb{R}^{p},\\theta^{\\top}z=1\\right\\} $, so that $\\beta=\\theta/T$ where $\\theta\\in\\Theta$ and $T>0$. We have\n",
"\n",
"\\begin{align*}\n",
"l\\left( \\theta,T\\right) =\\frac{N}{T}\\left\\{ \\hat{s}^{\\intercal}\n",
"\\Phi\\theta-T\\log\\sum_{y}\\exp\\left( \\frac{\\left( \\Phi\\theta\\right) _{y}}{T}\\right) \\right\\}\n",
"\\end{align*}\n",
"\n",
"and we define the *fixed temperature maximum likelihood estimator* by\n",
"\n",
"\\begin{align*}\n",
"\\theta\\left( T\\right) =\\arg\\max_{\\theta}l\\left( \\theta,T\\right)\n",
"\\end{align*}\n",
"\n",
" Note that $\\theta\\left( T\\right) =\\arg\\max_{\\theta\\in\\Theta}Tl\\left(\\theta,T\\right)$ where\n",
"\n",
"\\begin{align*}\n",
"Tl\\left( \\theta,T\\right) =N\\left\\{ \\hat{s}^{\\intercal}\\Phi\\theta-T\\log\\sum _{y}\\exp\\left( \\frac{\\left( \\Phi\\theta\\right) _{y}}{T}\\right) \\right\\}\n",
"\\end{align*}\n",
"\n",
"and we note that $Tl\\left( \\theta,T\\right) \\rightarrow N\\left\\{ \\hat{s}^{\\intercal}\\Phi\\theta-\\max_{y\\in\\mathcal{Y}}\\left\\{ \\left( \\Phi\\theta\\right)_{y}\\right\\} \\right\\} $ as $T\\rightarrow0$.\n",
"\n",
"We have\n",
"\n",
"\\begin{align*}\n",
"\\frac{Tl\\left( \\theta,T\\right) }{N}=\\hat{s}^{\\intercal}\\Phi\\theta-T\\log\\sum_{y}\\exp\\left( \\frac{\\left( \\Phi\\theta\\right) _{y}}{T}\\right)\n",
"\\end{align*}\n",
"\n",
"Let $\\theta\\left( 0\\right) =\\lim_{T\\rightarrow0}\\theta\\left(T\\right) $. Calling $m\\left( \\theta\\right) =\\max_{y\\in\\mathcal{Y}}\\left\\{\\left( \\Phi\\theta\\right) _{y}\\right\\} $, we have\n",
"\n",
"\\begin{align*}\n",
"\\theta\\left( 0\\right) \\in\\arg\\max_{\\theta}\\left\\{ \\hat{s}^{\\intercal}\\Phi\\theta-m\\left( \\theta\\right) \\right\\},\n",
"\\end{align*}\n",
"\n",
"or\n",
"\n",
"\\begin{align*}\n",
"\\theta\\left( 0\\right) \\in\\arg\\min_{\\theta}\\left\\{ m\\left( \\theta\\right)-\\hat{s}^{\\intercal}\\Phi\\theta\\right\\},\n",
"\\end{align*}\n",
"\n",
"Calling $m\\left( \\theta\\right) =\\max_{y\\in\\mathcal{Y}}\\left\\{ \\left(\\Phi\\theta\\right) _{y}\\right\\} $, one has \n",
"\n",
"\\begin{align*}\n",
"\\theta\\left( T\\right) \\in\\arg\\max\\left\\{ \\hat{s}^{\\intercal}\\Phi\\theta-m\\left( \\theta\\right) -T\\log\\sum_{y}\\exp\\left( \\frac{\\left(\\Phi\\theta\\right) _{y}-m\\left( \\theta\\right) }{T}\\right) \\right\\}\n",
"\\end{align*}\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Minimax-regret estimation\n",
"\n",
"Note that\n",
"\n",
"\\begin{align*}\n",
"\\theta\\left( 0\\right) \\in\\arg\\max\\left\\{ \\hat{s}^{\\intercal}\\Phi\\theta\n",
"-m\\left( \\theta\\right) \\right\\} .\n",
"\\end{align*}\n",
"\n",
"Define $R_{i}\\left( \\theta,y\\right) =\\left( \\Phi\\theta\\right)_{y}-\\left( \\Phi\\theta\\right) _{y_{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( \\theta,y\\right) =\\max_{y\\in\\mathcal{Y}}\\left\\{ \\left( \\Phi\\theta\\right)_{y}-\\left( \\Phi\\theta\\right)_{y_{i}}\\right\\}\n",
"\\end{align*}\n",
"\n",
"and the max-regret associated with the sample is $\\frac{1}{N}\\sum\\max_{y\\in\\mathcal{Y}}\\left\\{ R_{i}\\left( \\theta,y\\right) \\right\\} $, that is $\\max_{y\\in\\mathcal{Y}}\\left\\{ \\left( \\Phi\\theta\\right) _{y}\\right\\} - \\hat{s}^{\\intercal}X\\theta$.\n",
"\n",
"The minimax regret estimator\n",
"\n",
"\\begin{align*}\n",
"\\hat{\\theta}^{MMR}=\\min_{\\theta}\\left\\{ m\\left( \\theta\\right) -\\hat\n",
"{s}^{\\intercal}\\Phi\\theta\\right\\}\n",
"\\end{align*}\n",
"\n",
"which has a linear programming fomulation\n",
"\n",
"\\begin{align*}\n",
"& \\min_{m,\\theta}m-\\hat{s}^{\\intercal}\\Phi\\theta\\\\\n",
"s.t.~ & m-\\left( \\Phi\\theta\\right) _{y}\\geq\\forall y\\in\\mathcal{Y}\n",
"\\end{align*}"
]
},
{
"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_{\\theta,m}/\\max_{\\theta,m} \\theta^{\\intercal}d\\\\\n",
"s.t.~ & m-\\hat{s}^{\\intercal}X\\theta=V\\\\\n",
"& m\\geq\\left( \\Phi\\theta\\right)_{y}, \\quad \\forall y\\in\\mathcal{Y}%\n",
"\\end{align*}"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Link with exponential families and GLM\n",
"\n",
"See class notes"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Estimation with observed heterogeneity\n",
"\n",
"We now 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\n",
"$$U = \\Phi \\beta,$$ 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",
"See class notes."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Application\n",
"\n",
"Back to the dataset:"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"Phi_iy_k = np.column_stack((np.kron(np.identity(4)[0:4,1:4],np.repeat(1, ninds).reshape(ninds,1)), - travelmode['travel'].to_numpy(), - (travelmode['travel']*travelmode['income']).to_numpy(), - travelmode['gcost'].to_numpy()))"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"nbK = Phi_iy_k.shape[1]\n",
"phi_mean = Phi_iy_k.mean(axis = 0)\n",
"phi_stdev = Phi_iy_k.std(axis = 0, ddof = 1)\n",
"Phi_iy_k = ((Phi_iy_k - phi_mean).T/phi_stdev[:,None]).T"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"def log_likelihood(theta):\n",
" nbK = np.asarray(theta).shape[0]\n",
" Xtheta = Phi_iy_k.dot(theta)/sigma\n",
" Xthetamat_iy = Xtheta.reshape(nbchoices, ninds).T\n",
" max_i = np.amax(Xthetamat_iy, axis = 1)\n",
" expPhi_iy = np.exp((Xthetamat_iy.T -max_i).T)\n",
" d_i = np.sum(expPhi_iy, axis = 1)\n",
" \n",
" val = np.sum(Xtheta * muhat_iy) - np.sum(max_i) - sigma * np.sum(np.log(d_i))\n",
"\n",
" return -val"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"def grad_log_likelihood(theta):\n",
" nbK = np.asarray(theta).shape[0]\n",
" Xtheta = Phi_iy_k.dot(theta)/sigma\n",
" Xthetamat_iy = Xtheta.reshape(nbchoices, ninds).T\n",
" max_i = np.amax(Xthetamat_iy, axis = 1)\n",
" expPhi_iy = np.exp((Xthetamat_iy.T -max_i).T)\n",
" d_i = np.sum(expPhi_iy, axis = 1)\n",
" \n",
" temp_mat = (Phi_iy_k.T * expPhi_iy.T.flatten()).T\n",
" list_temp = []\n",
" for i in range(nbchoices):\n",
" list_temp.append(temp_mat[i*ninds:(i+1)*ninds,])\n",
" n_i_k = np.sum(list_temp,axis = 0)\n",
" \n",
" thegrad = muhat_iy.reshape(1,nbchoices*ninds).dot(Phi_iy_k).flatten() - np.sum(n_i_k.T/d_i, axis = 1)\n",
"\n",
" return -thegrad"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [],
"source": [
"theta0 = np.repeat(0,nbK)\n",
"sigma = 1\n",
"outcome = optimize.minimize(log_likelihood,method = 'CG',jac = grad_log_likelihood, x0 = theta0)"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
" fun: 280.9103677250631\n",
" jac: array([-1.56394686e-06, -3.55409114e-07, -1.80343209e-06, 1.80183436e-07,\n",
" 4.66669048e-07, 2.80306917e-07])\n",
" message: 'Optimization terminated successfully.'\n",
" nfev: 31\n",
" nit: 16\n",
" njev: 31\n",
" status: 0\n",
" success: True\n",
" x: array([-0.02769993, -0.35847009, 0.01487751, 0.01416175, 0.11784738,\n",
" -0.25344193])"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"outcome"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"-3.9456769581766116\n",
"[ 0.10929498 1.41440717 -0.05870185 -0.05587768 -0.46498768 1. ]\n"
]
}
],
"source": [
"temp_mle = 1 / outcome['x'][nbK - 1]\n",
"theta_mle = outcome['x']*temp_mle\n",
"print(temp_mle)\n",
"print(theta_mle)"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Academic license - for non-commercial use only - expires 2022-11-14\n",
"Using license file c:\\gurobi911\\gurobi.lic\n",
"Gurobi Optimizer version 9.1.1 build v9.1.1rc0 (win64)\n",
"Thread count: 6 physical cores, 12 logical processors, using up to 12 threads\n",
"Optimize a model with 841 rows, 216 columns and 5881 nonzeros\n",
"Model fingerprint: 0xece1c49f\n",
"Coefficient statistics:\n",
" Matrix range [3e-03, 5e+00]\n",
" Objective range [2e-01, 5e+01]\n",
" Bounds range [0e+00, 0e+00]\n",
" RHS range [1e+00, 1e+00]\n",
"Presolve removed 211 rows and 1 columns\n",
"Presolve time: 0.01s\n",
"Presolved: 630 rows, 215 columns, 2520 nonzeros\n",
"\n",
"Iteration Objective Primal Inf. Dual Inf. Time\n",
" 0 handle free variables 0s\n",
" 380 -1.3829488e+02 0.000000e+00 0.000000e+00 0s\n",
"\n",
"Solved in 380 iterations and 0.02 seconds\n",
"Optimal objective -1.382948769e+02\n",
"Value of the problem (Gurobi) = -138.29487687593996\n"
]
}
],
"source": [
"lenobj = nbK+ninds\n",
"c = np.concatenate((muhat_iy.reshape(1,nbchoices*ninds).dot(Phi_iy_k).flatten(),np.repeat(-1,ninds)))\n",
"\n",
"m = grb.Model('lp')\n",
"x = m.addMVar(lenobj, name='x', lb=-grb.GRB.INFINITY)\n",
"m.setObjective(c @ x, grb.GRB.MAXIMIZE)\n",
"cstMat = spr.hstack((spr.csr_matrix(-Phi_iy_k), spr.kron(two_d(np.ones(nbchoices)),spr.identity(ninds))))\n",
"rhs = np.zeros(ninds*nbchoices)\n",
"m.addConstr(cstMat @ x >= rhs)\n",
"nbCstr = cstMat.shape[0]\n",
"const_2 = np.array([0]*(nbK - 1))\n",
"const_2 = np.append(const_2, 1)\n",
"const_2 = np.append(const_2 ,[0]*ninds)\n",
"m.addConstr(const_2 @ x == 1)\n",
"m.optimize()\n",
"if m.status == grb.GRB.Status.OPTIMAL:\n",
" print(\"Value of the problem (Gurobi) =\", m.objval)\n",
" opt_x = m.getAttr('x')"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[ 0.13949426 0.05032141 0.03395604 -0.72749334 -0.02506117 1. ]\n",
"[ 0.10929498 1.41440717 -0.05870185 -0.05587768 -0.46498768 1. ]\n"
]
}
],
"source": [
"theta_lp = np.array(opt_x[:nbK])\n",
"print(theta_lp)\n",
"print(theta_mle)"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [],
"source": [
"indMax=100\n",
"tempMax=temp_mle\n",
"outcomemat = np.zeros((indMax+1,nbK-1))"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [],
"source": [
"def log_likelihood_fixedtemp(subsetoftheta, *temp):\n",
" val = log_likelihood(np.append(subsetoftheta, 1/temp[0]))\n",
" \n",
" return val"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [],
"source": [
"def grad_log_likelihood_fixedtemp(subsetoftheta, *temp):\n",
" val = grad_log_likelihood(np.append(subsetoftheta, 1/temp[0]))[:-1]\n",
" \n",
" return val"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [],
"source": [
"outcomemat[0,:] = theta_lp[:-1]\n",
"iterMax = indMax+1\n",
"for k in range(2,iterMax+1,1):\n",
" thetemp = tempMax * (k-1)/indMax\n",
" outcomeFixedTemp = optimize.minimize(log_likelihood_fixedtemp,method = 'CG',jac = grad_log_likelihood_fixedtemp, args = (thetemp,), x0 = theta0[:-1])\n",
" outcomemat[k-1,:] = outcomeFixedTemp['x']*thetemp"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([[ 0.13949426, 0.05032141, 0.03395604, -0.72749334, -0.02506117],\n",
" [ 0.14292345, 0.42333439, -0.06371649, -0.51242318, -0.1314772 ],\n",
" [ 0.14482454, 0.42218714, -0.06321331, -0.50753462, -0.12936482],\n",
" [ 0.14587841, 0.42033993, -0.06257653, -0.5042307 , -0.13058508],\n",
" [ 0.14624026, 0.41848991, -0.06210426, -0.50100175, -0.13399352],\n",
" [ 0.14642665, 0.41756089, -0.06139318, -0.49794246, -0.13827628],\n",
" [ 0.14665808, 0.41780041, -0.06035416, -0.49509681, -0.14277605],\n",
" [ 0.14695082, 0.41912092, -0.05907464, -0.4924095 , -0.14720936],\n",
" [ 0.1472558 , 0.42138082, -0.05767761, -0.48979281, -0.15145913],\n",
" [ 0.1475228 , 0.42447024, -0.05626612, -0.48716609, -0.15548361],\n",
" [ 0.14771955, 0.42831363, -0.05491026, -0.48446871, -0.15928106],\n",
" [ 0.147832 , 0.43285517, -0.05365016, -0.48166091, -0.16287184],\n",
" [ 0.14785899, 0.43804697, -0.05250383, -0.47872002, -0.16628727],\n",
" [ 0.14780664, 0.44384285, -0.05147502, -0.4756363 , -0.16956138],\n",
" [ 0.14768441, 0.4501967 , -0.05055895, -0.47240827, -0.1727272 ],\n",
" [ 0.14750252, 0.45706218, -0.049747 , -0.46904023, -0.17581347],\n",
" [ 0.14727078, 0.46439392, -0.04902871, -0.46553968, -0.17884437],\n",
" [ 0.14699769, 0.47214821, -0.04839391, -0.46191531, -0.18183989],\n",
" [ 0.14669068, 0.48028388, -0.04783308, -0.45817657, -0.18481554],\n",
" [ 0.14635581, 0.4887628 , -0.04733761, -0.454333 , -0.18778329],\n",
" [ 0.14599821, 0.49755006, -0.04690026, -0.45039344, -0.19075232],\n",
" [ 0.14562192, 0.50661408, -0.04651459, -0.44636624, -0.19372942],\n",
" [ 0.14523021, 0.51592631, -0.04617522, -0.44225927, -0.1967194 ],\n",
" [ 0.14482593, 0.52546154, -0.04587746, -0.4380792 , -0.19972597],\n",
" [ 0.14441128, 0.53519698, -0.04561735, -0.43383239, -0.20275147],\n",
" [ 0.14398796, 0.54511243, -0.04539156, -0.4295247 , -0.2057973 ],\n",
" [ 0.1435575 , 0.55519007, -0.0451971 , -0.4251611 , -0.20886437],\n",
" [ 0.1431211 , 0.56541397, -0.0450314 , -0.42074623, -0.21195309],\n",
" [ 0.14267966, 0.57577012, -0.0448923 , -0.41628419, -0.21506341],\n",
" [ 0.14223423, 0.58624609, -0.04477745, -0.41177899, -0.21819488],\n",
" [ 0.14178533, 0.59683058, -0.04468552, -0.40723351, -0.22134761],\n",
" [ 0.14133354, 0.60751392, -0.04461472, -0.40265106, -0.2245206 ],\n",
" [ 0.14087921, 0.61828703, -0.04456364, -0.3980343 , -0.22771335],\n",
" [ 0.14042282, 0.62914224, -0.0445311 , -0.3933857 , -0.23092522],\n",
" [ 0.13996489, 0.64007282, -0.04451558, -0.38870766, -0.23415508],\n",
" [ 0.13950548, 0.65107205, -0.04451636, -0.38400201, -0.23740251],\n",
" [ 0.13904476, 0.66213441, -0.0445324 , -0.37927088, -0.2406666 ],\n",
" [ 0.13858306, 0.67325491, -0.04456273, -0.37451579, -0.24394668],\n",
" [ 0.13812065, 0.68442903, -0.04460643, -0.36973808, -0.24724221],\n",
" [ 0.13765741, 0.69565258, -0.04466292, -0.36493967, -0.25055202],\n",
" [ 0.13719355, 0.70692191, -0.04473144, -0.36012167, -0.25387558],\n",
" [ 0.13672926, 0.71823353, -0.04481122, -0.35528524, -0.25721237],\n",
" [ 0.13626446, 0.7295844 , -0.04490188, -0.35043145, -0.26056178],\n",
" [ 0.13579942, 0.74097186, -0.04500259, -0.34556132, -0.26392307],\n",
" [ 0.13533403, 0.7523934 , -0.04511293, -0.34067596, -0.26729566],\n",
" [ 0.13486843, 0.76384653, -0.04523245, -0.33577616, -0.27067909],\n",
" [ 0.13440265, 0.77532938, -0.04536059, -0.3308628 , -0.27407258],\n",
" [ 0.13393689, 0.78684001, -0.04549687, -0.32593626, -0.27747614],\n",
" [ 0.13347077, 0.7983763 , -0.04564119, -0.32099779, -0.28088885],\n",
" [ 0.1330048 , 0.809937 , -0.04579276, -0.31604763, -0.28431068],\n",
" [ 0.13253869, 0.82152068, -0.04595141, -0.31108663, -0.28774075],\n",
" [ 0.13207258, 0.83312551, -0.04611692, -0.30611517, -0.2911789 ],\n",
" [ 0.13160643, 0.84475064, -0.04628883, -0.30113389, -0.29462481],\n",
" [ 0.13114011, 0.85639463, -0.04646697, -0.29614321, -0.29807827],\n",
" [ 0.13067406, 0.86805665, -0.04665084, -0.29114378, -0.30153838],\n",
" [ 0.13020788, 0.87973545, -0.04684038, -0.28613587, -0.30500532],\n",
" [ 0.12974168, 0.89143024, -0.04703527, -0.28112007, -0.30847847],\n",
" [ 0.12927578, 0.90314039, -0.04723521, -0.2760961 , -0.3119581 ],\n",
" [ 0.12880981, 0.91486466, -0.04744003, -0.27106525, -0.31544328],\n",
" [ 0.12834383, 0.92660242, -0.04764964, -0.26602702, -0.31893441],\n",
" [ 0.12787795, 0.93835309, -0.04786373, -0.26098232, -0.32243073],\n",
" [ 0.1274121 , 0.95011595, -0.04808207, -0.25593121, -0.3259322 ],\n",
" [ 0.12694638, 0.9618904 , -0.04830442, -0.25087406, -0.32943869],\n",
" [ 0.12648065, 0.97367593, -0.04853098, -0.24581086, -0.3329498 ],\n",
" [ 0.12601492, 0.98547194, -0.04876116, -0.2407422 , -0.33646546],\n",
" [ 0.12554925, 0.99727762, -0.04899523, -0.23566864, -0.33998514],\n",
" [ 0.12508389, 1.00909342, -0.04923243, -0.230589 , -0.34350961],\n",
" [ 0.12461853, 1.02091806, -0.0494732 , -0.22550493, -0.34703787],\n",
" [ 0.12415302, 1.03275119, -0.04971724, -0.22041607, -0.35057007],\n",
" [ 0.12368792, 1.04459299, -0.04996427, -0.21532268, -0.35410574],\n",
" [ 0.12322275, 1.05644253, -0.05021434, -0.2102249 , -0.35764511],\n",
" [ 0.12275764, 1.06829964, -0.05046731, -0.20512296, -0.3611879 ],\n",
" [ 0.1222923 , 1.08016371, -0.05072307, -0.20001745, -0.36473351],\n",
" [ 0.12182733, 1.092035 , -0.05098155, -0.19490767, -0.36828271],\n",
" [ 0.12136239, 1.10391293, -0.05124256, -0.18979413, -0.37183499],\n",
" [ 0.12089771, 1.11579742, -0.05150596, -0.18467661, -0.37539059],\n",
" [ 0.12043258, 1.12768748, -0.05177211, -0.1795563 , -0.37894824],\n",
" [ 0.11996796, 1.13958399, -0.05204018, -0.17443247, -0.38250907],\n",
" [ 0.11950335, 1.15148596, -0.05231066, -0.16930524, -0.38607249],\n",
" [ 0.11903875, 1.16339335, -0.05258328, -0.16417497, -0.3896385 ],\n",
" [ 0.1185742 , 1.17530596, -0.05285799, -0.15904171, -0.39320697],\n",
" [ 0.11810972, 1.18722358, -0.05313473, -0.15390554, -0.39677782],\n",
" [ 0.1176453 , 1.19914604, -0.05341341, -0.14876658, -0.40035098],\n",
" [ 0.11718094, 1.21107314, -0.05369399, -0.14362492, -0.40392636],\n",
" [ 0.11671664, 1.22300473, -0.05397639, -0.13848065, -0.40750389],\n",
" [ 0.1162524 , 1.23494063, -0.05426056, -0.13333385, -0.41108351],\n",
" [ 0.11578821, 1.2468807 , -0.05454644, -0.12818461, -0.41466514],\n",
" [ 0.11532409, 1.25882477, -0.05483398, -0.123033 , -0.41824872],\n",
" [ 0.11486001, 1.27077273, -0.05512312, -0.1178791 , -0.42183419],\n",
" [ 0.114396 , 1.28272443, -0.05541383, -0.11272299, -0.42542148],\n",
" [ 0.11393203, 1.29467973, -0.05570604, -0.10756472, -0.42901055],\n",
" [ 0.11346804, 1.30663857, -0.05599965, -0.10240454, -0.43260135],\n",
" [ 0.11300418, 1.31860074, -0.05629474, -0.0972422 , -0.4361938 ],\n",
" [ 0.11254036, 1.33056618, -0.05659122, -0.09207789, -0.43978786],\n",
" [ 0.11207659, 1.34253477, -0.05688903, -0.08691169, -0.44338348],\n",
" [ 0.11161287, 1.35450642, -0.05718813, -0.08174365, -0.44698063],\n",
" [ 0.11114919, 1.36648102, -0.0574885 , -0.07657383, -0.45057925],\n",
" [ 0.11068556, 1.37845848, -0.05779009, -0.07140228, -0.4541793 ],\n",
" [ 0.11022198, 1.39043872, -0.05809287, -0.06622904, -0.45778073],\n",
" [ 0.10975844, 1.40242164, -0.05839681, -0.06105416, -0.46138352],\n",
" [ 0.10929494, 1.41440717, -0.05870187, -0.0558777 , -0.46498762]])"
]
},
"execution_count": 17,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"outcomemat"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The zero-temperature estimator is:"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[ 0.14292345 0.42333439 -0.06371649 -0.51242318 -0.1314772 ]\n"
]
}
],
"source": [
"print(outcomemat[1,:])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The mle estimator is:"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[ 0.10929494 1.41440717 -0.05870187 -0.0558777 -0.46498762]\n"
]
}
],
"source": [
"print(outcomemat[indMax,])"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Gurobi Optimizer version 9.1.1 build v9.1.1rc0 (win64)\n",
"Thread count: 6 physical cores, 12 logical processors, using up to 12 threads\n",
"Optimize a model with 1680 rows, 426 columns and 11760 nonzeros\n",
"Model fingerprint: 0x93e351cb\n",
"Coefficient statistics:\n",
" Matrix range [3e-03, 5e+00]\n",
" Objective range [2e-01, 5e+01]\n",
" Bounds range [0e+00, 0e+00]\n",
" RHS range [5e-04, 7e+00]\n",
"Presolve time: 0.01s\n",
"Presolved: 1680 rows, 426 columns, 11760 nonzeros\n",
"\n",
"Iteration Objective Primal Inf. Dual Inf. Time\n",
" 0 handle free variables 0s\n",
" 1123 -2.9966123e+02 0.000000e+00 0.000000e+00 0s\n",
"\n",
"Solved in 1123 iterations and 0.06 seconds\n",
"Optimal objective -2.996612347e+02\n"
]
}
],
"source": [
"nbB = 2\n",
"thetemp = 1\n",
"epsilon_biy = special.digamma(1) -np.log(-np.log(np.random.uniform(0,1,ninds*nbchoices*nbB)))\n",
"lenobj = ninds*nbB+nbK\n",
"\n",
"newc = np.concatenate((muhat_iy.reshape(1,nbchoices*ninds).dot(Phi_iy_k).flatten(),np.repeat(-1/nbB,ninds*nbB)))\n",
"newm = grb.Model('new_lp')\n",
"x = newm.addMVar(lenobj, name='x', lb=-grb.GRB.INFINITY)\n",
"newm.setObjective(newc @ x, grb.GRB.MAXIMIZE)\n",
"mat1 = spr.kron(-Phi_iy_k, two_d(np.repeat(1,nbB)))\n",
"mat2 = spr.kron(two_d(np.repeat(1,nbchoices)),spr.identity(ninds*nbB))\n",
"newcstMat = spr.hstack((mat1, mat2))\n",
"rhs = epsilon_biy\n",
"newm.addConstr(newcstMat @ x >= rhs)\n",
"newm.optimize()"
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Value of the problem (Gurobi) = -299.661234703907\n"
]
}
],
"source": [
"if m.status == grb.GRB.Status.OPTIMAL:\n",
" print(\"Value of the problem (Gurobi) =\", newm.objval)\n",
" opt_x = np.array(newm.getAttr('x'))\n",
"newtheta_lp = opt_x[:nbK] / opt_x[nbK-1]"
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[ 0.10929498 1.41440717 -0.05870185 -0.05587768 -0.46498768 1. ]\n",
"[-0.01345289 1.2917597 -0.22401021 -0.00305825 -0.50896292 1. ]\n"
]
}
],
"source": [
"print(theta_mle)\n",
"print(newtheta_lp)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Finally probit"
]
},
{
"cell_type": "code",
"execution_count": 24,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Gurobi Optimizer version 9.1.1 build v9.1.1rc0 (win64)\n",
"Thread count: 6 physical cores, 12 logical processors, using up to 12 threads\n",
"Optimize a model with 1680 rows, 426 columns and 11760 nonzeros\n",
"Model fingerprint: 0x5e8fff22\n",
"Coefficient statistics:\n",
" Matrix range [3e-03, 5e+00]\n",
" Objective range [2e-01, 5e+01]\n",
" Bounds range [0e+00, 0e+00]\n",
" RHS range [2e+03, 2e+03]\n",
"Presolve time: 0.01s\n",
"Presolved: 1680 rows, 426 columns, 11760 nonzeros\n",
"\n",
"Iteration Objective Primal Inf. Dual Inf. Time\n",
" 0 handle free variables 0s\n",
" 750 -3.5266815e+05 0.000000e+00 0.000000e+00 0s\n",
"\n",
"Solved in 750 iterations and 0.04 seconds\n",
"Optimal objective -3.526681518e+05\n"
]
}
],
"source": [
"nbB = 2\n",
"thetemp = 1\n",
"epsilon_biy = np.random.normal(nbB*ninds*nbchoices)\n",
"lenobj = ninds*nbB+nbK\n",
"\n",
"newc = np.concatenate((muhat_iy.reshape(1,nbchoices*ninds).dot(Phi_iy_k).flatten(),np.repeat(-1/nbB,ninds*nbB)))\n",
"newm = grb.Model('new_lp')\n",
"x = newm.addMVar(lenobj, name='x', lb=-grb.GRB.INFINITY)\n",
"newm.setObjective(newc @ x, grb.GRB.MAXIMIZE)\n",
"mat1 = spr.kron(-Phi_iy_k, two_d(np.repeat(1,nbB)))\n",
"mat2 = spr.kron(two_d(np.repeat(1,nbchoices)),spr.identity(ninds*nbB))\n",
"newcstMat = spr.hstack((mat1, mat2))\n",
"rhs = epsilon_biy\n",
"newm.addConstr(newcstMat @ x >= rhs)\n",
"newm.optimize()"
]
}
],
"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.7.4"
}
},
"nbformat": 4,
"nbformat_minor": 2
}