Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
62,472 characters · 16 sections · 81 citation commands
Automatically Differentiable Random Coefficient Logistic Demand Estimation
\lstset{ language=Python, basicstyle=, numbers=left, numberstyle=\tiny\color{gray}, stepnumber=1, morekeywords={@fp}, numbersep=5pt, backgroundcolor=\color{white}, showspaces=false, showstringspaces=false, showtabs=false, frame=single, rulecolor=\color{black}, tabsize=2, captionpos=b, breaklines=true, breakatwhitespace=false, title=\lstname, keywordstyle=\color{codepurple}, commentstyle=\color{codegreen}, stringstyle=\color{black}, }
The random coefficient logistic demand model of berry1995automobile (henceforth BLP) has been a workhorse of the New Empirical Industrial Organization literature, allowing for varied substitution patterns across products, and accounting for endogeneity of price. The reliability of its estimation has been the subject of rigorous debate nevo2000practitioner,conlon2020best,knittel2014estimation, and the estimator itself has been the study of many proposed advances in econometric techniques as a sophisticated yet widely used structural model hong2020blp,forneron2020inference.
The most common implementation of the BLP estimator involves the use of a nested fixed point (NFP) as an inner loop within an outer loop of GMM estimation, although we acknowledge the Mathematical Programming with Equilibrium Constraints (MPEC) approach of dube2012improving, which is beyond the scope of this paper. dube2012improving,conlon2020best find that derivative-free optimization algorithms such as the Nelder-Meade or simplex algorithms often fail to converge or converge to the wrong solution. As such, the literature has settled on the use of analytical derivatives with a derivative-based optimization algorithm such as L-BFGS. nevo2000practitioner provides the analytical derivative for demand-only (DO) BLP in detail, and conlon2020best indicate that the same is possible for demand-and-supply (DS) BLP, although it involves tensor products.
Notwithstanding the ongoing debate knittel2014estimation,conlon2020best about whether the BLP estimator can be reliably implemented, it is clear that gradient-based estimation with tight tolerances in the NFP are required for any hope of reliability for the non-linear, non-convex optimization problem. However, the reliance on analytical derivatives poses challenges for the robust and reliable estimation of variants on the BLP model such as dynamic demand decisions involving NFPs within NFPs hendel2006measuring or sequential NFPs lee2013vertical. Quasi-Bayesian estimation using the Laplace Transform Estimator (LTE) of chernozhukov2003mcmc with the latest gradient-based Markov Chain Monte Carlo (MCMC) samplers such as Hamiltonian Monte Carlo (HMC) is also difficult in practice. By contrast, automatic differentiation offers easy extensibility, since automatically differentiable modifications preserve the automatic differentiability of the resulting objective function, and straightforward integration into the existing software ecosystem around the latest MCMC sampling techniques.
Our contribution is fourfold. First, we derive the BLP objective function (both DO and DS variants) using block diagonal matrix multiplication in place of for-loops so that the resulting objective function is reverse mode automatically differentiable. Second, we briefly outline what automatic differentiation is, and how the BLP objective function can be phrased as an easily extensible composition of automatically differentiable building blocks. Third, we point out the algorithmic improvements in our reference implementation, the admest\footnote{Automatically Differentiable M-ESTimation.} Python package, including where we depart from standard practice in the literature. Fourth, we evaluate these algorithmic improvements in Monte Carlo simulations examining the relative performance of 2S-GMM, CUE, and LTE estimators with automatic differentiation, with 2S-GMM using analytical gradients as a benchmark.
Preliminary simulation results confirm previous theory hansen1996finite,newey2004higher,newey2009generalized,donald2000jackknife that the CUE has a lower bias but higher variance compared to the 2S-GMM estimator given the weak instrument setting of the Monte Carlo simulation. On the performance of the LTE, we find that lifting the credible intervals from the MCMC sampler directly achieves good empirical coverage for the parameters estimated in the outer GMM loop. This is in contrast to hong2020blp, who found it necessary to correct for simulation error to achieve good empirical coverage. Using credible intervals directly results in poor coverage of the concentrated out linear parameters, but using the LTE posterior mean as a point estimate for frequentist uncertainty quantification on these parameters restores good empirical coverage compared to frequentist estimators. Our gradient-based No-U-Turn Sampler variant of HMC (NUTS-HMC) delivers higher effective sample sizes than a Random Walk Metropolis-Hastings (RWMH) implementation with significantly fewer draws forneron2020inference.
In Section (ref) we review the BLP model and show how to construct the the Continuously Updating Estimator (CUE) of hansen1996finite. Section (ref) gives a brief overview of automatic differentiation, and how it can be applied to the BLP model. Section (ref) details the algorithmic improvements of our choice of frequentist and Bayesian estimation techniques. Section (ref) documents the set up of our Monte Carlo simulation exercise, while Section (ref) contains the result of our empirical exercise. Section (ref) concludes. Appendix (ref) derives the 2S-GMM objective function and the analytical gradient without the use of for-loops.
As in berry1995automobile, we start with some individual $i$ in market $t$ that faces a choices of $J_t$ products, including the outside alternative denoted $j=0$. The utility of consumption for each product $u_{ijt}$ is given by
$$ u_{ijt} = \sum_{k\in K_1}x_k\beta_{k,r} - \alpha_r p_{jt} + \xi_{jt} + \varepsilon_{ijt} $$
where the parameter vector $(\beta_r, -\alpha_r)'$ has a mean component $(\beta, -\alpha)'$ and covariance matrix with Cholesky decomposition $\Sigma$, $\varepsilon_{ijt}$ is a Type 1 Extreme Value error, and $\xi_{jt}$ is an unobserved structural error. Under the berry1995automobile formulation, $\Sigma$ is assumed to be a diagonal as opposed to a lower triangular matrix, with the non-zero components denoted $\sigma_k$. Define $\theta_1 = (\beta, -\alpha)$, the full vector of features as $X_{1,jt}'$, and the subset of features with non-zero entries in the covariance matrix as $X_{2,jt}'$, and the standardized vector of individual-level taste heterogeneity as $\nu_i$, then we have
Consumers select the product that gives them the most utility, where the value of the outside good is normalized to zero so $u_{i0t} = 0\ \forall\ i,t$. Let $\mathcal{N}(t)$ denote the set of products in market $t$, then the observable market shares are given by
$$ s_{jt} = \int \frac{\exp(\delta_{jt} + \mu_{ijt})}{1 + \sum_{k \in \mathcal{N}(t)}\exp(\delta_{jt} + \mu_{kjt})} dF(\nu)$$
In practice, this integral is computed numerically using simulation under a distributional assumption for $\nu_i$. We assume independent normality (across products and individuals), as is common in the literature. berry1995automobile show that we can invert the simulated and observed shares to recover $\delta_{jt}$ using a contraction mapping. Once $\delta_{jt}$ is recovered, we can retrieve the structural error $\xi_{jt}$. Denote all observable data by $D$. With demand-side instruments $Z^D$, we have the moment conditions $g(\theta; D_{jt}) = Z^{D}_{jt}\xi_{jt}$, and $\mathbb{E}\left[g(\theta; D_{jt})\right] = 0$, which we can then estimate using GMM.
The following derivation is novel in that it eschews the use of for-loops while being numerically equivalent, avoiding the criticism conlon2020best of previous vectorization methods nevo2000practitioner,knittel2014estimation. Since each iteration of the for-loop across markets does not depend on computations from other iterations, the loop constitutes a finite-dimensional linear operator which definitionally has an equivalent expression as a series of matrix operations -- in this case, appropriately chosen block diagonal matrices.
In considering numerical integration, we make the subjective choice of assuming overlapping simulation draws hong2020blp, although there are potential efficiency gains from using separate draws for each market kristensen2017higher,lee1995asymptotic. We further use quasi-Monte Carlo nodes with equal weight rather than quadrature nodes, although the derivations easily extend to the case of quadrature nodes.
\paragraph{Share Inversion}
The contraction map given by berry1995automobile is
$$ \delta_{n+1}(\delta_n; \hat{s}(\delta_n; \theta_2), \mathcal{S}) = \delta_n + \log(\mathcal{S}) - \log(\hat{s}(\delta_n; \theta_2)) $$
where $\mathcal{S}$ is the observed share (and its dependence is henceforth suppressed for notational simplicity), and $\hat{s}$ is the simulated shares. For a given product $j$ in market $t$, with each simulated individual denoted by $i$, the simulated share is given by
$$ \hat{s}_{ji}(\delta_n; \theta_2) = \frac{\exp(\delta_n + \mu_{ji}(\theta_2))}{1 + \sum_j \exp(\delta_n + \mu_{ji}(\theta_2))} $$ $$ \mu_{ji}(\theta_2) = \sum_kx_{j,k}\nu_{i,k}\theta_{2,k} $$
At each iteration of the fixed point the term $\mu_i(\theta_2)$ does not change (this value is fixed in the outer loop of the non-linear GMM algorithm). Therefore, it makes sense to cache the values of $\mu_{ji}$ for re-use across each fixed point iteration, rather than computing it at every iteration. Note also that we can express the $\mu$ matrix for all goods and individuals as
$$ \mu(\theta_2) = X_2 \Omega_{\theta_2} \nu' \in \mathbb{R}^{J \times R}$$
where $\nu \in \mathbb{R}^{R\times K_2}$ is the matrix of simulation draws and $\Omega_{\theta_2}$ is a diagonal matrix with diagonal entries corresponding to $\theta_2$. The compute time of this $\mu$ matrix scales with the number of simulation draws $R$, so caching it yields significant savings in compute time especially when we increase $R$ for more precision in numerical integration.
To think about how to recover the shares across multiple markets in a single matrix operation, we first examine the case of a single market with number of products $J_t$, where we want to recover the matrix $\hat{s}(\delta_n;\theta_2) \in \mathbb{R}^{J_t \times R}$. We already have the matrix $\mu(\theta_2) \in \mathbb{R}^{J_t\times R}$, and we can perform column-wise addition (denoted $+^c$) of $\delta_n$ on it and apply the exponential function element-wise to get $\exp(\delta_n + \mu_{ji}(\theta_2))$, so the only term that remains is $\sum_j \exp(\delta_n + \mu_{ji}(\theta_2))$, which we want to compute at the individual level (so each entry within the same column has the same value, and this value should differ across the $R$ columns). Taking $\iota_{J_t} \in \mathbb{R}^{J_t\times J_t}$ (the “summation matrix”) to be the matrix of all ones, we can use the following formula
$$ \iota_{J_t} \exp(\delta_n +^c \mu(\theta_2)) \in \mathbb{R}^{J_t\times R} $$
This gives a single expression for the matrix of simulated shares as
$$ \hat{s}(\delta_n;\mu(\theta_2)) = \exp(\delta_n +^c \mu(\theta_2)) \oslash \left[1 + \iota_{J_t}\exp(\delta_n +^c \mu(\theta_2))\right] $$
where the addition of 1 is element-wise and $\oslash$ denotes element-wise division (as opposed to the Hadamard element-wise product denoted $\odot$). Moving to the many markets scenario, define instead $\iota_J$ to be the block diagonal with $\iota_{J_t}$ for each market $t$ in each block\footnote{The matrix $\iota_J$ is also constant across iterations of the fixed point, so it can be pre-computed and cached. In fact, it is invariant to the choice of $\theta_2$, so it should be passed in from outside of the outer non-linear GMM loop to save even more time.}, so we have
$$ \iota_J =
$$
Then the share matrix across all products, markets, and individuals can be computed by
$$ \mathbb{R}^{J \times R} \ni \hat{s}(\delta_n;\theta_2) = \exp(\delta_n +^c \mu(\theta_2)) \oslash \left[1 + \iota_J\exp(\delta_n +^c \mu(\theta_2))\right] $$
\paragraph{Numerical Safeguards}
conlon2020best note that in situations where $\delta_n + \mu_{ji}(\theta_2) \gg 0$, $\exp(\delta_n + \mu_{ji(\theta_2)})$ will overflow, leading to numerical issues when trying to divide infinity by infinity. They propose choosing (at the market level) $m_i = \max\{0, \max_{j}\delta_n + \mu_{ji}(\theta_2)\}$, and then note that
$$ \frac{\exp(\delta_n + \mu_{ji}(\theta_2))}{1 + \sum_j \exp(\delta_n + \mu_{ji}(\theta_2))} = \frac{\exp(\delta_n + \mu_{ji}(\theta_2) - m_i)}{\exp(-m_i) + \sum_j \exp(\delta_n + \mu_{ji}(\theta_2) - m_i)} $$
The choice of maximum provides the best protection against numerical overflow, but it is not expressible as a matrix operation. We propose re-centering using the mean instead. This guards against overflow and underflow, is expressible as a matrix operation, and is likely to be sufficient as long as $\|\theta_2\|_\infty$ is of reasonable size. The failure mode when $\|\theta_2\|_\infty \gg 0$ is unlikely to be encountered in practice, as this would lead to other numerical issues in the objective function.\footnote{conlon2020best recommend using bounds on the value of $\theta_2$ to ensure a well-behaved objective function. Within reasonable bounds, the performance of the mean is likely to be at least equivalent to using the max.} This operation is in fact similar to the summation within markets at the individual level that we have already detailed. Define the “averaging matrix”\footnote{One can retrieve $\iota_J$ from $\iota_{J/J}$ by performing a boolean check for non-zero elements (encoding true as 1 and false as 0), so in implementation one should construct $\iota_{J/J}$ first, and then perform the boolean check to retrieve $\iota_J$, with both operations happening outside the BLP algorithm.}
$$ \iota_{J/J} =
$$
then the matrix given by $\iota_{J/J} (\delta_n +^c \mu(\theta_2))$ gives the market level mean of $\delta_n +^c \mu(\theta_2)$ for a given individual. Incorporating this guard against numerical overflow, our final expression for the shares is given by
$$ \hat{s}(\delta_n;\theta_2) = \exp(\delta_n +^c \mu(\theta_2) - \iota_{J/J} [\delta_n +^c \mu(\theta_2)]) \oslash \left[1 + \iota_J\exp(\delta_n +^c \mu(\theta_2) - \iota_{J/J} [\delta_n +^c \mu(\theta_2)])\right] $$
To take the average over individuals, we can just take the row sums of this matrix.\footnote{It is advisable to cache $\hat{s}(\delta_n;\theta_2)$ before row-summation, as it will be needed in this form for other computations.} We are now ready to construct the CUE objective function for the BLP problem.
\paragraph{Demand-Only CUE}
We partition the parameter space $\theta$ into two components. $\theta_1 \in \mathbb{R}^{K_1}$ is the vector of linear parameters with the corresponding design matrix $X_1$, while $\theta_2 \in \mathbb{R}^{K_2}$ is the vector of non-linear parameters governing heterogeneous tastes with the corresponding design matrix $X_2$. We denote the fixed point solution of the contraction mapping $\lim_{n\to\infty}\delta_n(\delta_0; \mu(\theta_2)) := \lim_{n\to\infty}\delta_n(\delta_{n-1}; \mu(\theta_2))$, and the matrix of instruments by $Z\equiv Z^D$. The Nested Fixed Point (NFP) algorithm of berry1995automobile implemented as a CUE objective function is given by Algorithm (ref).\footnote{For comparison and completeness of documentation, the 2S-GMM version of demand-only NFP BLP along with its gradient without for-loops is given in Appendix (ref).}\footnote{We scale the objective function by $\frac{1}{2}$ to be consistent with the LTE formulation of chernozhukov2003mcmc.}
There are two remarks to be made about the formulation. First, while the rest of the literature uses a one-step linear GMM to compute the concentrated out parameters $\theta_1$, we use a 2-Stage GMM (2S-GMM) instead, since we do not have the weight matrix as a parameter.\footnote{The weight matrix is initialized using the 2-Stage Least Squares (2SLS) weight matrix $(Z'Z)^{-1}$.} Second, for the computation of the weight matrix within the CUE objective function, standard errors, and computing out concentrated out parameters we need an estimator for $Var(Z'\xi)$. Under the assumption of heteroskedascity, one popular choice is $\hat{Var}(Z'\xi) = Z'\Omega_\xi Z$ where $\Omega_\xi$ is a diagonal matrix with the squared residuals $\xi_{jt}^2$ as its diagonal entries. We find better performance using the following formulation as proposed by gallant2007statistical
$$\hat{Var}(Z'\xi) = Z'\Omega_\xi Z - Z'\xi\xi'Z $$
We use the standard GMM standard errors\footnote{We acknowledge the work of newey2009generalized in showing that the usual GMM standard errors are biased downwards for the CUE.} for $\hat{\theta}$ denoted $\Sigma_{\hat{\theta}}$, which berry1995automobile derive to be
$$ \Sigma_{\hat{\theta}} = \left[\frac{\partial \xi(\theta_2)'Z}{\partial \theta} \left[\sum_{i=1}^3 V_i\right]^{-1} \frac{\partial Z'\xi(\theta_2)}{\partial \theta}\right]^{-1} $$
where $V_1,\ V_2,\ V_3$ represent the variation arising from the process generating the product characteristics, from the consumer sampling process, and from the simulation process respectively. We have $V_1 = Var(Z'\xi)$, and with a large number of underlying consumers, $V_2$ is negligible. With a sufficient number of simulation draws relative to the number of products, $V_3$ is also negligible, so we use as an empirical estimator
$$ \Sigma_{\hat{\theta}} = \left[\left(\frac{\partial \xi(\theta_2)'Z}{\partial \theta}\right) \left(Z'\Omega_\xi Z - Z'\xi\xi'Z\right)^{-1} \frac{\partial Z'\xi(\theta_2)}{\partial \theta}\right]^{-1} $$
evaluated at the point estimate $\hat{\theta}$, where $\frac{\partial Z'\xi(\theta_2)}{\partial \theta}$ is computed using automatic differentiation.
conlon2020best show that as long as the supply side is correctly specified, its inclusion can increase the efficiency of estimation, especially for situations where the excluded demand-side instruments are relatively weak. Using the first order condition of firms under the assumption of Bertrand competition, we can generate an additional set of moment conditions. At the market $t$ level, we have the multi-product Bertrand markup $\eta_t(p_t,s_t,\theta_2)$ estimated by
$$ \eta_t(p_t, s_t, \theta_2) \equiv \Delta_t(p_t)^{-1}s_t = p_t - c_t $$
where $\Delta_t(p_t,s_t,\theta_2)$ is the intra-firm demand derivative given by
$$ \Delta_t(p_t) \equiv -\mathcal{H}_t \odot \frac{\partial s_t}{\partial p_t}(p_t) $$
$\mathcal{H}_t$ refers to the ownership matrix with each $(i,j)$ entry indicating whether the same firm produces products $i$ and $j$. With this, we can estimate marginal costs by taking $c_t = p_t - \eta_t(\theta_2)$, and parametrize it (across all markets and products) as
$$ f_{MC}(p - \eta(\theta_2)) = X_3\gamma + \omega$$
where $\omega$ is a structural error term, $\gamma$ is the vector of coefficients on the linear supply-side features, and $X_3$ is a design matrix of features that enter supply linearly. $f_{MC}$ is the functional form assumed for marginal costs, typically the identity function or the logarithm. This allows the construction of moment conditions $\mathbb{E}\left[\omega_{jt}Z^S_{jt}\right] = 0$. We can then estimate the full demand and supply system with multiple equation GMM.
\paragraph{Empirically Estimating Marginal Cost}
We focus on an empirical estimate for $\frac{\partial s_t}{\partial p_t}$. Denote the coefficient on price as $\alpha$.\footnote{Our notation here differs from the rest of the literature in that $\alpha < 0$ is the coefficient on price rather than its size.}, and $\alpha_i = \alpha + \theta_{2,\text{price}} \nu_{i,\text{price}}\ \forall\ i = 1,\ldots, R$ to be the vector of (simulated) individual level coefficients on price, and the matrix $\Omega_\alpha$ to be a diagonal matrix with $\alpha_i$ as its diagonal elements. In the case where price is not a non-linear feature, set $\theta_{2,\text{price}} = 0$ and this reduces to $\alpha_i = \alpha\ \forall\ i$. Then at the market $t$ level we have
$$ \frac{\partial s_j}{\partial p_k} =
$$
Following morrow2011fixed, we can decompose this into a diagonal matrix $\Lambda_t$ and a dense matrix $\Gamma_t$ so we have
To move to the many markets scenario without for-loops, we note that we can simply substitute the problem-level matrices for $\Lambda$. For $\Gamma$, we want a block diagonal matrix rather than a dense matrix. Hence, we “over-provision” a dense matrix and use Hadamard multiplication by $\iota_J$ to “slice” out the entries that we want\footnote{See Appendix (ref) for a similar approach taken for the analytical gradient, and associated trade-offs.}, so we have
$$ \Gamma = \hat{s}\Omega_{\alpha}\hat{s}' \odot \iota_j / R$$
We then estimate $\hat{\eta}$ with\footnote{We use the Moore-Penrose Pseudo-Inverse denoted $\dagger$ when inverting the matrix of intra-firm derivatives because the matrix can be computationally singular when the observed outside share is large and predicted shares are small.}
$$ \hat{\eta} = \left(-\mathcal{H}\odot \frac{\partial \hat{s}}{\partial p}\right)^\dagger \mathcal{S} $$
\paragraph{Demand-Supply CUE}
We partition the parameter space $\theta = (\theta_1, \theta_2, \theta_3)$, where $\theta_1 \in \mathbb{R}^{K_1}$ are the parameters on the features that enter demand linearly including the price parameter $\alpha$ (which is negative), $\theta_2 \in \mathbb{R}^{K_2}$ are the parameters on the features that enter demand non-linearly, and $\theta_3 \in \mathbb{R}^{K_3}$ are the parameters on the features that enter supply linearly. We present the DS CUE objective function where the outer loop only comprises the non-linear parameters $\theta_2$. This is a departure from conlon2020best, who include the linear parameter on price $\alpha$ in the outer GMM loop. They do this because it allows for the sharing of the weight matrix between the concentrated out parameters and the outer GMM loop. Our simulation results suggest that unless one imposes parameter bounds, this tends to have poor performance when applied to the CUE. Instead, we estimate the demand-side of BLP, using a 2-Stage linear IV GMM procedure to recover the concentrated out linear parameters including $\hat{\alpha}$, and then estimate the supply-side using $\hat{\alpha}$. Doing so forces $\hat{\alpha}$ to be in a feasible region (and in particular, negative) without the use of any parameter bounds.
On top of the data used for demand-only BLP, we need $\mathcal{H}$ the conduct matrix, $f_{MC}$ the functional form for marginal costs, $X_3$ the matrix of features that enter supply linearly, and $Z^S$ supply-side instruments, the last of which is used to construct $Z \equiv
$ the full block matrix of instruments. Denote the vector of firm indices to be $\vec{firm}$, then one trick to computing $\mathcal{H}$ is to compute the boolean matrix
$$ \mathcal{H}: \left[\vec{firm}(1/\vec{firm})'\right] \overset{?}{=} 1 $$
encoding True to be 1 and False to be 0. The CUE objective function is given by Algorithm (ref), and uncertainty quantification proceeds in similar fashion to demand-only BLP. Denote $\varepsilon = (\xi, \omega)$, then we have
$$ \Sigma_{\theta} = \left[\frac{\partial \varepsilon(\theta)'Z}{\partial \theta} \left(Z'\Omega_\varepsilon Z - Z'\varepsilon\varepsilon'Z\right)^{-1} \frac{\partial Z'\varepsilon(\theta)}{\partial \theta}\right]^{-1} $$
where $\frac{\partial Z'\varepsilon(\theta_2)}{\partial \theta}$ is computed using automatic differentiation.
Traditionally, gradient-based estimation of the BLP objective function has relied on deriving the analytical gradient. However, conlon2020best point out that the mere inclusion of a supply-side greatly complicates this endeavour, requiring the use of tensor products. Automatic differentiation offers a way to estimate gradients up to numerical precision (avoiding the error from a finite differences approximation) without the need for extensive derivation.
This section is not a substitute for a more extensive treatment of automatic differentiation found in rall1996introduction,jax2018github. Readers interested in the history of automatic differentiation and its use in machine learning can refer to iri1991history and baydin2018automatic respectively. Given a function $f: \mathbb{R}^n \to \mathbb{R}^m$, we think of the Jacobian evaluated at $x \in \mathbb{R}^n$ denoted $\partial f(x)$ to be an $m\times n$ matrix representing the differential of $f$ at every point where $f$ is differentiable. For some $y \in \mathbb{R}^m$, the Jacobian-vector product $\partial f(x)\cdot y$ represents the best linear approximation to the change $y-x$ in the neighborhood around $x$. Forward Mode automatic differentiation computes the Jacobian-vector product. The Jacobian of the function is build one column at a time by keeping track of the partial derivative of each intermediate step in the function, or the tangent of the primal. This is the most efficient for the case where $f: \mathbb{R} \to \mathbb{R}^m$, where the entire Jacobian matrix is computed in a single pass. Reverse Mode automatic differentiation computes the vector-Jacobian product, which is the transpose of the Jacobian vector product. The Jacobian matrix is computed one row at a time in two sweeps. The first forward sweep keeps track of the dependencies that intermediate values have on the original input, and then the reverse sweep computes the derivative by propagating the adjoints in reverse. The special case where $f: \mathbb{R}^n \to \mathbb{R}$, where the Jacobian is computed in a single pass, is known as back-propagation. Given a scalar-valued objective function, back-propagation is the most efficient of the two automatic differentiation modes, explaining its popularity in the machine learning literature and why it is our mode of choice.
Similar to the chain rule for analytical derivatives, compositions of automatically differentiable functions preserve automatic differentiability. This means that we get gradients of any variant of the original BLP algorithm “for free”, provided the variation is phrased as an automatically differentiable function. The simplest example of this is that any automatically differentiable moment function admits an automatically differentiable CUE objective function, since the weight matrix estimation step only involves matrix operations.
With modern machine learning libraries, compositions of the following building blocks are automatically differentiable:
By contrast, automatic differentiation does not work well with loop constructs, both for-loops and while-loops. Forward-mode differentiation is possible with loop-constructs by “unrolling” the loops at every evaluation of the function, a costly operation, but reverse-mode differentiation, and therefore back-propagation, is not possible. It is clear that the derivation of both the DO and DS BLP objective functions in Section (ref) only comprises compositions of building blocks 1.--4., and is thus automatically differentiable. The derivative of a fixed-point solution relies on the implicit function theorem, and is typically not supported out of the box by most major automatic differentiation libraries. Hence, we detail our implementation of the reverse-mode automatic derivative of a fixed-point solution.
\paragraph{Reverse Mode Automatic Derivative of Fixed-Point Solutions}
Numerical methods such as fixed-point iteration that effectively rely on a while-loop are not amenable to back-propagation, and are slow to automatically differentiate using forward mode. However, nevo2000practitioner already observed that we can use the implicit function theorem to retrieve the derivative of $\delta$ with respect to $\theta_2$ given by
$$ \frac{\partial\delta}{\partial\theta_2} = -\left[\frac{\partial \hat{s}}{\partial\delta}\right]^{-1} \frac{\partial \hat{s}}{\partial\theta_2} $$
As pointed out by jax2018github,jeon2021differentiable, this can be applied more generally. Start with some contraction mapping $x = f(\theta, x)$ which implies a solution function $\theta \mapsto x^*(\theta)$, where we want to retrieve the reverse-mode derivative (vector-Jacobian product) $v' \mapsto v'\partial x^*(\theta)$. Assume that the implicit function theorem holds, and that $x^*(\theta) = f(\theta,x^*(\theta))\ \forall\ \theta$ in some neighborhood around $\theta_0$ which we want to differentiate. Differentiating both sides of the equation, we have
$$ \partial x^*(\theta_0) = B + A\partial x^*(\theta_0) = (I-A)^{-1}B $$
where $A = \frac{\partial f(\theta_0, x^*(\theta_0))}{\partial x^*(\theta_0)},\ B = \frac{\partial f(\theta_0, x^*(\theta_0))}{\partial \theta_0}$. Therefore, we can evaluate vector-Jacobian products by
$$ v'\partial x^*(\theta_0) = v'(I-A)^{-1}B = w'B $$
where $w' = v'(I-A)^{-1} = v' + w'A$. We can equivalently think of $w'$ as the fixed point of the map $u' \mapsto v' + u'A$, so we can write the reverse mode derivative of a contraction mapping as the solution to another contraction mapping, avoiding the expensive operation of matrix inversion.
\paragraph{Quasi-Newton Methods}
Most of the BLP literature use quasi-Newton methods for frequentist optimization and estimation. When the objective function is convex, quasi-Newton methods that utilize an approximation to second-order Hessian information such as BFGS and its variants will dominate a gradient-descent based approach that does not. However, in the case of highly non-convex objective functions that can be ill-behaved outside a feasible region such as the BLP objective function, quasi-Newton methods can have issues with numerical stability and fail to converge knittel2014estimation. conlon2020best resolve this by noting that imposing parameter bounds helps with numerical stability, both in replications of berry1995automobile,nevo2001measuring, as well as simulation exercises. We are able to confirm this finding. However, the choice of parameter bounds is a subjective hyper-parameter that can be difficult to justify in empirical work, motivating the search for an optimization method that functions well even without parameter bounds.
\paragraph{Modern Gradient Descent Methods}
For ill-conditioned problems, gradient descent methods which only utilize first-order derivative information controlled by the learning rate offers more consistent results without parameter bounds across a range of reasonable starting values. The challenge here is to have an optimization method that is relatively robust to local stationary points, but still achieves good convergence. There are two main branches of modern gradient descent algorithms, Stochastic Gradient Descent (SGD) and adaptive learning rate or “momentum” based algorithms dogo2018comparative.
In SGD (or mini-batch SGD mSGD) schaul2013unit, each update is computed using the gradient with respect to a single (small batch of) observation(s). This helps with local stationary points because the stationary point of the objective function with respect to each (batch of) observation(s) is likely to be different, and this difference will allow the optimizer to “escape” local minima. We do not test or employ SGD in our implementation because it is not amenable to Just-In-Time-compilation (JIT-compilation). The only principled way to batch the data within the BLP setting is to group them by market. However, doing so will change the dimension of the design matrices at each iteration of the objective function, defeating JIT-compilation. SGD also relies on annealing the learning rate over iterations to control the over-shoot from using batched gradient descent, which is hard to calibrate in practice.
Adaptive learning rate optimization algorithms such as RMSProp hinton2012neural, AdaGrad duchi2011adaptive, and Adam kingma2014adam,reddi2019convergence adjust the learning rate on a per-parameter basis based on a memory of past gradients. The intuition here is that if there is a sudden and large change in the direction of the optimization update, then we should be cautious about moving in that direction. Hence, momentum-based optimizers keep a “memory” of previous gradients, and generates an update that attenuates sudden, large changes and amplifies movements in consistent directions. This momentum simultaneously addresses the exploding gradient problem and allows the optimizer to escape local optima. We use the AdaBelief algorithm zhuang2020adabelief implemented in Optax optax2020github. It is an Adam variant which incorporates curvature information without approximating the Hessian by utilizing deviations between predicted and actual gradients. We find that with adaptive learning rate gradient descent optimizers, imposing parameter bounds is no longer necessary, removing another source of hyper-parameter arbitrariness in estimation.
\paragraph{Choice of Sampler}
chernozhukov2003mcmc use a RWMH in their original paper on LTEs. However, forneron2020inference find that RWMH takes a long time to converge for the BLP problem, resulting in a low effective sample size even with a long Markov chain. hong2020blp propose using Hamiltonian Monte Carlo (HMC) duane1987hybrid,neal2011mcmc instead, which utilizes gradient information in sampling.\footnote{betancourt2017conceptual provides a comprehensive introduction to HMC and its variants.} We can think of the proposal process of RWMH as a blind-folded person throwing a ball with random direction and momentum. Wherever the ball lands, we check the value of the objective function relative to our current position, and assign a high acceptance probability if there is a large decrease in the objective function value. In HMC, we can think of the ground as having gravity wells, where regions of low objective function value have a higher gravity than regions of high objective function value. The same ball with random direction and momentum no longer travels in a straight line, but is instead pulled towards regions of low objective function value, speeding convergence. Because the ball is thrown with random momentum, there is still a chance that the sampler will escape the pull of the gravity wells and therefore explore beyond the closest local minimum. The curved path of the ball represents an integral that is hard to compute exactly. In practice, a simplectic integrator based on Hamiltonian dynamics is used, breaking up the curved path into $L$ steps of length $\varepsilon$, which are the hyper-parameters of the sampler.
We extend the work of hong2020blp by using a variant of Hamiltonian Monte Carlo that is less sensitive to the choice of hyper-parameters. HMC's hyper-parameters are $(L, \varepsilon, M)$, which are the number of steps in the leap-frog symplectic integrator, the size of each step in the integrator, and the mass matrix or the covariance matrix of the proposal distribution respectively\footnote{Correspondence with hong2020blp authors confirm their choice of parameters as $L=10,\ \varepsilon=0.1,\ M = I$, as in neal2011mcmc.}. The tuning of these hyper-parameters is artisanal and problem-specific, and hoffman2014no outline the effects of a poorly chosen set of $(L, \varepsilon, M)$. An $L$ that is too large results in wasted computation, and in the converse results in a poor exploration of the parameter space. $\varepsilon$ controls the granularity of exploration, and a well-chosen $M$ that approximates the posterior covariance matrix increases the rate of convergence of the Markov Chain. They propose a No-U-Turn Sampler (NUTS) variant of HMC that automatically tunes the hyper-parameters. Instead of a single ball thrown in one direction, NUTS uses two balls thrown in opposite directions with the same random momentum. It grows the trajectory of both balls using the simplectic integrator (choosing which to advance at random) until one of the paths turns on itself, indicating that remainder of the path is likely to backtrack. It then samples a point along this bi-directional path as the proposed candidate. This is an auto-tuning of the number of steps parameter $L$ on a per-sample basis. $\varepsilon$ and $M$ are tuned in a set of warm-up draws using prime-dual averaging for a target acceptance rate, and the empirical covariance matrix of the warm-up draws respectively. With minor modifications, NUTS-HMC has become the standard in the machine learning literature carpenter2017stan, phan2019composable,salvatier2016pymc3, and thus our MCMC sampler of choice as implemented in NumPyro phan2019composable.
\paragraph{Sampler Settings and Diagnostics}
Rather than a single long chain of MCMC samples that is common in the economics literature hong2020blp,chernozhukov2003mcmc,forneron2020inference, we use many smaller chains so that we can compute MCMC mixing diagnostics that are common in the computer science literature. The use of a single Markov chain can lead to a failure to diagnose non-convergence, and this has been noted as far back as gelman1992single. In particular, the seeming convergence of a Markov chain to a single mode may fail to detect multi-modality that is far apart -- which in economics implies a potential identification failure. Similarly, the performance of a single chain may fail to detect that the chain is stuck in a region of high-curvature, with a step-size that is too large to allow for acceptable proposals to “escape” -- which can imply the brittleness of frequentist estimation through optimization. We use the rank-normalized split-$\hat{R}$ of vehtari2021rank as a diagnostic, which checks for whether the variance of the first half of a chain fails to match the second half, and whether the variance of multiple chains together is greater than each individual chain. An affirmative of either condition shows up as a high $\hat{R}$, and indicates a failure to converge. The $\hat{R}$ is an oft-recommended convergence diagnostic used across a range of MCMC statistical packages salvatier2016pymc3,carpenter2017stan. carpenter2017stan suggests using at least four Markov chains in estimation using a split-$\hat{R}$ cut-off of 1.05 to indicate a lack of convergence.
MCMC sampling can also produce draws which are auto-correlated, which reduces the effective region explored by the sampler. This is clear by construction with a Gibbs sampler, but persists even with RWMH. Therefore, one common diagnostic is to compute the effective sample size $N_{eff}$ of the sampler, accounting for this auto-correlation forneron2020inference. If $N_{eff} \ll N$, where $N$ is the number of samples to be drawn, then this could indicate numerical issues with the MCMC sampler.\footnote{HMC-NUTS attempts to perform anti-correlated draws, so in fact $N_{eff} > N$ is possible.} We use the rank-normalized version of $N_{eff}$ vehtari2021rank, which abstracts away potential issues with heavy-tailed distributions and reparametrization.
Across all empirical exericse, we use the $R_d$ sequence of roberts2018unreasonable for numerical integration. It has good performance while eschewing the need to tune any hyper-parameters such as co-primes for Halton sequences. All empirical exercises were run GPU-accelerated on consumer hardware, inside the Docker container provided by admest.\footnote{The actual system comprises an Intel i9-10850K CPU and an Nvidia RTX 2070 Super GPU running Fedora 34 as the host operating system. The Docker container builds on an Ubuntu 18.04 image provided by Nvidia for CUDA/CUDNN development.} For the 2-Stage GMM models, we use the weight matrix
$$ W = \left[Z'\Omega_{\xi}Z\right]^{-1} $$
instead of the CUE weight matrix specified in Section (ref) to prevent numerical issues when computing the analytical standard errors.
Our simulation set-up largely mirrors that of conlon2020best. We randomly draw utilities and costs, but solve for equilibrium prices and shares. Over 100 different synthetic datasets, we set $T=20$ markets and choose a random number of firms from 2 to 10, and then have each firm produce between 3 to 5 products. Similar to conlon2020best, this procedure yields $200 < N < 600$ on average.
In order to make the choice of weight matrix consequential, we introduce heteroskedasticity by construction into the structural erorrs $(\xi_{jt},\omega_{jt})$. we first draw $(\psi^1_{jt}, \psi^2_{jt})$ from independent normal distributions with variance 0.2, and $\phi_{jt}^1,\ \phi_{jt}^2 \sim Unif(0.5, 2)$. We then taking $\xi_{jt} = \psi^1_{jt} \times \phi_{jt}^1,\ \omega_{jt} = \psi^2_{jt} \times \phi_{jt}^2$ as the structural errors for the simulation. Linear demand characteristics and supply characteristics follow the same formulation as conlon2020best with $[1,x_{jt},p_{jt}]$ and $[1,x_{jt},w_{jt}]$ respectively. The two exogenous characteristics are drawn from the standard uniform distribution, and we assume that the coefficient on $x_{jt}$ is random. We assume Bertrand competition for our conduct matrix, and a linear function on marginal costs so $c_{jt} = [1,x_{jt},w_{jt}]\theta_3 + \omega_{jt}$. The endogenous prices and shares are computed using the $\zeta$-markup approach of morrow2011fixed. In simulating the data, we use 1,000 draws from the $R_d$ sequence to perform numerical integration, to prevent fluctuations across simulations arising from integration error. We use similar true values to conlon2020best, with $\theta_1 = [-7,6,1],\ \theta_2 = [3],\ \theta_3 = [2,1,0.5]$. This generates $Corr(p_{jt},w_{jt}) \approx 0.2$, implying relatively weak cost-shifting instruments, and outside shares that average around 0.9.
We estimate both Demand-Only and Demand-Supply BLP models, drawing starting values from a uniform distribution (common to each simulation run) with support 50% above and below the true value. We use 100 $R_d$ sequence draws for numerical integration within the contraction mapping. For Demand-Only BLP, we estimate four models, given by
For Demand-Supply BLP, we estimate models b) and c) under similar settings. For the LTE estimator, we report uncertainty quantification based on the credible intervals of the sample in addition to computing the frequentist standard errors from a point estimate such as the posterior mean.
\paragraph{Relative Performance of Different Models}\\ For each model, we report the following statistics across all demand-related parameters (including all linear parameters, which are often of interest also): the mean bias, the Median Absolute Error (MAE), as well as the empirical coverage.
Our simulation results for the CUE confirm previous theoretical results from hansen1996finite,newey2009generalized,newey2004higher,donald2000jackknife. The CUE and LTE (which uses the CUE objective function) estimators have a lower bias but higher MAE (recall that $Corr(p_{jt},w_{jt}) \approx 0.2$ implies relatively weak instruments) compared to 2S-GMM. This is especially true with the inclusion of a supply-side. Of particular note is the substantially lower bias for $\theta_{1,price},\ \theta_{2,x}$ for the CUE estimator over the 2S-GMM estimator, which are often the key parameters of interest. This could be due not just to the inherently lower bias of the CUE, but also because the CUE formulation allows for the linear parameter on price $\alpha = \theta_{1,price}$ to be done within the objective function itself, which forces it into a feasible region.
Similar to hong2020blp, we find that the 2S-GMM estimators can have empirical coverage that is different from the nominal level. Our results differ in that we find under-coverage of the non-linear parameter $\theta_{2,x}$ rather than over-coverage. Similar to hong2020blp, using the LTE brings the empirical coverage close to the nominal level, but we find this using the credible interval rather than using the LTE as a point-estimate for frequentist uncertainty quantification. We also find that the CUE estimator produces similar results or slightly better results in terms of coverage, which suggests that the improvements in coverage could be due to the use of the CUE objective function rather than the use of LTE estimation. We find that using credible intervals from the sampler directly for non-linear parameters alongside frequentist uncertainty quantification for the linear parameters delivers the best results for LTE uncertainty quantification.
\paragraph{MCMC Diagnostics}\\ For each varaint of the simulation exercise, we report the average $\hat{R}$, the number of runs with $\hat{R} > 1.05$, the average $N_{eff}$, and the number of divergent chains.
Contrary to forneron2020inference, we find that the LTE estimated using NUTS-HMC when applied to the BLP problem has good MCMC diagnostics. Our $N_{eff}$ is high despite modest chain sizes, and the average $\hat{R}$ indicates a lack of evidence that our Markov chains fail to mix. We also have little to no divergent chains. The difference is likely due to our use of HMC-NUTS rather than RWMH for sampling.
We have shown how to phrase the BLP problem as an automatically differentiable moment function that is numerically equivalent to implementations in the literature based on for-loops. The use of automatic differentiaton allows us to implement more theoretically appealing estimators such as the CUE. Within the BLP setting, it also allows us to compute the linear parameter on price $\alpha$ “online” within the objective function itself since the weight matrix is no longer shared and we do not need to use multiple-equation GMM. This has resulted in better performance without the use of parameter bounds since $\alpha$ is forced into a feasible region. We anticipate that these techniques should be broadly applicable to other econometric models, particularly in situations where we have concentrated out parameters or sophisticated structural models.
The use of automatic differentiation has also enabled us to draw on the latest in machine learning estimation techniques. We provide intuition for why adaptive learning rate optimization methods are better suited to the BLP problem, and how NUTS-HMC sampling is superior to traditional RWMH for LTE estimation. Since automatically differentiable modifications preserve automatic differentiability, we envision that the robust estimation of BLP extensions and variants will now be much easier. For example, one can assume a dense covariance matrix on the random coefficients as opposed to a diagonal matrix as in the original formulation, or extend BLP to dynamic models of demand.
The implementation in admest outlined by this paper brings algorithmic improvements that better leverage modern hardware such as GPUs and TPUs for econometric estimation. We also incorporate hitherto neglected best practices for MCMC sampling that are applicable for all LTE-type estimation, and quasi- or full Bayesian approaches to econometric estimation more generally. The inclusion of adapative learning rate first order gradient descent methods should help with numerical stability for poorly behaved objective functions where traditional quasi-Newton methods tend to diverge. Once again, these advances should prove useful in broad classes of econometric estimation problems beyond BLP.
Our simulation exercises extend the work of conlon2020best,hong2020blp, and allow us to examine past theoretical work on the CUE and many weak moment conditions more generally hansen1996finite,newey2009generalized,newey2004higher,donald2000jackknife within the context of a non-linear IV-GMM with concentrated out parameters. Simulation results affirm that CUE-based frequentist and LTE estimation has a lower bias but higher MAE compared to traditional 2S-GMM methods. We also find that asymptotic normality is likely a poor approximation to the “medium sample” distribution of the BLP estimator, motivating the choice of the LTE over traditional frequentist uncertainty quantification, at least for the non-linear parameters that are directly sampled.
Future work can incorporate the bias correction of newey2004higher,chen2019mastering, which will not only improve the bias of the BLP estimator, but will also likely lead to empirical coverage closer to the nominal level. With automatic differentiation, the implementation of analytical corrections should be straightforward, since derivatives are “free”; and with advances in computation the time-costs of implementing empirical corrections should also decrease. The framework of automatically differentiable M-estimation can also be applied to structural estimation more generally, both in the field of industrial organization and beyond.
\printbibliography