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.
103,574 characters · 20 sections · 129 citation commands
Penalized GMM Framework for Inference on Functionals of Nonparametric Instrumental Variable Estimators
\addtocontents{toc}{\setcounter{tocdepth}{-1}}
Instrumental variables (IV) methods are widely used in applied research for estimation and inference in models containing endogenous regressors. In many cases, economic theory does not impose any functional form restrictions motivating nonparametric instrumental variables (NPIV) methods, where the function of interest is not assumed to be known up to a finite-dimensional parameter. In many cases, structural parameters of economic interest appear as functionals of that underlying unknown function. Examples are policy effects, average (weighted) partial effects, consumer surplus, measures of substitution patterns, and various counterfactuals from structural models. It is quite common for the estimation problem to be high-dimensional. There might be many control variables which we want to include in a flexible way along with the endogenous regressor, or a structural model may depend on many variables, e.g. in the demand for differentiated goods framework, the demand function depends on the vector of prices and product characteristics of all products in the market. In this paper, we are interested in estimation and inference on structural economic objects in presence of endogeneity when the dimensionality of the problem is (moderately) high.
Machine learning (ML) literature provides a collection of modern statistical tools for flexible estimation of various statistical objects that are especially powerful in high-dimensional settings. However, standard ML estimators, such as Lasso, boosting, or Neural Networks are unable to pick up causal relationships when endogenous regressors are present hartford2017deepIV. On the other hand, there is a new line of research in machine learning and computer science communities that offers a series of new algorithms that both addresses endogeneity and can be applied in high-dimensional environments, we refer to them as MLIV estimators. These algorithms are data-driven and exploit various forms of regularization to ameliorate the ill-posedness of the problem while maintaining the functional form flexibility. Examples include the DeepIV estimator hartford2017deepIV and its regularized version li2024regdeepiv, the Kernel IV regression singh2019kiv, the Dual IV regression muandet2020dual, the DeepGMM estimator bennett2019deepGMM, the Double Lasso estimator of gold2020, a series of estimators constructed using the minimax framework of dikkala2020minimax, the deep feature estimator xu2020learning, the boostIV estimator bakhitov2021boostIV, the minimax IV regression bennett2023minimax, and the two-stage ML estimator bruns-smith2025tsml. The goal of this paper is to use these novel methods to estimate and perform inference on various economic objects of interest that appear as functionals of the underlying structural function under endogeneity.
As standard ML algorithms, MLIV estimators produce inherently biased estimates. The main source of bias is regularization and/or model selection needed to balance out squared bias and variance to obtain overall small mean squared errors. In the NPIV context, regularization is particularly important as it plays a dual role. First, it allows to deal with the curse of dimensionality, as in the case of standard ML estimators. Second, it is necessary to solve the ill-posed problem. As a result, regularization and/or model selection bias leads to poor coverage unless it is corrected for. Furthermore, the bias term will propagate into the functional estimate if we simply plug-in an MLIV estimator into the functional formula. As chernozhukov2022adml point out, squared bias of plug-in estimators can shrink slower than the variance, leading to extremely poor confidence interval coverage.
In this paper, we provide an approach for performing valid asymptotic inference on functionals of MLIV estimators. Building on the automatic debiased machine learning framework of chernozhukov2022adml, hereafter CNS, we construct Neyman-orthogonal moment functions by adding the influence function adjustment term for the NPIV estimator derived by ichimura2022influence to the identifying moment conditions. The resulting debiased moment function has a zero derivative with respect to the MLIV estimator, ensuring insensitivity to local perturbations around the true value of the estimated function and allowing plug-in of noisy MLIV estimates without strongly violating the moment condition. We focus on regular functionals with a finite semiparametric asymptotic variance bound necessary for root-$n$ estimability, allowing for both linear and nonlinear functionals, though the conditions for root-$n$ rate are much tighter for the nonlinear case.
The influence function adjustment depends on the Riesz representer (RR) for the identifying moment function in the linear case, or the derivative of the identifying moment condition in the nonlinear case. In the NPIV framework, the closed-form RR is typically either very complicated to derive or even unknown. Our main methodological contribution is a penalized GMM (PGMM) framework that estimates the RR directly from the identifying moment conditions by exploiting their orthogonality, without requiring knowledge of the closed-form expression. The debiasing is therefore automatic in the sense that it depends only on the form of the identifying moment function but not on the form of the bias correction term. The PGMM estimator is novel and, to the best of our knowledge, is the only non-minimax automatic estimator of the RR in the NPIV framework, generalizing the Lasso minimum distance estimator of CNS to allow for a more general form of the influence function.
We derive the convergence rate for the PGMM estimator and provide conditions for root-$n$ consistency and asymptotic normality of the debiased MLIV estimator of the functional of interest. To accommodate for a large variety of MLIV estimators, we only require certain mean square consistency and convergence rates for MLIV estimators. The required conditions differ quite drastically for linear and nonlinear functionals. For linear functionals it is sufficient to require the MLIV estimator to converge at some positive rate in the projected mean square norm. It is well-known that NPIV estimators exhibit much faster convergence rates in the projected norm rather than in the standard mean square norm due to ill-posedness blundell2007, ChenPouzo2012, ChenPouzo2015QLR. However, for nonlinear functionals it is necessary to account for the linearization bias which requires the convergence rate to be faster than $n^{-1/4}$, which is a standard condition in the semiparametric literature newey1994. Moreover, the presence of nonlinearities in the identifying moment function results in the convergence rate condition in the standard mean square norm rather than the projected norm, which makes it harder to satisfy in practice.
In Monte Carlo experiments based on the average derivative functional with multiple endogenous regressors, we demonstrate that the PGMM-estimated Riesz representer performs on par with the known analytical expression. The plug-in estimator exhibits severely deteriorating coverage as the sample size grows, falling from 64% at $n=100$ to below 5% at $n=10000$. In contrast, the automatic debiased estimator achieves near-nominal coverage of 90--96% throughout, matching the performance of the analytical debiased estimator that uses the closed-form Riesz representer chen2023ann_npiv. Moreover, the automatic debiasing procedure exhibits greater numerical stability in small samples and high-dimensional settings, where the analytical approach can suffer from instability in the intermediate estimation steps required to evaluate the closed-form expression.
We apply our approach to estimate the mean own-price elasticity functional in the nonparametric demand for differentiated products framework berryhaile2016,compiani2022market that has been gaining popularity in the last years as an alternative to the standard parametric procedure of BerryLevinsohnPakes1995, hereafter, BLP. Building on the dimensionality reduction idea of gandhi_houde2019, we derive a semiparametric inverse demand estimation equation that remains tractable as the number of products grows. We then construct own-price elasticities as nonlinear functionals of the estimated inverse demand system using the implicit function theorem. To our knowledge, this is the first application of automatic debiased machine learning to nonlinear functionals in the endogenous setting.
In our Monte Carlo experiments, we demonstrate that the plug-in estimator exhibits small bias but severely undercovers, with empirical coverage ranging from 31% to 55%, reflecting that its standard errors dramatically understate the true sampling variability. In contrast, the debiased estimator achieves near-nominal coverage of 91--95% across all configurations, demonstrating that the automatic debiasing procedure delivers valid inference even for nonlinear functionals in moderately sized samples.
We use the Monte Carlo results as a basis for our empirical application where we estimate own-price elasticities using scanner data. We show that semiparametric elasticity estimates are $\approx 20\%$ larger in absolute value across products compared to the logit benchmark, reflecting richer substitution patterns. Moreover, the debiasing corrections vary across products in both magnitude and direction, with several products exhibiting corrections that substantially exceed the associated standard errors, underscoring the practical importance of the debiasing procedure for valid inference in empirical applications.
This paper connects several strands of literature. First, since the focus of the paper is functionals of nonparametric quantities, our methodology relates to the literature on semiparametric statistical theory van_der_vaart1991,bickel1993efficient,newey1994,robins_rotnitzky1995,van_der_vaart2000. These papers focus on functionals of densities or regressions in low dimensional settings, while in our paper we focus on functionals of MLIV estimators over domains that may include low, moderate, and high dimensional objects. A more recent work by chernozhukov2022locally generalizes and extends the insights from the classical theory by constructing Neyman-orthogonal moment conditions allowing for a wide range of ML estimators\footnote{chernozhukov2022locally provide high-level conditions for inference on functionals for conditional moment restriction models that nest the NPIV problem (see Theorem 19). Our results are complementary as we provide an estimator of the RR and give low-level conditions to derive its convergence rate.}. We follow chernozhukov2022locally and use Neyman-orthogonal moment functions with the influence function adjustment term for the NPIV estimator from ichimura2022influence.
Riesz representers are important objects in semiparametric theory as they appear in calculations of the asymptotic variance of functionals of nonparametric quantities ichimura2022influence,chernozhukov2022locally. For the same reason they appear in the influence function calculations, which makes estimation of RRs a cornerstone of the debiased machine learning literature. chernozhukov2022dml_gl and CNS propose Lasso and Dantzig minimum distance estimators of the RR based on the sparse approximation assumption. While the latter provides asymptotic results for regular functionals, the former provides finite sample analysis and also allows for irregular functionals. A recent paper by chernozhukov2022riesznet proposes to use neural networks and random forests to estimate the RR. On the other hand, chernozhukov2025adversarial take a different approach and allow for a more general estimator of the RR based on the minimax framework of dikkala2020minimax. ghassami2022minimax extend the doubly robust framework to settings where nuisance functions solve integral equations and propose a minimax kernel machine learning approach for their estimation, with applications to proximal causal inference. bennett2025inference use penalized minimax estimators for both the primary and debiasing nuisance functions and derive conditions for asymptotic normality of linear functionals that do not require closedness or strong identification of $\gamma$. Our work is complementary: while we focus on the NPIV framework, we provide an automatic PGMM-based approach to estimate the influence function adjustment and give conditions for valid asymptotic inference on debiased estimated for both linear and nonlinear regular functionals.
This work also contributes to the literature on estimation and inference on conditional restrictions models which nest the NPIV regression problem as a special case. Several NPIV estimators are now available including kernel-based estimators hall_horowitz2005,darolles2011npiv and series or sieve estimators newey_powell2003,blundell2007,ChenPouzo2012,chen2023ann_npiv. There are several papers focusing on linear regular functionals of NPIV estimators, see e.g., AiandChen2003, santos2011_IV, and severini_tripathi2012 among others. ChenPouzo2015QLR and chen_christensen2018optimal give conditions for pointwise and uniform asymptotic normality, respectively, of possibly nonlinear functionals of the sieve NPIV estimator. The results presented in this paper are complementary to the results on inference on functionals of NPIV estimators.
The paper also touches on the growing literature on flexible demand estimation in differentiated product markets. compiani2022market follows the nonparametric identification arguments of berry2014identification and demonstrates the performance of the NPIV estimator in a very simple case of two products with two characteristics. He uses Bernstein polynomials along with shape restrictions to alleviate the curse of dimensionality and nonparametrically estimate the inverse demand function. In that regard, our modeling approach is complementary: we leverage the dimensionality reduction idea of gandhi_houde2019 to construct a semiparametric inverse demand estimation equation that combined with our debiasing procedure allows the practitioner to apply it to more realistic settings. lu2023semi consider a similar framework, but they focus on applications with large amounts of products instead of large amount of markets. A closely related paper by singh2023choice models choice probabilities directly instead of demand inversion and applies the control function approach to account for price endogeneity. fosgerau2020inverse and monardo2021measuring consider a different class of inverse product differentiation models which generalize the inverse demand function of the nested logit model.
The remainder of the paper is organized as follows. Section (ref) briefly introduces the NPIV framework, discusses practical issues and MLIV estimators. In Section (ref), we describe the objects of interest and provide several economic examples. We also illustrate how to construct the debiased estimator and the estimator of its asymptotic variance. Finally, we introduce the PGMM estimator of the RR. Section (ref) gives conditions necessary to derive a convergence rate for the PGMM estimator. Section (ref) gives conditions for root-$n$ consistency and asymptotic normality of the debiased estimator for linear functionals. In Section (ref) we introduce additional conditions necessary to extend our results to nonlinear functionals. Section (ref) examines the performance of the debiased estimator in a simple Monte Carlo exercise. In Section (ref) we introduce the semiparametric demand estimation framework and estimate the mean own-price derivative functional using both simulated and real data. Section (ref) gives conclusions and provides possible extensions. All additional details and proofs are left for the Appendix.
Notation: For a vector $x \in {\mathbb R}^{n}$, let $|x|_{1}$, $||x||$, and $||x||_{\infty}$ denote its $\ell_{1}$-, $\ell_{2}$-, and $\ell_{\infty}$-norms respectively. For an $m \times n$ matrix $A$, we define $||A||_{\infty} = \max_{j,k}|A_{jk}|$. Let $||A||_{\ell_{\infty}} = \max_{i}\sum_{j=1}^{n}|A_{ij}|$ denote the induced $\ell_{\infty}$-norm of A. For $S \subseteq \{1,\dots,\,n\}$ let $x_S$ be the modification of $x$ that places zeros in all entries of $x$ whose index does not belong to $S$. For a random variable $X$, let $L_{2}(X)$ denote a space of all measurable and square integrable functions.
We start by briefly discussing the NPIV framework, consequences of ill-posedness for practitioners, and motivating the use of MLIV estimators.
Consider the nonparametric instrumental variables framework of newey_powell2003,
where $Y$ is an explanatory variable, $X$ is a vector of potentially endogenous regressors, $Z$ is a vector of instruments, and $\varepsilon$ is an error term. Suppose that $\gamma_{0}$ is identified and the completeness condition holds, i.e.\ for all measurable real functions $\delta$ with finite expectation,
Intuitively, this condition implies that there is enough variation in the instruments to explain the variation in the endogenous covariates. In the linear model, completeness reduces to the usual rank condition.
The unknown function $\gamma_{0}$ solves the integral equation
where $f$ denotes the conditional pdf of $X$ given $Z$. Solving for $\gamma$ directly is an ill-posed problem as it involves inverting linear compact operators kress1989. Ill-posedness implies that the solution to (ref) is not continuous in ${\mathbb E}[Y|Z]$ and $f(x|z)$: one cannot construct an estimator of $\gamma$ by simply plugging in consistent estimators of ${\mathbb E}[Y|Z]$ and $f(x|z)$ and solving for $\gamma$.
A standard solution is regularization, which means constructing an estimator of $\gamma_{0}$ so that ill-posedness does not affect consistency. In essence, regularization avoids estimation of higher-order terms that drive up variance. Traditional approaches include replacing $\gamma$ with a finite-dimensional approximation kress1989, Tikhonov regularization hall_horowitz2005,CFR2007regNPIV, and functional form restrictions chetverikov2017npiv.
Ill-posedness slows convergence rates of NPIV estimators relative to standard nonparametric regression. To quantify this, let $T:L_{2}(X)\mapsto L_{2}(Z)$ denote the conditional expectation operator,
and define the measure of ill-posedness as
where $\Gamma \subseteq L_{2}(X)$ and $||T(\gamma - \gamma_{0})||=\sqrt{{\mathbb E}\{{\mathbb E}[\gamma - \gamma_{0}|Z]\}^{2}}$ is the projected mean square norm. Assuming $\tau$ is bounded,
so convergence in the mean square norm is always slower than in the projected norm. Importantly, fast rates in the projected norm are achievable even when mean square rates are slow, since the projected norm sidesteps ill-posedness blundell2007,ChenPouzo2012,dikkala2020minimax.
Standard NPIV methods provide flexible approaches to nonparametric estimation under endogeneity. Nevertheless, ill-posedness poses several challenges. From the practitioner's standpoint, it limits what can be learned about $\gamma_{0}$: horowitz2011applied notes that only low-order approximation terms can be estimated with desirable precision, reflecting a fundamental characteristic of the estimation problem rather than a limitation of any particular method. Using a Gaussian example, newey2013nonparametric illustrates the connection between ill-posedness and instrument strength, showing that stronger instruments yield more precise estimates of higher-order terms.
The curse of dimensionality, which affects all nonparametric estimators, becomes more acute in the NPIV context. In severely ill-posed cases, convergence rates may be only polynomial in $\log n$ rather than $n$ blundell2007,darolles2011npiv,ChenPouzo2012. As a result, variance of NPIV estimators can be substantially higher than that of standard nonparametric regression estimators, even in moderate dimensions. bakhitov2025mliv provides Monte Carlo evidence demonstrating that variance of series NPIV estimators grows rapidly with the number of endogenous regressors, while the problem is much less severe under exogeneity.
One promising solution is to leverage machine learning algorithms designed specifically for the IV setting. Standard ML estimators, such as Lasso, boosting, and neural networks, fit the conditional expectation ${\mathbb E}[Y|X]$ rather than the structural function $\gamma$, and thus fail to capture causal relationships when endogeneity is present. However, recent literature has developed MLIV estimators that combine sophisticated regularization with the IV moment conditions to solve the ill-posed problem while maintaining functional form flexibility. Examples include DeepIV hartford2017deepIV,li2024regdeepiv, Kernel IV singh2019kiv, Dual IV muandet2020dual, DeepGMM bennett2019deepGMM, Double Lasso gold2020, minimax estimators dikkala2020minimax,bennett2023minimax, iterative estimators xu2020learning,bakhitov2021boostIV, and the two-stage ML estimator bruns-smith2025tsml. We refer the reader to bakhitov2025mliv for a detailed comparison of some of these methods and empirical evidence on their finite-sample performance.
In practice, the structural function itself is rarely the object of interest; rather, practitioners seek economically meaningful quantities such as average partial effects. Consider, for example, a demand estimation problem. The demand level itself does not bear a lot of economic meaning while objects like partial effects of demand shifters, consumer surplus, price and income elasticities or diversion ratios are potential objects of interest. These are functionals of the structural function. The goal of this paper is to provide valid inference on such functionals when $\gamma$ is estimated using MLIV methods.
This paper focuses on estimation and inference on functionals of a flexible (i.e. nonparametric) structural function $\gamma_{0}$ in presence of endogenous regressors, i.e. within the framework of the nonparametric instrumental variables model. Let $W_{i} \equiv (Y_{i}, X_{i}, Z_{i})$ be a data observation. Let $m(W, \,\gamma)$ denote a functional of $\gamma$ that depends on an observation $W$. We consider parameters of interest of the form
For expositional convenience, in this Section we will focus on functionals that depend linearly on $\gamma$. In Section (ref) we extend our results to nonlinear functionals. The object of interest $\theta_{0}$ is an expectation of some functional $m(W,\,\gamma_{0})$ over the data distribution. Hence, we are interested in mean effects, which restricts a set of possible functionals of interest, such as, for example, a simple evaluation functional $\theta_{0} = \gamma_{0}(\bar{X})$, where $\bar{X} \in \text{supp}(X)$. However, our framework is still general enough and covers a wide range of economically important objects.
Below, we give several examples of the types of objects under consideration, including both linear and nonlinear functionals.
Suppose that we are given $\hat\gamma$, an MLIV estimator of $\gamma_{0}$. A natural approach to estimate $\theta_{0}$ is to simply plug-in $\hat\gamma$ into $m$ and replace the expectation with the sample average,
That said, the plug-in estimator will not be root-$n$ consistent if the first-order bias does not vanish at root-$n$ rate, which is the case when $\hat\gamma$ involves regularization and/or model selection chernozhukov2022locally. In the NPIV model, regularization is essential to dealing with ill-posedness rendering all NPIV/MLIV estimators regularized estimators.
The plug-in estimator suffers from the first-order bias because the moment condition the moment condition defining $\theta_{0}$ is not orthogonal to local perturbations of $\gamma$ around $\gamma_{0}$. Namely, let $\delta$ be a local perturbation around $\gamma_{0}$, then the Gateaux derivative in the direction $\delta$ is
Thus, obtaining an orthogonal moment condition is a crucial step for establishing our results.
We consider functionals $m(W,\,\gamma)$ such that there exists a function $\alpha_{0}(Z)$ with ${\mathbb E}[\alpha^{2}_{0}(Z)] < \infty$ and
As discussed in ichimura2022influence, if there exists $v(X)$ with ${\mathbb E}[v^{2}(X)] < \infty$ and ${\mathbb E}[m(w,\,\gamma)] = {\mathbb E}[v(X)\gamma(X)]$, then the existence of $\alpha_{0}(Z)$ requires $v(X) = {\mathbb E}[\alpha_{0}(Z)|X]$. As pointed out in severini_tripathi2012, this is a necessary condition for root-$n$ estimability of $\theta_{0}$. Moreover, by the Riesz representation theorem, the existence of such $\alpha_{0}(Z)$ is equivalent to ${\mathbb E}[m(W,\,\gamma)]$ being a mean square continuous functional of $\gamma$. Henceforth, we refer to $\alpha_{0}(Z)$ as a Riesz representer. newey1994 shows that mean square continuity of ${\mathbb E}[m(W,\,\gamma)]$ is equivalent to the semiparametric efficiency bound of $\theta_{0}$ being finite. Thus, our approach focuses on regular functionals. Similar uses of the Riesz representation theorem can be found in ai_chen2007estimation, ackerberg2014asymptotic, hirshberg2020debiased, and CNS among others.
ichimura2022influence establish the form of the orthogonal moment function for NPIV estimators
where $\alpha(Z)[Y - \gamma(X)]$ is the influence function. Note that the moment function in (ref) is Neyman-orthogonal to local perturbations $(\delta,\,\beta)$ of $(\gamma_{0},\,\alpha_{0})$ such that
where the first two terms cancel out by the Riesz representation theorem and the last term is zero by the exogeneity condition. This property makes the orthogonal moment condition an excellent basis for constructing a debiased estimator of $\theta_{0}$ in the NPIV setting where estimators are typically regularized. Similar uses of the Neyman-orthogonal moment condition can be found in chen2023ann_npiv for NPIV sieve estimators and in gautier2011high for the high-dimensional linear IV regression.
Moreover, the exogeneity condition and iterated expectations imply
for any $\alpha(Z)$, meaning that the expectation of the influence function is zero regardless of $\alpha$. This implies
which allows us to use (ref) to estimate $\theta_{0}$. The debiased estimator $\hat\theta$ can be constructed by plugging $\hat\gamma$ and $\hat\alpha$ into the moment function $\psi(W, \theta, \gamma, \alpha)$ in place of $\gamma$ and $\alpha$ and solving for $\hat{\theta}$ from setting the sample moment $\psi(W, \theta, \hat{\gamma}, \hat{\alpha})$ to zero.
Note that the debiased estimator $\hat\theta$ requires an estimator of $\alpha_{0}$. Typically in the NPIV setting, the form of $\alpha_{0}$ is very complicated to derive or even unknown. Consider the weighted average derivative example from above. The RR is a solution to the following integral equation
where $f_{0}(X)$ is the marginal pdf of $X$. As a result, it is desirable to have a flexible approach for automatic estimation of the RR. The next subsection describes how to construct such an estimator.
chernozhukov2022locally show that we can exploit the orthogonality of the debiased moment function $\psi(W, \theta, \gamma, \alpha)$ to estimate $\alpha_{0}$. The Gateaux derivative of $\psi(W, \theta, \gamma, \alpha)$ in the direction $\delta$ is
where the last equality comes from $m(W,\gamma)$ being linear in $\gamma$. This can be thought of as a population moment condition for $\alpha_{0}$.
We assume that the Riesz representer estimator takes the form $\hat\alpha = b(Z)'\hat\rho$, where $b(Z)$ is a $p$-dimensional dictionary of basis functions with $p$ being possibly much larger than $n$. Let $d(X)$ be a $q$-dimensional dictionary of basis functions that represent deviations from $\gamma_{0}$. Using $d(X)$, we can construct a vector of moment conditions to estimate $\rho$. Let $d_{j}(X)$ be an element of $d(X)$, then we can form a sample moment condition corresponding to the population moment condition (ref) by replacing the expectation with a sample average and $\alpha_{0}(Z)$ with $b(Z)'\rho$ to obtain
Note that we require $q\geq p$ to ensure identification and estimability of $\rho$.
To allow for a high-dimensional $\alpha$ specification, we follow caner_kock2018 and use the penalized GMM (PGMM) framework. Let $\hat{\psi}(\rho) = (\hat{\psi}_{1}(W,\,\rho), \,\dots,\,\hat{\psi}_{q}(W,\,\rho))'$ where $\hat{\psi}_{j}(W,\,\rho)$ is defined in (ref). Then a solution to the PGMM problem takes the form
where $\hat{\Omega}_q = \hat{\Omega}/q$, $\hat\Omega$ is a $q \times q$ positive semi-definite matrix, and $2\lambda_{n}|\rho|_{1}$ is a penalty term. This framework allows for $q \geq p > n$ and can be viewed as a Lasso extension of the standard GMM.
Let $\hat{G} = \frac{1}{n}\sum_{i =1}^{n} d(X_{i})b'(Z_{i})$ and $\hat{M} = \frac{1}{n}\sum_{i =1}^{n} m(W_{i},\,d)$ be unbiased estimators of $G = {\mathbb E}[d(X)b'(Z)]$ and $M={\mathbb E}[m(W,\,d)]$, respectively. Then we can rewrite (ref) in matrix form as
The estimator $\hat\rho_{L}$ can be interpreted as a minimum distance version of the high-dimensional GMM estimator of caner_kock2018. Implementation details can be found in Appendix (ref).
The natural choice is to use the optimal weight matrix $\hat\Omega^{\text{opt}} = \left(\frac{1}{n}\sum_{i=1}^{n}\hat\psi_i(\tilde{\rho})\hat\psi_i(\tilde{\rho})'\right)^{-1}$, where $\tilde\rho$ is the preliminary PGMM estimate based on the identity weight matrix. However, in high-dimensional scenarios $(q > n)$ this is not an appropriate choice as $\hat\Omega^{\text{opt}}$ will be rank-deficient. Instead, we propose to use a diagonal weight matrix of the form
where $\hat\sigma^{2}_{j}(\tilde\rho) = \frac{1}{n}\Sigma_{i=1}^{n}\psi_{j}^{2}(W,\,\tilde\rho)$ is the sample variance of the $j$-th moment evaluated at the preliminary PGMM estimate. The only downside of $\hat\Omega^{d}$ is a slight efficiency loss when moments are correlated, but it still allows to handle conditional heteroskedasticity.
The estimation procedure can be summarized in the following pseudo-algorithm:
Next, we informally discuss the key conditions behind the asymptotic normality result. Since $\hat\theta$ is constructed by plugging-in $\hat\gamma$ and $\hat\alpha$ in the orthogonal moment condition, asymptotic properties of $\hat\theta$ depend on the asymptotic behavior of $\hat\gamma$ and $\hat\alpha$. First, to allow for a wide range of MLIV estimators, we assume that $\hat\gamma$ satisfies some projected mean square convergence rate condition as an estimator of $\gamma_{0}$. Specifically, we require
where $\kappa_{n}^{\gamma}$ can be slower than root-$n$ rate\footnote{The result also holds for the standard mean square rate condition, i.e. $||\hat\gamma - \gamma_{0}|| = O_p(\kappa_{n}^{\gamma})$, however, for NPIV/MLIV estimators this rate is slower due to ill-posedness.}. As pointed out in Section (ref), it is possible to obtain a fast rate under the projected mean square norm. Hence, it is a weak high-level assumption that can be satisfied by a variety of MLIV estimators such as Double Lasso gold2020, Kernel IV singh2019kiv, and a series of estimators constructed using the minimax framework of dikkala2020minimax among others.
The second condition is the mean square convergence rate of $\hat\alpha$. For the ease of exposition, assume that $\hat\alpha$ satisfies the following mean square convergence rate condition,
We derive an exact expression for $\kappa_{n}^{\alpha}$ in Section (ref).
Finally, under quite standard regularity conditions asymptotic normality can be established provided that
which is satisfied when $\sqrt{n}\kappa_{n}^{\gamma}\kappa_{n}^{\alpha} \rightarrow 0$. Hence, there is a trade-off between the convergence rates of $\hat\gamma$ and $\hat\alpha$. It is possible to allow for a slower convergence rate of $\hat\gamma$ at the expense of a faster convergence rate of $\hat\alpha$ and vice versa.
In this section we provide the mean square convergence rate for the PGMM estimator $\hat{\alpha}$ which is necessary for the asymptotic analysis of $\hat\theta$. We start by introducing some conditions.
The first part of Assumption (ref) is pretty standard and requires a consistent estimate of the weight matrix. The rate $O_p(\sqrt{\log q / n})$ is the natural high-dimensional extension of the classic $O_p(n^{-1/2})$ rate bickel2008covariance. The second part of the assumption, as discussed in caner_kock2018, might be restrictive as it requires a high-dimensional matrix to be uniformly bounded in $\ell_{\infty}$- norm, but for the notational convenience we keep it. Assumption (ref) is trivially satisfied by the identity weight matrix. We also demonstrate that the main result still follows though under the diagonal weight matrix and relaxed assumptions.
Note that the convergence rate of the PGMM estimator defined in (ref) depends on the convergence rates of $\hat\Omega$, $\hat G$, and $\hat M$. Assumption (ref) ensures that $\hat\Omega$ is consistent. To obtain a convergence rate for $\hat G$, we impose the following condition.
This condition implies
Unlike the standard Lasso, the second moment matrix convergence rate depends on the number of moments, i.e. the number of elements in $d(X)$, rather than the number of elements in $b(Z)$.
Next, we impose a sparse approximation condition for $\alpha_{0}$.
Intuitively, this assumption controls the squared approximation error from using the linear combination $b'\bar\rho$ to approximate $\alpha_{0}$. Note that Assumption (ref) does not necessarily require $\alpha_{0}$ to be equal to the linear combination of $\bar{s}$ terms, it states that there exists a sparse $\bar\rho$ with $\bar{s}$ non-zero elements such that the approximation error is bounded by $C\bar{s}\varepsilon^{2}_{n}$. In other words, Assumption (ref) is general enough to accommodate both exact and approximate sparsity of $\alpha_{0}$. Approximate sparsity allows for a large number of potential regressors (possibly much larger than the sample size) when relatively few important regressors give a good approximation but the identity of those few is not known, which is different from a standard series approximation where typically the first $\bar{s}$ regressors are assumed to achieve a good approximation bradic2021minimax. Thus, very sparse approximations allow to keep $\bar{s}$ relatively small which results in faster convergence rates. For a more detailed discussion of approximation bias conditions we refer the reader to CNS.
Let $S = \{1,\dots,\,p\}$, $S_{\rho}$ be a subset of $S$ with $\rho_{j} \neq 0$, and $S_{\rho}^{c}$ be the complement of $S_{\rho}$ in $S$. Let $\rho_{L}$ be the population coefficients, i.e.
The PGMM estimator $\hat\rho_{L}$ estimates the population coefficients $\rho_{L}$, which in turn might be different from the approximation coefficients $\bar\rho$. The following condition is essential to derive the oracle inequality for $\hat\rho_{L}$, and hence, the convergence rate for $\hat\alpha_{L}=b'\hat\rho_{L}$.
Assumption (ref) is the modified population restricted eigenvalue condition as in caner_kock2018. To accommodate for the PGMM estimator the condition is imposed on $G'\Omega_{q} G$ rather than ${\mathbb E}[b(Z)b'(Z)]$ as in the classic restricted eigenvalue condition of bickel2009. Showing that its empirical counterpart is bounded uniformly away from zero will be used to put a bound on the estimation error of $\hat\alpha_{L}$.
This condition is needed to put a bound on $||M||_{\infty}$ which is necessary to establish the oracle inequality for $\hat\rho_{L}$, and hence, the convergence rate for $\hat\alpha_{L}$. Moreover, note that by Assumption (ref), $\varepsilon_{n} = \varepsilon^{M}_{n} = \varepsilon^{G}_{n} = \sqrt{\log q/n}$. This simplifies the analysis, but is not necessary for establishing the results below. Also, let $|\bar\rho|_{1} \leq \bar{A} < \infty$. We can allow for the norm to grow with $n$ at a certain rate, however, it does not change the main results, hence, for simplicity we put a bound on $|\bar\rho|_{1}$.
The presence of endogeneity results in a slower rate of convergence for the RR estimator compared to the exogenous counterpart in CNS. The MD Lasso estimator of CNS converges at the $\sqrt{\bar{s}}\lambda_{n}$ rate, while the PGMM estimator is slower by a factor of $\sqrt{\bar{s}}$. Note that the convergence rate only depends on the number of approximation elements $\bar{s}$, but is independent of the number of relevant moments.
We relax the convergence rate condition in Assumption (ref) by allowing for a diagonal weight matrix. As a result, we have to increase the penalty to ensure $\bar{s}^{2}\lambda_n \rightarrow 0$. Otherwise, the main result remains the same.
In this Section, we provide conditions ensuring root-$n$ consistency and asymptotic normality of the debiased estimator $\hat\theta$. Under the specified conditions, we can do inference in a standard way. First, we focus on linear functionals and then provide additional conditions to extend the results to nonlinear functionals in Section (ref).
We impose the following conditions.
Assumption (ref) is purely technical, and we maintain it for simplicity. Assumption (ref) allows for estimators $\hat\gamma$ that are mean square consistent. Assumption (ref) requires $\hat\gamma$ to converge to $\gamma_{0}$ in the projected norm at a rate equal to $\kappa_{n}^{\gamma}$ which is typically slower than root-$n$. Note that this condition is weaker than convergence in standard mean square norm (see Section (ref)). This specification is general enough and allows for various MLIV estimators.
This condition is sufficient to guarantee $\sqrt{n}||\hat\alpha_{L}-\alpha_{0}||\,||T(\hat\gamma - \gamma_{0})|| \xrightarrow{p} 0$, leading to asymptotic normality of $\hat\theta$.
It is possible to extend the results from Section (ref) to allow for estimation of $\theta_{0} = {\mathbb E}[m(W,\,\gamma_{0})]$ for nonlinear $m(W,\,\gamma)$. The estimator is similar to the linear case except we estimate the RR of the linearization of $m(W,\,\gamma)$ leading to a different $\hat M$ needed. In this section, we show how to construct such an estimator and provide additional conditions that are sufficient for valid asymptotic inference for nonlinear functionals. As we mentioned in the introduction, due to nonlinearity of $m(W,\,\gamma)$, we have to impose restrictions on the convergence rate of $\hat\gamma$ in terms of the standard mean square norm, not the projected norm as in the linear case. We provide more details below.
To account for nonlinearity of $m(W,\,\gamma)$ in $\gamma$, we assume linearity of the Gateaux derivative of a nonlinear functional. To be more precise, let $\zeta$ be a deviation from $\gamma$. We assume that $m(W,\,\gamma)$ is Gateaux differentiable with the derivative $D(W,\,\gamma,\,\zeta)$, meaning that
for a scalar $\tau$, and that $D(W,\,\gamma,\,\zeta)$ is linear in $\zeta$. Moreover, assume that $\alpha_{0}(Z)$ satisfies
In other words, Equation (ref) implies that $D(W,\,\gamma,\,\zeta)$ is a mean-square continuous functional of $\zeta$, which corresponds to Assumption 3 of ichimura2022influence, meaning that $\alpha_{0}(Z)$ is a Riesz representer of the Gateaux derivative of $m(W,\gamma)$ with respect to $\gamma$ evaluated at $\gamma=\gamma_{0}$. Thus, by the Riesz representation theorem, for $D(W, \gamma_{0}, d) = (D(W, \gamma_{0}, d_{1}),\dots,\,D(W, \gamma_{0}, d_{q}))'$,
We can construct an estimator $\hat\theta$ exactly like in Equation (ref) except we need a different estimator of $\alpha_{0}(Z)$ based on (ref). Despite $\gamma$ enters $m(W,\,\gamma)$ nonlinearly, the estimator will still have zero first-order bias and be root-$n$ consistent and asymptotically normal under suficient regularity conditions. See newey1994, ichimura2022influence, and chernozhukov2022locally for more details.
An estimator $\hat{\alpha}_{\ell}$ can be constructed exactly as described in Section (ref) except being based on a different $\hat M_{\ell}$, where it is convenient to bring back the $\ell$ subscript. Let $\hat{\gamma}_{\ell,\ell'}$ be based on observations not in either $I_{\ell}$ or $I_{\ell'}$, then the unbiased estimator $\hat{M}_{\ell}$ is given by
where $\hat{M}_{\ell j}$ is the Gateaux derivative of the moment function with respect to $\gamma$ in the direction of the $j^{\text{th}}$ dictionary function. This estimator uses further sample splitting where $\hat{M}$ is constructed by averaging over observations that are not used in $\hat{\gamma}_{\ell, \ell'}$. This additional sample splitting allows $\hat{M}_{\ell}$ to depend on an estimator of $\gamma$ as required when $m(W,\,\gamma)$ is nonlinear in $\gamma$.
To establish the convergence rate for $\hat M_{\ell}$, we impose the following condition.
As CNS point out, the presence of the initial estimator $\hat\gamma_{\ell,\ell'}$ in $\hat M_{\ell}$ makes the convergence rate of $||\hat{M}_{\ell} - M_{\ell}||_{\infty}$ slower, $\kappa_{n}^{\gamma}$ instead of $\sqrt{\log q/n}$. Thus, $\varepsilon_{n} = \varepsilon^{M}_{n} = \kappa_{n}^{\gamma}$, which requires $\lambda_{n}$ to converge to zero slightly slower than $\kappa_{n}^{\gamma}$. Though it affects the convergence rate of $\hat\alpha$, the main results of Section (ref) still hold (see Appendix (ref)).
This condition controls the size of the linearization remainder in a linearization using the Gateaux derivative. It implies that ${\mathbb E}[m(W,\,\gamma)]$ is Frechet-differentiable in $||\gamma - \gamma_{0}||$ at $\gamma_{0}$ with derivative ${\mathbb E}[D(W,\,\gamma_{0},\,\gamma-\gamma_{0})]$.
It is a standard assumption to accommodate for nonlinearity of $m(W,\,\gamma)$. This might be a very tight restriction to satisfy given overall slow convergence rates of NPIV estimators, especially in the severely ill-posed case. However, as discussed in CNS, it is not known whether it is possible to weaken the $n^{-1/4}$ condition for nonlinear functionals, which goes back to newey1994.
This section presents Monte Carlo evidence illustrating finite sample performance of our automatic debiasing procedure. We focus on the average derivative functional and compare three approaches: (i) the plug-in estimator, (ii) debiased ML using the analytical Riesz representer, and (iii) our adaptive debiased ML using the PGMM-estimated Riesz representer.
Our design builds on newey_powell2003, santos2012inference, and ChenPouzo2015QLR, modified to allow for multiple regressors and instruments. We generate i.i.d.\ draws
where $k$ denotes the number of regressors and instruments. The response variable is
The composite error structure implies that endogeneity weakens as $k$ increases, since each regressor $X_j$ is correlated with only one component of $v_i$. We therefore restrict attention to $k \leq 10$.
Following bakhitov2025mliv, we specify the structural function as
where $X_{-1} = (X_2, \ldots, X_k)'$. This specification ensures that the parameter of interest,
is constant across simulations.
We compare three estimators. The plug-in (PI) estimator averages the estimated derivative without bias correction. The debiased ML (DML) estimator employs the analytical Riesz representer, specifically, the identity score estimator of chen2023ann_npiv. The automatic debiased ML (ADML) estimator uses our PGMM framework to learn the Riesz representer directly from the orthogonality conditions.
For all methods, we construct dictionaries $b(Z)$ and $d(X)$ using cubic polynomials with interactions, yielding $p = q$ basis functions. The structural function $\gamma$ is estimated using the Double Lasso estimator of gold2020. We report results from $1000$ replications for $k \in \{2, 5, 10\}$ and $n \in \{100, 500, 1000, 10000\}$, using five-fold cross-fitting throughout. Table (ref) reports absolute bias, median standard error, and coverage of a nominal 95% confidence interval.
The plug-in estimator exhibits substantial bias across all specifications. More strikingly, the coverage probability deteriorates as the sample size increases: for $k=2$, coverage falls from $64.3\%$ at $n=100$ to $4.6\%$ at $n=10000$; for $k=10$, it drops from $54.4\%$ to $0.4\%$. This reflects the well-known phenomenon that the squared bias of plug-in estimators shrinks slower than the variance chernozhukov2022locally, causing the standard errors to increasingly understate the true uncertainty.
The analytical debiasing procedure performs well in large samples but exhibits instability when $n$ is small. At $n=100$, the absolute bias reaches $1.020$ for $k=2$ and $0.242$ for $k=10$, accompanied by severely inflated standard errors (the median SE reaches $2.957$ for $k=5$ and $39.717$ for $k=10$). This occurs because the identity score estimator requires basis function expansions for intermediate estimation steps\footnote{See Section 4.2 in chen2023ann_npiv.}, which become numerically unstable when the effective sample size is small relative to the dimension. Despite the elevated bias and variance, DML maintains close-to-nominal coverage ($93$--$97\%$) as the inflated standard errors appropriately capture the increased uncertainty, consistent with the slight conservatism documented in chen2023ann_npiv.
The automatic debiasing procedure demonstrates substantially greater stability in small samples. At $n=100$, ADML achieves bias of $0.178$, $0.039$, and $0.001$ for $k=2$, $5$, and $10$, with moderate standard errors ($0.198$, $0.162$, $0.165$) that avoid the inflation seen in the analytical approach. ADML correctly quantifies its uncertainty, yielding coverage probabilities close to the nominal level. The procedure exhibits slight undercoverage ($90$--$93.5\%$ at $n=100$), contrasting with DML's conservatism and reflecting different bias-variance tradeoffs. As sample size increases, both procedures converge to similar performance, with ADML delivering comparable or lower bias throughout.
In this Section, we introduce a new framework for demand estimation that combines the nonparametric identification arguments of berry2014identification with the dimensionality reduction techniques of gandhi_houde2019, making it applicable to real data sets with more than two products unlike compiani2022market whose approach suffers from the curse of dimensionality.
We follow berry2014identification and present a general model of demand first, later on we will impose additional restrictions on the form of the indirect utility function. In market $t$, $t=1,\dots,\,T$, there is a continuum of consumers choosing from a set of products $\mathcal{J} = \{0,\,1,\dots,\,J\}$ which includes the outside option. The choice set in market $t$ is characterized by a set of product characteristics $\chi_{t}$ partitioned as follows:
where $x_{t} \equiv (x_{1t}, \dots, x_{Jt})$ is a vector of exogenous observable characteristics (e.g. exogenous product characteristics or market-level income), $p_{t} \equiv (p_{1t}, \dots, p_{Jt})$ are observable endogenous characteristics (typically, market prices) and $\xi_{t} \equiv (\xi_{1t}, \dots, \xi_{Jt})$ represent unobservables potentially correlated with $p_{t}$ (e.g. unobserved product quality). Let $\mathcal{X}$ denote the support of $\chi_{t}$. Then the structural demand system is given by
where $\Delta^{J}$ is a unit $J$-simplex. The function $\sigma$ gives, for every market $t$, the vector $s_{t}$ of shares for the $J$ goods.
Following berry2014identification, we partition the exogenous characteristics as $x_{t} = \left(x_{t}^{(1)}, x_{t}^{(2)} \right)$, where $x_{t}^{(1)} \equiv \left(x_{1t}^{(1)}, \dots, x_{Jt}^{(1)}\right)$, $x_{jt} \in \mathbb{R}$ for $j \in \mathcal{J} \backslash \{0\}$, and define the linear indices
and let $\delta_{t} \equiv (\delta_{1t}, \dots, \delta_{Jt})$. Without loss of generality, we can normalize $\beta_{j} = 1$ for all $j$ as unobserved characteristics $\xi_{jt}$ have no natural scale (see berry2014identification for more details). Given the definition of the demand system, for every market $t$,
Following berry2013connected and berry2014identification, under the connected substitutes assumption, there exists at most one vector $\delta_{t}$ such that $s_{t} = \sigma\left(\delta_{t}, p_{t}, x_{t}^{(2)}\right)$, meaning that for every inside good we can write
We can rewrite (ref) in a more convenient form to get the following estimation equation
Note that in (ref) the inverse demand is indexed by $j$, meaning that we have to estimate $J$ inverse demand functions, which is exactly why the approach of compiani2022market suffers from the curse of dimensionality.
To circumvent this problem, we exploit the symmetry properties of the demand system. Under the linear in characteristics utility specification, gandhi_houde2019 show that the inverse demand function admits a symmetric representation. Specifically, define the state vector for product $j$ in market $t$ as
where $\Delta_{jkt} = \tilde{x}_{jt} - \tilde{x}_{kt}$ denotes characteristic differences with $\tilde{x}_{jt} \equiv \left(p_{jt}, x_{jt}^{(2)}\right)$, and index $k$ runs over all products except $j$, including the outside good $k = 0$. Note that $\omega_{jt}$ contains (i) the shares of all rival products (from which the own share $s_{jt}$ can be recovered), and (ii) the differences in characteristics with respect to all alternatives, including the outside good, so that $\Delta_{j0t} = (p_{jt}, x_{jt}^{(2)})$ captures own-product levels.
Building on Proposition 1 of gandhi_houde2019, we establish the following result.
The key insight underlying Theorem (ref) is that the standard outside good normalization ($\delta_{0t} = 0$) pins down any market-specific constant that arises in the symmetric representation. Since the state vector $\omega_{jt}$ contains the difference $\Delta_{j0t} = (p_{jt}, x_{jt}^{(2)})$ with respect to the outside good, it encodes sufficient information to recover the outside good's state $\omega_{0t}$, thereby absorbing the market effect into $\gamma$.
Equation (ref) has several important features. First, the function $\gamma$ is not indexed by $j$, the same function applies to all products, with only the inputs $\omega_{jt}$ varying across products. This allows us to pool observations across products and markets for estimation. Second, $\gamma$ is symmetric in rival products: permuting the labels of competing goods $k \neq j$ leaves $\gamma(\omega_{jt})$ unchanged. This symmetry dramatically reduces the effective dimensionality of the estimation problem. Third, under the multinomial logit model, we have $\gamma(\omega_{jt}) = \beta_p p_{jt} + x_{jt}^{(2)\prime}\beta_x$, so that $\gamma$ can be seen as a generalization of the linear functional form in the logit model.\footnote{This estimation equation is inherently connected to the Independence of Irrelevant Alternatives (IIA) property. If we specify $\gamma(\omega_{jt}) = \beta_p p_{jt} + x_{jt}^{(2)\prime}\beta_x + g(\{\Delta_{jkt}\}_{j \neq k})$, we will get the IIA-test form Section 3.3 of gandhi_houde2019, who propose testing $g = 0$ as a diagnostic for IIA violations and instrument relevance.} Thus, (ref) is a semiparametric demand model in the spirit of lu2023semi, though derived from different primitives. While lu2023semi exploit the inclusive value structure to summarize cross-product information in a market-level aggregate, our approach retains product-level rival information in $\omega_{jt}$ and leverages symmetry for dimensionality reduction. Additionally, they treat the nonparametric component as having a random true value (induced by the dependence of the inclusive value on unobserved quality), whereas in our framework $\gamma$ is a fixed function evaluated at endogenous arguments, leading to a standard NPIV estimation problem.
Let $y_{jt} \equiv \log(s_{jt} / s_{0t}) - x_{jt}^{(1)}$, then we can rewrite equation (ref) in a more convenient form
Equation (ref) is the main structural equation where $\gamma$ is a complex nonparametric function characterizing the relationship between the inverse demand and product attributes and shares. Dimensionality of the input vector $\omega_{jt}$ depends on both the dimensionality of the characteristics space and the number of products in the market, thus, $\omega_{jt}$ is potentially high-dimensional. This will always be the case if we want to augment standard datasets with unstructured data such as product reviews, package images, etc. Since both the market shares $s_t$ and prices $p_t$ depend on the unobservable characteristics $\xi_t$, $\mathbb{E}[\xi_{jt}|\omega_{jt}] \neq 0$, and hence, $\omega_{jt}$ is endogenous.
In order to estimate $\gamma$, we need to construct a vector of instruments $z_{jt}$. BerryLevinsohnPakes1995 argue that the vector of product characteristics $x_{jt}$ is exogenous with respect to the structural error term $\xi_{jt}$, i.e. ${\mathbb E}[\xi_{jt}|x_{jt}] = 0$. This exogeneity condition can be used to construct demand side instruments $z_{jt}$. Instrument construction is a well-known problem in demand estimation, since it can lead to weak identification and distorted inference. We refer the reader to reynaert_verboven2014 and gandhi_houde2019 for a more detailed discussion.
To construct demand side instruments, we follow gandhi_houde2019 and use the transformed characteristics space $z_{jt} = (\{\Delta^{x}_{jkt}\}_{j \neq k})$, where $\Delta^{x}_{jkt} = x_{jt} - x_{kt}$, such that ${\mathbb E}[\xi_{jt}|z_{jt}] = 0$. Note that since $\omega_{jt}$ includes $z_{jt}$, it enforces strong correlation between endogenous inputs and instruments. If data permit, one can augment the instrument space with supply side instruments, such as cost shifters. Let $c_{jt}$ be a cost shifter for product $j$ in market $t$, then the instrument space becomes $z_{jt} = \left(\left\{\Delta^{x}_{jkt},\,\Delta^{c}_{jkt}\right\}_{j\neq k}\right)$, where $\Delta^{c}_{jkt} = c_{jt} - c_{kt}$.
One of the main primitives in demand estimation is substitution patterns, which allow the researcher to investigate the responsiveness of consumer choices to changes in market structure and, thus, understand the nature of competition between firms. A conventional measure of substitution is price elasticities, reliable estimation and inference on which is therefore a first-order concern for applied researchers. Traditional approaches rely on parametric demand models such as logit, nested logit, BLP BerryLevinsohnPakes1995 that impose strong functional form restrictions on the mapping from structural parameters to elasticities. The semiparametric framework of Section (ref) removes these restrictions, but introduces a new challenge: the price elasticity is a nonlinear functional of $\gamma$, and therefore requires the machinery developed in Section (ref). For expositional clarity, we focus on own-price elasticities below; the analysis extends to cross-price elasticities naturally.
The own-price elasticity for product $j$ in market $t$ is defined as the percentage change in market share in response to a percentage change in own price,
By the implicit function theorem applied to equation (ref), the share derivative matrix satisfies
where $L \in {\mathbb R}^{J \times J}$ is the log-share Jacobian with $L_{jk} = \partial\log(s_{jt}/s_{0t})/\partial s_{kt}$, $\Gamma^p \in {\mathbb R}^{J \times J}$ collects the price derivatives $\Gamma^p_{jk} = \partial\gamma(\omega_{jt})/\partial p_{kt}$, and $\Gamma^s \in {\mathbb R}^{J \times J}$ collects the share derivatives $\Gamma^s_{jk} = \partial\gamma(\omega_{jt})/\partial s_{kt}$. Setting $A(\gamma) \equiv L - \Gamma^s(\gamma)$, the own-price elasticity becomes
The matrix $A$ captures the net effect of shares on the inverse demand system: $L$ accounts for the mechanical relationship between log-share ratios and share levels, while $\Gamma^s$ reflects how $\gamma$ responds to share movements. The interplay between $A^{-1}$ and $\Gamma^p$ encodes the full substitution pattern implied by the model, including indirect effects through equilibrium share responses. See Appendix (ref) for a detailed derivation.
Under the logit model, $\gamma(\omega_{jt}) = \beta_p p_{jt} + x_{jt}^{(2)\prime}\beta_x$, so $\Gamma^s = 0$ and $\Gamma^p = \beta_p I$. In this case $A = L$ and the elasticity simplifies to $\varepsilon_{jj} = \beta_p\,p_{jt}(1 - s_{jt})$, recovering the familiar closed-form logit formula. More generally, departures from the logit model introduce nonzero share derivatives $\Gamma^s \neq 0$, which generate richer and more flexible substitution patterns through the matrix inverse.
The functional of interest is the average own-price elasticity,
where the expectation is taken over the distribution of markets. Since $\varepsilon_{jj}(\gamma)$ depends nonlinearly on $\gamma$ through the matrix inverse $A^{-1}(\gamma)$, the average own-price elasticity (ref) is a nonlinear functional of $\gamma$. To construct a debiased estimator with valid inference, we appeal to the framework of Section (ref).
The key step is to verify that the Gateaux derivative of $\varepsilon_{jj}(\gamma)$ with respect to $\gamma$ is linear in the perturbation direction, so that the Riesz representer is well-defined and the PGMM estimator of Section (ref) can be applied. Proposition (ref) in Appendix (ref) establishes that for a perturbation $\zeta$\footnote{For notational convenience, we write $D_\gamma \varepsilon_{jj}[\zeta]$ for the Gateaux derivative $D(W, \gamma, \zeta)$ specialized to the elasticity functional.},
where $Z^p_{jk} = \partial\zeta(\omega_{jt})/\partial p_{kt}$ and $Z^s_{jk} = \partial\zeta(\omega_{jt})/\partial s_{kt}$ are the price and share derivative matrices of the perturbation $\zeta$. Expression (ref) consists of two terms: the first captures how a perturbation of $\gamma$ affects the price derivative matrix $\Gamma^p$, while the second captures how it affects the share derivative matrix $\Gamma^s$ and, consequently, the matrix inverse $A^{-1}$. Proposition (ref) shows that $D_\gamma\varepsilon_{jj}[\zeta]$ is linear in $\zeta$ since the matrices $A$ and $\Gamma^p$ depend on the current function $\gamma$ at which the derivative is evaluated, not on the perturbation direction $\zeta$, while $Z^p$ and $Z^s$ are linear in $\zeta$ by the linearity of differentiation. Therefore, Theorem (ref) and Corollary (ref) are applicable, ensuring that the debiased estimator is asymptotically normal under the conditions stated therein.
A critical distinction from linear functionals is the role of double cross-fitting. For a linear functional, the target vector $\hat{M}_\ell$ in the PGMM problem does not depend on an estimate of $\gamma$ and can be computed directly from the basis functions. For the elasticity functional, however, the Gateaux derivative (ref) involves $\Gamma^p(\gamma)$, $\Gamma^s(\gamma)$, and $A^{-1}(\gamma)$, all of which depend on $\gamma$. Consequently, $\hat{M}_\ell$ must be constructed using an estimate $\hat\gamma$, and additional sample splitting is required to prevent the estimation error in $\hat\gamma$ from contaminating the Riesz representer estimate (see Appendix (ref) for details).
The asymptotic theory of Sections (ref)--(ref) applies with the market as the unit of observation which is a natural asymptotic regime for scanner data applications nevo2001measuring. Markets are assumed to be independent, so cross-fitting partitions markets rather than individual product-level observations. Let $W_{jt} \equiv (y_{jt},\,\omega_{jt},\,z_{jt})$ be a data tuple, and let ${\mathcal T}_\ell$, $\ell = 1, \ldots, L$, be a partition of the observation index set $\{1, \ldots, T\}$ into $L$ distinct subsets of approximately equal size. Then for test fold $\ell$, the target vector is constructed as
where $\hat\gamma_{\ell,\ell'}$ is trained on all markets outside both folds $\ell$ and $\ell'$, and $d_j$ denotes the $j$-th basis function. This double cross-fitting ensures that each summand in (ref) evaluates the Gateaux derivative at markets from fold $\ell'$ using a $\hat\gamma_{\ell,\ell'}$ that was not fitted on those markets. At the same time, the entire vector $\hat{M}_\ell$ uses no data from the test fold $\ell$, preserving the independence required for valid inference on the Riesz representer $\hat\alpha_{jj,\ell}$. As discussed in Section (ref), the presence of $\hat\gamma_{\ell,\ell'}$ in $\hat{M}_\ell$ slows its convergence rate from $\sqrt{\log q / T}$ to $\kappa_T^\gamma$, the convergence rate of the MLIV estimator in the standard mean square norm, which in turn requires a slightly stronger rate condition on $\hat\gamma$ (Assumption (ref)).
Given the estimates $\hat\gamma_{\ell}$ and $\hat\alpha_{jj,\ell}$ for each fold $\ell = 1,\dots,\,L$, the debiased estimator of the average own-price elasticity takes the form
where $\varepsilon_{jj}\!\left(W_t,\,\hat\gamma_\ell\right)$ is the plug-in own-price elasticity for product $j$ in market $t$ and $\hat\alpha_{jj,\ell}(z_{jt})$ is the PGMM estimate of the Riesz representer of the Gateaux derivative of $\varepsilon_{jj}$. The variance can be consistently estimated by
The full estimation procedure is described in Appendix (ref).
We now evaluate the finite sample performance of the debiased estimator for the average own-price elasticity functional $\theta_{jj} = {\mathbb E}[\varepsilon_{jj}(\gamma)]$. This is a nonlinear functional of $\gamma$, requiring the double cross-fitting procedure and the stronger convergence conditions described in Section (ref). Unlike the average derivative studied in Section (ref), the elasticity depends on $\gamma$ through the matrix inverse $(L - \Gamma^s)^{-1}$, making both the functional evaluation and the Riesz representer estimation more challenging.
We simulate data from a simple logit model. The mean valuation in the logit model is given by $\delta_{jt} = \beta_{p} p_{jt} + x_{jt}'\beta_{x} + \xi_{jt}$. Product shares can be calculated using the following formulae, for $j=1,\dots,\,J$ and $t=1,\dots,\,T$,
We set the total number of product characteristics besides the price to be equal to 4, i.e. $x^{(1)}_{jt}$ is a scalar and $x^{(2)}_{jt}$ is a three-dimensional vector. We draw all the observed product characteristics from ${\mathcal U}(0,1)$, while the unobserved characteristics, $\xi_{jt}$, are distributed as ${\mathcal N}(1,\,0.15^2)$ for all $j$ and $t$. The price is
where $c_{jt} \sim {\mathcal U}(0,1)$ is a cost shifter, and $e_{jt} \sim {\mathcal U}(0, 0.1)$ is idiosyncratic noise. The price coefficient is $\beta_{p} = -2$ and the coefficients on product characteristics are $\beta_{x} = (1,\,-0.5,\,0.5,\,1)'$.
To estimate $\gamma$, we choose KIV over the Double Lasso estimator used in Section (ref) due to its substantially lower computational cost in the nonlinear functional setting which requires double cross-fitting. We construct dictionaries $d(\omega_{jt})$ and $b(z_{jt})$ using quadratic polynomials with and without interactions, respectively. Under the specified DGP, $\omega_{jt} = \left(\{s_{kt}, \Delta_{jkt}\}_{j \neq k}\right)$ and $z_{jt} = \left(\left\{\Delta^{x}_{jkt},\,\Delta^{c}_{jkt}\right\}_{j\neq k}\right)$, and hence, $dim(\omega_{jt}) = dim(z_{jt})$ and $p>q$. We run 500 replications for $J \in \{2, 5\}$\footnote{Note that the estimation problem quickly becomes high-dimensional. Under $J=2$, $dim(\omega_{jt}) = 15$, which grows up to $30$ when $J=5$.} products and $T \in \{100, 200, 400, 800\}$ markets. We use five-fold cross-fitting, $L = 5$. Without loss of generality, we focus on the elasticity of the first product ($j = 1$).
Under logit, the true own-price elasticity for product $j$ in market $t$ has the closed-form expression $\varepsilon_{jj,t} = \beta_p\, p_{jt}(1 - s_{jt})$. The population parameter of interest is $\theta_0 = {\mathbb E}[\varepsilon_{jj,t}]$, which we approximate using a large pre-simulation draw of $T = 100{,}000$ markets. The resulting values are approximately $-4.22$ for $J = 2$ and $-4.28$ for $J = 5$.
Table (ref) reports absolute bias, median estimated standard error, and empirical coverage of nominal 95% confidence intervals for the plug-in and debiased estimators. Unlike the average derivative application in Section (ref), the analytical form of the Riesz representer for the elasticity functional is intractable, hence, we omit the analytical debiasing approach from the comparison.
The plug-in estimator severely undercovers across all specifications. While its absolute bias is small, its estimated standard errors dramatically understate the true sampling variability. In contrast, the debiased estimator achieves near-nominal coverage throughout, ranging from 91.2% to 95.0%. These results align with and extend the evidence from Section (ref) on linear functionals. The key qualitative finding that plug-in inference fails due to standard error underestimation rather than excessive bias carries over from average derivatives to own-price elasticities.
Nevertheless, the nonlinear setting introduces additional challenges causing the coverage for the debiased estimator to be slightly lower than in the linear case, particularly at small $T$ with $J = 2$. This behavior is consistent with the additional estimation noise introduced by double cross-fitting and the PGMM estimation of the $\gamma$-dependent Riesz representer. Despite these additional difficulties, the automatic debiasing procedure reliably delivers valid inference for the nonlinear elasticity functional across all sample sizes considered.
We apply our estimator to retail scanner data from the IRI Academic Database bronnenberg2008IRI. It records unit sales by UPC code, store, and week for a sample of supermarkets over 2001--2012, together with product characteristics. We focus on one year span of 2003 and top ten most sold products which account for approximately $69\%$ of all sales. To exploit variation in product attributes, we do not aggregate to the brand level; instead, products are defined by the combination of a brand name and a set of observable characteristics, e.g. Coke Classic or Caffeine-free Diet Pepsi.
Carbonated beverages are sold in different packages and package sizes. We restrict our attention to cans and define a product unit as a 12 oz can. Hence, we construct market shares based on the total amount of cans sold. Prices are defined as the ratio of total revenue to total number of units sold. We aggregate the data to geographic region-month level resulting in an unbalanced panel with 5,993 observations at the product-region-month level.
Data on product characteristics includes beverage flavor, sugar, caffeine and calorie levels. All characteristics are represented by categorical variables, thus, for computational reasons we aggregate product attributes in larger groups (see Appendix (ref) for more details). After aggregation and dropping collinear characteristics, we are left with five product attributes we use for estimation. We have Caffeine and Sugar dummy variables indicating whether a product contains caffeine and sugar, respectively. The remaining three variables represent different flavor categories: Cola, Lemonade, and Pepper, with the base category representing other flavors.
We begin with a logit model as a parametric benchmark. In the absence of external cost-shifter data, we use BLP-style instruments (sums of rivals' attributes) to identify the price coefficient. Given the binary nature of product characteristics, most rival sums are nearly collinear with own attributes, leaving only two instruments with independent variation: sums of rivals' sugar and caffeine levels. One popular approach to achieve stronger price coefficient identification without external data is to use Hausman instruments hausman1996, however, in our case they capture endogenous demand shifts rather than track down exogenous cost shocks\footnote{See Appendix (ref) for more details.}, which is a common critique of this approach nevo2001measuring. We therefore use BLP instruments throughout.
Before estimating the semiparametric model, we present evidence that the data reject the logit specification. Following gandhi_houde2019, we conduct an Independence of Irrelevant Alternatives (IIA) test using local differentiation instruments.\footnote{See Appendix (ref) for details on the test construction.} The test strongly rejects IIA ($F = 68.5$, $p < 0.0001$), indicating that the data contain richer substitution patterns than logit can rationalize. This motivates the use of a semiparametric approach that does not impose parametric restrictions on demand.
Similarly to Section (ref), we estimate $\gamma$ using KIV and construct debiased estimates using five-fold cross-fitting. Since the dimensionality of $\omega_{jt}$ and $z_{jt}$ grows with the number of products, we exploit the symmetry of the inverse demand function to construct symmetric basis functions\footnote{Using symmetric dictionaries implicitly restricts the Riesz representer $\alpha_0(z_{jt})$ to be a symmetric function of rival characteristics. This restriction is justified by the structure of the moment conditions that define $\alpha_0$. First, the own-price elasticity $\varepsilon_{jj}$ is invariant to permutations of rivals, so the Gateaux derivative $D_\gamma\varepsilon_{jj}[d_k]$ inherits this symmetry whenever $\gamma$ and the basis function $d_k$ are symmetric (Proposition (ref)). Then, since the instruments $z_{jt}$ are also constructed symmetrically, the defining moment condition ${\mathbb E}[D_\gamma\varepsilon_{jj}[d_k]] = {\mathbb E}[\alpha_0(z_{jt})\,d_k(\omega_{jt})]$ is symmetric, and hence $\alpha_0$ must be as well.} for dictionaries $d(\omega_{jt})$ and $b(z_{jt})$ (see Appendix (ref) for details). It allows us to reduce the dimensionality of the PGMM problem leading to substantial computational gains.
Table (ref) reports own-price elasticity estimates for each of the ten products under three approaches: the logit benchmark, the plug-in KIV estimator (PI), and the debiased estimator (ADML).
First, the semiparametric estimators yield substantially larger elasticities in absolute value than the logit. The mean own-price elasticity increases from $-3.09$ under the logit to $-3.75$ (plug-in) and $-3.73$ (ADML). This is consistent with the IIA rejection: by constraining all cross-price elasticities to be proportional to market shares, the logit restricts substitution patterns in a way that attenuates own-price elasticities. While the magnitudes differ, the qualitative ranking of products is stable across methods. More popular products like Coke Classic and Pepsi Classic have smaller elasticities in absolute value, while niche products such as Caffeine-Free Diet Pepsi and Mountain Dew have larger elasticities. This implies that consumers of mainstream products are less sensitive to their price changes compared to more niche products, which is quite intuitive.
Second, the debiasing corrections vary substantially across products, both in magnitude and direction. For several products (Dr.\ Pepper, Mountain Dew Classic, and Mountain Dew Other) the two estimates are virtually identical, suggesting that the plug-in bias is negligible for these products. For all Coke products and Sprite, the debiased estimates are less elastic (positive debiasing term) compared to the plug-in ones, while all Pepsi products exhibit more elastic own-price elasticities compared to the plug-in ones (negative debiasing term).
Third, both plug-in and debiased estimates have similar standard errors, with the exception of Diet Coke that also has the largest correction. We observe that debiasing plays a big role in this application resulting in half the products (Diet Coke, Caffeine-free Diet Coke, Sprite, Diet Pepsi, and Caffeine-free Diet Pepsi) having debiased estimates differ from plug-in estimates by much more than the associated standard errors. These results are also consistent with the CNS demand application.
In this paper, we have given an automatic method of debiasing functionals of machine learners under endogeneity. We have shown how to use a PGMM minimum distance estimator to perform debiasing using only the form of the object of interest, without knowing the form of the bias correction term. We allow for a wide range of MLIV estimators that satisfy certain convergence rate conditions. We have shown root-$n$ consistency and asymptotic normality and given a consistent asymptotic variance estimator for both linear and nonlinear functionals. For linear functionals we require MLIV estimators to converge fast enough in the projected mean square norm, while for nonlinear functionals we require fast enough convergence in the standard mean square norm, which is a more stringent requirement due to ill-posedness. Relaxing the convergence rate condition for nonlinear functionals as well as extending the approach to irregular functionals are promising directions for future research.
Finally, we have applied our debiasing procedure to estimate the mean own-price elasticity functional in the nonparametric demand for differentiated products framework. We have obtained evidence from both simulated and real scanner data that plug-in estimates are biased and have smaller variance compared to the debiased estimates, which reflects the bias-variance trade-off occurring due to regularization. Using scanner data, we have demonstrated that debiasing corrections vary substantially across products, both in magnitude and direction. In particular, half the products exhibit profound debiasing effects highlighting the importance of debiasing in our application.