{ "cells": [ { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "# Foundations of Computational Economics #19\n", "\n", "by Fedor Iskhakov, ANU\n", "\n", "" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "## Measuring the volume of illegal trade with linear programming\n", "\n", "" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "\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", "" ] }, { "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']" ] }, "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": [ "" ] }, "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": [ "
" ] }, "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": [ "
" ] }, "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": [ "
" ] }, "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 }