Estimation of BLP models with high-dimensional controls

Estimation of BLP models with high-dimensional controls

\doublespacing

\thispagestyle{empty}

\setcounter{page}{1}

Introduction

The estimation of demand functions for differentiated products is a primary focus within Industrial Organization (IO) and has become increasingly influential across various fields of empirical microeconomics. A seminal work by berry1995automobile established a powerful framework for this purpose, demonstrating its utility in applications such as analyzing automobile markets. The BLP model has since been widely adopted in empirical industrial organization, demonstrating considerable versatility across a range of applications. For example, nevo2000mergers,nevo2001measuring applies the BLP framework to study demand in the ready-to-eat cereal industry. petrin2002quantifying uses the model to quantify the welfare effects of the introduction of the minivan in the U.S. automobile market. berry1999voluntary evaluate the impact of voluntary export restraints on the U.S. automobile industry, while goldberg2001evolution analyze price dispersion in the European car market.

A significant limitation in the practical application of the standard BLP model is the lack of clear guidance on how to select product characteristics in the demand specification. Researchers must decide which product characteristics to include, but there is no generally accepted, theory-driven criterion for doing so. In practice, the literature adopts a variety of approaches. Some studies rely on data-driven methods: for example, gillen2019blp use LASSO to select covariates in their analysis of Mexican elections. Others reduce dimensionality through statistical techniques, such as principal component analysis, as in backus2021common for the cereal industry. Still others choose a small set of characteristics based on ad hoc judgment, as in chidmi2007brand. While these strategies differ in implementation, they share a common feature: none provides a principled, theory-based criterion for variable selection.

This absence of guidance creates a fundamental trade-off. On the one hand, including many product characteristics risks triggering the curse of dimensionality, which can lead to imprecise estimates or even estimation failure. This concern is empirically relevant: gillen2019blp begin with over 200 demographic and social variables, while backus2021common consider 23 product characteristics in their study. On the other hand, restricting the model to a small set of variables—whether by ad hoc choice or aggressive dimensionality reduction—risks omitting important confounders, potentially introducing bias. The problem is further exacerbated in modern high-dimensional settings, where the number of potential covariates can exceed the number of observations. In such cases, traditional estimation methods become infeasible, and informal selection strategies become increasingly unreliable. For example, beauchamp2011molecular study a setting with over 360,000 variables but only 7,500 observations.

We present an approach to estimating and performing inference on the BLP model where the dimension of product characteristics is large—potentially much larger—than the number of markets. In the BLP model, the consumer utility from choosing product $j$ in market $t$ is given by

The utility equals the value a consumer gets from a product’s characteristics, $x_{jt}$, $p_{jt}$, and $\xi_{jt}$, and an individual-specific taste shock, $\epsilon_{ijt}$. A significant computational and statistical challenge arises when the dimension of the product characteristics vector, $d_x$, is large relative to the number of market observations, $T$. This challenge is further exacerbated when the linear term, $x_{jt}'\beta_0$, is generalized to a nonlinear function $f_{u0}(x_{jt})$. In such high-dimensional settings, traditional estimation methods, such as the GMM estimator of Berry, Levinsohn, and Pakes (2004), become impractical or infeasible, necessitating additional structure to enable informative inference. To resolve this, we construct a Neyman orthogonal estimator and then apply it to a high-dimensional BLP framework under approximate sparsity. Specifically, we posit that the nuisance parameters—including the unknown function $f_{u0}(x_{jt})$—can be represented by a sparse linear component $x_{jt}'\beta_0$ plus a negligible approximation error. The coefficient vector $\beta_0$ is assumed to be sparse, containing only a small number of non-zero elements relative to the sample size $T$, i.e. $\lVert\beta_0\rVert_0\ll T$. This assumption reflects the realistic empirical scenario where researchers have access to a vast set of potential product characteristics but lack prior knowledge of the small subset that truly explains consumer heterogeneity and market shares. Crucially, approximate sparsity ensures that the nuisance parameters can be estimated at a sufficiently fast rate, which is a prerequisite for establishing the asymptotic properties of the main parameters of interest. The first contribution of this paper is to develop a novel estimation and inference method for the BLP framework that accommodates potentially high-dimensional product characteristics, $x_{jt}$. Our approach integrates Neyman orthogonalization—developed by belloni2018high and chernozhukov2018double—into the BLP setting, allowing for estimation in the presence of high-dimensional nuisance parameters. The second contribution is to adapt this approach to settings with approximate sparsity. Specifically, we combine machine learning methods with orthogonal estimation in a two-step procedure. The estimation procedure is stated as follows:

We establish theoretical results showing that the proposed estimator is consistent and asymptotically normal. Importantly, our approach allows inference on any pre-specified, low-dimensional subset of parameters and maintains flexibility regarding the choice of machine learning estimators used in the first step. Finally, we provide evidence from simulation studies demonstrating strong finite-sample performance. Compared to naive plug-in estimators, our approach exhibits lower bias and improved rejection rates.

Literature Review. First, this paper contributes to the extensive literature on high/infinite-dimensional BLP models. berry2014identification establish the foundational conditions for the identification of the BLP model. Building on this framework, dunker2023nonparametric extend these results by providing conditions for the identification of the densities of the random coefficients. A significant strand of this literature treats the random coefficient density as an unknown function within an infinite-dimensional space. Notable contributions include lu2023semi and wang2023sieve, which employ sieve estimation, a widely used nonparametric method, to approximate the density of random coefficients. Similarly, compiani2022market characterizes the inverse demand function itself as an unknown function in an infinite-dimensional space and employs sieve estimation for its approximation. These works enable inference on functionals or coefficients in the linear utility component. Other recent developments, such as rafi2024nonparametric and singh2024choice, leverage nonparametric methods like Neural Networks to estimate functionals such as price elasticities; however, they often rely on simplified models that omit either random coefficients or unobserved product characteristics.

Several studies have specifically attempted to integrate high-dimensional characteristics into the BLP framework (e.g., liu2021double, rafi2024nonparametric, gillen2014demand, gillen2019blp, and sawada2020estimating). Nevertheless, these approaches frequently necessitate simplifying assumptions that may limit model realism—such as the omission of unobserved characteristics $\xi_{jt}$. More importantly, many of these methods (gillen2014demand, gillen2019blp, and sawada2020estimating) lack a formal theoretical foundation for conducting valid inference on the parameters of interest in a high-dimensional setting.

Our work is among the first to address the BLP model explicitly incorporating high-dimensional controls. We differ from this existing literature by maintaining parametric specifications for the random coefficient density while allowing for a nonlinear, unknown function of high-dimensional characteristics within the utility function. Under these assumptions, we develop an estimation and inference procedure that yields valid inference while including random coefficients and unobserved product characteristics.

Second, this paper connects to the literature on high-dimensional inference, notably the work on Neyman Orthogonalization and regularized methods by chernozhukov2018double and belloni2018high. The existing literature has primarily established the theoretical properties of these methods in linear or partially linear settings. Our primary contribution here is the successful application and extension of the Neyman orthogonalization and regularization approach to a highly nonlinear econometric setting, the BLP model. This demonstrates the flexibility and robustness of these techniques and resolves the key challenge of providing valid inference after model selection in a context where it was previously infeasible.

Organization of the Paper. The remainder of this paper is structured as follows. Section 2 discusses the theoretical underpinnings of the standard BLP model and its estimation methodology. Section 3.1-3.4 presents a high-dimensional extension of the BLP model along with the associated Neyman orthogonalization approach. Theoretical results regarding the consistency and asymptotic normality are provided. Section 3.5-3.6 applies the general theory to a specific BLP model under approximate sparsity, detailing the estimation procedure for the high-dimensional nuisance parameters and the construction of the Neyman orthogonal estimator for the parameters of interest. In Section 4, we conduct a simulation study to evaluate the performance of the estimator, and Section 5 concludes. Appendices provide detailed proofs for the theoretical results established in Section 3.

Notation. We use the following empirical process notations. Cadinality of a set $S$ is denoted by $|S|$. Given a vector $\beta$, and a set of indices $S$, $\beta_S$ is the vector that has the same value as $\beta$ on the indices in $S$ and zero elsewhere, i.e., $\beta_{S,j}=\beta_j$ if $j\in S$ and $\beta_{S,j}=0$ if $j\notin S$. $\mathbb{E}_J:=\frac{1}{J}\sum_{j=1}^J$. $\mathbb{E}_T:=\frac{1}{T}\sum_{t=1}^T$, $\mathbb{E}_{JT}:=\frac{1}{JT}\sum_{j=1}^J\sum_{t=1}^T$, and $\mathbb{E}_{T,L}:=\frac{1}{T}\sum_{l=1}^L\sum_{t\in I_l}$, where $I_l$ is the index subset of the $\{1,\dots,T\}$ and $T_l=|I_l|$. The $l_p$ norm of a vector $\left\lVerta\right\rVert_p=(\sum_i |a_i|^p)^{1/p}$. The $l_p$ norm of a matrix $\lVertX\rVert_p = (\sum_{i,j} |X_{ij}|^p)^{1/p}$. The operator norm of a symmetric matrix $X$ is $\lVertX\rVert_{op}=\sup_{\left\lVerta\right\rVert_2=1}\lVertXa\rVert_2$. $a_t\lesssim b_t$ means there exists a constant $C$ such that $a_t\leq C b_t$ for all $t$. $a_t\lesssim b_t\quad w.p.a.1$ means there exists a constant $C$ such that $P(a_t\leq C b_t)\to 1$ as $t\to\infty$. At last, for a function $f(w,\hat{\eta}(w))$, where $\hat{\eta}$ is an estimated function using observations independent of $w$, define $\lVertf(w,\hat{\eta}(w))\rVert_{L^2}=\left(\int \lVertf(w,\hat{\eta}(w))\rVert_2^2 F(dw)\right)^{1/2}$.

Standard BLP Model

Framework

Berry, Levinsohn, and Pakes (1995) begin with the following problem. Consider a market with $J$ competing products and an outside good, denoted as $0$. The vector of product characteristics of product $j$ in market $t$ will be denoted by $(\xi_{jt},x_{jt}',p_{jt})$. $\xi_{jt}\in\mathbb{R}$ represents the product characteristic which is not observed by the econometrician whereas $x_{jt}\in\mathbb{R}^{d_x}$ are observed. $p_{jt}\in\mathbb{R}$ is assumed to be an endogenous variable, e.g., price. Assume that the conditional mean of the unobserved characteristics is zero, i.e.,

In addition to those exogenous variables $x_{jt}$, we allow the existance of endogenous variables $p_{jt}$ in the sense of being realted to $\xi_{jt}$. This thus requires the use of instruments for identification. An individual $i$ in market $t$ chooses a product from a set of $J$ products plus an outside option. The utility of individual $i$ from choosing product $j$ in market $t$ is given by equation (ref). The utility of the outside option is normalized to

The idiosyncratic error term $\epsilon_{ijt}$ are independent and identically distributed (i.i.d.) extreme value random variables. $\alpha_i$ is a random variable as consumers are assumed to have different preferences. We assume that $\alpha_i=\alpha_0+b_i$, where $b_i$ follows a distribution $F(b_i,\sigma_0)$ that is determined by a parameter $\sigma_0$. So the model defines a map from the parameters $\theta_0=(\sigma_0,\alpha_0,\beta_0')\in\Theta$, where $\Theta$ is a compact space, and the vector of product characteristics $(\xi_{jt},x_{jt}',p_{jt})$ to the market shares $s_{t}\in\mathbb{R}^{J}$. Furthermore, we maintain mutual independence assumptions: $\alpha_i$, $\epsilon_{ijt}$, and $(\xi_{jt},x_{jt},p_{jt})$ are mutually independent. Consequently the utility can be rewritten as

We define the linear component of the utility $y_{jt}=x_{jt}'\beta_{0} + p_{jt}\alpha_0+ \xi_{jt}$. For each market $t$, we then define the corresponding market-level vectors: $y_t=(y_{1t},\dots,y_{Jt})'$, the price vector $p_t=(p_{1t},\dots,p_{Jt})'$, the unobserved characteristics vector $\xi_t=(\xi_{1t},\dots,\xi_{Jt})'$. The parameters to be estimated are $\theta_0$. Choice variable

The individual market share of product $j$ in market $t$, $s_{jt}$ and the corresponding vector of all shares in the market, $s_t$, are given by

This is sometimes referred to as random coefficient logit models. Berry, Levinsohn, and Pakes (1995) show that given fixed $\sigma$, $p_t$ and $s_t$, there is a unique $y\in\mathbb{R}^{J}$ that solves the following equation \footnote{Precisely, they show that given $p_t$, $s_t$, $x_t=(x_{1t},\cdots,x_{Jt})'$, $\beta,\alpha, \sigma$, there is a unique $\xi\in\mathbb{R}^J$ that solves the equation \[ s_t = f_s(p_t,x_t'\beta+p_t\alpha+\xi,\sigma). \] It is similar to show the uniqueness of $y$, as it can be seen as a special case where $\beta,\alpha=0$. } \[ s_t = f_s(p_t,y,\sigma) \] In other words, $f_s(p_t,\cdot,\sigma)$ is invertible and we can impute $y$ and consequently $\xi$ for each $\sigma$. We denote the inverse function as $f_s^{-1}(s_t,p_t,\sigma):\mathbb{R}\to\mathbb{R}^{J}$ and define imputed values of $y_t$ and $\xi_t$ as follows:

where $\theta=(\sigma,\alpha,\beta')$. Let $y_{jt}(\sigma)$ and $\xi_{jt}(\theta)$ be the $j$-th element of $y_t(\sigma)$ and $\xi_t(\theta)$, respectively. Obviously, these mappings map the true value $\sigma_0$ to the true values of linear component of the utility, $y_t$ and unobserved characteristics $\xi_t$.

Estimation

We here present the standard GMM estimator for the BLP model as proposed in berry2004limit. The data are observed at the market level. One typically observes $(x_t,p_t,s_t,z_t)$ for $t=1,\dots,T$. To address endogeneity, $z_t=(z_{1t},\dots,z_{Jt})'$ represents instrumental variables in market $t$. Then we define $\tilde{z}_{jt}$ as a function of $x_{jt}$ and $z_{jt}$. For example, a natural choice is $\tilde{z}_{jt}=(x_{jt}',z_{jt})'$. Then the estimator of $\theta_0$ is defined as \[ \hat{\theta}=arg\min\limits_{\theta}\hat{g}(\theta)'\hat{W}\hat{g}(\theta) \] where $\hat{g}(\theta)=\mathbb{E}_{JT}\xi_{jt}(\theta)\tilde{z}_{jt}$ and $\hat{W}$ is an estimator of a positive definite weighting matrix. In a two-step GMM estimation, for example, $\hat{W}$ is derived using the first step estimator, i.e., $\hat{W}=\mathbb{E}_T\left[\mathbb{E}_{J}\xi_{jt}(\tilde{\theta})\tilde{z}_{jt}\mathbb{E}_{J}\xi_{jt}(\tilde{\theta})\tilde{z}_{jt}'\right]$, where the preliminary estimator $\tilde{\theta}$ is derived from an initial unweighted GMM minimization.

The estimator developed by Berry, Levinsohn, and Pakes (2004) works poorly or even fails when $d_x$ is large relative to the sample size as the simluation studies show. Our proposed estimator addresses this issue and achieves the desired asymptotic normality.

BLP models with high-dimensional controls

Framework

Having many product characteristics is common in practice. For example, they may include the nutritional content of food products or the features of electronic devices, which can be numerous. This, however, creates a challenge for estimation and inference. In this framwork we generalize the utility function to the following form:

A key assumption that makes it possible to perform estimation and inference in such cases is the following:

We define the parameters $\theta_0=(\sigma_0,\alpha_0,f_{z0},f_{u0})$. Note that we allow the functions $f_{z0}$, $f_{u0}$, the distribution of the variables, and the parameters $\alpha_0$, $\sigma_0$ to vary with the sample size $T$, though we suppress the dependence on $T$ for notational simplicity. Combining equations (ref) and (ref), we can see that the model is a nonlinear extension to the high-dimensional linear IV models as studied by chernozhukov2018double. Researchers can apply any machine learning methods to estimate the functions $f_{z0}$ depending on the context.

In the following sections, we treat $(f_{z0},f_{u0})$ as nuisance parameters, while considering $\theta_{10}=(\sigma_0,\alpha_0)$ as structural parameters of interest. However, this designation is flexible: under the specification $f_{u0}(x_{jt})=x_{jt}\beta_0$, any subset of $\beta_0$ can be included as parameters of interest, with the remaining elements treated as nuisance parameters.

Estimation procedure

The estimation strategy adheres to a standard Neyman orthogonal approach and proceeds in two steps:

In the case where $\sigma_0$ is known, our model collapses to the high-dimensional linear instrumental variable (IV) framework described by Eqs. 4.5–4.6 chernozhukov2018double, with the Neyman orthogonal moment function correspondingly reducing to their Equation (4.7). \footnote{ They consider the partially linear IV model

and define the Neyman orthogonal moment function as \[ \psi(w,\theta,g,m) = (Z-m(X))\cdot(Y-Z\theta-g(X)) \] where $w=(Y,Z,X)$, $g$ and $m$ are the nuisance functions. } While our approach differs slightly from the 'partialling-out' Neyman orthogonal moment functions—such as those defined in Eq. 4.8 chernozhukov2018double \footnote{ Equation 4.8 defines the Neyman orthogonal moment function \[ \psi(w,\theta,l,m,r) = (D-r(X))\cdot (Y-l(X)-\theta(Z-m(X))) \] where $l_0(X) = \mathbb{E}[Y|X]$, $m_0(X)=\mathbb{E}[Z|X]$ and $r_0(X)=\mathbb{E}[D|X]$ are the nuisance functions. } or the procedure in Chernozhukov et al. (2015, Algorithm 1 chernozhukov2015post)—it maintains the Neyman orthogonality condition. We adopt this specific formulation due to its superior finite-sample performance in simulations and its flexibility in accommodating various choices of parameters of interest.

Neyman Orthogonal Estimation of $\theta_{10}$

Suppose we have obtained preliminary estimators $\hat{f}_z$, $\hat{f}_u$ that satisfy the desired preliminary convergence rates and we are interested in $\theta_{10}=(\sigma_0,\alpha_0)$, a subset of all parameters $\theta_0$. Let $w_t=(x_t,z_t,p_t,s_t)$ represents observations in market $t$. $(w_t)_{t=1}^{T}$ is modelled as independent and identically distributed. A natural approach to estimation of $\theta_{10}$ would be, for example, simplely plugging in the preliminary estimator $\hat{f}_u$ into the moment function and then minimizing the quadratic form of the empirical moment function, i.e., \[ \check{\theta}_1^{(1)}=\arg\min_{\theta_1}\hat{g}(\theta_1,\hat{f}_u)'W\hat{g}(\theta_1,\hat{f}_u) \] where $\hat{g}(\theta_1,\hat{f}_u)=\mathbb{E}_{JT}[\mathbb{E}_J\tilde{z}_{jt}\cdot(y_{jt}(\sigma)-\alpha p_{jt}-\hat{f}_u(x_{jt}))]$, $\tilde{z}_{jt}$ is a function of $x_{jt}$ and instruments $z_{jt}$, and $W$ is a positive definite weighting matrix. The estimator $\check{\theta}_1^{(1)}$ will generally have a slower than $1/\sqrt{T}$ convergence rate, as the argument in belloni2018high shows. The underlying reason is that in the proof of the GMM estimator's asymptotic normality, the remainder terms contain first-order effect (bias) introduced by the plug-in of the preliminary estimator $\hat{f}_u$. This issue can be addressed using the Neyman orthogonal approach, which is often combined with sample splitting to ensure that the remaining error terms vanish in probability. In this problem, we find that the remainder terms when using the Neyman orthogonal estimator also contain \[ \sqrt{T}\mathbb{E}_T[\mathbb{E}_J(\hat{f}_z(x_{jt})-f_{z0}(x_{jt}))\cdot\xi_{jt}]+\sqrt{T}\mathbb{E}_T[\mathbb{E}_Jv_{jt}\cdot (\hat{f}_u(x_{jt})-f_{u0}(x_{jt}))] \] The use of sample splitting along with the i.i.d. assumption allows simple and tight control of such terms. Now define the Neyman orthogonal moment function as

where $f_z$ and $f_u$ are the nuisance functions. The true values of the nuisance functions are $f_{z0}(x_{jt})$ and $f_{u0}(x_{jt})$. This Neyman orthogonal moment function is equivalent to the "partialling-out" moment function defined as

where $f_z,f_y,f_p$ are the nuisance functions. The true values of these nuisance functions are $f_{z0}(x_{jt})$, $f_{y0}(x_{jt})=\mathbb{E}[y_{jt}(\sigma_0)|x_{jt}]$ and $f_{p0}(x_{jt})=\mathbb{E}[p_{jt}|x_{jt}]$. It can be seen that $f_{y0}(x_{jt})-\alpha_0 f_{p0}(x_{jt})=f_{u0}(x_{jt})$. Furthermore it can be proved that the estimators derived from these two moment functions are asymptotically equivalent, in the sense that they yield the same asymptotic variance. We note that the Neyman orthogonal moment function $\psi(w_t,\theta_1,f_z,f_u)$ implicitly depends on the sample size $T$, as the underlying data-generating process evolves with $T$. For notational brevity, we suppress this dependence throughout the paper.

Partition the observation indices $\{1,\ldots,T\}$ into $L$ groups $I_l\; (l=1,...,L)$. For each $l$, construct estimators $\hat{f}_z^{l}$ and $\hat{f}_u^{l}$ using all observations not in $I_l$. Then plug in the nuisance parameter estimators $\hat{f}_z^{l}$ and $\hat{f}_u^{l}$ into the Neyman orthogonal moment function and derive the empirical Neyman orthogonal moment function as follows: \[ \psi_t(\theta_1,\hat{f}_z^{l},\hat{f}_u^{l})=\psi(w_t,\theta_1,\hat{f}_z^{l},\hat{f}_u^{l}) \] and \[ \mathbb{E}_{T,L}[\psi_t(\theta_1,\hat{f}_z^{l},\hat{f}_u^{l})]=\frac{1}{T}\sum_{l = 1}^{L} \sum_{t\in I_l}\psi_t(\theta_1,\hat{f}_z^{l},\hat{f}_u^{l}). \] Then we have the estimator for $\theta_1$ as

where $\hat{W}$ can be $I_2$ or estimated weighting matrix, $\Theta_\alpha, \Theta_{\sigma}$ are compact sets. Given that $\hat{f}_z^l$ and $\hat{f}_u^l$ have the convergence rate $\lVert\hat{f}_z^l(x_{jt})-f_{z0}(x_{jt})\rVert_{L^2}=o_p(T^{-1/4})$, and $\lVert\hat{f}_u^l(x_{jt})-f_{u0}(x_{jt})\rVert_{L^2}=o_p(T^{-1/4})$ and then following the techniques developed in Chernozhukov et al. (2018), the Neyman orthogonal estimator of $\theta_1$ can be proved to be asymptotically normally distributed. However, one may think of a more direct estimator below given equation ((ref))

Although this moment function satisfies Neyman orthogonality, this estimator does not incorporate sample splitting. Furthermore, similar to LASSO method in a linear model yields a convergence rate slower than $n^{-1/2}$, the estimator $\check{\theta}_1^{(2)}$ exhibits a slower rate of convergence than $\check{\theta}_1$.

Theory for the asymptotic normality of $\check{\theta}_1$

In this subsection, we provide a set of conditions, the nuisance parameter convergence rates and regularity requirements, that are sufficient for the asymptotic normality of Neyman orthogonal estimator $\check{\theta}_1$.

\noindentComment. Condition (ref), (ref), (ref), are standard regularity conditions in the BLP literature. Condition (ref) is a regularity condition on the density function of the random coefficients, ensuring the smoothness of the inverse function $y_{jt}(\sigma)$.

\noindentComment. Condition (ref) bounds the parameters and random variables in the model. One may replace the boundedness condition with conditions on moments, as stated in Condition SE and SM in belloni2014inference. Here we choose the boundedness condition for simplicity.

\noindentComment. Condition (ref) is the standard identification condition. Freyberger (2015) assumes a slightly weaker condition, but there is no significant difference. \footnote{Berry, Levinsohn, and Pakes (2004) and Freyberger (2015) assume a slightly weaker condition: For all $\delta>0$, $\exists C(\delta)$ s.t. \[ \lim P(\inf_{\theta\notin\mathcal{N}_{\theta_0}(\delta)}\left\lVert\mathbb{E}_Tg(\theta)-\mathbb{E}_Tg(\theta_0)\right\rVert_2\geq C(\delta))=1 \] However there is no significant difference. } Condition (ref) is the standard weighting matrix condition. Condition (ref) is equivalent to the statement that $\lVert\mathbb{E}[\psi(w_t,\theta_1,f_{z0},f_{u0})]\rVert_2\to 0 \implies \lVert\theta_1-\theta_{10}\rVert_2\to 0$. This is the high-dimensional extension of the standard identification condition in a fixed-dimensional setting and is commonly assumed in the high-dimensional literature, such as beyhum2024high. \footnote{ Specifically, beyhum2024high assumes in Assumption 1 that $\forall$ $\eta>0$, there exists $\epsilon>0$, s.t. $\inf_{\left\lVert\theta-\theta_0\right\rVert_2>\eta}R(\theta_0)-R(\theta)\geq\epsilon$, where $R$ is the population objective function. This is equivalent to condition (ref). } This condition is trivially satisfied when the distribution of the data does not vary with $T$ and the parameter space is compact.

\noindentComment. Unlike identification in fixed-dimensional contexts, Condition (ref) imposes the inequality uniformly over $T$. This requirement is crucial as both the underlying data distribution and the model parameters are allowed to vary with $T$. Condition (ref) (ref) are commonly assumed in the existing BLP literature. For example berry2004limit and freyberger2015asymptotic assume that $\Gamma=\frac{\partial}{\partial \theta}\mathbb{E}[g(\theta)]|_{\theta=\theta_0}$ has full rank, which is the same as condition (ref). freyberger2015asymptotic also assumes that the matrix $\mathbb{E} z_{jt}z_{jt}'\xi_{jt}^2$ is positive definite, which is the same as condition (ref). Condition (ref) depicts the local property of the Neyman orthogonal moment function around the true parameter value and is crucial for the asymptotic normality of $\check{\theta}_1$, whereas condition (ref) and (ref) depict the global property of the moment function and are crucial for the consistency of $\check{\theta}_1$.

\noindentComment. The convergence rates of the preliminary estimators are crucial for the asymptotic normality of $\check{\theta}_1$. This convergence rate requirement is standard in Neyman orthogonal estimation literature, such as chernozhukov2018double.

\noindentComment. Theorem (ref) establishes the consistency of the Neyman orthogonal estimator. This result is crucial for the estimator to be asymptotically normal. The consistency proof builds upon the standard framework of newey1994large, with slight adaptations to account for the existance of the first step estimators. Theorem (ref) establishes the asymptotic normality of the Neyman orthogonal estimator, which is the main result of this article. The asymptotic normality proof strategy builds upon the standard framework with slight adaptations to account for the existance of the nuisance parameter estimators.

Estimation of nuisance parameters $f_{z0}$, $f_{u0}$

The motivation for estimating $f_{z0}$, $f_{u0}$ is to mitigate the first-order bias that arises when the estimator $\hat{f}_u$ is plugged into a non-orthogonal moment function. Note that one can choose any machine learning method to estimate $f_{z0}$ and $f_{u0}$ depending on underlying data assumptions. In this subsection, we demonstrate a specific estimation approach for $f_{z0}$ and $f_{u0}$ utilizing $l_1$-regularization. To achieve the desired convergence rate, we impose approximate sparsity assumption on functions, $f_{u0}$ and $f_{z0}$ as well as the auxiliary function $f_{p0}$

To illustrate the estimation strategy under this approximate sparsity assumption, we first define parameter space, moment function and empirical moment function as follows:

Then we present the following approximate sparsity assumption on the functions $f_{z0}$, $f_{u0}$ and $f_{p0}$:

where $r_{jt}^p:=f_{p0}(x_{jt},z_{jt})-x_{jt}'\beta_{px0}-z_{jt}'\beta_{pz0}$, $r_{jt}^u:=f_{u0}(x_{jt})-x_{jt}'\beta_0$ and $r_{jt}^z:=f_{z0}(x_{jt})-\Pi_0x_{jt}$ are the approximation errors. This motivates the use of LASSO method to estimate $f_{p0}$ and $f_{z0}$, as studied in belloni2014high. In the food example, it is reasonable to assume that the instrument, such as own-cost-shifters in backus2021common, and the price can be well approximated by a small subset of the product nutritions.

Relying on the estimator $\hat{f}_p$, with a convergence rate of $T^{-1/4}$, and the approximate sparsity on $f_{u0}$, we can combine two separate $l_1$-penalized minimizations to obtain a preliminary estimator $\hat{f}_u$ with a convergence rate of $T^{-1/4}$.

Estimation of $f_{z0}$ and $f_{p0}$

Note equations ((ref)) and ((ref)) imply that instruments $z_{jt,i}$, where $i=1,\dots,d_z$, can be represented as a linear combination of product characteristics $x_{jt}$ plus the approximation error and the noise term. Similary, price $p_{jt}$ can be represented as a linear combination of product characteristics $x_{jt}$ and instruments $z_{jt}$ plus the approximation error and the noise term.

where $\Pi_i$ is the i-th row of $\Pi_0$, i.e., $\Pi_0=(\Pi_1,\dots,\Pi_{d_z})'$. We apply the approach described in belloni2014inference geared for non-Gaussian cases to each equation. With estimated $\hat{\Pi}_i$ and $(\hat{\beta}_{px}', \hat{\beta}_{pz}')$, we can construct the estimators of $f_{z0}$ and $f_{p0}$. More formally, for each $l=1,\dots,L$, consider the sample $I_l^c$ after the sample splitting, we define

and

The penalty term $\lambda_\Pi$ and $\lambda_{\beta_{px},\beta_{pz}}$ are chosen by the user. The simulation applies cross-validation to select $\lambda_\Pi$ and $\lambda_{\beta_{px},\beta_{pz}}$. Alternatively, the researcher may adopt the theoretically suggested values $\lambda_\Pi=C_{\Pi}\sqrt{d_0\log(d_x\vee T)/T}$ and $\lambda_{\beta_{px},\beta_{pz}}=C_{\beta_{px},\beta_{pz}}\sqrt{\log(d_x\vee T)/T}$, where $C_{\Pi}$ and $C_{\beta_{px},\beta_{pz}}$ are constants. Then define $\hat{\Pi}=(\hat{\Pi}_1,\dots,\hat{\Pi}_{d_z})'$. Consequently we have $\hat{f}_z^l(x_{jt})=\hat{\Pi}^lx_{jt}$ and $\hat{f}_p^l(x_{jt},z_{jt})=x_{jt}'\hat{\beta}_{px}^l+z_{jt}'\hat{\beta}_{pz}^l$.

Estimation of $f_{u0}$

This part combines the $l_1$-penalized minimization that motivates from the Restricted Minimum Distance (RMD) estimator in belloni2018high and LASSO method. The difficulty of this problem lies in the high-dimensionality and the nature of nonlinearity and noncovexity of the moment function $g(\theta)$ itself. The issue of high-dimensionality and nonconvexity are addressed by the penalized minimizer, while the nonlinearity issue is addressed by the LASSO method. The idea is to first obtain a well-performed estimator of $\sigma$, the "nonlinear parameter" in BLP literature, and then transforming this problem into an approximately linear one, which is addressed using LASSO.

To be more precise, in step A, the $l_\infty$ norm is used since it outperforms the $l_2$ norm in high dimensions due to several reasons. First, the $l_2$ norm aggregates squared deviations across all $d_{\tilde{z}}$ dimensions, causing the total noise floor to accumulate at a rate proportional to the dimensionality. In contrast, the $l_\infty$ norm focuses exclusively on the largest absolute deviation. Under mild conditions (such as sub-Gaussian noise), this maximum error grows only at a rate of $\sqrt{\log d_{\tilde{z}}}$, allowing for significantly tighter control of the estimation error as the dimension increases and thus enabling a faster convergence rate for the estimator. Second, the $l_2$ norm acts as an "averaging" operator, which allows errors from many irrelevant or noisy dimensions to be "smeared" across the entire parameter vector, potentially biasing the final estimate. The $l_\infty$ norm, however, acts as a uniform constraint. By ensuring that the worst-case component is minimized, it prevents any single dimension from being significantly off-track.

Since $p_{jt}$ is endogenous, we need to include instruments in $\tilde{z}_{jt}$ to ensure the identification. In practice, one can choose higher order terms of $z_{jt}$ and $x_{jt}$ to be included in $\tilde{z}_{jt}$, which is common in the BLP literature. Step B is intuitive as it is a LASSO regression with the "predicted price" $\hat{f}_{p}(x_{jt},z_{jt})$ as the replacement of the actual price, $p_{jt}$. More formally, for each $l=1,\dots,L$, consider the sample $I_l^c$ after the sample splitting. The estimation procedure for $f_{u0}$ is as follows:

The tuning parameters $\lambda_{\tilde{\theta}}$ and $\lambda_{\beta}$ are chosen by cross validation in the simulation. Alternatively, the researcher may choose the theoretically suggested values $\lambda_{\tilde{\theta}}=C_{\tilde{\theta}}\sqrt{\log(d_{\tilde{z}}\vee T)/T}$ and $\lambda_{\beta}=C_{\beta}\sqrt{d_0\log(d_x\vee T)/T}$, where $C_{\tilde{\theta}}$ and $C_{\beta}$ are constants.

Theory for the convergence rates of nuisance parameter estimators

In this section, we provide regularity conditions that are sufficient for the desired convergence rate of the nuisance parameter estimators, $f_{z0}$ and $f_{u0}$. Note that $\tilde{z}_{jt}$ is the vector function of $x_{jt}$ and $z_{jt}$ used in the moment function $g(\theta)$.

\noindentComment. Condition (ref) and (ref) impose boundedness on the eigenvalues, parameters and random variables to ensure the convergence of the constructed estimator of $f_{u0}$. One may replace the boundedness condition with conditions on moments, as stated in Condition SE and SM in belloni2014inference. Here we choose the boundedness condition for simplicity.

\noindentComment. Condition (ref) is the key assumption ensuring that the $l_1$-penalized estimators $\hat{f}_z$, $\hat{f}_u$, and $\hat{f}_p$ converge at the desired rate. Specifically, we require that the approximation errors vanish at the rate of $\sqrt{d_0/T}$, which matches the oracle convergence rate of the estimated coefficients that is achievable if the identities of the relevant controls were known a priori. This rate ensures that approximation errors do not dominate the overall estimation error when employing LASSO to estimate $f_{z0}$ and $f_p$. Such conditions are standard in the high-dimensional econometrics literature (see, e.g., belloni2011high, belloni2014inference).

The next condition concerns the behaviour of the Gram matrices \\ $\frac{1}{J|I_l^c|}\sum_{j=1}^J\sum_{t\in I_l^c}(f_p(x_{jt},z_{jt}),x_{jt}')'(f_p(x_{jt},z_{jt}),x_{jt}')$, $\frac{1}{JT_l}\sum_{j=1}^J\sum_{t\in I_l}(z_{jt}',x_{jt}')'(z_{jt}',x_{jt}')$,\\ and $\frac{1}{J|I_l^c|}\sum_{j=1}^J\sum_{t\in I_l^c}x_{jt}x_{jt}'$ for $l=1,\dots,L$, which are crucial for the convergence of the LASSO estimators. We say a semi-definite matrix $M$ satisfies the restricted eigenvalue condition over $S$ with parameters $(\kappa,\nu)$ if \[ \Delta'M\Delta \geq \kappa\left\lVert\Delta\right\rVert_2^2\quad \forall \Delta\in\mathbb{C}_{\nu}(S):=\{\Delta\in\mathbb{R}^p|\left\lVert\Delta_{S^c}\right\rVert_1\leq \nu\left\lVert\Delta_S\right\rVert_1\} \] Define the support indices for the parameters in the LASSO estimation as follows. Let $S_1$ be the support of the vector $(\alpha_0,\beta_0')$. Let $S_2$ be the union of the support of the vectors $(\beta_{pz0}',\beta_{px0}')$ and $(\beta_{pz0}',\Pi_i')$, $i=1,\dots,d_z$, i.e., $S_2 = \{j: \beta_{pz0,j}\neq 0\}\cup\{j+d_z: \beta_{px0,j}\neq 0\}\cup\left(\cup_{i=1}^{d_z} \{j+d_z: \Pi_{ij}\neq 0\}\right)$.

\noindentComment Condition (ref) is a standard assumption in LASSO literature (wainwright2019high). Note that Condition (ref) implies that $\frac{1}{J|I_l^c|}\sum_{j=1}^J\sum_{t\in I_l^c}x_{jt}x_{jt}'$ also satisfies the restricted eigenvalue condition over a new set $S_3$ with parameters $(\kappa_2,\nu_2)$ \footnote{ Specifically, $S_3=\{j: \beta_{px0,j}\neq 0\}\cup\left(\cup_{i=1}^{d_z} \{j: \Pi_{ij}\neq 0\}\right)$. }. If it is assumed that $f_{p0}(x_{jt},z_{jt})$ is linear in $z_{jt}$ and $x_{jt}$, then condition (ref) implies condition (ref) under mild conditions.\footnote{ For example, if $\lVert\beta_{px0}\rVert_1\ll \lVert\beta_{pz0}\rVert_1$. }

\noindentComment. Condition (ref) ensures the identifiability of $\sigma_0$, a property that holds in standard BLP models. Specifically, in a standard framework with fixed $d_x$ and linear $f_{u0}$, it can be shown that under mild regularity conditions (e.g., freyberger2015asymptotic), this condition is satisfied with the norm $\lVert\cdot\rVert_{\Theta_{\alpha,f_u}}=\lVert\cdot\rVert_{2}$. Furthermore, in a special case where $\left\lVert(\alpha,f_u)-(\alpha_0,f_{u0})\right\rVert_{\Theta_{\alpha,f_u}}\geq \lVert\beta-\beta_0\rVert_{2}$, Step B and estimation of $f_p$ may be omitted. In other words, $\tilde{\beta}$ already achieves the desired convergence rate. In addition, condition (ref) ensures there is a fixed $c$ such that the inequality holds regardless of the parameter space's dimensionality and thus $\sigma_0$ is identified, while condition (ref) guarantees the validity of the instrument $\epsilon_{jt}^z$ for identification of $\theta_{10}$, which includes $\alpha$. Consequently, neither condition implies the other.

The following theorem establishes the desired convergence rate of the preliminary estimators of the nuisance parameters under the suitable conditions.

\noindentComment. Theorem (ref) asserts that the preliminary estimators $\hat{f}_z$ and $\hat{f}_u$ converge at desired rates. This result plays the same role as the lemma 1 does in belloni2014high and is crucial for the Neyman orthogonal estimator to achieve the desired asymptotic normality. We set the tuning parameter $\lambda_{\Pi}=C_{\Pi}\sqrt{\frac{d_0\log(d_x\vee T)}{T}}$ which converges to 0 at a slower rate than $\lambda_{\beta_{pz},\beta_{px}}=C_{\Pi}\sqrt{\frac{\log(d_x\vee T)}{T}}$. The reason is the following.

When applying LASSO method with tuning parameter $\lambda_{\Pi}$ to estimate $\Pi_{i}$ in equation ((ref)), we must bound the convergence rate of the parameter estimation error, denoted by $\lVert\hat{\Pi}_{i}-\Pi_{i}\rVert_2$. In contrast, when applying LASSO method with tuning parameter $\lambda_{\beta_{px},\beta_{pz}}$ to estimate $\beta_{px0}$ and $\beta_{pz0}$ in equation ((ref)), we only need to control the convergence of the prediction error $\hat{f}_p-f_p$. Parameter recovery is a more stringent requirement because it demands that the estimated coefficients converge to the true values. Prediction, in contrast, can be accurate even if the model is misspecified (i.e., the true functional form differs from the assumed one), as long as the fitted values are close to the true values. The total error in parameter estimation arises from two sources: the noise term and the approximation error resulting from the model's functional form. To ensure that $\lVert\hat{\Pi}_{i}-\Pi_{i}\rVert_2$ converges to zero, the tuning parameter $\lambda_{\Pi}$ must be chosen large enough to dominate both sources of error. In particular, it must be sufficiently large to “kill” the approximation error—otherwise, it will not vanish, and the parameter estimates will not be consistent. If only prediction error matters, the condition is much weaker. A model can yield excellent predictions even when it is misspecified. Consequently the tuning parameter $\lambda_{\beta_{px},\beta_{pz}}$ can be chosen smaller without harming predictive performance.

Monte Carlo Simulation

The purpose of this section is to provide a comprehensive examination of the finite-sample performance of the proposed Neyman Orthogonal estimator. We achieve this through a Monte Carlo simulation, which is designed to evaluate key statistical properties, such as bias, variance, and coverage rates of confidence intervals. We emphasize that while our theoretical results establish the estimator's asymptotic properties, these simulations are crucial for validating its practical utility in the realistic setting of limited sample sizes. The data-generating process is based on the core model outlined in Section 2 and the parameter settings are inspired by the Monte Carlo studies in belloni2014high

Monte Carlo examples

The data generating process is based on the model:

Parameters of interest, $\theta_1=(\sigma_0,\alpha_0)$, is set to be $(1,-1)$. Nuisance parameters $\beta_0=2\cdot(1,(1/2)^2,\dots,(1/4)^2,0,\dots,0)$, a quartic decay pattern is assumed. $z_{jt}$ is set to be a linear function of $x_{jt}$, with $\Pi_{0,i}=(\underbrace{0,\dots,0}_{i},1,(1/2)^2,\dots,(1/5)^2,0,\dots,0)\cdot\frac{1}{2}$ also following a quartic decay pattern, where $\Pi_{0,i}$ is the $i$-th row of $\Pi_0$. $\beta_{p0}=(1.2,\underbrace{0,\dots,0}_{dim(x_{jt})-5},1.2,1.2,1.2,1.2)'$, also with a quartic decay pattern. Note that $(\beta_{px0}',\beta_{pz0}')'=\mathbb{E}[(x_{jt}',z_{jt}')'(x_{jt}',z_{jt}')]^{-1}\mathbb{E} p_{jt}\cdot (x_{jt}',z_{jt}')$. The exact values of $\beta_{px0}$ and $\beta_{pz0}$ are not crucial for the simulation and thus we do not specify them here. Assume $dim(x_{jt})=200$. The first element of $x_{jt}$, $x_{jt,1}$, is set to be 1, and the rest of the elements follow a uniform distribution, i.e., $x_{jt,k}\sim \sqrt{3}\cdot U[-1,1],\quad k\geq 2$ so that the variance of each element of $x_{jt,k}$ is 1. $\xi_{jt}\sim U[-1,1].$ We assume additionally that there are two variables $\eta_{jt,1}$ and $\eta_{jt,2}$ that link both $\epsilon_{jt}^p$ and $\epsilon_{jt}^z$. Suppose $\eta_{jt,1}\sim U[-1,1]$ and $\eta_{jt,2}\sim U[-1,1]$. Let $dim(z_{jt})=4$, i.e., there are four instruments. $(\epsilon^z_{jt,1},\epsilon^z_{jt,2})= (((\eta^z_{jt,1})^2,(\eta^z_{jt,2})^2)-\frac{1}{3})\cdot 1.34 $. $(\epsilon^z_{jt,3},\epsilon^z_{jt,4})= (\eta^z_{jt,1},\eta^z_{jt,2})\cdot 0.86 $. $u^p_{jt}=\xi_{jt} + \frac{1}{10}\cdot(x_{jt,2}^2-1)+\frac{1}{5}\cdot(x_{jt,3}^2-1)+\eta_{jt,1}+\frac{1}{2}\eta_{jt,2}+\frac{1}{5}(\exp(\eta_{jt,1})-\mathbb{E}\exp(\eta_{jt,1})) +\frac{1}{5}(\exp(\eta_{jt,2})-\mathbb{E}\exp(\eta_{jt,1}))$, where $\mathbb{E}\exp(\eta_{jt,1})$ equals $1.175201$. We set the number of goods in each market $J=4$ and the number of markets $T=50$. Define $\tilde{z}_{jt} = (x_{jt}',z_{jt}')'$. Tuning parameters are chosen automatically by cross-validation. Inference results are based on conventional $t$-tests with standard errors calculated using standard plug-in methods. The estimation employs sample splitting and the market indices are splitted into $6$ groups, i.e., $L=6$. Specifically, we split the $50$ market indices into $6$ different sets, $(1,\dots,8)$, $(9,\dots,16)$, $(17,\dots,24)$, $(25,\dots,32)$, $(33,\dots,40)$, and $(41,\dots,50)$.

We report results for six different estimators, which can be categorized as either infeasible benchmarks or feasible procedures. Note that in this high-dimensional simulation setting, directly minimizing the GMM objective function is computationally infeasible. The estimators are based on the empirical moment condition $\psi(\theta,w_{t},\Pi,\beta)=\mathbb{E}_J\phi_j(\theta,w_{t},\Pi,\beta)$, where $w_t$ encompasses all data in market $t$. Two of the procedures are infeasible oracle estimators: Oracle 1 and Oracle 2, which assume knowledge of the true nuisance parameters and are thus unavailable in practice. Oracle 1 uses the true $\beta_0$ and include a subset of $x_{jt}$ whose coefficients in $\beta_0$ are nonzero as instruments. Oracle 2 uses the true $\beta_0$ and $\Pi_0$, and considers the residual term $z_{jt}-\Pi_0 x_{jt}$, i.e., $\epsilon_{jt}^z$, as instruments. These estimators represent the optimal performance attainable within the Neyman Orthogonal framework and serve as our primary benchmarks. The other four procedures we consider are feasible. The preliminary estimator is an estimator obtained by minimizing the infinity norm of the empirical moment function, $\left\lVert\hat{g}(\theta)\right\rVert_\infty$. The non-Neyman orthogonal estimator is a non-orthogonalized estimator derived without applying the Neyman orthogonality condition. It serves to illustrate the performance improvement achieved through the Neyman orthogonal approach by directly plugging in the preliminary estimator of $\beta$ into the moment function and minimizing the $l_2$ norm. The Neyman orthogonal estimator and optimal weighting Neyman orthogonal estimator are the two feasible procedures proposed in this paper. We plug in the preliminary estimators of $\Pi$ and $\beta$ into the Neyman orthogonal moment function and minimize the $l_2$ norm. The optimal weighting Neyman orthogonal estimator further incorporates an estimated optimal weighting matrix into the estimation process. The detailed definitions of these estimators are provided below, and the results are summarized in Tables 1 and 2.

We summarize the simulation results for estimating parameter $\sigma$ in Table 1, reporting bias, standard error, root-mean-square-error (RMSE), and rejection rate of 5% level tests (Rej. Rate). As expected, the Oracle estimators, which represent infeasible benchmarks, demonstrate strong performance with relatively low RMSE (0.165, 0.377 for $\sigma$ and 0.155, 0.316 for $\alpha$) and rejection rates close to the nominal level (0.040, 0.030 for $\sigma$ and 0.060, 0.045 for $\alpha$). The Neyman Orthogonal procedures achieve the rejection rate of 0.045 for $\sigma$ and 0.07 for $\alpha$, indicating satisfied performance, though this comes at the cost of higher variance as reflected in their larger standard errors and RMSE values. Notably, the optimized Neyman Orthogonal variant performs similarly to the standard version. In contrast, the procedures without Neyman orthogonalization—both the Preliminary and No Neyman Orthogonal estimators—exhibit substantial bias with rejection rates of 0.610 and 0.425 respectively for $\sigma$ and 0.965, 0.915 for $\alpha$, significantly exceeding the nominal 5% level. This pronounced over-rejection problem in conventional approaches underscores the critical importance of employing Neyman orthogonalization methods when conducting inference in high-dimensional settings with nuisance parameters, even when those parameters exhibit sparsity and quartic decay patterns.

Conclusion

This study advances the estimation of demand in differentiated product markets by addressing the challenges of high-dimensional product characteristics within the Berry, Levinsohn, and Pakes (1995) framework. By integrating machine learning techniques into a Neyman orthogonal estimation framework, we provide a solution for settings where the number of potential characteristics may exceed the number of market observations. Our theoretical results establish that this approach achieves $\sqrt{T}$-asymptotic normality for key parameters of interest, such as price coefficients and heterogeneity, effectively insulating them from the slower convergence rates inherent in high-dimensional nuisance parameter estimation.

A central contribution of this research is the application of these techniques under the approximate sparsity condition. From an economic perspective, this condition aligns with the concept of bounded rationality, suggesting that while the space of product attributes is vast, consumers' purchasing decisions are driven by a sparse subset of key characteristics. Our Monte Carlo simulations confirm the practical utility of this framework, demonstrating its efficacy in finite samples compared to traditional methods that struggle in high-dimensional environments.

By bridging modern de-biased machine learning with structural Industrial Organization, this methodology offers a rigorous and adaptable toolkit for empirical researchers. Future work could extend this orthogonal estimation framework to multi-market or time-series contexts, further exploring how data-driven covariate selection can clarify consumer behavior in increasingly complex digital and physical marketplaces.

\addcontentsline{toc}{chapter}{Appendices}