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.
104,183 characters · 22 sections · 56 citation commands
Nested Pseudo-GMM Estimation of Demand for Differentiated Products
\thispagestyle{empty}
\setcounter{page}{1}
\setstretch{1.6}
The berry_levinshon_1995 (hereafter BLP) random coefficients discrete-choice model has become a workhorse of demand modeling in the empirical literature, as it accommodates rich substitution patterns between products and avoids the independence of irrelevant alternatives property inherent in standard multinomial logit models. The original estimation approach relies on a nested fixed-point algorithm, which combines the minimization of a generalized method of moments (GMM) criterion function in the outer loop with a fixed-point iteration in the inner loop to invert the demand system and recover mean utilities. This method has been widely adopted in empirical research, including nevo2000practitioner, petrin2002quantifying, bayer2007unified, and fan2013ownership, among many others.\footnote{See ackerberg2007econometric for a comprehensive survey, and conlon2020best for a discussion of best practices in implementing this algorithm.}
It is well known that the nested fixed-point algorithm is computationally inefficient for estimating the GMM parameters in the BLP model. While the outer-loop minimization steps are relatively inexpensive, the inner-loop fixed-point iterations are computationally intensive, as they involve inversion of a high-dimensional demand system. The primary source of inefficiency is that the algorithm fully solves the costly fixed-point problem for every candidate parameter vector, even when those guesses are far from the final estimate.
Several studies have proposed alternative methods to reduce the computational burden of estimating the BLP model. dube2012improving introduce a mathematical program with equilibrium constraints (MPEC) approach, which minimizes the GMM objective function subject to the market share equations treated as constraints. By jointly optimizing the Lagrangian over both the structural parameters and the mean utilities, the MPEC method requires solving the market share inversion only once—at the optimum. ABLP2015 (ABLP2015) propose a method that uses a linear approximation to analytically invert the market share equations, leading to a pseudo-GMM minimization problem over the model parameters. We refer to this approach as the ABLP estimator.
Inspired by the Nested Pseudo Likelihood (NPL) estimator of aguirregabiria_mira_2002,aguirregabiria_mira_2007 in the context of dynamic discrete choice models, we propose a Nested Pseudo-GMM (NP-GMM) method for estimating the BLP demand model. Our approach incorporates several features that substantially reduce the computational burden of GMM estimation. NP-GMM is a computationally motivated estimator that targets the same structural parameters under the same identification conditions as GMM. It trades a small asymptotic efficiency loss for large speed gains.
First, similar to ABLP2015 and in contrast to the standard nested fixed-point algorithm, our method swaps the order of the fixed-point iteration and the GMM optimization. This avoids the repeated—and computationally intensive—inversion of the demand system, at the cost of performing more, but less expensive, iterations in the minimization of a pseudo-GMM criterion. Second, in the inner loop, we treat the vector of individual outside-option choice probabilities as incidental objects (i.e., as sufficient statistics for competitive effects). Conditional on these probabilities, mean utilities can be recovered product by product using a closed-form expression: that is, the market-share–mean-utility inversion is both closed-form and, crucially, separable across products. Third, while the gradient of the full GMM objective function lacks an analytical expression, the gradient of our pseudo-GMM criterion function is simple and explicit. This facilitates efficient computation of both the objective function and its derivatives and makes our algorithm particularly well-suited for parallel processing and multithreading. Together, these innovations lead to a significantly faster estimation procedure. Importantly, while the defining feature in ABLP2015's ABLP method is a Taylor approximation / linearization of the fixed point mapping, our method instead uses an exact identity, not an approximation.
Our NP-GMM method outperforms alternative Nested Pseudo (NP) algorithms that more directly adapt the insights from aguirregabiria_mira_2002, aguirregabiria_mira_2007 to the BLP framework. Our approach extends these alternative methods in several important ways. First, rather than fixing the mean utilities as in a direct adaptation, we fix individual choice probabilities for the outside option to implement the pseudo-GMM step in the inner loop. This yields simple, closed-form expressions that are well-suited for parallelization and multithreading. Second, we show that our estimation problem can be reformulated as one in which the set of parameters is augmented with an auxiliary set, which is held fixed during the pseudo-GMM step. This recasts the problem into a structure similar to that of the original Nested Pseudo estimators.
We establish consistency and asymptotic normality of our NP-GMM estimator, providing a new asymptotic framework for NP estimation in high- and infinite-dimensional settings, where existing proofs do not apply. In many applications of NP methods—such as those in lin2024endogeneity and dearing2025efficient—it is standard to assume a finite number of states, which underpins the asymptotic arguments in aguirregabiria_mira_2002, aguirregabiria_mira_2007. However, the BLP model differs fundamentally from standard dynamic discrete choice models: the dimensionality of the mean utility vector grows with the sample size and tends to infinity. This renders conventional asymptotic arguments for NP estimators inapplicable. In contrast, our reformulation offers a novel framework for establishing the statistical properties of NP-GMM methods in such settings.
We conduct a series of Monte Carlo experiments to evaluate the computational efficiency and finite-sample performance of the proposed method. Our comparison focuses on the ABLP estimator, which is the most closely related approach and, to date, the fastest available method for estimating random-coefficients demand models using aggregate data. As shown by ABLP2015, the ABLP estimator delivers substantial computational gains relative to traditional nested fixed-point and MPEC methods. We then illustrate the implementation of our methodology using wine demand data from the Liquor Control Board of Ontario (LCBO), demonstrating the scalability and feasibility of estimating demand models with a large number of products.
A key advantage of our method lies in its compatibility with multithreading. While our algorithm benefits significantly from parallelization when computing the GMM objective function and its gradients, the ABLP estimator performs poorly under multithreading, particularly when evaluating gradients—a major component of its overall computational cost. As a result, ABLP gains little from parallelization. When comparing both methods using their respective optimal number of threads, we find that our algorithm evaluates the GMM objective function and gradients approximately five times faster than ABLP. Even for a moderate number of products, our proposed method is approximately twice as fast as the fastest currently available estimator, and its computational advantage grows more than proportionally with the number of products in the demand system.
Our paper also relates to a broader literature that develops computationally efficient estimators for variants of the BLP model. Assuming a discretization of consumer tastes and a reduced-form pricing function, fox2011simple propose a linear regression estimator subject to linear inequality constraints. kalouptsidi2012market also considers discrete consumer types and leverages a duality between market shares and consumer heterogeneity. More recently, lu2023semi introduce a semi-nonparametric approach and propose a two-step estimator based on sieve GMM. See also wang2023sieve and compiani2022market for related sieve-based methods for estimating nonparametric distributions of random coefficients in the BLP framework. In addition, several papers examine the numerical performance and implementation best practices of the standard BLP estimator and its alternatives; see, for example, knittel2014estimation and conlon2020best.
The rest of the paper is organized as follows. Section (ref) introduces the BLP demand model, the GMM estimator, and three estimation algorithms: the Nested Fixed Point method, the ABLP approach, and our proposed NP-GMM algorithm. Section (ref) presents Monte Carlo evidence comparing the performance of the ABLP and NP-GMM estimators. Section (ref) illustrates the implementation of the proposed algorithm and demonstrates its computational performance using wine demand data from the LCBO. Section (ref) summarizes and concludes. The code implementing the different estimators and the Monte Carlo experiments is available on GitHub at \url{https://github.com/luo8yao/np-gmm-blp}.
We consider a set of markets indexed by \( t \in \{1,2,\ldots,T\} \). Each market is populated by a mass \( M_t \) of consumers and offers a set of products \( \mathcal{J} = \{1,2,\ldots,J\} \). In each market, consumers choose one product \( j \in \mathcal{J} \) or opt not to purchase, which we denote by \( j = 0 \). Product \( j \) in market \( t \) is characterized by observed attributes \( \boldsymbol{x}_{jt} \) and an unobserved quality component \( \xi_{jt} \). The vector \( \boldsymbol{x}_{jt} \) contains \( K \) product characteristics, including price, and we let \( \boldsymbol{x}_t \) denote the collection \( (\boldsymbol{x}_{1t},\ldots,\boldsymbol{x}_{Jt}) \). The scalar \( \xi_{jt} \) captures market- and product-specific utility components observed by consumers but unobserved by the researcher; we collect these terms in the vector \( \boldsymbol{\xi}_t = (\xi_{1t},\ldots,\xi_{Jt}) \). Consumer $i$ in market $t$ obtains the utility from purchasing product $j$:
The utility of the outside good, the "no-purchase" option, is $u_{i0t} = \epsilon_{i0t}$. The parameter vector $\boldsymbol{\beta}_{it} = (\beta^{1}_{it}, \beta^{2}_{it}, ..., \beta^{K}_{it})^{\prime}$ contains consumer $i$'s (in market $t$) tastes for the $K$ characteristics; $\epsilon_{ijt}$ is an additional idiosyncratic shock that follows the Type I extreme value distribution.
For $k \in \{1,2,...,K\}$, each random coefficient $\beta^{k}_{it}$ has the following structure:
where $\beta^k$ and $\sigma^{k}$ are fixed parameters; $(\nu_{it}^{1}, \ldots, \nu_{it}^{K})$ are independent random draws from the standard normal distribution; $\boldsymbol{d}_{it}$ is a vector of observable demographics of consumer $i$ in market $t$, which may be measured at the consumer or market level; and $\boldsymbol{\pi}^{k}$ is a vector of parameters.
Let $\boldsymbol{\theta}$ be the vector with all the model parameters. It is convenient to distinguish two groups of parameters: $\boldsymbol{\theta}^{\prime} = \left( \boldsymbol{\beta}^{\prime}, \, \boldsymbol{\sigma}^{\prime} \right)$, with $\boldsymbol{\beta}^{\prime} = \left( \beta^{1}, \, ..., \, \beta^{K} \right)$ and $\boldsymbol{\sigma}^{\prime} = \left( \sigma^{k}, \, \boldsymbol{\pi}^{k}: \, k=1,2, ..., K \right)$. Let $\delta_{jt}$ be the average utility of product $j$ in market $t$: that is, $\delta_{jt} = \boldsymbol{x}_{jt}^{\prime} \, \boldsymbol{\beta} + \xi_{jt}$. Define the vector of average utilities for all the products in market $t$ as $\boldsymbol{\delta}_{t}$. The market share of product $j$ in market $t$ is represented as $s_{jt}$. We use function $\mathscr{s}_{j}(\boldsymbol{\delta}_{t}, \boldsymbol{x}_t, \boldsymbol{\sigma})$ to represent the prediction of the model about market share $s_{jt}$. That is:
where $F_{t}(\cdot)$ represents the distribution function of $(\boldsymbol{\nu}_{it}, \boldsymbol{d}_{it})$ in market $t$. This integral is typically evaluated via Monte Carlo simulation, using $N$ draws from the distribution $F_t$. Throughout the paper, we assume that market shares are computed using a finite number of such draws.\footnote{This can be interpreted either as assuming a finite population or that $N$ is sufficiently large for the simulation error to be negligible—an assumption that is nearly universal in this literature.} Accordingly, we treat the set of draws $\{ \boldsymbol{\nu}_{it}, \boldsymbol{d}_{it}: i = 1, 2, \dots, N \}$ as fixed and known to the researcher, and we use the following expression for the model-predicted market shares:
Define the vector of market shares $\boldsymbol{s}_{t} \equiv \left( s_{1t}, \, ..., \, s_{Jt} \right)$ and the corresponding vector-valued mapping $\boldsymbol{\mathscr{s}}(\boldsymbol{\delta}_{t}, \boldsymbol{x}_t, \boldsymbol{\sigma}) \equiv \left( \mathscr{s}_{1}(\boldsymbol{\delta}_{t}, \boldsymbol{x}_t, \boldsymbol{\sigma}), \, ..., \, \mathscr{s}_{J}(\boldsymbol{\delta}_{t}, \boldsymbol{x}_t, \boldsymbol{\sigma}) \right)$. The demand system, in vector form, is:
By Berry's Inversion Theorem (berry1994, berry1994), this mapping is invertible in $\boldsymbol{\delta}_{t}$. For any value of $(\boldsymbol{x}_t, \boldsymbol{\sigma})$, there exists an inverse function $\boldsymbol{\mathscr{s}}^{-1}(\boldsymbol{s}_{t}, \boldsymbol{x}_t, \boldsymbol{\sigma})$ such that $\boldsymbol{\delta}_{t} = \boldsymbol{\mathscr{s}}^{-1}(\boldsymbol{s}_{t}, \boldsymbol{x}_t, \boldsymbol{\sigma})$, or at the product level,
While the inverse function has a closed-form expression in the standard logit and nested logit models, this is not the case in the random coefficients logit model. In this case, the inversion must be computed numerically, which is computationally intensive—particularly in models with a large number of products.
Based on the inversion of the demand system, we can write this demand model as a partially linear regression model:
Price is included in \( \boldsymbol{x}_{jt} \) and is potentially correlated with \( \xi_{jt} \), so consistent estimation requires instrumental variables. The econometric model is completed by specifying the validity and identification power of instrumental variables. Specifically, we assume the existence of a vector of $q$ observable variables $\boldsymbol{z}_{jt}$ that satisfies the moment conditions $\mathbb{E}\left( \boldsymbol{z}_{jt} \, \xi_{jt} \right) = \boldsymbol{0}$. These conditions identify the parameter vector $\boldsymbol{\theta}$. A common choice for the instruments includes observable product characteristics other than the price of product $j$, along with the same characteristics for products that compete with product $j$.
Let $m(\boldsymbol{\theta})$ be the vector of $q$ sample moments that correspond to the population moments $\mathbb{E}\left( \boldsymbol{z}_{jt} \, \xi_{jt} \right)$.
The GMM estimator of $\boldsymbol{\theta}$ is defined as:
where $\boldsymbol{W}$ is a positive definite $q \times q$ weighting matrix.
The GMM estimator $\widehat{\boldsymbol{\theta}}$ satisfies the first-order conditions: $\frac{\partial G\widehat{\boldsymbol{\theta}})}{\partial \boldsymbol{\theta}} = 0$. berry_levinshon_1995 propose the Nested Fixed-Point (NFXP) algorithm that computes this estimator by applying Newton’s method to solve this system of nonlinear equations (ref). Specifically, it performs the following iterative procedure until convergence. At iteration \(n+1 \geq 1\):
At every Newton iteration, and for every market $t$ in the data, the NFXP algorithm fully solves the inversion problem that defines the vector of average utilities $\boldsymbol{\delta}_{t}$ as an implicit function of $(\boldsymbol{s}_{t}, \boldsymbol{x}_t, \boldsymbol{\sigma})$. Specifically, for a fixed value of $\boldsymbol{\sigma}$ (the one from Newton's iteration $n$), it applies the following fixed point iterations until convergence. At iteration \(\tau +1 \geq 1\):
Upon convergence of the fixed-point algorithm, the solution $\boldsymbol{\delta}_t$ represents the evaluation of the inverse function $\boldsymbol{\mathscr{s}}^{-1}(\boldsymbol{s}_t, \boldsymbol{x}_t, \boldsymbol{\sigma})$. This result can then be used to compute the criterion function $G(\boldsymbol{\theta}^{n})$ and its derivatives, enabling a new Newton iteration that updates the parameter vector to a new value $\boldsymbol{\theta}^{n+1}$.
In empirically common situations --especially when the outside good share is small-- the modulus of this contraction mapping is often close to one, which means the fixed-point algorithm may require many iterations to converge (dube2012improving, dube2012improving, conlon2020best, conlon2020best). The NFXP algorithm must solve $T$ such fixed-point problems at each trial value of $\boldsymbol{\theta}$ during the Newton iterations. These characteristics make the NFXP algorithm computationally intensive.
Before introducing our estimator, we begin by describing a direct application of the Nested Pseudo-Likelihood or GMM method from aguirregabiria_mira_2002, aguirregabiria_mira_2007 to the BLP model. This section serves two purposes. First, it introduces the core ideas shared by a broad class of Nested Pseudo-GMM (NP-GMM) estimators. Second, it demonstrates that a straightforward, off-the-shelf application of the standard NP-GMM approach—commonly used in dynamic discrete choice models—does not work in the context of the BLP model. These limitations motivate the modifications to this method that we present in Section (ref).
Let $\boldsymbol{\delta} = (\boldsymbol{\delta}_{t}: t=1, 2, ..., T)$ be the vector of mean utilities across all the markets and products in the sample. The method treats these mean utilities as incidental parameters to be estimated together with the structural parameters in $\boldsymbol{\theta}$. Specifically, we define an extended or pseudo-GMM criterion function that includes $(\boldsymbol{\delta}, \boldsymbol{\theta})$ as arguments:
A key feature concerns the algorithm used to compute the estimator $(\widehat{\boldsymbol{\delta}}, \widehat{\boldsymbol{\theta}})$ that satisfies conditions AM-1 and AM-2. This algorithm shares several characteristics with the NFXP algorithm: it is a nested procedure that alternates between minimizing the pseudo-GMM criterion with respect to $\boldsymbol{\theta}$ and performing iterations to solve for the mean utilities $\boldsymbol{\delta}_t$ implied by the demand system. However, unlike the NFXP algorithm, the roles of the inner and outer loops are reversed.
The algorithm proceeds as follows. It starts from an initial value for the vectors of means utilities: $\widehat{\boldsymbol{\delta}}^{0} = (\widehat{\boldsymbol{\delta}}^{0}_{t}: \, t \,= \,1, \, 2, \, ..., \,T)$. A natural way of initializing the mean utilities is using the ones that we would have if the model did not have random coefficients (i.e., $\boldsymbol{\sigma} = 0$) such that we had a standard logit model:
Then, at every (outer) iteration $n \geq 1$ we perform two tasks:
The main computational advantage of this algorithm over NFXP is that it avoids repeatedly solving the $T$ fixed-point problems. Instead, the NP-GMM algorithm solves these $T$ fixed-point problems only once -—at convergence. This efficiency comes at the cost of solving multiple pseudo-GMM minimization problems. However, these minimizations are computationally inexpensive because the pseudo-GMM criterion function is quadratic in $\boldsymbol{\theta}$ such that the pseudo-GMM estimator has a simple closed-form solution.
However, the approach described above suffers from a fundamental limitation. The pseudo-GMM criterion function in equation (ref) depends on the incidental parameters $\boldsymbol{\delta}$ and a subset of the structural parameters, $\boldsymbol{\beta}$, but not on the remaining structural parameters, $\boldsymbol{\sigma}$. As a result, these parameters cannot be identified or estimated within this framework.
In the next section, we introduce an alternative pseudo-GMM criterion function that depends on a set of incidental parameters distinct from $\boldsymbol{\delta}$. This new criterion enables the estimation of all structural parameters within a NP-GMM framework. Crucially, it retains the computational simplicity of the original pseudo-GMM approach.
Let $s_{0it}$ denote the market share of the outside good for consumer $i$ in market $t$. This consumer-level market share can be expressed as a function of the mean utilities $\boldsymbol{\delta}_{t}$ and $\boldsymbol{\sigma}$:
Combining the definition of $s_{0it}$ with equation (ref) for the model-implied market shares, we have:
That implies the following regression-like expression:
where $\boldsymbol{\lambda}_{t}$ represents the vector of consumer-level market shares for the outside good, $(s_{0it}: \, i=1, 2, \dots N)$, and:\footnote{Note that when $\boldsymbol{\sigma} =0$, the function $h\left( \boldsymbol{\lambda}_{t}, \boldsymbol{x}_{jt}, \boldsymbol{\sigma} \right)$ simplifies to $\ln(s_{0t})$, and the regression equation (ref) corresponds to the standard logit model.}
An interesting property of the vector $\boldsymbol{\lambda}_{t}$ is that—although it does not vary across products $j$—it captures all the information from the market shares of products other than $j$ that is required to compute the inverse demand function $\delta_{jt} = \mathscr{s}^{-1}_{j}(\boldsymbol{s}_{t}, \boldsymbol{x}_{t}, \boldsymbol{\sigma})$. Specifically, the model implies:
This property has an important implication for constructing the Pseudo-GMM criterion function used in our estimator: once we treat $\boldsymbol{\lambda}_{t}$ as a vector of incidental parameters, we no longer need to include $\boldsymbol{\delta}_{t}$ explicitly. Thus, conditional on $\boldsymbol{\lambda}_{t}$, recovering mean utilities $\boldsymbol{\delta}_{t}$ requires only evaluating the closed form expression $h\left( \boldsymbol{\lambda}_{t}, \boldsymbol{x}_{jt}, \boldsymbol{\sigma} \right)$ separately for each product, and does not require solving a $J$-dimensional fixed point.
Given the regression equation (ref) and treating $\boldsymbol{\lambda} =(\boldsymbol{\lambda}_{t}: t=1, 2, \dots T)$ as incidental parameters, we can define the following pseudo-GMM criterion function: {
}
In the remainder of this section, we present two main results. First, we explain why the NP-GMM estimator differs from the standard GMM estimator. Second, we establish the consistency and asymptotic normality of the NP-GMM estimator under identification assumptions that are no stronger than those required for the consistency of the GMM estimator.
To illustrate the difference between the GMM and NP-GMM estimators, it is useful to define two mappings that relate the incidental parameters $\boldsymbol{\delta}$ and $\boldsymbol{\lambda}$ to the structural parameters $\boldsymbol{\theta}$. We begin by defining Berry's fixed-point mapping $\Psi$:
This mapping is continuously differentiable in both $\boldsymbol{\delta}$ and $\boldsymbol{\sigma}$. By Berry's Inversion Theorem berry1994, for any fixed value of $\boldsymbol{\sigma}$, the mapping $\Psi\left( \cdot , \boldsymbol{\sigma} \right)$ is a contraction. Consequently, there exists a unique solution $\boldsymbol{\delta}$ to the fixed-point problem for each $\boldsymbol{\sigma}$. We denote this solution by $\boldsymbol{\delta}^\ast(\boldsymbol{\sigma})$. Given $\boldsymbol{\sigma}$ and the associated vector of mean utilities $\boldsymbol{\delta}^\ast(\boldsymbol{\sigma})$, we can then define the mapping $\boldsymbol{\lambda}^\ast(\boldsymbol{\sigma})$, which yields the consumer-level market shares for the outside good:
Taking into account the mappings $\delta^{\ast}$ and $\lambda^{\ast}$, we can establish a key relationship between the GMM and the NP-GMM criterion functions. For any vector $\boldsymbol{\sigma}$, we have:
This identity implies the following expression for the gradient of the GMM criterion function:
It is important to note that neither $\partial Q/\partial \boldsymbol{\lambda}^{\prime}$ nor $\partial \lambda^{\ast}/\partial \boldsymbol{\sigma}^{\prime}$ are generally equal to zero. As a result, the condition $\partial G/\partial \boldsymbol{\sigma}^{\prime} = 0$ does not imply $\partial Q/\partial \boldsymbol{\sigma}^{\prime} = 0$, and vice versa.
We are now in a position to formalize the distinction between the GMM and the NP-GMM estimators, denoted by $\widehat{\boldsymbol{\theta}}_{gmm}$ and $\widehat{\boldsymbol{\theta}}_{np}$, respectively. By construction, these estimators satisfy the following first-order optimality conditions:
Using equation (ref), it follows that the NP-GMM estimator does not satisfy the first-order conditions defining the GMM estimator, and vice versa. Therefore, the two methods yield distinct estimators.
Throughout the paper, we use the hat notation $\widehat{\cdot}$ to denote statistics and functions that involve sampling error, and the subscript $_{0}$ to indicate their population counterparts. Let $\boldsymbol{\theta}_{0}$ be the true value the structural parameters in the population, and let $\boldsymbol{\lambda}_{0}$ be the corresponding value for $\boldsymbol{\lambda}$, i.e., $\boldsymbol{\lambda}_{0} = \lambda^{\ast}(\boldsymbol{\theta}_{0})$. Remember that the mapping $\lambda^{\ast}(\boldsymbol{\theta})$ is deterministic, i.e., it does not incorporate sampling/estimation error. Define the population counterpart of the sample criterion functions, for GMM and NP-GMM estimation:
The following assumptions summarize the identification and regularity conditions needed for the consistency and asymptotic normality of the NP-GMM estimator. We study asymptotic properties of the estimators as either the number of products or the number of markets grows large, that is, as $JT \rightarrow \infty$.
Assumption (ref) imposes standard conditions ensuring correct model specification and parameter identification. At the true parameter values, the moment conditions are satisfied, and the corresponding criterion function attains its minimum at zero. For any parameter value different from the true one, the moment conditions are violated, and the criterion function remains strictly positive. Importantly, this identification assumption holds equivalently under both standard GMM estimation, based on $G_{0}$, and the NP-GMM framework, based on $Q_{0}$.
Assumption (ref) presents standard regularity conditions which are used to derive the asymptotic distribution of the NP-GMM estimator.
Proof: See Appendix (ref).
Proposition (ref) establishes the asymptotic distribution of the NP-GMM estimator. To present this result, it is useful to express the first-order conditions for the NP-GMM estimator as follows:
Here, $\xi_{jt}(\boldsymbol{\lambda}, \boldsymbol{\theta})$ represents $\ln \left( s_{jt} \right) - \boldsymbol{x}_{jt}^{\prime} \boldsymbol{\beta} - h\left( \boldsymbol{\lambda}_{t}, \boldsymbol{x}_{jt}, \boldsymbol{\sigma} \right)$, and $\boldsymbol{z}_{jt}^{*}(\boldsymbol{\lambda}, \boldsymbol{\theta})$ denotes the effective vector of instruments—equal in dimension to the number of parameters in $\boldsymbol{\theta}$—obtained by multiplying the original instrument vector $\boldsymbol{z}_{jt}$ by the weighting matrix $\boldsymbol{W}$ and the gradient vector of $\xi_{jt}(\boldsymbol{\lambda}, \boldsymbol{\theta})$ with respect to $\boldsymbol{\theta}$. We also define the following matrices:
Proof: See Appendix (ref).
The asymptotic variance of the GMM estimator is given by $\Omega_{\theta \theta}^{-1}$. In this model, the Jacobian matrix $\Lambda_{\theta}$ is nonzero, and it is straightforward to show that this implies that matrix $\boldsymbol{V}_{np} - \Omega_{\theta \theta}^{-1}$ is positive definite. As a result, the NP-GMM estimator is asymptotically less efficient than the standard GMM estimator. Nevertheless, as shown in our Monte Carlo experiments, the resulting loss in efficiency is negligible in practice.
In the existing literature on Nested Pseudo estimation methods—including aguirregabiria_mira_2002, aguirregabiria_mira_2007, as well as more recent applications such as lin2024endogeneity and dearing2025efficient—the asymptotic analysis is conducted under the maintained assumption that the number of incidental parameters is finite. This typically follows from a finite state space, which implies a fixed and finite number of conditional choice probabilities that does not grow with the sample size. In sharp contrast, the BLP model departs fundamentally from this setting: the dimension of the mean utility vector (the collection of $\delta$’s) increases with the number of products and therefore diverges as the sample size grows. As a result, existing proofs of consistency and asymptotic normality for NP-GMM or NPL estimators do not apply in this environment. This paper, specifically Propositions (ref) and (ref), fills this gap by establishing consistency and asymptotic normality of the NP-GMM estimator in a setting where the number of incidental parameters grows without bound as the number of products tends to infinity, thereby providing a new asymptotic framework for NP-GMM estimation in high- and infinite-dimensional settings.
This section presents the algorithm we use to compute our NP-GMM estimator. The algorithm is similar to the one described above for the Aguirregabiria-Mira method, but with two key differences. First, the pseudo-GMM criterion function $Q$ now depends on the full set of structural parameters. Second, before each pseudo-GMM iteration, we update the value of $\widehat{\boldsymbol{\lambda}}$ using its closed-form expression as a function of $(\widehat{\boldsymbol{\delta}}, \widehat{\boldsymbol{\sigma}})$.
NP-GMM algorithm:
Upon convergence, the algorithm provides the NP-GMM estimator defined by conditions NP-1 to NP-3.
Remark 1: The main computational advantage of this algorithm over NFXP is that it eliminates the need to repeatedly solve the $T$ fixed-point problems. Instead, the NP-GMM algorithm solves these problems only once, at the point of convergence. This efficiency comes at the cost of solving multiple pseudo-GMM minimization problems. However, these minimizations are relatively inexpensive, as the pseudo-GMM criterion has a simple closed-form expression and is globally concave in $\boldsymbol{\theta}$.\footnote{As with any existing method for estimating the BLP model, our algorithm yields a local optimum upon convergence, which may not be the global optimum if multiple local optima exist. The standard practice applies: the algorithm should be run from multiple starting values to assess the presence of multiple local minima. In such cases, the estimator is defined as the local optimum that achieves the lowest value of the GMM criterion function.}
In contrast to the Approximate BLP method proposed by ABLP2015—which we describe in Section (ref)—our criterion function does not rely on a linear approximation of the nonlinear demand system. Instead, it leverages the structure of the model and treats certain endogenous objects as incidental parameters. As our Monte Carlo experiments demonstrate, this difference has important implications for the relative performance of the two methods.
Remark 2: Step 1, which updates $\widehat{\boldsymbol{\lambda}}$, is relatively inexpensive and can be fully parallelized across $NT$ threads. The pseudo-GMM iteration in Step 2 is slightly more computationally demanding than in the standard Aguirregabiria-Mira version of the NP-GMM. This is because the criterion function is no longer quadratic, and we no longer have a closed-form expression for the estimator. However, the criterion function still has a simple analytical form and can be concentrated in $\boldsymbol{\sigma}$ by exploiting its quadratic structure in $\boldsymbol{\beta}$. Most importantly, the function is easy to evaluate, globally convex in $\boldsymbol{\theta}$, and has a closed-form expression for its gradient. As a result, standard gradient-based optimization methods can be used to reliably obtain the global minimum, i.e., the pseudo-GMM estimator.
The Newton–Kantorovich outer iteration in Step 3 can be replaced by a standard fixed-point iteration based on Berry’s mapping. We implemented this alternative version of the algorithm as well. In our experiments, the version using Newton iterations performs markedly better: it requires far fewer outer iterations to converge and reduces the incidence of non-convergence to essentially zero.
ABLP2015 proposed an estimation procedure for the BLP model parameters based on a first-order approximation of the log-demand system around an arbitrary vector of mean utilities, denoted by $\boldsymbol{\delta}_{0}$. For any value $(\boldsymbol{\delta}, \boldsymbol{\sigma})$:
where $\approx$ denotes a first-order Taylor series expansion, $\nabla_{\delta'} \equiv \partial/\partial \delta'$, and $ \nabla_{\delta^{\prime}} \ln\boldsymbol{\mathscr{s}}(\boldsymbol{\delta}_{0}, \boldsymbol{\sigma})$ is a $J \times J$ Jacobian matrix.
Thanks to this linear approximation, the inversion of the system of demand equations to obtain the vector of average utilities $\boldsymbol{\delta}$ can be approximated using matrix inversion or linear projection. We denote this approximation to the inversion problem using the operator $\Psi^{ablp}(\boldsymbol{\delta}_{0}, \boldsymbol{\sigma})$. That is:
The Approximate-BLP (ABLP) estimator is based on the following criterion function:
where $\Psi^{ablp}_{jt}(\boldsymbol{\delta}_{0}, \boldsymbol{\sigma})$ represents element $jt$ of the vector $\Psi^{ablp}(\boldsymbol{\delta}_{0}, \boldsymbol{\sigma})$.
Approximate BLP algorithm: ABLP2015 proposed the following algorithm to compute this ABLP estimator.
Upon convergence, the algorithm provides the ABLP estimator defined by conditions ABLP-1 and ABLP-2.
Remark 3: There is a fundamental difference between the Pseudo-GMM estimation in the ABLP algorithm and that in our NP approach. Each Pseudo-GMM step in the ABLP procedure requires repeatedly evaluating $\Psi^{ablp}_{jt}(\boldsymbol{\delta}_{0}, \boldsymbol{\sigma})$ —once for each trial value of $\boldsymbol{\sigma}$ —which, in turn, involves computing the matrix inverse $\left[ \nabla_{\delta^{\prime}} \ln \boldsymbol{\mathscr{s}}(\boldsymbol{\delta}_{0}, \boldsymbol{\sigma}) \right]^{-1}$ for each of those values. By contrast, the Pseudo-GMM step in our NP-GMM algorithm only requires evaluating simple, closed-form expressions for each product individually when computing both the criterion function and its gradient.
Remark 4: The distinction between ABLP and our NP-GMM algorithm becomes even more significant when parallel computing is introduced. A key advantage of our approach is the straightforward computation of the sample counterparts of the error terms $\xi_{jt}$, which can be efficiently parallelized across multiple threads. As demonstrated in our numerical experiments in Section (ref), this leads to substantial reductions in CPU time during the Pseudo-GMM estimation step.\footnote{Appendix (ref) derives the gradients of our Pseudo-GMM objective function, which have a notably simple form and are straightforward to vectorize.} In contrast, parallel computing yields negligible time savings for the Pseudo-GMM step in the ABLP method. This is because evaluating ABLP's criterion function and its gradients requires inverting matrices that depend on all products simultaneously, limiting the effectiveness of parallelization across products.
Table (ref) summarizes key features of the main algorithms used for the estimation of BLP demand models, and highlights the novel features of the NP-GMM method.
In this section, we use simulated datasets to compare the computational speed of our approach with that of the ABLP method. ABLP2015 show that the ABLP algorithm is significantly faster than both the NFXP and MPEC algorithms, especially in large-sample settings.
We adopt the same Monte Carlo simulation setting as in dube2012improving, which is also used by ABLP2015 to compare the ABLP and MPEC approaches.\footnote{All experiments are conducted on a system equipped with an Intel(R) Core(TM) i7-8850H CPU @ 2.60GHz, 16GB of RAM, running Windows 10 (64-bit) and MATLAB R2023a. For additional speed gains, we also provide Julia code implementing our method.}
We conduct a series of Monte Carlo experiments using samples in which the number of markets $T$ takes values $50$, $100$, or $500$, and the number of products $J$ takes values $25$, $50$, $100$, $200$, $400$, $800$, $1600$, and $3200$. In each market $t$, all $J$ products are present.
Let the deterministic component of utility be defined as: $\boldsymbol{x}_{jt}^{\prime} \boldsymbol{\beta} = \beta_0 + \beta_1 x_{1,j} + \beta_2 x_{2,j} + \beta_3 x_{3,j} + \beta_p p_{jt}$, where the vector of coefficients is given by $\boldsymbol{\beta}^{\prime} = (\beta_0, \beta_1, \beta_2, \beta_3, \beta_p) = (0, 1.5, 1.5, 0.5, -3)$. The three exogenous product characteristics -—$x_{1,j}$, $x_{2,j}$, and $x_{3,j}$-— are constant across markets $t$ and are drawn from a multivariate normal distribution given by:
In each market $t$, we simulate $N = 1,000$ individual-level taste draws, denoted by $\boldsymbol{\nu}_{it}^{\prime}$ $= (\nu_{0,it}, \nu_{1,it}, \nu_{2,it}, \nu_{3,it}, \nu_{p,it})$, which are independently drawn from a standard normal distribution. The scale parameters associated with these consumer-level random coefficients are $\boldsymbol{\sigma}^{\prime} =$ $(\sigma_0, \sigma_1, \sigma_2, \sigma_3, \sigma_p)$ $= (\sqrt{0.5}, \sqrt{0.5}, \sqrt{0.5}, \sqrt{0.5}, \sqrt{0.2})$. Note that this DGP abstracts from observed demographics, $\boldsymbol{d}_{i}$.
Each product $j$ has a market-specific vertical characteristic $\xi_{jt}$, which is independently drawn across products and markets from a standard normal distribution. Product prices vary across markets and are determined by the following equation:
where $\omega_{jt}$ is a marginal cost shock, independently drawn from a uniform distribution $U(0,1)$ across products and markets. Note that price $p_{jt}$ is an endogenous regressor in the demand estimation due to its correlation with the unobserved demand shock $\xi_{jt}$.
We address price endogeneity using instrumental variables derived from observable marginal cost shifters. Specifically, the DGP includes six observable cost shifters per product, denoted by ${w_{k,jt}}$, which are generated according to the following equation:
where $e_{k,jt} \sim i.i.d. \, U(0,1)$. By construction, each instrument $w_{k,jt}$ is strongly correlated with price $p_{jt}$ and independent of the unobserved component $\xi$.
The full vector of instrumental variables $\boldsymbol{z}_{jt}$ includes 42 variables, comprising: a constant term (1 variable); the product's own $x$ variables in levels, squares, and cubes (9 variables); the product's own $w$ variables in levels, squares, and cubes (18 variables); the product of all own $x$ variables (1); the product of all own $w$ variables (1); the cross-products between $x_{1,j}$ and each $w$ variable (6); and the cross-products between $x_{2,j}$ and each $w$ variable (6).
We begin by comparing the two methods in terms of both the computational speed of the algorithms and the statistical properties of the estimators. For a fair comparison, we apply the same convergence criterion and identical initial values to both algorithms. Convergence is assessed using the infinity norm of the differences between successive iterations of the vectors $(\widehat{\boldsymbol{\delta}}^{n} - \widehat{\boldsymbol{\delta}}^{n-1})$ and $(\widehat{\boldsymbol{\theta}}^{n} - \widehat{\boldsymbol{\theta}}^{n-1})$. Specifically, an algorithm is considered to have converged when
Since both NP-GMM and ABLP may yield multiple fixed points, we implement each algorithm using five distinct starting values—identical across the two methods—for each simulated dataset. For each algorithm, we select the fixed point that minimizes the value of the criterion function. To generate the initial values $\widehat{\boldsymbol{\delta}}^{0}$ and $\widehat{\boldsymbol{\sigma}}^{0}$, we follow the approach proposed by dube2012improving (dube2012improving) to construct five random starting points.\footnote{To generate a random initial value for $\widehat{\boldsymbol{\sigma}}^{0}$, dube2012improving (dube2012improving) use the formula: $\widehat{\sigma}_{k}^{0} = 0.5 \cdot |\widehat{\beta}_{k}^{logit}| \cdot \widetilde{U}_{k}$, where $\widehat{\boldsymbol{\beta}}^{logit}$ is the IV estimate of $\boldsymbol{\beta}$ from a standard logit model with instruments $\boldsymbol{z}_{jt}$. $\widetilde{U}_{k}$ is a random draw from a U(0,1). Given $\widehat{\boldsymbol{\sigma}}^{0}$, the initial value $\widehat{\boldsymbol{\delta}}^{0}$ is obtained by solving the fixed-point problem that defines the vector of mean utilities.} For each sample and algorithm, we define a successful estimation as one in which at least one of the five starting points converges.
Table (ref) and Figure (ref) compare the performance of the ABLP and NP-GMM methods using simulated samples with $T = 100$ markets and the numbers of products ranging from $J=25$ to $J=400$. Each method is implemented with its optimal number of computational threads.\footnote{See Section (ref) on the calculation of the optimal number of threads.} Across all values of $J$, NP-GMM achieves significantly lower wall-clock times than ABLP (Panel A).\footnote{Wall-clock time is the standard computer-science measure of the total time required to complete a task, as recorded by a real clock. It includes all delays, such as input/output operations, system interruptions, and waiting time. This contrasts with CPU time, which measures only the periods during which the processor is actively executing instructions.} This advantage exists even though NP-GMM requires more outer iterations (Panel D), more inner iterations (Panel E), and more criterion-function evaluations (Panel F). The source of this advantage is clearly shown in Panel B of Table (ref). ABLP incurs a substantially higher computational cost in each inner-loop GMM minimization because, at every trial value of $\sigma$, it repeatedly computes and inverts a Jacobian matrix. By conditioning on the incidental parameters $\boldsymbol{\lambda}$, the NP-GMM algorithm avoids this expensive step, leading to a far lower per-iteration cost and ultimately much faster overall computation despite requiring more iterations.
Figure (ref) shows how the performance gap between ABLP and NP-GMM methods widens more than proportionally as $J$ increases. An OLS regression of log–Time-per-Inner-Iteration on log-$J$ yields a slope of $0.96$ for the NP-GMM algorithm —indicating roughly linear growth— and a slope of $1.32$ for the ABLP algorithm. Consequently, the relative speed advantage of NP-GMM becomes increasingly pronounced as the size of the problem grows. For example, extrapolating these rates to $J = 12{,}800$ products yields predicted times per-inner-iteration of approximately $1076$ seconds for ABLP versus only $33$ seconds for NP-GMM.
Figure (ref) reports the average wall-clock time required by the NP-GMM and ABLP algorithms to evaluate their respective GMM criterion functions and gradients. The averages are computed over 1,000 simulated samples with $J=25$, $T=500$, and $N=1{,}000$. All evaluations are conducted using parallel computing, with the number of threads varying from 1 to 6. For the evaluation of the GMM objective function, NP-GMM and ABLP differ only in their computation of $\boldsymbol{\delta}$. NP-GMM relies on Equation (ref), which provides a simple closed-form mapping in $\boldsymbol{\sigma}$ for a given $\boldsymbol{\lambda}$. In contrast, ABLP uses Equation (ref), a mapping in $\boldsymbol{\sigma}$ that crucially requires matrix inversion for its calculation. Because $\boldsymbol{\delta}$ is a simpler closed-form function of $\boldsymbol{\sigma}$ under NP-GMM, our estimator is substantially faster than ABLP in computing both the GMM objective function and its associated gradients.
A further advantage of our market-share inversion is its ability to exploit parallel computing more effectively. The ABLP2015 MATLAB code uses the -singleCompThread flag or option to disable multithreading, in which case CPU time and wall-clock time are comparable across algorithms. When multiple threads are allowed, however, we observe substantial reductions in wall-clock time for evaluating the GMM objective function and its gradients under NP-GMM. By contrast, ABLP’s gradient evaluation does not benefit from multithreading and in fact deteriorates as the number of threads increases, as shown in Figure (ref).
To understand why the optimal number of threads in parallel computing is finite—and why this optimum is substantially larger for NP-GMM than for ABLP—it is useful to consider the following simple decomposition. Let $\ell$ denote the number of threads, $T(\ell)$ the total execution time using $\ell$ threads, $T_{\text{parallel}}$ the time required to perform tasks that can be parallelized, $T_{\text{serial}}$ the time required for inherently serial tasks, and $T_{\text{overhead}}(\ell)$ the overhead introduced by parallelization, which is increasing in $\ell$. Total execution time can then be written as
As the number of threads increases, parallel execution incurs overhead costs $T_{\text{overhead}}(\ell)$ that are absent in serial computation, including thread creation and management, synchronization, inter-thread communication, and operating-system scheduling. These overhead costs grow with the number of threads and eventually offset the gains from dividing the parallel workload across more processors. For example, if the overhead is linear in $\ell$, so that $T_{\text{overhead}}(\ell)=k_{0}+k_{1} \, \ell$, the optimal number of threads is finite and given by $\ell^{*} =\sqrt{T_{\text{parallel}}/k_{1}}$.
This framework also clarifies why the optimal number of threads is substantially larger for NP-GMM than for ABLP. Under NP-GMM, the computation of the $\delta$’s and the associated gradients relies on the closed-form expression in Equation (ref). These calculations are fully parallelizable and therefore contribute almost entirely to $T_{\text{parallel}}$. In contrast, under ABLP the computation of the $\delta$’s and gradients is based on Equation (ref), which requires matrix inversion. Such operations are only partially parallelizable and involve substantial serial components, implying that a larger fraction of the computation falls into $T_{\text{serial}}$. As a result, NP-GMM can efficiently exploit a larger number of threads before overhead dominates, leading to a higher optimal degree of parallelization than in ABLP.
Given this, we determine the optimal number of threads for each estimator to minimize overall evaluation time and apply it to ensure a fair comparison. Under the model setting $J = 25$, $T = 500$, and $N = 1,000$, the optimal number of threads is 2 for ABLP and 6 for NP-GMM. See Appendix (ref) for further details on optimal thread configurations under different model settings. Using these optimized settings, we employ MATLAB's profiler to track the wall-clock time for 1,000 evaluations of the GMM objective function and associated gradients, as shown in Figure (ref). NP-GMM demonstrates significant speed advantages in computing both the GMM objective function and its gradients. Overall, our NP-GMM algorithm is approximately five times faster than ABLP in evaluating the objective function and gradients. See Appendix (ref) for more details.
Tables (ref) and (ref) summarize the performance of the two estimators and algorithms in two Monte Carlo experiments: one with $J = 25$ and $T = 500$, and another with $J = 100$ and $T = 100$. Additional experiments with comparable sample sizes $J \times T$ but different combinations of $J$ and $T$ yield similar statistical properties for both estimators.\footnote{In our DGP, and thus in all our experiments, the number of random draws for consumer-level random coefficients is fixed at $N = 1,000$.} Each experiment is based on 20 simulated datasets. For each dataset, estimation is initialized from five random starting values, resulting in a total of 100 estimations per method. The maximum number of outer iterations is capped at 100, and the algorithm terminates once this limit is reached.
Table (ref) summarizes the computational and statistical performance of the two algorithms in these experiment. Both algorithms exhibit excellent convergence behavior. For each experiment and dataset, at least one run converged for each estimator. In the experiment with $J=25$ and $T=500$, out of 100 estimations per algorithm (20 datasets × 5 starting values), ABLP failed to converge in 2 cases (from a single dataset), while NP-GMM failed in 4 cases (from two datasets). In the experiment with $J=100$ and $T=100$, each algorithm failed to converge in one case.
Consistent with the results reported above, the wall-clock time required to implement NP-GMM is between two and three times lower than that of ABLP. This difference appears systematically across all datasets and initial parameter values.
The bottom panel of Table (ref) reports the statistical properties of the estimators. The results show that both estimators exhibit very similar finite-sample performance in terms of bias, variance, and mean squared error. The NP-GMM estimator performs slightly better in terms of RMSE in the two experiments.
Table (ref) reports detailed Monte Carlo results on the means and standard deviations of the parameter estimates under the two methods. For each parameter, the two methods exhibit very similar finite-sample biases and variances.
In this section, we use wine demand data from the LCBO to illustrate the implementation of the proposed methodology and to assess its computational performance relative to ABLP. The empirical results are intended to demonstrate the scalability and feasibility of estimating demand models with a large number of products. The LCBO is a provincially owned near monopoly that oversees the wholesale and retail distribution of alcoholic beverages in Ontario. A detailed description of the industry and the raw data is provided in aguirregabiria2023decentralized (aguirregabiria2023decentralized).
The dataset covers all 634 LCBO stores in Ontario and contains monthly sales data from October 2011 to October 2013, for a total of 25 months. LCBO classifies stores into six size categories --AAA, AA, A, B, C, and D-- from largest to smallest. We focus on the five AAA stores, which are the highest-volume outlets and together account for a substantial share of total LCBO wine sales. Restricting attention to these stores allows us to study demand in large, stable markets with a broad and consistently available product assortment.
We abstract from consumers’ store choice and treat each store-month observation as a separate market. Over the sample period, the LCBO sold 14,364 distinct wine products. Many of these products are offered only temporarily or in a limited number of stores, which creates substantial product turnover and sparsity at the store-month level. To avoid unbalanced product availability and to simplify the interpretation and estimation of the demand system, we restrict attention to the 76 store–month markets (out of 125) in which at least 3,000 products record positive sales. Within each selected market, we designate the products with the highest sales as inside goods and aggregate sales of all remaining wine products into the outside good. We fix the number of inside goods at $J=2900$, thereby trimming products with extremely small market shares. Our full working sample therefore consists of the top 2,900 wine products in each of the 76 markets, yielding a total of 220,400 product–market observations.
Product-level variables include standardized monthly sales (measured in 750-ml bottle equivalents), standardized price per 750-ml bottle, alcohol content, sugar content, and indicators for red and white wine, with other wine types serving as the omitted category. We define market size as twice the maximum monthly sales observed at each store over the sample period. Table (ref) presents summary statistics for our working sample and the variables in our model.\footnote{The LCBO also sells non-wine spirits which we exclude from our sample. We rely on the LCBO’s own classification of products as “wine”. Ten out of the 14,364 products that LCBO classifies as wines have alcohol content above 30%. One of these—Zwack Unicum Slivovitz 3-Year-Old—is sufficiently popular to appear in our estimation sample, which explains why the maximum observed alcohol content is 47%. All other wines in the sample have alcohol levels below 23%.}
In our empirical model, the product attribute vector $\boldsymbol{x}_{jt}$ includes a constant, alcohol content, sugar content, price, and indicators for red and white wine. The deterministic component of utility is:
We allow for random coefficients on the constant, alcohol content, sugar content, and price. In each market $t$, individual taste heterogeneity is captured by draws $\boldsymbol{\nu}_{it} = (\nu_{0,it}, \nu_{\text{alc},it}, \nu_{\text{sugar},it}, \nu_{p,it})$, which are independently drawn from a standard normal distribution. We simulate 1,000 consumers in each market. The corresponding scale parameters are denoted by $\boldsymbol{\sigma} = (\sigma_0, \sigma_{\text{alc}}, \sigma_{\text{sugar}}, \sigma_p)$. As in our Monte Carlo experiments, we abstract from observed demographics.
We treat price as endogenous. In addition to the product characteristics summarized in Table (ref), the data include information on each wine’s country of origin and a broad classification of wine style (e.g., dry or sweet). We construct BLP-style instruments using the average alcohol and sugar content of other wines from the same country of origin, as well as those of other wines within the same style category.
These instruments are motivated by cost-side and technological similarities among wines that share a country of origin or style. Alcohol and sugar content are closely linked to grape characteristics and fermentation technology, reflecting exogenous production conditions that systematically affect wine quality and pricing outcomes combris1997estimation,ashenfelter2008predicting. Accordingly, variation in the average alcohol and sugar content of competing wines within narrowly defined competitive groups acts as a cost shifter that predicts equilibrium prices, consistent with standard differentiated-product pricing models in which firms’ optimal prices respond to rivals’ cost conditions berry_levinshon_1995,nevo2001measuring. Empirically, these instruments are highly predictive of price, with a first-stage $F$-statistic of 1{,}903.98, suggesting that weak-instrument concerns are unlikely.
To study how the relative performance of the algorithms varies with the number of products, we estimate the model using different values of $J$, holding the set of markets $T$ fixed at the 76 identified above. We report results for the full sample with $J = 2{,}900$ products, as well as for smaller subsamples with $J = 800$ and $J = 1{,}600$ products, drawn at random from the full product set.
Both ABLP and NP-GMM benefit from parallelization as $J$ increases. Table (ref) in Appendix (ref) reports mean wall-clock time per evaluation of the objective function and its gradients using one to six threads. For NP-GMM, six threads consistently yield the fastest evaluation. For ABLP, five threads are optimal when $J = 800$, and six threads are optimal when $J = 1{,}600$ or $2{,}900$. These thread settings are used in all subsequent estimations.
Estimation is conducted with five random starting points, with both estimators initialized from the same starting values and updating the mean utility vector $\boldsymbol{\delta}$ via Newton's method in the outer loops. Table (ref) reports parameter estimates obtained using the two methods for each of the three samples. Across all samples and parameters, except for $\beta_0$ and $\sigma_0$, the estimates produced by the two methods are extremely close, with only negligible differences between them. Parameters $\beta_0$ and $\sigma_0$ are typically weakly identified in this setting due to limited variation in baseline preferences once observed characteristics and unobserved product quality are accounted for (see, e.g., dube2012improving and ABLP2015).
The near-identical point estimates produced by ABLP-GMM and NP-GMM across all configurations underscore that the results are not driven by algorithmic choices. Overall, the model delivers a plausible and robust characterization of wine demand at the LCBO. Across all product sets, the estimates paint a very coherent picture of consumer preferences in the LCBO wine market that is both economically intuitive and remarkably stable as the number of products grows.
Both alcohol and sugar content enter average utility with negative coefficients, suggesting that consumers on average dislike higher-alcohol and sweeter wines. Price has an economically meaningful negative effect on utility, generating substantial price sensitivity. The magnitude is stable across product counts. The positive coefficients on the red and white wine dummies indicate that, relative to the omitted “other wine” category, both red and white wines enjoy a significant utility premium. The random-coefficient estimates reveal important heterogeneity in tastes. There is meaningful heterogeneity in preferences for alcohol content and price, indicating that while the average consumer dislikes high alcohol and high prices, there is a non-trivial subset of consumers who are less price sensitive and who positively value stronger wines. By contrast, the estimated dispersion in sugar preferences is small, suggesting relatively homogeneous attitudes toward sweetness.
Figure (ref) shows total wall-clock times, illustrating that the performance gap between ABLP and NP-GMM widens disproportionately as $J$ increases.
Table (ref) reports computational performance metrics, highlighting the speed advantage of NP-GMM over ABLP. Although NP-GMM requires more objective function evaluations, its per-evaluation cost is substantially lower. As $J$ doubles, NP-GMM's per-evaluation time roughly doubles, whereas ABLP's increases by a factor of three. Across the three values of $J$, evaluations of the ABLP objective function and gradients are between 12 and 41 times slower than those of NP-GMM. Overall, NP-GMM is between 7 and 14 times faster than ABLP in this application.
This paper proposes a new computational approach to estimating the random-coefficients logit demand model that substantially reduces the computational burden of existing methods while preserving their desirable large-sample properties. The key insight is that, conditional on consumer-level probabilities of choosing the outside option, the BLP share inversion admits a closed-form representation that decouples the evaluation of unobserved product characteristics across products. Exploiting this representation allows us to reformulate estimation as a nested pseudo-GMM problem in which the most computationally demanding component of standard BLP estimation—the repeated numerical solution of the share inversion—is no longer required at each parameter trial.
This reformulation leads to an estimation algorithm with three important features. First, it replaces repeated fixed-point iterations or constrained optimization over mean utilities with simple, regression-like calculations that are analytically tractable and admit closed-form gradients. Second, it shifts computational effort from the inner loop to the outer loop of the algorithm, yielding a criterion function that is considerably cheaper to evaluate even if a larger number of outer iterations is required. Third, because the resulting objective and gradient evaluations are naturally separable across products and simulation draws, the algorithm scales particularly well with parallel computing.
Monte Carlo experiments and an empirical application show that these features translate into substantial computational gains at no cost in terms of finite-sample performance. Even in moderately sized demand systems, the proposed estimator is approximately twice as fast as the fastest currently available alternative, and its relative advantage grows more than proportionally with the number of products. Importantly, these gains do not come at the expense of statistical reliability: across all designs considered, the pseudo-GMM estimator exhibits sampling variability and finite-sample bias that are very close to those of fully iterated GMM estimators.
Beyond its immediate computational benefits, the proposed approach opens several avenues for future research. First, the idea of conditioning on auxiliary statistics that summarize competitive interactions may be useful in other models of differentiated-product demand, including extensions with richer heterogeneity, alternative choice sets, or dynamic consumer behavior. Second, the algorithm’s favorable scaling properties make it well suited for applications involving large product spaces, high-frequency data, or repeated estimation in simulation-based policy analysis.
\baselineskip14pt