“Black Market Performance: Illegal Trade in Beijing License Plates”
by Øystein Daljord, Mandy Hu, Guillaume Pouliot, Junji Xiao
From abstract:
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.
Let $ \mathbb{P}_0 $ be the distribution of car sales prices from pre-lottery time, and $ \mathbb{P}_1 $ the analogous distribution post-lottery.
Under assumptions
the sales distributions should not change from the pre- to the post lottery period, i.e. $ \mathbb{P}_0 = \mathbb{P}_1 $
Data on manufacturer suggested retail prices (MSRP) of the registered vehicles.
import pandas as pd
dt = pd.read_stata('_static/data/beijin_data.dta')
dt.dropna(inplace=True) # drop rows with nan
print('Data has %d observations and %d variables'%tuple(dt.shape)) # print expects tuple
print(dt.head(n=10))
Data has 243677 observations and 3 variables year month MSRP 0 2010 9 153.313139 1 2010 9 44.543519 2 2011 2 88.812069 3 2010 11 210.732564 4 2011 4 101.591900 5 2010 12 56.428979 6 2011 1 140.571004 7 2010 8 170.066283 8 2011 11 111.935614 9 2010 7 191.988099
print(dt['MSRP'].describe())
q99 = dt['MSRP'].quantile(0.99)
dt = dt[dt['MSRP']<q99]
print(dt['MSRP'].describe())
count 243677.000000 mean 166.374743 std 114.990191 min 19.728006 25% 91.780464 50% 134.771217 75% 209.744757 max 1130.669488 Name: MSRP, dtype: float64 count 241240.000000 mean 161.169630 std 102.551081 min 19.728006 25% 91.229416 50% 134.283479 75% 207.808600 max 617.902636 Name: MSRP, dtype: float64
import matplotlib.pyplot as plt
%matplotlib inline
plt.rcParams['figure.figsize'] = [12, 8]
def plot2hist(d1,d2,bins=10,labels=['1','2']):
'''Plots two overlapping histograms'''
plt.hist(d1,bins=bins,density=True,histtype='step',label=labels[0])
plt.hist(d2,bins=bins,density=True,histtype='step',label=labels[1])
plt.legend()
plt.show()
dt10 = dt[dt['year']==2010]['MSRP']
dt11 = dt[dt['year']==2011]['MSRP']
plot2hist(dt10,dt11,labels=['2010','2011'])
Linear programming problem solved by scipy.optimize.linprog
(equality constraints automatically converted)
https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.linprog.html
# Code up the model
#
#
#
#
#
#
#
#
#
#
#
#
#
#
#
#
#
#
#
#
#
#
# Answer below
N = 5 # number of bins to represent distribution
dt['gr'] = pd.qcut(dt.MSRP,q=N,labels=False) # N quantiles
gr10 = dt[dt.year==2010].groupby('gr')
gr11 = dt[dt.year==2011].groupby('gr')
d10 = gr10.MSRP.count()/dt[dt.year==2010].MSRP.count()
d11 = gr11.MSRP.count()/dt[dt.year==2011].MSRP.count()
print(d10,d11,sep='\n\n')
gr 0 0.231451 1 0.211331 2 0.205582 3 0.189273 4 0.162363 Name: MSRP, dtype: float64 gr 0 0.123699 1 0.172511 2 0.186457 3 0.226024 4 0.291310 Name: MSRP, dtype: float64
import numpy as np
# Set up transportation problem
costs = np.ones((N,N)) - np.eye(N) # costs matrix
origins = np.array(d10) # origins
destinations = np.array(d11) # destinations
plt.rcParams['figure.figsize'] = [5, 5]
plt.spy(costs)
<matplotlib.image.AxesImage at 0x7fbd40c5e350>
# convert to linear programming problem
C = costs.reshape(N*N)
A1 = np.kron(np.eye(N),np.ones((1,N))) # sums of x for each origin
A2 = np.kron(np.ones((1,N)),np.eye(N)) # sums of x for each destination
A = np.vstack((A1,A2)) # concatenate vertically
plt.spy(A)
b = np.concatenate((origins,destinations))
# Solve the transportation problem
from scipy.optimize import linprog
res = linprog(c=C,A_eq=A[:-1],b_eq=b[:-1],bounds=(0,None),method='simplex')
print(res.message)
X = res.x.reshape((N,N)) # reshape back to X_ij
plt.spy(X)
print(X)
black_market_estim = 1 - np.diag(X).sum() # do not count the stationary diagonal
print('With N=%d the lower bound on black market share is %1.5f'%(N,black_market_estim))
Optimization terminated successfully. [[0.12369875 0. 0. 0. 0.10775178] [0. 0.17251076 0. 0.01762504 0.02119497] [0. 0. 0.18645705 0.01912521 0. ] [0. 0. 0. 0.18927336 0. ] [0. 0. 0. 0. 0.16236309]] With N=5 the lower bound on black market share is 0.16570