19_black_market.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Foundations of Computational Economics #19\n",
    "\n",
    "by Fedor Iskhakov, ANU\n",
    "\n",
    "<img src=\"_static/img/dag3logo.png\" style=\"width:256px;\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "## Measuring the volume of illegal trade with linear programming\n",
    "\n",
    "<img src=\"_static/img/lab.png\" style=\"width:64px;\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "<img src=\"_static/img/youtube.png\" style=\"width:65px;\">\n",
    "\n",
    "[https://youtu.be/4z7MU73cx0M](https://youtu.be/4z7MU73cx0M)\n",
    "\n",
    "Description: Application of the optional transport problem."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Application: measuring illegal trade\n",
    "\n",
    "**“Black Market Performance: Illegal Trade in Beijing License Plates”**\n",
    "\n",
    "by\n",
    "Øystein Daljord, Mandy Hu, Guillaume Pouliot, Junji Xiao\n",
    "\n",
    "*From abstract:*\n",
    "\n",
    "We estimate the incentives to trade in the black market for license plates that emerged following the recent rationing of new car sales in Beijing by lottery. Under weak assumptions on car preferences, we use optimal transport methods and comprehensive data on car sales to estimate that at least 12% of the quota is illegally traded.\n",
    "\n",
    "[PDF for the paper (right-click and Save as…)](_static/pdf/DaljordHuPouliotXiao2019.pdf)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Øystein Daljord (1979-2020)\n",
    "\n",
    "<img src=\"_static/img/chicago-booth-oeystein-daljord.jpg\" style=\"width:400px;\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Black market of license plates\n",
    "\n",
    "- Measure the size of black market for license plates  \n",
    "- Case of Beijing license plates regulation  \n",
    "- Allocation by random lottery should have no effect on car sales  \n",
    "- In reality, there is sizable shift in distribution of cars  \n",
    "- Optimal transportation method is ideal to compute the lower bound on the volume of illegal trade of license plates  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Beijing license plate lottery\n",
    "\n",
    "- Cars driving in Beijing are required to have Beijing license plates  \n",
    "- From Jan 2011 license plates are rationed to a quota of about 35% of the previous year’s sales  \n",
    "- License plates are allocated by a lottery with simple application  \n",
    "- A Beijing household needs a license plate before it can register a new car  \n",
    "- License plates are non-transferable  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Material shift in distribution of cars\n",
    "\n",
    "- From cheaper to more expensive car models  \n",
    "- Hard to explain if lottery is a truly random allocation of license plates to the car purchasers  \n",
    "- No similar shifts in sales in comparable cities without rationing policy, in the same time period  \n",
    "- No supply side responses to the rationing policy  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Modeling framework\n",
    "\n",
    "Let $ \\mathbb{P}_0 $ be the distribution of car sales prices from pre-lottery time,\n",
    "and $ \\mathbb{P}_1 $ the analogous distribution post-lottery.\n",
    "\n",
    "Under assumptions\n",
    "\n",
    "1. Pricing policy did not change between 2010 and 2011  \n",
    "1. Demand structure did not change between 2010 and 2011  \n",
    "1. Lottery is uniform  \n",
    "\n",
    "\n",
    "the sales distributions should not change from the pre- to the post lottery period, i.e. $ \\mathbb{P}_0 = \\mathbb{P}_1 $"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Data\n",
    "\n",
    "Data on manufacturer suggested retail prices (MSRP) of the registered vehicles."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "hide-output": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Data has 243677 observations and 3 variables\n",
      "   year  month        MSRP\n",
      "0  2010      9  153.313139\n",
      "1  2010      9   44.543519\n",
      "2  2011      2   88.812069\n",
      "3  2010     11  210.732564\n",
      "4  2011      4  101.591900\n",
      "5  2010     12   56.428979\n",
      "6  2011      1  140.571004\n",
      "7  2010      8  170.066283\n",
      "8  2011     11  111.935614\n",
      "9  2010      7  191.988099\n"
     ]
    }
   ],
   "source": [
    "import pandas as pd\n",
    "dt = pd.read_stata('_static/data/beijin_data.dta')\n",
    "dt.dropna(inplace=True)  # drop rows with nan\n",
    "print('Data has %d observations and %d variables'%tuple(dt.shape))  # print expects tuple\n",
    "print(dt.head(n=10))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "hide-output": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "count    243677.000000\n",
      "mean        166.374743\n",
      "std         114.990191\n",
      "min          19.728006\n",
      "25%          91.780464\n",
      "50%         134.771217\n",
      "75%         209.744757\n",
      "max        1130.669488\n",
      "Name: MSRP, dtype: float64\n",
      "count    241240.000000\n",
      "mean        161.169630\n",
      "std         102.551081\n",
      "min          19.728006\n",
      "25%          91.229416\n",
      "50%         134.283479\n",
      "75%         207.808600\n",
      "max         617.902636\n",
      "Name: MSRP, dtype: float64\n"
     ]
    }
   ],
   "source": [
    "print(dt['MSRP'].describe())\n",
    "q99 = dt['MSRP'].quantile(0.99)\n",
    "dt = dt[dt['MSRP']<q99]\n",
    "print(dt['MSRP'].describe())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "hide-output": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "plt.rcParams['figure.figsize'] = [12, 8]\n",
    "\n",
    "def plot2hist(d1,d2,bins=10,labels=['1','2']):\n",
    "    '''Plots two overlapping histograms'''\n",
    "    plt.hist(d1,bins=bins,density=True,histtype='step',label=labels[0])\n",
    "    plt.hist(d2,bins=bins,density=True,histtype='step',label=labels[1])\n",
    "    plt.legend()\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "hide-output": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAs8AAAHSCAYAAAAT0iZvAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3df7BfZZ0n+PdHgoTqgbDGRGMu2cQi4yS0KyKilJalMBlQusiWaHdYZ5tuqeIf2HV0tqYu84ebtpqupEoFLLWrqMYxw/aINDNKqqF0GGmqd6ltMA7OIL+qs5CFiyg0YJRe+TnP/nEP9OVykzwx997v5d7XqyqV73nO85z7OXmKy/ue+5xzqrUWAADg0N4w6gIAAOD1QngGAIBOwjMAAHQSngEAoJPwDAAAnYRnAADotGzUBRyON7/5zW39+vWjLgMAgEXsRz/60d+11lbNtO91FZ7Xr1+fPXv2jLoMAAAWsar6fw+0z7INAADoJDwDAEAn4RkAADq9rtY8AwAwd1544YVMTEzk2WefHXUp82L58uUZGxvL0Ucf3T1GeAYAIEkyMTGR4447LuvXr09VjbqcOdVay5NPPpmJiYls2LChe5xlGwAAJEmeffbZrFy5ctEH5ySpqqxcufKwr7ILzwAAvGIpBOeX/SbnKjwDALBgPPLII/nIRz6STZs25eSTT85VV12VJHnqqaeyZcuWbNy4MVu2bMnTTz+dJLn//vtzxhln5JhjjskXv/jFVx3re9/7Xt7xjnfkpJNOyo4dO2alPmueAQCY0Qd23JpHf/HrWTve2hOOze3jZx60z7Jly/KlL30pp556an71q1/lPe95T7Zs2ZJvfvObOeusszI+Pp4dO3Zkx44d2blzZ970pjflK1/5Sr773e++6jgvvfRSLrnkktxyyy0ZGxvLe9/73px33nnZvHnzEZ2D8AwAwIwe/cWvs2/HubN2vPXjNx2yz5o1a7JmzZokyXHHHZdNmzbl0UcfzY033pjbbrstSXLhhRfmwx/+cHbu3JnVq1dn9erVuemmVx/7zjvvzEknnZS3v/3tSZJt27blxhtvPOLwbNkGAAAL0r59+3LXXXflfe97X37+85+/EqrXrFmTxx9//KBjH3300Zx44omvbI+NjeXRRx894pqEZwAAFpxnnnkm559/fq688socf/zxhz2+tfaattm4GVJ4BgBgQXnhhRdy/vnn51Of+lQ+/vGPJ0ne8pa35LHHHkuSPPbYY1m9evVBjzE2NpZHHnnkle2JiYm87W1vO+LahGcAABaM1louuuiibNq0KZ/73OdeaT/vvPOya9euJMmuXbuydevWgx7nve99b/72b/82Dz30UJ5//vlcd911Oe+88464PjcMAgCwYNx+++259tpr8853vjOnnHJKkuRP/uRPMj4+nt/93d/NNddck3Xr1uUv/uIvkiQ/+9nPctppp+WXv/xl3vCGN+TKK6/Mvffem+OPPz5f/epXc/bZZ+ell17Kpz/96Zx88slHXJ/wDADAjNaecGzXEzIO53iH8sEPfnDG9cpJ8oMf/OA1bW9961szMTExY/+Pfexj+djHPnZ4RR6C8AwAwIwO9UzmpciaZwAA6CQ8AwBAJ8s2+I3M9us6F7qe14kCAIuf8MxvZLZf17nQzebNEgDA65dlGwAA0El4BgBgwXjkkUfykY98JJs2bcrJJ5+cq666Kkny1FNPZcuWLdm4cWO2bNmSp59+Okly//3354wzzsgxxxyTL37xi6861qc//emsXr06v/3bvz1r9Vm2AQDAzK54Z7L/4dk73op1yWfvPmiXZcuW5Utf+lJOPfXU/OpXv8p73vOebNmyJd/85jdz1llnZXx8PDt27MiOHTuyc+fOvOlNb8pXvvKVfPe7333Nsf7gD/4gl156aX7/939/1k5BeAYAYGb7H06275+9421fccgua9asyZo1a5Ikxx13XDZt2pRHH300N954Y2677bYkyYUXXpgPf/jD2blzZ1avXp3Vq1fnpptee3/Shz70oezbt2/26o9lGwAALFD79u3LXXfdlfe97335+c9//kqoXrNmTR5//PGR1CQ8AwCw4DzzzDM5//zzc+WVV+b4448fdTmvEJ4BAFhQXnjhhZx//vn51Kc+lY9//ONJkre85S157LHHkiSPPfZYVq9ePZLausJzVZ1TVQ9U1d6qGp9h/zFV9e1h/x1VtX7KvsuG9geq6uwp7SdU1Q1VdX9V3VdVZ8zGCQEA8PrVWstFF12UTZs25XOf+9wr7eedd1527dqVJNm1a1e2bt06kvoOGZ6r6qgkX0vy0SSbk1xQVZundbsoydOttZOSXJFk5zB2c5JtSU5Ock6Srw/HS5KrknyvtfZPkrwryX1HfjoAALye3X777bn22mtz66235pRTTskpp5ySm2++OePj47nllluycePG3HLLLRkfn7ye+7Of/SxjY2P58pe/nD/+4z/O2NhYfvnLXyZJLrjggpxxxhl54IEHMjY2lmuuueaI6+t52sbpSfa21h5Mkqq6LsnWJPdO6bM1yfbh8w1JvlpVNbRf11p7LslDVbU3yelVdU+SDyX5gyRprT2f5PkjPhsAAGbPinVdT8g4rOMdwgc/+MG01mbc94Mf/OA1bW9961szMTExY/9vfetbh1dfh57wvDbJI1O2J5K870B9WmsvVtX+JCuH9r+ZNnZtkl8neSLJv6mqdyX5UZLPtNb+/jc5CQAA5sAhnsm8FPWsea4Z2qb/OHCgPgdqX5bk1CR/2lp7d5K/T/KatdRJUlUXV9WeqtrzxBNPdJQLAABzoyc8TyQ5ccr2WJKfHqhPVS1LsiLJUwcZO5FkorV2x9B+QybD9Gu01q5urZ3WWjtt1apVHeUCAMDc6AnPP0yysao2VNUbM3kD4O5pfXYnuXD4/Ikkt7bJxSq7k2wbnsaxIcnGJHe21n6W5JGqescw5qy8eg01AAAjcKD1xovRb3Kuh1zzPKxhvjTJ95McleQbrbV7quoLSfa01nYnuSbJtcMNgU9lMmBn6Hd9JoPxi0kuaa29NBz6f0ny50MgfzDJHx529QAAzJrly5fnySefzMqVKzP57IfFq7WWJ598MsuXLz+scT03DKa1dnOSm6e1fX7K52eTfPIAYy9PcvkM7T9OctrhFAsAwNwZGxvLxMRElsp9ZsuXL8/Y2NhhjekKzwAALH5HH310NmzYMOoyFjSv5wYAgE7CMwAAdBKeAQCgk/AMAACdhGcAAOgkPAMAQCfhGQAAOgnPAADQSXgGAIBOwjMAAHQSngEAoJPwDAAAnYRnAADoJDwDAEAn4RkAADoJzwAA0El4BgCATsIzAAB0Ep4BAKCT8AwAAJ2EZwAA6CQ8AwBAJ+EZAAA6Cc8AANBJeAYAgE7CMwAAdBKeAQCgk/AMAACdhGcAAOgkPAMAQCfhGQAAOgnPAADQSXgGAIBOwjMAAHQSngEAoJPwDAAAnYRnAADoJDwDAEAn4RkAADoJzwAA0El4BgCATsIzAAB0Ep4BAKCT8AwAAJ2EZwAA6CQ8AwBAJ+EZAAA6Cc8AANBJeAYAgE7CMwAAdBKeAQCgk/AMAACdhGcAAOgkPAMAQCfhGQAAOgnPAADQSXgGAIBOXeG5qs6pqgeqam9Vjc+w/5iq+vaw/46qWj9l32VD+wNVdfaU9n1VdXdV/biq9szGyQAAwFxadqgOVXVUkq8l2ZJkIskPq2p3a+3eKd0uSvJ0a+2kqtqWZGeS36uqzUm2JTk5yduS/Keq+settZeGcR9prf3dLJ4PAADMmZ4rz6cn2dtae7C19nyS65JsndZna5Jdw+cbkpxVVTW0X9dae6619lCSvcPxAADgdacnPK9N8siU7YmhbcY+rbUXk+xPsvIQY1uS/1hVP6qqiw+/dAAAmF+HXLaRpGZoa519Djb2A621n1bV6iS3VNX9rbW/fs0XnwzWFyfJunXrOsoFAIC50XPleSLJiVO2x5L89EB9qmpZkhVJnjrY2Nbay38/nuQ7OcByjtba1a2101prp61ataqjXAAAmBs94fmHSTZW1YaqemMmbwDcPa3P7iQXDp8/keTW1lob2rcNT+PYkGRjkjur6req6rgkqarfSvLPkvzkyE8HAADmziGXbbTWXqyqS5N8P8lRSb7RWrunqr6QZE9rbXeSa5JcW1V7M3nFedsw9p6quj7JvUleTHJJa+2lqnpLku9M3lOYZUn+XWvte3NwfgAAMGt61jyntXZzkpuntX1+yudnk3zyAGMvT3L5tLYHk7zrcIsFAIBR6grPLEFXvDPZ//ABd+9bnmT7vFUzP1asSz5796irAAAWMOGZme1/ONm+/4C714/flH07zp3HgubB9hWjrgAAWOC6Xs8NAAAIzwAA0E14BgCATsIzAAB0Ep4BAKCT8AwAAJ2EZwAA6CQ8AwBAJ+EZAAA6Cc8AANBJeAYAgE7CMwAAdBKeAQCgk/AMAACdhGcAAOgkPAMAQCfhGQAAOgnPAADQSXgGAIBOwjMAAHQSngEAoJPwDAAAnYRnAADoJDwDAEAn4RkAADoJzwAA0El4BgCATsIzAAB0Ep4BAKCT8AwAAJ2EZwAA6CQ8AwBAJ+EZAAA6Cc8AANBJeAYAgE7CMwAAdBKeAQCgk/AMAACdhGcAAOgkPAMAQCfhGQAAOgnPAADQSXgGAIBOwjMAAHRaNuoCFosP7Lg1j/7i16MuY9bsW56sH7/pgPvXnnDsPFYDALAwCM+z5NFf/Dr7dpw76jJmz/YsrvMBAJgFlm0AAEAn4RkAADoJzwAA0El4BgCATsIzAAB0Ep4BAKCT8AwAAJ2EZwAA6CQ8AwBAJ+EZAAA6dYXnqjqnqh6oqr1VNT7D/mOq6tvD/juqav2UfZcN7Q9U1dnTxh1VVXdV1V8e6YkAAMBcO2R4rqqjknwtyUeTbE5yQVVtntbtoiRPt9ZOSnJFkp3D2M1JtiU5Ock5Sb4+HO9ln0ly35GeBAAAzIeeK8+nJ9nbWnuwtfZ8kuuSbJ3WZ2uSXcPnG5KcVVU1tF/XWnuutfZQkr3D8VJVY0nOTfJnR34aAAAw93rC89okj0zZnhjaZuzTWnsxyf4kKw8x9sok/yrJfzvYF6+qi6tqT1XteeKJJzrKBQCAudETnmuGttbZZ8b2qvqdJI+31n50qC/eWru6tXZaa+20VatWHbpaAACYIz3heSLJiVO2x5L89EB9qmpZkhVJnjrI2A8kOa+q9mVyGciZVfV//Ab1AwDAvOkJzz9MsrGqNlTVGzN5A+DuaX12J7lw+PyJJLe21trQvm14GseGJBuT3Nlau6y1NtZaWz8c79bW2j+fhfMBAIA5s+xQHVprL1bVpUm+n+SoJN9ord1TVV9Isqe1tjvJNUmuraq9mbzivG0Ye09VXZ/k3iQvJrmktfbSHJ0LAADMqUOG5yRprd2c5OZpbZ+f8vnZJJ88wNjLk1x+kGPfluS2njoAAGCUvGEQAAA6Cc8AANBJeAYAgE7CMwAAdBKeAQCgk/AMAACdhGcAAOgkPAMAQCfhGQAAOgnPAADQSXgGAIBOwjMAAHQSngEAoJPwDAAAnYRnAADoJDwDAEAn4RkAADoJzwAA0El4BgCATsIzAAB0Ep4BAKCT8AwAAJ2EZwAA6CQ8AwBAJ+EZAAA6Cc8AANBJeAYAgE7CMwAAdBKeAQCgk/AMAACdhGcAAOi0bNQFwIKxYl2yfcWMu/YtT7J9XquZHyvWJZ+9e9RVAMDrhvAMLztIiFw/flP27Th3HouZJwf4YQEAmJllGwAA0El4BgCATsIzAAB0Ep4BAKCT8AwAAJ08bQOWsoM8nm9R8mg+AI6Q8AxL2VILkkvpBwUA5oRlGwAA0El4BgCATsIzAAB0Ep4BAKCT8AwAAJ2EZwAA6CQ8AwBAJ+EZAAA6Cc8AANBJeAYAgE7CMwAAdBKeAQCgk/AMAACdhGcAAOgkPAMAQCfhGQAAOgnPAADQSXgGAIBOXeG5qs6pqgeqam9Vjc+w/5iq+vaw/46qWj9l32VD+wNVdfbQtryq7qyq/1JV91TVH83WCQEAwFw5ZHiuqqOSfC3JR5NsTnJBVW2e1u2iJE+31k5KckWSncPYzUm2JTk5yTlJvj4c77kkZ7bW3pXklCTnVNX7Z+eUAABgbvRceT49yd7W2oOtteeTXJdk67Q+W5PsGj7fkOSsqqqh/brW2nOttYeS7E1yepv0zND/6OFPO8JzAQCAObWso8/aJI9M2Z5I8r4D9WmtvVhV+5OsHNr/ZtrYtckrV7R/lOSkJF9rrd0x0xevqouTXJwk69at6ygXZt/aE47N+vGbRl3GvFl7wrG5ffzMUZcBAAtOT3iuGdqmXyU+UJ8Djm2tvZTklKo6Icl3quq3W2s/eU3n1q5OcnWSnHbaaa5OMxJLLUgupR8UAOBw9CzbmEhy4pTtsSQ/PVCfqlqWZEWSp3rGttZ+keS2TK6JBgCABasnPP8wycaq2lBVb8zkDYC7p/XZneTC4fMnktzaWmtD+7bhaRwbkmxMcmdVrRquOKeqjk3yT5Pcf+SnAwAAc+eQyzaGNcyXJvl+kqOSfKO1dk9VfSHJntba7iTXJLm2qvZm8orztmHsPVV1fZJ7k7yY5JLW2ktVtSbJrmHd8xuSXN9a+8u5OEEAAJgtPWue01q7OcnN09o+P+Xzs0k+eYCxlye5fFrbf03y7sMtFgAARskbBgEAoJPwDAAAnYRnAADoJDwDAEAn4RkAADoJzwAA0El4BgCATsIzAAB0Ep4BAKCT8AwAAJ2EZwAA6CQ8AwBAJ+EZAAA6Cc8AANBJeAYAgE7CMwAAdBKeAQCgk/AMAACdhGcAAOgkPAMAQCfhGQAAOgnPAADQSXgGAIBOwjMAAHQSngEAoJPwDAAAnYRnAADoJDwDAEAn4RkAADoJzwAA0El4BgCATsIzAAB0Ep4BAKCT8AwAAJ2EZwAA6CQ8AwBAJ+EZAAA6Cc8AANBJeAYAgE7CMwAAdBKeAQCgk/AMAACdhGcAAOgkPAMAQCfhGQAAOgnPAADQSXgGAIBOwjMAAHQSngEAoJPwDAAAnYRnAADoJDwDAECnZaMuAFh41p5wbNaP3zTqMmbdvuWZ8bzWnnBsbh8/cwQVAfB6IzwDr7Fog+T2ZN+Oc1/TvBh/UABgbli2AQAAnYRnAADoJDwDAECnrvBcVedU1QNVtbeqxmfYf0xVfXvYf0dVrZ+y77Kh/YGqOntoO7Gq/qqq7quqe6rqM7N1QgAAMFcOGZ6r6qgkX0vy0SSbk1xQVZundbsoydOttZOSXJFk5zB2c5JtSU5Ock6Srw/HezHJv2ytbUry/iSXzHBMAABYUHquPJ+eZG9r7cHW2vNJrkuydVqfrUl2DZ9vSHJWVdXQfl1r7bnW2kNJ9iY5vbX2WGvtPydJa+1XSe5LsvbITwcAAOZOz6Pq1iZ5ZMr2RJL3HahPa+3FqtqfZOXQ/jfTxr4qJA9LPN6d5I7DqHt+XfHOZP/DB+2yb3mS7fNSzfxYsW7UFQAALDg94blmaGudfQ46tqr+UZJ/n+RftNZ+OeMXr7o4ycVJsm7diALd/oeT7fsP2mX9+E0zPj8WAIDFo2fZxkSSE6dsjyX56YH6VNWyJCuSPHWwsVV1dCaD85+31v7Dgb54a+3q1tpprbXTVq1a1VEuAADMjZ7w/MMkG6tqQ1W9MZM3AO6e1md3kguHz59IcmtrrQ3t24ancWxIsjHJncN66GuS3Nda+/JsnAgAAMy1Qy7bGNYwX5rk+0mOSvKN1to9VfWFJHtaa7szGYSvraq9mbzivG0Ye09VXZ/k3kw+YeOS1tpLVfXBJP9zkrur6sfDl/rXrbWbZ/sEAQBgtvSsec4Qam+e1vb5KZ+fTfLJA4y9PMnl09r+r8y8HhoAABYsbxgEAIBOwjMAAHQSngEAoJPwDAAAnYRnAADoJDwDAEAn4RkAADoJzwAA0El4BgCATsIzAAB0Ep4BAKCT8AwAAJ2EZwAA6CQ8AwBAJ+EZAAA6Cc8AANBJeAYAgE7CMwAAdBKeAQCgk/AMAACdhGcAAOgkPAMAQCfhGQAAOgnPAADQSXgGAIBOwjMAAHQSngEAoJPwDAAAnYRnAADoJDwDAEAn4RkAADoJzwAA0El4BgCATsIzAAB0Ep4BAKCT8AwAAJ2EZwAA6CQ8AwBAJ+EZAAA6Cc8AANBJeAYAgE7CMwAAdBKeAQCgk/AMAACdhGcAAOgkPAMAQCfhGQAAOgnPAADQadmoCwCYNyvWJdtXvKZ53/Ik2+e9mvmxYl3y2btHXQXAoiE8A0vHAULk+vGbsm/HufNczDyZ4YcFAH5zlm0AAEAn4RkAADoJzwAA0El4BgCATsIzAAB0Ep4BAKCT8AwAAJ2EZwAA6NQVnqvqnKp6oKr2VtX4DPuPqapvD/vvqKr1U/ZdNrQ/UFVnT2n/RlU9XlU/mY0TAQCAuXbI8FxVRyX5WpKPJtmc5IKq2jyt20VJnm6tnZTkiiQ7h7Gbk2xLcnKSc5J8fTheknxzaAMAgNeFnivPpyfZ21p7sLX2fJLrkmyd1mdrkl3D5xuSnFVVNbRf11p7rrX2UJK9w/HSWvvrJE/NwjkAAMC86AnPa5M8MmV7YmibsU9r7cUk+5Os7BwLAACvCz3huWZoa519esYe/ItXXVxVe6pqzxNPPHE4QwEAYFb1hOeJJCdO2R5L8tMD9amqZUlWZHJJRs/Yg2qtXd1aO621dtqqVasOZygAAMyqnvD8wyQbq2pDVb0xkzcA7p7WZ3eSC4fPn0hya2utDe3bhqdxbEiyMcmds1M6AADMr0OG52EN86VJvp/kviTXt9buqaovVNV5Q7drkqysqr1JPpdkfBh7T5Lrk9yb5HtJLmmtvZQkVfWtJP93kndU1URVXTS7pwYAALNrWU+n1trNSW6e1vb5KZ+fTfLJA4y9PMnlM7RfcFiVAgDAiHnDIAAAdBKeAQCgk/AMAACdhGcAAOgkPAMAQCfhGQAAOgnPAADQSXgGAIBOwjMAAHQSngEAoJPwDAAAnYRnAADoJDwDAEAn4RkAADoJzwAA0El4BgCATsIzAAB0Ep4BAKCT8AwAAJ2EZwAA6CQ8AwBAJ+EZAAA6LRt1AQDAEbrincn+h0ddxfxasS757N2jroIlSHgGgNe7/Q8n2/ePuor5tX3FqCtgibJsAwAAOgnPAADQSXgGAIBOwjMAAHRywyCw5K094disH79p1GXMiX3L85pzW3vCsbl9/MwRVQTw+iY8A0veog6S25N9O859VdNi/UEBYD5YtgEAAJ1ceQZYzFase83zcPctT7J9JNXMDy/PAOaQ8AywmM0QIteP3/SapRyLipdnAHPIsg0AAOgkPAMAQCfhGQAAOgnPAADQSXgGAIBOwjMAAHTyqDoAFpcZnm296K1YN+oKYMkQngFYXLwgBZhDwjPAErP2hGOzfvymUZcxb9aecGxuHz9z1GUw25bqbxj8cDhywjPAErPUguRS+kFhSVmKIXKp/bCwQLlhEAAAOgnPAADQSXgGAIBO1jwDwCLygR235tFf/HrUZcwbN4Qy34RnAFhEHv3Fr7Nvx7mjLmPeuCGU+WbZBgAAdBKeAQCgk2UbACxqS/GlMMDcEZ4BWNTcTAbMJss2AACgk/AMAACdhGcAAOhkzTMA8Lq1lG4I3bd88iU41vGPlvAMALxuLakguT1L6u2RC5VlGwAA0KnrynNVnZPkqiRHJfmz1tqOafuPSfJvk7wnyZNJfq+1tm/Yd1mSi5K8lOR/ba19v+eYAABMsWJd9uV/SraPupD581hWZc32vaMu41UOGZ6r6qgkX0uyJclEkh9W1e7W2r1Tul2U5OnW2klVtS3JziS/V1Wbk2xLcnKStyX5T1X1j4cxhzomAAAv++zdo65g3q3ZvmLUJbxGz7KN05Psba092Fp7Psl1SbZO67M1ya7h8w1JzqqqGtqva60911p7KMne4Xg9xwQAgAWlJzyvTfLIlO2JoW3GPq21F5PsT7LyIGN7jgkAAAtKz5rnmqGtdfY5UPtMoX36MScPXHVxkouHzWeq6oED1DnVm5P8XUe/fn8006m8Wu2c1a+4GM3+vHCkzMnCZF4WJvOyMJmXhWn25qUjg82B//5AO3rC80SSE6dsjyX56QH6TFTVsiQrkjx1iLGHOmaSpLV2dZKrO+p8RVXtaa2ddjhjmHvmZeExJwuTeVmYzMvCZF4WpsU8Lz3LNn6YZGNVbaiqN2byBsDd0/rsTnLh8PkTSW5trbWhfVtVHVNVG5JsTHJn5zEBAGBBOeSV59bai1V1aZLvZ/Kxct9ord1TVV9Isqe1tjvJNUmuraq9mbzivG0Ye09VXZ/k3iQvJrmktfZSksx0zNk/PQAAmD01eYF4camqi4flHiwg5mXhMScLk3lZmMzLwmReFqbFPC+LMjwDAMBc8HpuAADotKjCc1WdU1UPVNXeqhofdT1LSVV9o6oer6qfTGl7U1XdUlV/O/z93w3tVVVfGebpv1bVqaOrfHGrqhOr6q+q6r6quqeqPjO0m5sRqqrlVXVnVf2XYV7+aGjfUFV3DPPy7eGG6gw3XX97mJc7qmr9KOtfzKrqqKq6q6r+ctg2JyNWVfuq6u6q+nFV7RnafA8bsao6oapuqKr7h//HnLFU5mXRhOf6h9eIfzTJ5iQX1OTrwZkf30xyzrS28SQ/aK1tTPKDYTuZnKONw5+Lk/zpPNW4FL2Y5F+21jYleX+SS4b/LszNaD2X5MzW2ruSnJLknKp6f5KdSa4Y5uXpJBcN/S9K8nRr7aQkVwz9mBufSXLflG1zsjB8pLV2ypRHn/keNnpXJflea+2fJHlXJv+7WRLzsmjCc7zye6Raa3+dySetTDX1te27kvyPU9r/bZv0N0lOqKo181Pp0tJae6y19p+Hz7/K5De3tTE3IzX8+z4zbB49/GlJzkxyw9A+fV5enq8bkpxVVSN5a8BiVlVjSc5N8mfDdsWcLFS+h41QVR2f5EOZfNpaWmvPt9Z+kSUyL4spPHvl98LzltbaY8lkiEuyemg3VyMw/Fr53UnuiLkZuWF5wI+TPJ7kliT/T5JftNZeHLpM/bd/ZV6G/fuTrDdcYA8AAAKUSURBVJzfipeEK5P8qyT/bdheGXOyELQk/7GqflSTbx1OfA8btbcneSLJvxmWOf1ZVf1Wlsi8LKbw3PMacRYGczXPquofJfn3Sf5Fa+2XB+s6Q5u5mQOttZdaa6dk8g2rpyfZNFO34W/zMseq6neSPN5a+9HU5hm6mpP594HW2qmZ/NX/JVX1oYP0NS/zY1mSU5P8aWvt3Un+Pv+wRGMmi2peFlN47nmNOPPr5y//Wmb4+/Gh3VzNo6o6OpPB+c9ba/9haDY3C8Twq87bMrkm/YSqevnlVVP/7V+Zl2H/irx2mRRH5gNJzquqfZlc9ndmJq9Em5MRa639dPj78STfyeQPm76HjdZEkonW2h3D9g2ZDNNLYl4WU3j2yu+FZ+pr2y9McuOU9t8f7r59f5L9L/+ah9k1rMG8Jsl9rbUvT9llbkaoqlZV1QnD52OT/NNMrkf/qySfGLpNn5eX5+sTSW5tHtI/q1prl7XWxlpr6zP5/49bW2ufijkZqar6rao67uXPSf5Zkp/E97CRaq39LMkjVfWOoemsTL5NeknMy6J6SUpVfSyTVwpefuX35SMuacmoqm8l+XCSNyf5eZL/Pcl3k1yfZF2Sh5N8srX21BDovprJp3P8f0n+sLW2ZxR1L3ZV9cEk/2eSu/MP6zj/dSbXPZubEamq/yGTN9MclcmLGNe31r5QVW/P5FXPNyW5K8k/b609V1XLk1ybyTXrTyXZ1lp7cDTVL35V9eEk/1tr7XfMyWgN//7fGTaXJfl3rbXLq2plfA8bqao6JZM3174xyYNJ/jDD97Ms8nlZVOEZAADm0mJatgEAAHNKeAYAgE7CMwAAdBKeAQCgk/AMAACdhGcAAOgkPAMAQCfhGQAAOv3/1HIsBSyLIXsAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 864x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "dt10 = dt[dt['year']==2010]['MSRP']\n",
    "dt11 = dt[dt['year']==2011]['MSRP']\n",
    "plot2hist(dt10,dt11,labels=['2010','2011'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Optimal transport problem\n",
    "\n",
    "$$\n",
    "\\min \\sum_{i=1}^{m}\\sum_{j=1}^{n} cost_{ij} x_{ij}, \\text{ subject to}\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\sum_{i=1}^{m} x_{ij} = origin_j, j \\in \\{1,\\dots,n\\},\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\sum_{j=1}^{n} x_{ij} = destination_i, i \\in \\{1,\\dots,m\\},\n",
    "$$\n",
    "\n",
    "$$\n",
    "x_{ij} \\ge 0 \\text{ for all } i,j\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### A linear programming problem\n",
    "\n",
    "Linear programming problem solved by `scipy.optimize.linprog` (equality constraints automatically converted)\n",
    "\n",
    "$$\n",
    "\\max(c \\cdot x) \\text{ subject to }\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\begin{array}{l}\n",
    "A_{ub}x \\le b_{ub} \\\\\n",
    "A_{eq}x = b_{eq} \\\\\n",
    "l \\le x \\le u\n",
    "\\end{array}\n",
    "$$\n",
    "\n",
    "- stack all $ x_{ij} $ into a single vector  \n",
    "- express equality constraints for origins and destinations as inequalities  \n",
    "\n",
    "\n",
    "[https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.linprog.html](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.linprog.html)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "hide-output": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "# Code up the model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "hide-output": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "#\n",
    "#\n",
    "#\n",
    "#\n",
    "#\n",
    "#\n",
    "#\n",
    "#\n",
    "#\n",
    "#\n",
    "#\n",
    "#\n",
    "#\n",
    "#\n",
    "#\n",
    "#\n",
    "#\n",
    "#\n",
    "#\n",
    "#\n",
    "#\n",
    "#\n",
    "# Answer below"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "hide-output": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "gr\n",
      "0    0.231451\n",
      "1    0.211331\n",
      "2    0.205582\n",
      "3    0.189273\n",
      "4    0.162363\n",
      "Name: MSRP, dtype: float64\n",
      "\n",
      "gr\n",
      "0    0.123699\n",
      "1    0.172511\n",
      "2    0.186457\n",
      "3    0.226024\n",
      "4    0.291310\n",
      "Name: MSRP, dtype: float64\n"
     ]
    }
   ],
   "source": [
    "N = 5  # number of bins to represent distribution\n",
    "dt['gr'] = pd.qcut(dt.MSRP,q=N,labels=False)  # N quantiles\n",
    "gr10 = dt[dt.year==2010].groupby('gr')\n",
    "gr11 = dt[dt.year==2011].groupby('gr')\n",
    "d10 = gr10.MSRP.count()/dt[dt.year==2010].MSRP.count()\n",
    "d11 = gr11.MSRP.count()/dt[dt.year==2011].MSRP.count()\n",
    "print(d10,d11,sep='\\n\\n')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "hide-output": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x7fbd40c5e350>"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAASsAAAEyCAYAAACiffbnAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAKiklEQVR4nO3dT4jnd33H8de7mxWlFjxkDyEbOh5EGgKNsAQhnoKH9Q/aYwR7EvZSYYWCaC/FQ2+lePESbFBQDIIeJAQkYFJbsNHdGFu3qxDE4qKwW0Kp8VCJefcwQ5ndzM7vt8nvN7/f+7ePBwzMb+Y333l/Z377nO/3O7/ZT3V3ALbdH216AIBliBUwglgBI4gVMIJYASOIFTDC2FhV1fmq+nlVvVxVn9v0PKtSVU9W1fWq+ummZ1mlqnqgqp6rqqtVdaWqLm56plWoqrdX1Q+r6icH+/WFTc+0alV1qqp+XFVPb3KOkbGqqlNJvpTkQ0keTPKJqnpws1OtzFeSnN/0EGvwWpK/7u4/S/L+JH+1I9+z/03yWHf/eZKHk5yvqvdveKZVu5jk6qaHGBmrJI8kebm7f9Hdv0/yVJKPb3imleju7yd5ZdNzrFp3/6a7Xzx4/bfZf/Dfv9mp3rre9+rBzdMHLzvzTOuqOpvkI0m+vOlZpsbq/iS/OnT7WnbggX+3qKq9JO9L8sJmJ1mNg9Okl5JcT/Jsd+/Efh34YpLPJnl904NMjVUd8bad+Wm2y6rqnUm+leQz3f0/m55nFbr7D939cJKzSR6pqoc2PdMqVNVHk1zv7subniWZG6trSR44dPtskl9vaBaWVFWnsx+qr3f3tzc9z6p1938neT67c83x0SQfq6pfZv9Sy2NV9bVNDTM1Vj9K8p6qendVvS3J40m+s+GZOEZVVZJ/THK1u/9h0/OsSlWdqap3Hbz+jiQfTPKzzU61Gt39+e4+29172f839r3u/uSm5hkZq+5+Lcmnk3w3+xdqv9ndVzY71WpU1TeS/CDJe6vqWlV9atMzrcijSf4y+z+dXzp4+fCmh1qB+5I8V1X/lv0fos9290Z/xb+ryn8RA0ww8sgKuPuIFTCCWAEjiBUwglgBI4yPVVVd2PQM62C/5tnVfduW/RofqyRb8YVcA/s1z67u21bs1y7ECrgLrOVJoffee2/v7e2tfLtHuXHjRs6cOXMin+vy5a34e07Yad191H9UkHvW8cn29vZy6dKldWx6o/b/vA3YBKeBwAhiBYwgVsAIYgWMIFbACGIFjCBWwAhiBYwgVsAIYgWMIFbACGIFjCBWwAhiBYwgVsAIYgWMIFbACGIFjCBWwAhiBYwgVsAIS8Wqqs5X1c+r6uWq+ty6hwK41cJYVdWpJF9K8qEkDyb5RFU9uO7BAA5b5sjqkSQvd/cvuvv3SZ5K8vH1jgVws2VidX+SXx26fe3gbQAnZplYHbUM8RvWnK+qC1V1qaou3bhx461PBnDIMrG6luSBQ7fPJvn1rXfq7ie6+1x3nztz5syq5gNIslysfpTkPVX17qp6W5LHk3xnvWMB3OyeRXfo7teq6tNJvpvkVJInu/vK2icDOGRhrJKku59J8syaZwG4Lc9gB0YQK2AEsQJGECtgBLECRhArYASxAkYQK2AEsQJGECtgBLECRhArYASxAkYQK2AEsQJGECtgBLECRhArYASxAkYQK2AEsQJGqO43LK781jdatfqNboF1fK22RdVRC2/DyevuIx+MjqyAEcQKGEGsgBHEChhBrIARxAoYQayAEcQKGEGsgBHEChhBrIARxAoYQayAEcQKGEGsgBHEChhBrIARxAoYQayAEcQKGEGsgBHEChhBrIARxAoYYWGsqurJqrpeVT89iYEAjrLMkdVXkpxf8xwAx1oYq+7+fpJXTmAWgNu6Z1UbqqoLSS6sansAh1V3L75T1V6Sp7v7oaU2WrV4owMt87Waqqo2PQIkSbr7yAej3wYCI4gVMMIyT134RpIfJHlvVV2rqk+tfyyAmy11zeqON+qa1TiuWbEtXLMCRhMrYASxAkYQK2AEsQJGECtgBLECRhArYASxAkYQK2AEsQJGECtgBLECRhArYASxAkYQK2AEsQJGECtgBLECRhArYASxAkZY2fLxd4NdXgFmV1fu2eXv2d3GkRUwglgBI4gVMIJYASOIFTCCWAEjiBUwglgBI4gVMIJYASOIFTCCWAEjiBUwglgBI4gVMIJYASOIFTCCWAEjiBUwglgBI4gVMIJYASOIFTDCwlhV1QNV9VxVXa2qK1V18SQGAzisFi1uWVX3Jbmvu1+sqj9JcjnJX3T3fxzzMbu5YuYOs8gp26K7j/ymLTyy6u7fdPeLB6//NsnVJPevdjyA493R8vFVtZfkfUleOOJ9F5JcWMlUALdYeBr4/3esemeSf0ryd9397QX33c1zih3mNJBt8aZPA5Okqk4n+VaSry8KFcA6LHOBvZJ8Nckr3f2ZpTbqyGocR1Zsi9sdWS0Tqw8k+eck/57k9YM3/013P3PMx+zmI3+HiRXb4k3H6s0Qq3nEim3xlq5ZAWyaWAEjiBUwglgBI4gVMIJYASOIFTCCWAEjiBUwglgBI4gVMIJYASOIFTCCWAEjiBUwglgBI4gVMIJYASOIFTCCWAEjiBUwwh0tH8/u2tVVYHZ11Z5kd79nt+PIChhBrIARxAoYQayAEcQKGEGsgBHEChhBrIARxAoYQayAEcQKGEGsgBHEChhBrIARxAoYQayAEcQKGEGsgBHEChhBrIARxAoYQayAEcQKGGFhrKrq7VX1w6r6SVVdqaovnMRgAIfVokUga38lxT/u7ler6nSSf0lysbv/9ZiP2d2VJRnFIqfzdPeRO7ZwRebe/26/enDz9MHL7j4CgK201DWrqjpVVS8luZ7k2e5+Yb1jAdxsqVh19x+6++EkZ5M8UlUP3XqfqrpQVZeq6tKqhwRYeM3qDR9Q9bdJftfdf3/MfZwmshVcs5rndteslvlt4JmqetfB6+9I8sEkP1vteADHW3iBPcl9Sb5aVaeyH7dvdvfT6x0L4GZ3fBq41EadBrIlnAbO86ZPAwG2gVgBI4gVMIJYASOIFTCCWAEjiBUwglgBI4gVMIJYASOIFTCCWAEjiBUwglgBI4gVMIJYASOIFTCCWAEjiBUwglgBI4gVMMIyS3HBWLu6Akyymyv3nDt37rbvc2QFjCBWwAhiBYwgVsAIYgWMIFbACGIFjCBWwAhiBYwgVsAIYgWMIFbACGIFjCBWwAhiBYwgVsAIYgWMIFbACGIFjCBWwAhiBYwgVsAIYgWMIFbACEvHqqpOVdWPq+rpdQ4EcJQ7ObK6mOTqugYBOM5Ssaqqs0k+kuTL6x0H4GjLHll9Mclnk7x+uztU1YWqulRVl1YyGcAhC2NVVR9Ncr27Lx93v+5+orvPdfe5lU0HcGCZI6tHk3ysqn6Z5Kkkj1XV19Y6FcAtFsaquz/f3We7ey/J40m+192fXPtkAId4nhUwwj13cufufj7J82uZBOAYjqyAEcQKGEGsgBHEChhBrIARxAoYQayAEcQKGEGsgBHEChhBrIARxAoYQayAEcQKGEGsgBHEChhBrIARxAoYQayAEcQKGEGsgBHuaHWbO/BfSf5zTdu+1b0Hn2/X2K95TnTfquqkPtVJ7tef3u4d1d0nNMN6VNWlXVyy3n7Ns6v7ti375TQQGEGsgBF2IVZPbHqANbFf8+zqvm3Ffo2/ZgXcHXbhyAq4C4gVMIJYASOIFTCCWAEj/B/lWx5/Cmh+rAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 360x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "# Set up transportation problem\n",
    "costs = np.ones((N,N)) - np.eye(N)  # costs matrix\n",
    "origins = np.array(d10)        # origins\n",
    "destinations = np.array(d11)   # destinations\n",
    "plt.rcParams['figure.figsize'] = [5, 5]\n",
    "plt.spy(costs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "hide-output": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAATIAAACSCAYAAADGrUafAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAIfElEQVR4nO3dX4hmdR3H8fenXaX8ExluUWqpIZF0oe4QliGWEWaRBQUGhnVjF1oaQlg3dtNdmV2EsKkpZEqolUSUYkp1szSzLek6RWKmq5s7EqR0Y+a3i+fZXGZnds8zM+d55nd8v2CZZ86emef73bPz4Zwz5/f7paqQpJa9btYFSNJ6GWSSmmeQSWqeQSapeQaZpOYZZJKaN5UgS3JRkr8keTzJddN4z2lK8mSSR5LsTjI/63rWK8mtSfYnefSgbW9O8kCSv44/njDLGtdqld6+meSZ8fHbneTiWda4HklOSfJQksUke5JcPd4+iOO3mt6DLMkW4PvAx4Azgc8lObPv952BD1XVWVU1N+tCNsBtwEXLtl0HPFhVZwAPjj9v0W0c2hvAd8fH76yq+uWUa9pILwPXVtV7gHOBK8c/b0M5fiuaxhnZ+4DHq+qJqnoJuAu4ZArvqzWqqt8C/1y2+RLg9vHr24FPTbWoDbJKb4NRVfuqatf49YvAInASAzl+q5lGkJ0EPH3Q53vH24akgPuTLCS5YtbF9OStVbUPRj8swFtmXM9GuyrJn8aXnoO47EpyKnA2sJOBH79pBFlW2Da0cVHnVdU5jC6fr0xy/qwL0kRuAt4FnAXsA74z23LWL8lxwD3ANVX1wqzr6ds0gmwvcMpBn58MPDuF952aqnp2/HE/8FNGl9ND81yStwGMP+6fcT0bpqqeq6r/VtUrwA9o/PglOYpRiN1RVfeONw/2+MF0guwPwBlJTktyNHApcN8U3ncqkhyb5PgDr4GPAo8e/quadB9w+fj15cDPZ1jLhjrwAz72aRo+fkkC3AIsVtUNB/3VYI8fQKYx+8X419k3AluAW6vqW72/6ZQkOZ3RWRjAVuDHrfeX5E7gAuBE4DngeuBnwE+AdwBPAZ+tquZumq/S2wWMLisLeBL40oH7Sa1J8kHgd8AjwCvjzd9gdJ+s+eO3mqkEmST1ySf7JTXPIJPUPINMUvMMMknNM8gkNW9qQTbgoTuA/bVuyP0NubcDpnlGNvR/TPtr25D7G3JvgJeWkgaglwdikwz6Kdvt27cfsm1paYlt27bNoJrpsL92Dam3hYWF56vqkGYMsjVwNIQ0G0kWVpq8tNOl5dCnqpbUtiMG2WtoqmpJjepyRuZU1ZI2tS5B1mmq6iRXJJkfwipCktqytcM+naaqrqodwA4Y/s1+SZtLlzOywU9VLaltXYJs0FNVS2rfES8tq+rlJFcBv+bVqar39F6ZJHXUywOxc3NzNT8/3Hv+o/UdhssHfrVZreuBWEnazAwySc0zyCQ1zyCT1DyDTFLzDDJJzTPIJDXPIJPUPINMUvMMMknN6zKNj5YZ+hCeIQ/BGvqxe63yjExS8wwySc3rsvjIKUkeSrKYZE+Sq6dRmCR11eUe2cvAtVW1K8nxwEKSB6rqsZ5rk6ROjnhGVlX7qmrX+PWLwCIrLD4iSbMy0T2yJKcCZwM7+yhGktaic5AlOQ64B7imql5Y4e//vxzc0tLSRtYoSYfVKciSHMUoxO6oqntX2qeqdlTVXFXNbdu2bSNrlKTD6vJbywC3AItVdUP/JUnSZLqckZ0HfB74cJLd4z8X91yXJHXWZTm437PyauOStCk41lKHGPJ4xCGPI4VhH7vDcYiSpOYZZJKaZ5BJap5BJql5Bpmk5hlkkppnkElqnkEmqXkGmaTmGWSSmucQJb2mDH0Iz9CHYK3GMzJJzZtkhtgtSf6Y5Bd9FiRJk5rkjOxqRguPSNKm0nWq65OBjwM391uOJE2u6xnZjcDXgFd6rEWS1qTLnP2fAPZX1cIR9nMVJUkz0XXO/k8meRK4i9Hc/T9avpOrKEmalS4rjX+9qk6uqlOBS4HfVNVlvVcmSR35HJmk5k30ZH9VPQw83EslkrRG6WPIRpLO37TFISOTDAOxv81l0iE8Q+6vtd4AkixU1dzy7V5aSmqeQSapeQaZpOYZZJKaZ5BJap5BJql5Bpmk5hlkkppnkElqnkEmqXm9rKK0fft25ufnO+3b4pCKSeqwv7V9375MWsOQ+2utt8PxjExS8wwySc3ruvjIm5LcneTPSRaTvL/vwiSpq673yL4H/KqqPpPkaOCYHmuSpIkcMciSvBE4H/gCQFW9BLzUb1mS1F2XS8vTgSXgh+OVxm9OcuzynVxFSdKsdAmyrcA5wE1VdTbwb+C65Tu5ipKkWekSZHuBvVW1c/z53YyCTZI2hS7Lwf0DeDrJu8ebLgQe67UqSZpA199afhm4Y/wbyyeAL/ZXkiRNplOQVdVu4JCVSzbC0IdU2N+r7G9t37cvffU26ffeCD7ZL6l5Bpmk5hlkkppnkElqnkEmqXkGmaTmGWSSmmeQSWqeQSapeQaZpOYZZJKa18tycH0Z0tiw9dbQ2rg+sL+DtdbfZl9GzzMySc3ruorSV5PsSfJokjuTvL7vwiSpqyMGWZKTgK8Ac1X1XmALcGnfhUlSV10vLbcCb0iyldFScM/2V5IkTabLVNfPAN8GngL2Af+qqvv7LkySuupyaXkCcAlwGvB24Ngkl62wn8vBSZqJLpeWHwH+VlVLVfUf4F7gA8t3cjk4SbPSJcieAs5NckxGD4dcCCz2W5YkddflHtlORmtZ7gIeGX/Njp7rkqTOuq6idD1wfc+1SNKaNDVEaRKbfUjFejlc61WtHTuwv4NN+v9zJQ5RktQ8g0xS8wwySc0zyCQ1zyCT1DyDTFLzDDJJzTPIJDXPIJPUPINMUvPSx5CGJEvA35dtPhF4fsPfbPOwv7YNub8h9fbOqjpknrBegmwlSearam4qbzYD9te2Ifc35N4O8NJSUvMMMknNm2aQDX0yRvtr25D7G3JvwBTvkUlSX7y0lNQ8g0xS8wwySc0zyCQ1zyCT1Lz/AZO34Z1Rkp+uAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 360x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# convert to linear programming problem\n",
    "C = costs.reshape(N*N)\n",
    "A1 = np.kron(np.eye(N),np.ones((1,N)))  # sums of x for each origin\n",
    "A2 = np.kron(np.ones((1,N)),np.eye(N))  # sums of x for each destination\n",
    "A = np.vstack((A1,A2))  # concatenate vertically\n",
    "plt.spy(A)\n",
    "b = np.concatenate((origins,destinations))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "hide-output": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimization terminated successfully.\n",
      "[[0.12369875 0.         0.         0.         0.10775178]\n",
      " [0.         0.17251076 0.         0.01762504 0.02119497]\n",
      " [0.         0.         0.18645705 0.01912521 0.        ]\n",
      " [0.         0.         0.         0.18927336 0.        ]\n",
      " [0.         0.         0.         0.         0.16236309]]\n",
      "With N=5 the lower bound on black market share is 0.16570\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAASsAAAEyCAYAAACiffbnAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAKX0lEQVR4nO3dwYukd53H8c93JyMKLnjYPkgmbDyIrAibME0QchtcmFXRawJ6EuayQgRB9Og/IF72MmhwQTEIepDAIgETRHBjemIUx9EliIuDwvQiorko0e8eutedmfRM16xV/fS35vWCgqqumurvb6bn3U899XQ/1d0BOO3+ZukBAFYhVsAIYgWMIFbACGIFjCBWwAhjY1VVF6vqZ1X1alV9eul51qWqnq6qG1X146VnWaeqeqiqnq+qa1V1taqeWnqmdaiqN1fV96vqh4fr+uzSM61bVZ2pqh9U1bNLzjEyVlV1Jsm/JvnnJO9O8mRVvXvZqdbmS0kuLj3EBrye5JPd/Q9J3pvkX7bk3+wPSS509z8meSTJxap678IzrdtTSa4tPcTIWCV5LMmr3f3z7v5jkmeSfHjhmdaiu7+T5DdLz7Fu3f3r7n758Prvc/DF/+CyU/31+sBrhzfPHl625kjrqjqX5ANJvrD0LFNj9WCSX950+3q24Av/flFVDyd5NMmLy06yHocvk15JciPJc929Fes69Pkkn0ry56UHmRqrOuJjW/PdbJtV1VuTfD3JJ7r7d0vPsw7d/afufiTJuSSPVdV7lp5pHarqg0ludPeVpWdJ5sbqepKHbrp9LsmvFpqFFVXV2RyE6ivd/Y2l51m37v5tkheyPfscH0/yoar6RQ52tVyoqi8vNczUWL2U5J1V9Y6qelOSJ5J8c+GZuIuqqiRfTHKtuz+39DzrUlU7VfW2w+tvSfK+JD9ddqr16O7PdPe57n44B//Hvt3dH1lqnpGx6u7Xk3w8ybdysKP2a919ddmp1qOqvprke0neVVXXq+pjS8+0Jo8n+WgOvju/cnh5/9JDrcHbkzxfVT/KwTfR57p70bf4t1X5FTHABCO3rID7j1gBI4gVMIJYASOIFTDC+FhV1aWlZ9gE65pnW9d2WtY1PlZJTsVf5AZY1zzburZTsa5tiBVwH9jIQaFVtZVHmp4/f/7EPtf+/n52dnZO7POdlG1dV3Lya7ty5VT8fPHadfdRv6hArO6Fo/05TQ5+3HL73ClWXgYCI4gVMIJYASOIFTCCWAEjiBUwglgBI4gVMIJYASOIFTCCWAEjiBUwglgBI4gVMIJYASOIFTCCWAEjiBUwglgBI4gVMIJYASOsFKuqulhVP6uqV6vq05seCuB2x56Kq6rOJPnPJP+U5HqSl5I82d0/ucuf2cpzVjkVF6eJU3G90WNJXu3un3f3H5M8k+TD6xwO4DirxOrBJL+86fb1w48BnJgHVnjMUZtkb3g9VFWXklz6qycCOMIqsbqe5KGbbp9L8qvbH9Tdl5NcTrZ3nxWwnFVeBr6U5J1V9Y6qelOSJ5J8c7NjAdzq2C2r7n69qj6e5FtJziR5uruvbnwygJsce+jC/+tJt/RloEMXOE0cugBwCokVMIJYASOIFTCCWAEjiBUwglgBI4gVMIJYASOIFTCCWAEjiBUwglgBI4gVMIJYASOIFTCCWAEjiBUwglgBI4gVMIJYASOscpLTe3b+/Pns7e1t4qkXta1nE0m298w92/xvdr+xZQWMIFbACGIFjCBWwAhiBYwgVsAIYgWMIFbACGIFjCBWwAhiBYwgVsAIYgWMIFbACGIFjCBWwAhiBYwgVsAIYgWMIFbACGIFjCBWwAhiBYwgVsAIx8aqqp6uqhtV9eOTGAjgKKtsWX0pycUNzwFwV8fGqru/k+Q3JzALwB2tbZ9VVV2qqr2q2tvf31/X0wIkWWOsuvtyd+929+7Ozs66nhYgiXcDgSHEChhhlUMXvprke0neVVXXq+pjmx8L4FYPHPeA7n7yJAYBuBsvA4ERxAoYQayAEcQKGEGsgBHEChhBrIARxAoYQayAEcQKGEGsgBHEChhBrIARxAoYQayAEcQKGEGsgBHEChhBrIARxAoYQayAEY49uw3/p7uXHmFjqmrpEbhH2/j1uLu7e8f7bFkBI4gVMIJYASOIFTCCWAEjiBUwglgBI4gVMIJYASOIFTCCWAEjiBUwglgBI4gVMIJYASOIFTCCWAEjiBUwglgBI4gVMIJYASOIFTCCWAEjHBurqnqoqp6vqmtVdbWqnjqJwQButspJTl9P8snufrmq/jbJlap6rrt/suHZAP7i2C2r7v51d798eP33Sa4leXDTgwHc7J72WVXVw0keTfLiEfddqqq9qtrb399fz3QAh1aOVVW9NcnXk3yiu393+/3dfbm7d7t7d2dnZ50zAqwWq6o6m4NQfaW7v7HZkQDeaJV3AyvJF5Nc6+7PbX4kgDdaZcvq8SQfTXKhql45vLx/w3MB3OLYQxe6+7tJ6gRmAbgjR7ADI4gVMIJYASOIFTCCWAEjiBUwglgBI4gVMIJYASOIFTCCWAEjiBUwglgBI4gVMIJYASOIFTCCWAEjiBUwglgBI4gVMIJYASMce3Yb7g/dvfQIG3Fw2ku2gS0rYASxAkYQK2AEsQJGECtgBLECRhArYASxAkYQK2AEsQJGECtgBLECRhArYASxAkYQK2AEsQJGECtgBLECRhArYASxAkYQK2AEsQJGECtghGNjVVVvrqrvV9UPq+pqVX32JAYDuNkqJzn9Q5IL3f1aVZ1N8t2q+vfu/o8NzwbwF8fGqg9O1fva4c2zh5ftPH0vcGqttM+qqs5U1StJbiR5rrtf3OxYALdaKVbd/afufiTJuSSPVdV7bn9MVV2qqr2q2tvf31/3nMB97p7eDezu3yZ5IcnFI+673N273b27s7OzpvEADqzybuBOVb3t8PpbkrwvyU83PRjAzVZ5N/DtSf6tqs7kIG5f6+5nNzsWwK1WeTfwR0kePYFZAO7IEezACGIFjCBWwAhiBYwgVsAIYgWMIFbACGIFjCBWwAhiBYwgVsAIYgWMIFbACGIFjCBWwAhiBYwgVsAIYgWMIFbACGIFjCBWwAirnIoLxurupUfYmKpaeoQTZcsKGEGsgBHEChhBrIARxAoYQayAEcQKGEGsgBHEChhBrIARxAoYQayAEcQKGEGsgBHEChhBrIARxAoYQayAEcQKGEGsgBHEChhBrIARxAoYQayAEVaOVVWdqaofVNWzmxwI4Cj3smX1VJJrmxoE4G5WilVVnUvygSRf2Ow4AEdbdcvq80k+leTPd3pAVV2qqr2q2tvf31/LcAD/69hYVdUHk9zo7it3e1x3X+7u3e7e3dnZWduAAMlqW1aPJ/lQVf0iyTNJLlTVlzc6FcBtjo1Vd3+mu89198NJnkjy7e7+yMYnA7iJ46yAER64lwd39wtJXtjIJAB3YcsKGEGsgBHEChhBrIARxAoYQayAEcQKGEGsgBHEChhBrIARxAoYQayAEcQKGEGsgBHEChhBrIARxAoYQayAEcQKGEGsgBHEChihunv9T1q1n+S/1v7ER/u7JP99Qp/rJFnXPNu6tpNc199395GndN9IrE5SVe119+7Sc6ybdc2zrWs7LevyMhAYQayAEbYhVpeXHmBDrGuebV3bqVjX+H1WwP1hG7asgPuAWAEjiBUwglgBI4gVMML/AHEJC1mbT8wIAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 360x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Solve the transportation problem\n",
    "from scipy.optimize import linprog\n",
    "res = linprog(c=C,A_eq=A[:-1],b_eq=b[:-1],bounds=(0,None),method='simplex')\n",
    "print(res.message)\n",
    "X = res.x.reshape((N,N)) # reshape back to X_ij\n",
    "plt.spy(X)\n",
    "print(X)\n",
    "black_market_estim = 1 - np.diag(X).sum() # do not count the stationary diagonal\n",
    "print('With N=%d the lower bound on black market share is %1.5f'%(N,black_market_estim))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Displacement of the car sales distributions\n",
    "\n",
    "- Main result: significant evidence for a sizable black market share!  \n",
    "- Computed the **lower boundary** on the fraction of illegal trade (why?)  \n",
    "- Grain of salt: this is one of possible mechanisms, need to eliminate other possible routes (see the paper)  \n",
    "- What robustness checks should be run? Technical parameter $ N $ clearly affects the numerical result  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Further learning resources\n",
    "\n",
    "- Full paper\n",
    "  [https://www.jstor.org/stable/j.ctt1q1xs9h](https://www.jstor.org/stable/j.ctt1q1xs9h)  \n",
    "- [PDF for the paper (right-click and Save as…)](_static/pdf/DaljordHuPouliotXiao2019.pdf)  "
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "date": 1612589585.5666718,
  "download_nb": false,
  "filename": "19_black_market.rst",
  "filename_with_path": "19_black_market",
  "kernelspec": {
   "display_name": "Python",
   "language": "python3",
   "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.6"
  },
  "title": "Foundations of Computational Economics #19"
 },
 "nbformat": 4,
 "nbformat_minor": 4
}