Author: Conghan Zheng
Last updated:20 DEC 2024
Compiled with: Python 3.11.11
This example performs a structural demand estimation using the algorithm described in Berry, Levinsohn & Pakes (1995).
Suppose there are M markets indexed by $m = 1,\dots, M$, and let $J_m$ be the number of options available to each consumer in market $m$. Following Nevo (2000), the (indirect) choice-specific utility that consumer $i$ in market $m$ receives from product $j$ (where $j=0$ represents the "no purchase" option, i.e. the external good) is given by
\begin{align*} \tilde{u}_{ijm} = \underbrace{x_{jm} \beta_i + \alpha_i (y_i - p_{jm})}_{\substack{\text{deterministic component}}} + \underbrace{\xi_{jm} + \varepsilon_{ijm}}_{\text{random component}} \end{align*}where:
Note that the utility from income, $\alpha_i y_i$, plays no role in the consumer's choice, so we can actually just drop it or transform the price $p_{jm}$ into $\frac{p_{jm}}{y_i}$. From now on, we consider the following monotonic transformation of the original utility specification: \begin{align} u_{ijm} = \underbrace{x_{jm} \beta_i - \alpha_i p_{jm} }_{ \equiv V_{ijm}} + \xi_{jm} + \varepsilon_{ijm} \tag{M.1} \end{align}
The key issue motivating the BLP estimation approach is the endogeneity of price. The price of any product depends on all its attributes. Some attributes are observed by the producer but not measured by the researcher and are included in $\xi$, but they affect the demand and/or the cost of the product. Therefore, the price $p_{jm}$ depends on $\xi_{jm}$.
The BLP approach to this problem is to move $\xi_{jm}$ to the observed part of utility. We split the deterministic component $V_{ijm}$ into two parts: $\bar{V}_{jm}$, the part that varies across products and markets but not consumers, and $\tilde{V}_{ijm}$, the part that varies across consumers as well as markets and products. So we can write \begin{align} (M.1): \quad u_{ijm} &= V_{ijm} + \xi_{jm} + \varepsilon_{ijm} \nonumber \\ &\equiv \bar{V}_{jm} + \tilde{V}_{ijm} + \xi_{jm} + \varepsilon_{ijm} \nonumber \\ &= \underbrace{\bar{V}_{jm} + \xi_{jm}}_{\equiv \delta_{jm}} +\tilde{V}_{ijm} + \varepsilon_{ijm} \nonumber \\ &\equiv \delta_{jm} + \tilde{V}_{ijm} + \varepsilon_{ijm} \tag{M.2} \end{align} where $\delta_{jm}$ is constant for each product in each market. A choice model based on this utility specification has no endogeneity, a constant included for each market-product level absorbs $x_{jm}$.
One difficulty in decomposing $V_{ijm}$ into $\bar{V}_{jm}$ and $\tilde{V}_{ijm}$ is that all terms in the indirect utility (M.1) except for $\xi_{jm}$ has a subscript $i$ (think of $\alpha_i p_{jm}$, $x_{jm} \beta_i$, and $\varepsilon_{ijm}$). To describe how consumer preferences vary as a function of the individual characteristics, we model the distribution of consumer taste parameters $\alpha_i$ and $\beta_i$. \begin{align} (\alpha_i,\beta_i)' &= (\alpha,\beta)' + \Pi \cdot \eta_i + \Sigma \cdot \nu_i , \quad \theta \equiv (\Pi, \Sigma) \tag{M.3} \end{align} where $\alpha$ and $\beta$ are the average levels of the two parameters, $\eta_i$ and $\nu_i$ are the individually varying part, and the second equality depends on $\Sigma$ being diagonal. According to Nevo (2000), given that no individual data is observed, the individual characteristics consist of two components: the demographics $\eta_i$ which is referred to as observed, and the additional characteristics $\nu_i$ which are referred to as unobserved.
Even though we do not observe individual data, demographic variables $\eta_i$ are still considered as "observed". We may not know the individual income of each agent, but we may know some feature of the income distribution. Finally, we might have a sample from the joint distribution of several demographic variables (e.g., income, age, family size, race, and education, etc. from census data).
The additional characteristics $\nu_i$ could include things like whether the individual owns a dog, which could be very important in the decision of which car to buy (suppose we are considering the automobile markets as in the BLP paper), but are unlikely to be collected by any data project.
Using (M.2) and (M.3), we can rewrite the original model (M.1) as \begin{align} u_{ijm} &= \delta_{jm} + \tilde{V}_{ijm}(\theta) + \varepsilon_{ijm} \nonumber \\ &= \delta_{jm} + (x_{jm},p_{jm}) \cdot (\Pi \eta_i + \Sigma \nu_i) + \varepsilon_{ijm} \tag{M.4} \\ \delta_{jm} &= \bar{V}_{jm}(\alpha,\beta) + \xi_{jm} \nonumber \\ &= x_{jm} \beta - \alpha p_{jm} + \xi_{jm} \tag{M.5} \end{align} $(\alpha,\beta)$ and $\theta$ are all the parameters of the model, usually vector $(\alpha, \beta)$ is called the linear parameters, and vector $\theta$ is called the nonlinear parameters (the reasons for the names will be discussed below).
The model (M.4) and (M.5) tells us that $(\alpha,\beta)$ (the linear parameters) and $\theta$ (the nonlinear parameters) are all the parameters we need to estimate in this model.
We proceed with the estimation of the choice model with the utility specification (M.4). Recall the additive random utility interpretation of multinomial models ($U_{j} = x_j \beta + \varepsilon_j$). We get a similar form in the current case, except that the parameters $(\eta_i, \nu_i)$ have randomness. A multi-choice logit model where the parameter $\beta$ in the deterministic component is allowed to be random is called a Random Parameters Logit (RPL, also called a Mixed Logit). RPL model is suitable for modelling the aggregation across heterogeneous consumers (in our context, random coefficients actually mean individual coefficients).
Given that the unobservables $\{\varepsilon_{ijm}\}$ are independent and type I extreme value distributed, the individual choice probability is \begin{align*} S_{ijm} &= \mathbf{P}(u_{ijm} \geq u_{ikm}, \forall k = 1,\dots,M) \\ &= \mathbf{P}\bigg(\varepsilon_{ijm} - \varepsilon_{ikm} \leq \left(\delta_{jm} +\tilde{V}_{ijm}(\theta)\right) - \left( \delta_{km} + \tilde{V}_{ikm}(\theta) \right), \forall k = 1,\dots,M \bigg) \\ &= \int_{\varepsilon_{ijm}: \, u_{ijm} \geq u_{ikm}, \forall k} \mathrm{d} F(\varepsilon_{ijm}) \\ &= \frac{e^{\delta_{jm} + \tilde{V}_{ijm}(\theta)}}{\sum_l e^{\delta_{lm} + \tilde{V}_{ijm}(\theta)}} \end{align*} We can integrate out the randomness (in our context, "individualness") from $\eta_i$ and $\nu_i$ analytically, and get the market share of product $j$ in market $m$: \begin{align} S_{jm} &\overset{(M.3)}{=} \int_{\nu_i} \int_{\eta_i} S_{ijm} \mathrm{d} F(\eta_i, \nu_i \vert \theta) \nonumber \\ &= \int_{\nu_i} \int_{\eta_i} \frac{e^{\delta_{jm} + \tilde{V}_{ijm}(\theta)}}{\sum_l e^{\delta_{lm} + \tilde{V}_{ilm}(\theta)}} \mathrm{d} F(\eta_i, \nu_i \vert \theta) \tag{M.6} \end{align}
Note that the integral in (M.6) is multidimensional. BLP assumed that the distribution $F(\cdot)$ is the product of $k+1$ independent normals (1 from $\alpha_i$, $k$ from $\beta_{i}$). The above integrals are typically evaluated by Monte Carlo simulation with $N_m$ (number of sampled consumers in market $m$) draws of $(\eta_i, \nu_i)$ from $F(\eta_i, \nu_i \vert \theta)$: \begin{align*} \hat{S}_{jm} &= \frac{1}{N_m} \sum_{i=1}^{N_m} \frac{e^{\delta_{jm} + \tilde{V}_{ijm}(\theta)}}{\sum_l e^{\delta_{lm} + \tilde{V}_{ilm}(\theta)}} \tag{M.7} \end{align*} In principle, other numerical integration methods could also be used.
Once the choice model is estimated, the estimated constants $\delta$ can be used as the dependent variable in a linear regression (M.5) to estimate $(\alpha,\beta)$. Since price is endogenous in this regression, it should be estimated by IV. BLP suggests two instruments: the average non-price attributes of other products of the same manufacturer; and the average non-price attributes of other firms' products.
But before we do that, we need to address the problem that there may be a very large number of constants $\delta_{jm}$ to estimate: $J \times M \times T$ (number of products times number of markets times number of years).
BLP provides an algorithm to estimate them quickly, within the iterative process for the other parameters. We already know that $\hat{S}_{jm}$ is a function of $\delta_{jm}$ and $\theta$ from (M.7). For a correctly specified model and a given value of $\theta$, the constants $\delta_{jm}$ determine the predicted market shares $\hat{S}_{jm}$ for each product, and can therefore be set so that the predicted shares equal actual shares: $\hat{S}_{jm}(\delta_{jm}) = S_{jm}$, at each trial value of $\theta$.
Instead of estimating the constants $\delta_{jm}$ by the usual gradient-based methods, the iterative procedure recalibrates the constants so that $\hat{S}_{jm}(\delta_{jm}) = S_{jm}$. The constants are iteratively adjusted by \begin{align} \delta_{jm}^{t+1} = \delta_{jm}^{t} + \ln \left( \frac{S_{jm}}{\hat{S}_{jm}(\delta^{t}_{jm})} \right) \tag{M.8} \end{align} Think about how the process (M.8) moves each constant in the "right" direction.
BLP shows that the iterative adjustment process is a contraction mapping, so it is guaranteed to converge to a unique set of constants $\delta$. For any given value of $\theta$, the calibrated constants $\delta$s are uniquely determined (by the iterative adjustment process), that is, $\hat{\delta} = \hat{\delta}(\theta)$. This implies that the choice probability $S_{ijm}$ is now a function of $\theta$ alone: \begin{align*} \hat{S}_{ijm}(\theta) &= \hat{S}_{ijm}(\delta(\theta), \theta) \end{align*} The optimization of the whole problem is only about $\theta$.
For any given $\theta$, the inner loop of the NFP (M.8) solves the share equations $\hat{S}_{jm} = S_{jm}$ until the successive iterates $\delta_{jm}^{t+1}$ and $\delta_{jm}^{t}$ are sufficiently close. \begin{align*} \quad & \hat{\delta}_{jm} (\theta) = \mathop{\mathrm{argmin}}\limits_{\delta_{jm}} \enspace \lVert \delta_{jm}^{t} - \delta_{jm}^{t+1} \rVert \end{align*}
In each iteration of the outer loop, (M.5) is solved by GMM. For the given $\theta$, we solve \begin{align} \left(\hat{\alpha}(\theta), \hat{\beta}(\theta) \right) = \mathop{\mathrm{argmin}}\limits_{(\alpha, \beta)} \enspace g' W^{-1} g \tag{M.9} \end{align} and we loop over all possible $\theta$ to minimize the GMM objective function.
The observation specific moments $g_{jm}$ in (M.9) are \begin{align} g_{jm} &= \xi_{jm} \cdot z_{jm} \overset{(M.5)}{=} \big[\delta_{jm}(\theta) - x_{jm} \beta + \alpha p_{jm}\big] \cdot z_{jm} \tag{M.10} \end{align} The optimal weight matrix $W$ is the asymptotic covariance of $g$: $W = g g'$.
Python: (best one to choose for now)
R:
Stata:
Decide the numbers of draws $N_m$, arbitrary starting values $\theta^{0}$ and $\delta^0$. Draw random values for $(\eta_i, \nu_i)$ across $i = 1,\dots, N_m$, from i.i.d multivariate normal (of length $k+1$) $F(\eta_i, \nu_i \vert \theta^{0})$. Approximate the integral $\hat{S}_{jm}$ using (M.7).
Use starting values $\theta^0$ and $\delta^0$ and the observed market shares $S$, solve $\hat{\delta}$ by the iterative provess (M.8).
Select arbitrary starting values $(\alpha^0, \beta^0)$, calculate the moment conditions $g(\alpha^0, \beta^0, \hat{\delta})$ using (M.10) and the starting weighting matrix $W^0 = g(\alpha^0, \beta^0, \hat{\delta}) g(\alpha^0, \beta^0, \hat{\delta})'$. Solve for the GMM estimates for $\alpha$ and $\beta$. Update the weighting matrix using the estimates, iterate until the GMM estimates converges. Take the converged weighting matrix and estimates as $\hat{W}$ and $(\hat{\alpha}, \hat{\beta})$.
Take the new values $(\hat{\delta}, \hat{\alpha}, \hat{\beta}, \hat{W})$, find a new value for $\theta$, repeat steps (1) to (3), until the moment expression $g' W^{-1} g$ from (M.9) is close enough to 0.
The estimates are consistent even if you don't iterate on the weighting matrix in step 3. While an optimal 2-step GMM or an efficient 1-step continuously updated GMM will improve the statistical efficiency (Nevo, 2000; Dubé et al., 2012).
I am trying to replicate Table 1 of Nevo (2000). For illustration purposes, I don't use the pyblp package.
"""
Preliminaries
"""
# !pip install pandas pyblp numba
import pandas as pd
import numpy as np
from scipy.optimize import minimize
import pyblp
from numba import jit, njit, prange
import time
import warnings
warnings.filterwarnings("ignore", category=Warning) ## Supress all warnings
"""
BLP Estimation
- x1: non-price product attributes (linear)
- x2: non-price product attributes (nonlinear)
- p: price
- v: monte carlo draws for nu
- demographic: consumer demographics
- J: vector of number of goods per market
- M: number of markets
- N: number of draws
- delta: mean utility
- sigma: Sigma for nonlinear parameters
- pi: Pi for nonlinear parameters
"""
class DeltaManager:
"""
Manages delta values throughout the estimation process
"""
def __init__(self, delta):
self.delta = np.array(delta)
class BLPEstimator:
def __init__(self):
self.time_start = time.time()
self._initialize_data()
self._setup_model()
def _initialize_data(self):
"""
Read in data
"""
## Product data
self.product = pd.read_csv(pyblp.data.NEVO_PRODUCTS_LOCATION)
self.product["cons"] = 1
## Consumer data
self.consumer = pd.read_csv(pyblp.data.NEVO_AGENTS_LOCATION)
self.demographic = self.consumer[["income", "income_squared", "age", "child"]].values
"""
Some data cleansing
"""
## Market structure
self.J = self.product.groupby("market_ids").sum(numeric_only=True).cons.values ## Number of goods per market
self.M = len(self.J) ## Number of markets
self.N = 20 ## Number of simulations per market (Accorigin to Nevo, 2000, page 538)
def _setup_model(self):
"""
We specify the regressors in $\bar{V}$ (linear) and $\tilde{V}$ (non-linear). The two sets of regressors are not necessarily the same, because we may not want random coefficients (the $\beta_i$s) over all the regressors.
"""
## The (linear) regressors in V_bar
self.X1 = np.hstack((self.product[["prices"]],
pd.get_dummies(self.product["product_ids"])))
## The (non-linear) regressors in V_tilde, over which we want random coefficients (so no fixed effects)
self.X2 = self.product[["cons", "prices", "sugar", "mushy"]].values
## Number of nonlinear parameters (k=3 because we have "prices", "sugar", "mushy")
self.k = self.X2.shape[1] - 1
## Insruments
iv_cols = [col for col in self.product.columns if 'demand' in col]
self.Z = np.hstack((self.product[iv_cols],
pd.get_dummies(self.product["product_ids"])))
"""
Monte carlo simulation for nu{i}
- k+1: k from x, 1 from p
- M: different draws for each market
- N: number of draws
"""
np.random.seed(0)
self.v = np.reshape(np.random.standard_normal((self.k + 1) * self.N * self.M), (self.M * self.N, self.k + 1))
"""
Delta
"""
## Shares of the outside good (base category) (See Nevo, 2000, page 520)
self.product["outside"] = self.product["shares"].groupby(
self.product["market_ids"]).transform("sum")
## Initial delta: log(share) - log(share(outside good)) (normalization)
delta_0 = np.log(self.product["shares"]) - np.log(self.product["outside"])
## Initialize a delta object using delta_0
self.d = DeltaManager(delta_0)
@staticmethod
@njit(parallel=True)
def _emax_iter(out, x2, v, demographic, delta, sigma, pi, J, M, N):
"""
Parallel computation of (deterministic) utility elements
"""
for i in prange(N):
mj = 0
for m in prange(M):
mktSize = J[m]
for j in prange(mktSize):
out[mj, i] = delta[mj] + \
x2[mj, 0] * (v[N * m + i, 0] * sigma[0] +
np.dot(pi[0,:], demographic[N * m + i,:])) + \
x2[mj, 1] * (v[N * m + i, 1] * sigma[1] +
np.dot(pi[1,:], demographic[N * m + i,:])) + \
x2[mj, 2] * (v[N * m + i, 2] * sigma[2] +
np.dot(pi[2,:], demographic[N * m + i,:])) + \
x2[mj, 3] * (v[N * m + i, 3] * sigma[3] +
np.dot(pi[3,:], demographic[N * m + i,:]))
mj += 1
return out
@staticmethod
@jit(forceobj=True)
def _indirect_utility(x2, v, demographic, delta, sigma, pi, J, M, N):
"""
Compute indirect utility matrix
"""
sigma = np.abs(sigma) ## positive Sigma
out = np.zeros((sum(J), N)) ## initialization
return BLPEstimator._emax_iter(out, x2, v, demographic, delta, sigma, pi, J, M, N)
@staticmethod
@jit(forceobj=True)
def _mkt_share(x2, v, demographic, delta, sigma, pi, J, M, N):
"""
Compute market shares
- q: quantities consumed
- s: market shares
"""
q = np.zeros((np.sum(J), N)) ## Quantities bought (initialization)
u = BLPEstimator._indirect_utility(x2, v, demographic, delta, sigma, pi, J, M, N) ## Indirect utilities
exp_u = np.exp(u) ## Exponentiate the utilities
first_good = 0 ## Outside good
for m in range(M):
mktSize = J[m] ## market size of market m
numer = exp_u[first_good:first_good + mktSize,:] ## numerator for shares
denom = 1 + numer.sum(axis=0) ## denominator for shares
q[first_good:first_good + mktSize,:] = numer/denom ## quantity(ijm)
first_good += mktSize
s = np.matmul(q, np.repeat(1/N, N)) ## Shares. We use equal weight for all goods.
return [q, s]
@staticmethod
@jit(forceobj=True)
def _solve_delta(s, x2, v, demographic, delta, sigma, pi, J, M, N, tol):
"""
Solve delta using the iterative process / contraction mapping
"""
delta_old = delta.copy()
diff = 1
while diff > tol:
sigma_jm = BLPEstimator._mkt_share(x2, v, demographic, delta_old, sigma, pi, J, M, N)[1] ## The estimated shares
delta_new = delta_old + np.log(s/sigma_jm) ## Contraction mapping
diff = np.max(np.abs(delta_new - delta_old)) ## Update the difference between iterates
delta_old = delta_new.copy()
return delta_old
def _objective(self, params, s, x1, x2, v, demographic, J, M, N, tol, Z, W):
"""
Calculate value of GMM objective function
"""
## Initialization: theta_0
sigma = params[0:4]
pi = params[4:].reshape((4,4))
pi[[0,2,3],1] = 0
pi[[0,2,3],3] = 0
pi[1,2] = 0
if np.min(sigma) < 0:
return 1e20
## Inner loop
self.d.delta = self._solve_delta(s, x2, v, demographic, self.d.delta, sigma, pi, J, M, N, tol)
delta_vector = self.d.delta.reshape((-1,1))
## Linear parameters (alpha, beta) given theta_0, FOC
b = np.linalg.inv(x1.T @ Z @ W @ Z.T @ x1) @ (x1.T @ Z @ W @ Z.T @ delta_vector)
xi = delta_vector - x1 @ b ## Error term \xi
g = Z.T @ xi / np.size(xi, axis=0)
return float(np.sum(J) ** 2 * g.T @ W @ g) ## np.sum(J): Number of obs, JxM
def estimate(self):
"""
Perform BLP estimation
- We choose the same starting values for theta and W as Nevo (2000)
"""
## Initial values for theta
sigma = [.377, 1.848, 0.004, 0.081]
pi1 = [3.09, 0, 1.186, 0]
pi2 = [16.6, -.66, 0, 11.6]
pi3 = [-0.193, 0, 0.03, 0]
pi4 = [1.468, 0, -1.5, 0]
params = np.hstack((sigma, pi1, pi2, pi3, pi4))
## Initial weighting matrix
W0 = np.linalg.inv(self.Z.T @ self.Z)
## Get the optimal weight matrix
params_init_wt = minimize(self._objective, params,
args=(self.product.shares.values, self.X1, self.X2,
self.v, self.demographic, self.J, self.M,
self.N, 1e-4, self.Z, W0),
method="Nelder-Mead")
params_2 = params_init_wt.x
sigma_2 = params_2[0:4]
pi_2 = params_2[4:].reshape((4,4))
self.d.delta = self._solve_delta(self.product.shares.values, self.X2,
self.v, self.demographic, self.d.delta,
sigma_2, pi_2, self.J, self.M, self.N, 1e-4) ## Inner loop
delta_vector = self.d.delta.reshape((-1,1))
b = np.linalg.inv(self.X1.T @ self.Z @ W0 @ self.Z.T @ self.X1) @ \
(self.X1.T @ self.Z @ W0 @ self.Z.T @ delta_vector)
xi = delta_vector - self.X1 @ b
## Update weight matrix
g_ind = self.Z * xi
vg = g_ind.T @ g_ind / np.sum(self.J)
## Optimal weighting matrix
Wo = np.linalg.inv(vg)
## NFXP with optimal weight matrix
params_optimal = minimize(self._objective, params,
args=(self.product.shares.values, self.X1, self.X2,
self.v, self.demographic, self.J, self.M,
self.N, 1e-4, self.Z, Wo),
method="Nelder-Mead")
params_estimated = params_optimal.x
## Calculate final results
final_obj = self._objective(params_estimated, self.product.shares.values,
self.X1, self.X2, self.v, self.demographic,
self.J, self.M, self.N, 1e-4, self.Z, Wo)
delta_vector = self.d.delta.reshape((-1,1))
b = np.linalg.inv(self.X1.T @ self.Z @ Wo @ self.Z.T @ self.X1) @ \
(self.X1.T @ self.Z @ Wo @ self.Z.T @ delta_vector)
return params_estimated, b, final_obj
def print_results(self, params_estimated, b, final_obj):
elapsed_time = time.time() - self.time_start
elapsed_minutes = int(elapsed_time / 60)
elapsed_seconds = int(elapsed_time % 60)
print("Elapsed time: {}m {}s \n".format(elapsed_minutes, elapsed_seconds))
print('-' * 20, " Replication results ", '-' * 20, '\n')
print("Context: \n Length of x vector (k) = ", self.k,
",\n Number of markets (M) = ", self.M,
", \n Number of draws (N) = ", self.N, '\n')
print("GMM objective value = ", round(final_obj/np.sum(self.J), 1), "; Nevo's result = 14.9 \n")
print("alpha(price) = ", round(b[0][0], 3), "; Nevo's result = -32.433 \n")
print("Sigma = \n", params_estimated[0:(self.k+1)], "\n")
print("Pi = \n", params_estimated[(self.k+1):].reshape(((self.k+1), (self.k+1))), "\n")
"""
Usage
"""
if __name__ == "__main__":
estimator = BLPEstimator()
params_estimated, b, final_obj = estimator.estimate()
estimator.print_results(params_estimated, b, final_obj)
Elapsed time: 0m 35s -------------------- Replication results -------------------- Context: Length of x vector (k) = 3 , Number of markets (M) = 94 , Number of draws (N) = 20 GMM objective value = 15.4 ; Nevo's result = 14.9 alpha(price) = -36.908 ; Nevo's result = -32.433 Sigma = [2.49897746e-04 9.33870292e+00 2.45505416e-06 3.07782234e-02] Pi = [[ 3.34507418 0. 1.58247038 0. ] [ 1.64308535 0.2381023 0. 11.13164642] [-0.27366374 0. 0.02471842 0. ] [ 2.20558565 0. -1.99445009 0. ]]
Berry, S., Levinsohn, J., & Pakes, A. (1995). Automobile prices in market equilibrium. Econometrica: Journal of the Econometric Society, 841-890.
Nevo, A. (2000). A practitioner's guide to estimation of random‐coefficients logit models of demand. Journal of economics & management strategy, 9(4), 513-548.
Berry, S., Levinsohn, J., & Pakes, A. (2004). Differentiated products demand systems from a combination of micro and macro data: The new car market. Journal of political Economy, 112(1), 68-105.
Rasmusen, E. (2007). The BLP Method of Demand Curve Estimation in Industrial Organization.
Train, K. E. (2009). Discrete choice methods with simulation. Cambridge university press.
Dubé, J. P., Fox, J. T., & Su, C. L. (2012). Improving the numerical performance of static and dynamic aggregate discrete choice random coefficients demand estimation. Econometrica, 80(5), 2231-2267.