A Practitioner’s Guide to Estimation of Random-Coeffcients Logit Models of Demand
import numpy as np
import pandas as pd
from scipy.optimize import minimize
from numba import jit, njit, prange
from numba import float64 as f64
from numba import int32 as i32
import pyblp
# index names
mid = 'market_ids'
pid = 'product_ids'
# product
prod = pd.read_csv(pyblp.data.NEVO_PRODUCTS_LOCATION)
prod = pd.concat([prod, pd.get_dummies(prod['product_ids'])], axis=1)
prod['con'] = 1
prod['delta0'] = np.log(prod['shares']) - np.log(1-prod.groupby('market_ids')['shares'].transform('sum'))
demand_inst = [mid, pid] + [col for col in prod.columns if 'demand' in col]
inst = pd.concat((prod[demand_inst], pd.get_dummies(prod["product_ids"])), axis=1)
# demographics
demog = pd.read_csv(pyblp.data.NEVO_AGENTS_LOCATION)
# column names
cn_s = [pid, 'shares'] # market share
cn_x1 = [pid, 'prices'] + prod['product_ids'].unique().tolist() # product characteristics (linear)
cn_x2 = [pid, 'prices', 'con', 'sugar', 'mushy'] # product characteristics (nonlinear)
cn_d = ['income', 'income_squared', 'age', 'child'] # demographics
cn_i = list(inst.columns)[1:]
cn_delta = [pid, 'delta0']
# tensors
S = np.array([prod[prod[mid]==m][cn_s].sort_values(pid)[cn_s[1:]].values for m in prod[mid].unique()])
log_S = np.log(S)
X1 = np.array([prod[prod[mid]==m][cn_x1].sort_values(pid)[cn_x1[1:]].values for m in prod[mid].unique()])
X2 = np.array([prod[prod[mid]==m][cn_x2].sort_values(pid)[cn_x2[1:]].values for m in prod[mid].unique()])
D = np.array([demog[demog[mid]==m][cn_d].values for m in demog[mid].unique()])
Z = np.array([inst[inst[mid]==m][cn_i].sort_values(pid)[cn_i[1:]].values for m in inst[mid].unique()])
delta0 = np.array([prod[prod[mid]==m][cn_delta].sort_values(pid)[cn_delta[1:]].values for m in prod[mid].unique()])
Nu = np.random.standard_normal((X1.shape[0], D.shape[1], X2.shape[-1]))
S.shape, delta0.shape, X1.shape, X2.shape, Nu.shape, D.shape, Z.shape
def gen_theta2(Nx2, Nd):
theta2 = np.zeros(int(Nx2*(Nx2+1)/2+Nd*Nx2))+1e-5
diag_idx, j = [], 0
for i in range(Nx2):
j += i + 1
diag_idx.append(j-1)
bounds = [(0, None) if i in diag_idx else (None, None) for i in range(len(theta2))]
return theta2, bounds
def get_L_Pi(theta2, Nx2, Nd):
Nl = int(Nx2*(Nx2+1)/2)
L = np.zeros((Nx2, Nx2))
L[np.tril_indices(Nx2)] = theta2[:Nl]
Pi = theta2[Nl:].reshape((Nx2, Nd))
return L, Pi
@jit(f64[:, :](f64[:, :], f64[:, :], i32, i32), nopython=True)
def m_share(delta, mu, Nc, Np):
u = delta + mu
s = np.zeros(Np)
for i in range(Nc):
u_max = max(0, np.max(u[:, i]))
u0 = np.exp(-u_max)
u[:, i] = np.exp(u[:, i]-u_max)
s = s + u[:, i] / (u0+u[:, i].sum())
s = s/Nc
for i in range(Np):
if s[i] == 0:
s[i] = 1e-99
return np.log(s).reshape(-1, 1)
@jit(f64[:, :](f64[:, :], f64[:, :], f64[:, :], i32, i32, f64), nopython=True)
def fp_iteration(log_S, delta0, mu, Nc, Np, epsilon):
err = 1
while err > epsilon:
log_s0 = m_share(delta0, mu, Nc, Np)
delta1 = delta0 + log_S - log_s0
err = np.abs(delta1-delta0).max()
delta0 = delta1
return delta0
@njit(parallel=True)
def solve_delta(log_S, delta0, mu, Nt, Nc, Np, epsilon):
delta = np.zeros(delta0.shape)
for i in prange(Nt):
delta[i] = fp_iteration(log_S[i], delta0[i], mu[i], Nc, Np, epsilon)
return delta
def gmm(X1, Z, delta):
# 2SLS
W0 = np.linalg.inv(Z.T @ Z)
theta1_2sls = np.linalg.inv(X1.T @ Z @ W0 @ Z.T @ X1) @ X1.T @ Z @ W0 @ Z.T @ delta
ksi_2sls = delta - X1 @ theta1_2sls
# GMM
W = np.linalg.inv((Z * ksi_2sls).T @ (Z * ksi_2sls))
theta1_gmm = np.linalg.inv(X1.T @ Z @ W @ Z.T @ X1) @ X1.T @ Z @ W @ Z.T @ delta
ksi_gmm = delta - X1 @ theta1_gmm
obj_gmm = ksi_gmm.T @ Z @ W @ Z.T @ ksi_gmm
return theta1_gmm, obj_gmm.item()
def search_theta2(theta2, X1, X2, Z, D, Nu, log_S, epsilon, report=1000):
# "hot start"
global delta, theta1
L, Pi = get_L_Pi(theta2, X2.shape[2], D.shape[2])
mu = X2 @ (D @ Pi.T + Nu @ L.T).transpose((0, 2, 1))
delta = solve_delta(log_S, delta, mu, X1.shape[0], D.shape[1], X1.shape[1], epsilon)
theta1, obj_gmm = gmm(np.vstack(X1), np.vstack(Z), np.vstack(delta))
# report # of outer loops that has been excuted
global O
O += 1
if O % report == 0:
print(f'# out: {O}; Loss: {obj_gmm}')
return obj_gmm
%%time
O = 0
theta2_init, bounds = gen_theta2(X2.shape[2], D.shape[2])
delta, theta1 = delta0.copy(), None
optimal_params = minimize(search_theta2,
theta2_init,
args=(X1, X2, Z, D, Nu, log_S, 1e-12, 2000),
method='L-BFGS-B',
#bounds=bounds,
options={'maxiter': 1e20})
optimal_params.x, theta1
get_L_Pi(optimal_params.x, X2.shape[2], D.shape[2])