© 2018-2023 by Alfred Galichon. Past and present support from NSF grant DMS-1716489, ERC grant CoG-866274 are acknowledged, as well as inputs from contributors listed here.
If you reuse material from this masterclass, please cite as:
Alfred Galichon, 'math+econ+code' masterclass series. https://github.com/math-econ-code/mec_equil
See slides D1a_getting-started.pdf.
For each $z \in \{0,...,n-1\}$, we define $h_z$ and $v_z$ the coordinates (horizontal and vertical) of $z$.
import numpy as np
import pandas as pd
import scipy.optimize as opt
from scipy.spatial.distance import cdist
from time import time
def vec(X):
return np.reshape(X,(-1, 1))
np.random.seed(777)
nbz = 3
z_df=pd.DataFrame ({'h':np.random.rand(nbz),'v':np.random.rand(nbz)}, columns = ['h','v'])
For each driver $i$, we shall generate the position (horizontal and vertical coordinates), the concavity of preferences $\tau_i$, and $\lambda_i$ the value of $i$'s time.
nbi = 200
np.random.seed(778)
i_df=pd.DataFrame ({'h':np.random.rand(nbi),
'v':np.random.rand(nbi),
'tau':np.random.rand(nbi),
'lambda':10 * np.random.rand(nbi) },
columns = ['h','v','tau','lambda'])
For each passenger $j$, we shall generate the position (horizontal and vertical coordinates), elasticity of price-time substitution $\sigma_j$, and the valuation of $j$'s time $\epsilon_j$.
nbj = 500
np.random.seed(779)
j_df=pd.DataFrame ({'h':np.random.rand(nbj),
'v':np.random.rand(nbj),
'sigma':1/np.random.rand(nbj),
'epsilon':20*np.random.rand(nbj),
'eta':np.random.rand(nbj) },
columns = ['h','v','sigma','epsilon','eta'])
The average speed walking is $4$ km/h and driving is $25$ km/h. Let's compute the time $T_{iz}$ that it takes driver each driver $i$ to drive to pickup at each $z$.
avge_speed_drive = 25
avge_speed_walk = 4
T_iz=cdist(i_df[['h','v']],z_df[['h','v']] ) / avge_speed_drive
T_jz=cdist(j_df[['h','v']],z_df[['h','v']] ) / avge_speed_walk
Now let's compute the supply. The utility that each driver $i$ have for each option $z$ is $$u_{iz}\left( p_z,\varepsilon \right) =p_{z}^{1-\tau _{i}}-\lambda _{i}T_{iz},$$ and the option that driver $i$ has for the exit option is $0$. We have $$s_z(p_z) = \sum_i 1{\{ z = argmax_{z^\prime} (u_{iz^\prime}(p_z)) \}}.$$
def s_z(p_z):
u_iz = (p_z[:,np.newaxis]**(1- i_df[['tau']].values.reshape(-1)) - T_iz.transpose()*i_df[['lambda']].values.reshape(-1)).transpose()
u_i0 = np.zeros((nbi,1))
return np.asarray([np.count_nonzero(np.argmax(np.concatenate((u_iz,u_i0),axis=1), axis=1) == z) for z in range(nbz)] )
We now compute the demand. The cost that each passenger $j$ has for each option $z$, which is given by $$c_{jz}\left( p_{z}\right) =\left( p_{z}^{1-\frac{1}{% \sigma _{j}}}+\left( \epsilon _{j}T_{jz}\right) ^{1-\frac{1}{\sigma _{j}}% }\right) ^{\frac{\sigma _{j}}{\sigma _{j}-1}},$$ and the cost associated with the exit option is $\eta_j$. Similar to above, we have $$d_z(p_z) = \sum_j 1{\{ z = argmin_{z^\prime} (v_{jz^\prime}(p_z)) \}}.$$
def d_z(p_z):
expos_j = 1-1/j_df[['sigma']].values.reshape(-1)
c_jz = ((p_z[:,np.newaxis]**expos_j + ( T_jz.transpose()*j_df[['epsilon']].values.reshape(-1) )**expos_j )**(1/expos_j)).transpose()
c_j0 = j_df[['eta']].values.reshape((-1,1))
return np.asarray([np.count_nonzero(np.argmin(np.concatenate((c_jz,c_j0),axis=1), axis=1) == z) for z in range(nbz)])
def Q_z(p_z):
return s_z(p_z)-d_z(p_z)
p_z=np.random.rand(nbz)
Q_z(p_z)
array([199, -31, -7])
Actually, we would like to have a smooth approximation of supply and demand. Thus set $$s^{smooth}_z(p_z,T) = \sum_i \frac {e^{\frac {u_{iz}} {T} }} {\sum_{z^\prime} e^{\frac {c_{iz^\prime}} {T} }},$$ and $$d^{smooth}_z(p_z,T) = \sum_j \frac {e^{\frac {-c_{jz}} {T} }} {\sum_{z^\prime} e^{\frac {-c_{jz^\prime}} {T} }}.$$ Note (math exercise!) that as $T \to 0$, we tend to the previous functions $s$ and $d$.
def ssmooth_z(p_z,T):
u_iz = (p_z[:,np.newaxis]**(1- i_df[['tau']].values.reshape(-1)) - T_iz.transpose()*i_df[['lambda']].values.reshape(-1)).transpose()
u_i0 = np.zeros((nbi,1))
max_i = np.max(np.concatenate((u_iz,u_i0),axis=1),axis=1)
utilde_i = max_i - T*np.log ((np.exp( (u_i0.reshape(-1) - max_i)/T) ).reshape(-1) + np.sum(np.exp( ( ( u_iz.transpose() - max_i).transpose())/T) ,axis=1))
return (np.sum( np.exp( ( ( u_iz.transpose() - utilde_i).transpose() )/T) ,axis = 0) )
def dsmooth_z(p_z,T):
expos_j = 1-1/j_df[['sigma']].values.reshape(-1)
c_jz = ((p_z[:,np.newaxis]**expos_j + ( T_jz.transpose()*j_df[['epsilon']].values.reshape(-1) )**expos_j )**(1/expos_j)).transpose()
c_j0 = j_df[['eta']].values.reshape((-1,1))
min_j = np.min(np.concatenate((c_jz,c_j0),axis=1),axis=1)
ctilde_j = min_j-T*np.log ((np.exp( (- c_j0.reshape(-1) + min_j)/T) ).reshape(-1) + np.sum(np.exp( ( (- c_jz.transpose() + min_j).transpose())/T) ,axis=1))
return (np.sum( np.exp( ( ( - c_jz.transpose() + ctilde_j).transpose() )/T) ,axis = 0) )
def Qsmooth_z(p_z,T):
return (ssmooth_z(p_z,T)-dsmooth_z(p_z,T))
EquilibriumProblem class¶We are now in a good place to create a Python class, with all the data that will be useful.
class EquilibriumProblem:
def __init__(self,nbz=0,q_z=np.array([]),p0_z=np.array([]),pmin=-1e10,pmax = 1e10):
if q_z.size == 0:
if nbz != 0:
q_z = np.zeros(nbz)
else:
raise Exception("Either nbz or q_z needs to be provided.")
elif nbz == 0:
nb_z = q_z.size
if p0_z.size == 0:
p0_z = np.zeros(nbz)
self.nbz = nbz
self.q_z = q_z
self.p0_z = p0_z
self.Q_z = Q_z
self.pmin = pmin
self.pmax = pmax
self.eq_nbsteps = []
self.eq_p_z = []
self.eq_diffp_z = []
self.eq_delta_z = []
self.eq_code = []
def Q_z(self,p_z):
return Q_z(p_z)
Q_z = lambda p_z: Qsmooth_z(p_z,0.001)
EqPb = EquilibriumProblem(nbz,pmin=0)
p_z = np.zeros(nbz)
EqPb.Q_z(p_z)
array([-83.12771217, -74.08853699, -15.33273283])
We now would like to construct a toolbox to solve the equilibrium problem $$Q(p)=q.$$ We will see two related basic methods, the Gauss-Seidel algorithm and the Jacobi algorithm. Both of them are based ont the notion of coordinate updating.
To start with, we need to introduce a coordinate update function, denoted $cu^z(p),$ such that $$Q_z(cu^z_z(p),p_{-z})=q_z\\ cu^z_{-z}(p)=p_{-z} $$ so that $cu^z_z$ finds the equilibrium price in the market for $z$, provided all the other are markets are at equilibrium.
def replace_entry(ind,vect,value):
np.put(vect,ind,value)
return vect
def cu(self,z,p_z):
oldpz = p_z[z]
newpz = opt.brentq(lambda thep : self.Q_z(replace_entry(z,p_z,thep))[z]-self.q_z[z],self.pmin,self.pmax)
p_z[z] = oldpz
return newpz, (newpz-oldpz)
Above, self is intended to be an object of EquilibriumProblem class. We can just add the method cu to that class.
EquilibriumProblem.cu = cu
# Try it:
p_z = np.zeros(nbz)
p_z[0],_=EqPb.cu(0,p_z)
print(p_z)
p_z[1],_=EqPb.cu(1,p_z)
print(p_z)
p_z[2],_=EqPb.cu(2,p_z)
print(p_z)
[0.00016693 0. 0. ] [0.00016693 0.00027716 0. ] [0.00016693 0.00027716 0.00028675]
p_z = np.zeros(nbz)
pp_z = np.zeros(nbz)
pp_z[0],_=EqPb.cu(0,p_z)
print(pp_z)
pp_z = np.zeros(nbz)
pp_z[1],_=EqPb.cu(1,p_z)
print(pp_z)
pp_z = np.zeros(nbz)
pp_z[2],_=EqPb.cu(2,p_z)
print(pp_z)
[0.00016693 0. 0. ] [0.00000000e+00 1.73776256e-05 0.00000000e+00] [0.000000e+00 0.000000e+00 1.892368e-12]
The Gauss-Seidel algorithm consists of setting $p^t$ such that
$ p^1: Q_1(p^1_1,p^0_{-1})=q_1 , p^1_{-1} = p^0_{-1} \\ p^2: Q_2(p^2_2,p^1_{-2})=q_2, p^2_{-2} = p^1_{-2} \\ ... \\ p^Z: Q_z (p^Z_Z,p^{Z-1}_{-Z})=q_Z, p^Z_{-Z} = p^{Z-1}_{-Z} \\ p^{Z+1}: Q_{1} (p^{Z+1}_1,p^{Z}_{-1}) = q_1, p^{Z+1}_{-1} = p^{Z+1}_{-1} \\ ... $
In other words, one full iteration of the Gauss-Seidel algorithm amounts to iterating the map
$f^{GS}(p)=(cu^1(p),cu^2\circ cu^1(p),...,cu^Z\circ cu^{Z-1}\circ ... \circ cu^1(p)),$
which we now implement.
def fGS_z(self,p_z):
pp_z = np.copy(p_z)
diffp_z = np.zeros(self.nbz)
for z in range(nbz):
pp_z[z],diffp_z[z] = self.cu(z,pp_z)
return(pp_z,diffp_z)
The Jacobi algorithm consists of setting $p^t$ such that
$ Q_1(p^1_1,p^0_{-1}) = q_1 \\ Q_2(p^1_2,p^0_{-2}) = q_2 \\ ... \\ Q_z (p^1_Z,p^0_{-Z}) = q_Z \\ Q_{1} (p^2_1,p^1_{-1}) = q_1\\ ... $
In other words, one full iteration of the Gauss-Seidel algorithm amounts to
$f^{J}(p)=(cu^1(p),cu^2(p),...,cu^Z(p)).$
Let's implement $f^J$:
def fJ_z(self,p_z):
pp_z = np.zeros(self.nbz)
diffp_z = np.zeros(self.nbz)
for z in range(nbz):
pp_z[z],diffp_z[z]=self.cu(z,p_z)
return(pp_z,diffp_z)
The first five steps of the Gauss-Seidel algorithm will therefore be:
p_z=np.zeros(nbz)
for i in range(5):
p_z,_ = fGS_z(EqPb,p_z)
print(p_z)
EqPb.Q_z(p_z) - EqPb.q_z
[0.00016693 0.00027716 0.00028675] [0.0005505 0.00088081 0.00090476] [0.00122043 0.00164713 0.00162068] [0.00204997 0.00205102 0.00202434] [0.00239362 0.00240875 0.00237726]
array([-6.97605609e+00, -1.45690873e+00, 4.84945417e-13])
The first five steps of the Jacobi algorithm will therefore be:
p_z=np.zeros(nbz)
for i in range(5):
p_z,_ = fJ_z(EqPb,p_z)
print(p_z)
EqPb.Q_z(p_z) - EqPb.q_z
[1.66926158e-04 1.73776256e-05 1.89236800e-12] [0.00024788 0.00027722 0.00014836] [0.00050723 0.00047984 0.00028993] [0.00074001 0.00082783 0.00051344] [0.00113698 0.00115414 0.00085994]
array([-14.92463311, -7.67069782, -9.16727684])
Compared with Gauss-Seidel above, Jacobi does not make use of all the relevant information (i.e. latest price update) at a given point in time. But its structure makes it more naturally suited for parallelization -- we will return to that.
Let's provide a solve method to EquilibriumProblem which computes both Gauss-Seidel and Jacobi.
def solve(self,maxit = 1000, method = 'GS', valtol=1e-5, steptol=1e-9,output=0):
start_time = time()
code = 0
delta_z = np.zeros(self.nbz)
p_z = np.copy(self.p0_z)
if method=='GS':
method_name = 'Gauss-Seidel'
elif method=='J':
method_name = 'Jacobi'
else:
raise Exception('Method code '+method+' is not available.')
for i in range(maxit):
if method=='GS':
p_z,diffp_z = self.fGS_z(p_z)
elif method=='J':
p_z,diffp_z = self.fJ_z(p_z)
else:
raise Exception('Method code '+method+' is not available.')
delta_z = self.Q_z(p_z) - self.q_z
if output > 1 :
print("p=",p_z)
if np.max(np.abs(delta_z) ) < valtol :
code = 0
break
if np.max(np.abs(diffp_z) ) < steptol :
code = 1
break
code = 2
comp_time = time() - start_time
if output > 0 :
print(method_name, 'method converged in', i, 'iterations.')
print('Value of p =',p_z)
print('Value of p\'_z - p_z =',diffp_z)
print('Value of e(p)-q =',delta_z)
print('Time elapsed =',comp_time,' s.')
print('Code =',code)
self.eq_p_z = p_z
self.eq_diffp_z = diffp_z
self.eq_delta_z = delta_z
self.comp_code = code
self.comp_nbsteps = i
self.comp_time = comp_time
return code
EquilibriumProblem.fGS_z = fGS_z
EquilibriumProblem.fJ_z = fJ_z
EquilibriumProblem.solve = solve
code = EqPb.solve(valtol=1e-5,steptol = 1e-12,method='GS',output=1)
Gauss-Seidel method converged in 26 iterations. Value of p = [0.00404359 0.00398981 0.00398592] Value of p'_z - p_z = [-4.25764874e-10 -3.68276937e-11 -6.54447174e-10] Value of e(p)-q = [-3.54566518e-06 9.11272851e-06 1.52858348e-09] Time elapsed = 11.164602994918823 s. Code = 0
code = EqPb.solve(valtol=1e-5,steptol = 1e-12,method='J',output=1)
Jacobi method converged in 64 iterations. Value of p = [0.00404359 0.00398981 0.00398593] Value of p'_z - p_z = [-4.68941770e-10 7.31033839e-10 1.58366912e-09] Value of e(p)-q = [3.16996350e-06 3.27073918e-06 9.37327900e-06] Time elapsed = 47.91896080970764 s. Code = 0
See the basics in slides D1.1_mathematical-results.pdf.
More on convergence will require Perron-Froebenius theory, which will be discussed in Advanced Lecture AL1.
Parallel implementation on the Cloud will be addressed in Advanced Lecture AL5.