Automatic Debiased Machine Learning of Structural Parameters with General Conditional Moments

Automatic Debiased Machine Learning of Structural Parameters with General Conditional Moments

\newgeometry{height=9in,left=1in,right=0.75in,bottom=1in}

Introduction

Models defined by conditional moment restrictions (CMRs) are ubiquitous in economics, appearing in a wide range of settings such as regressions, quantile models, dynamic discrete choice models, simultaneous equations models, demand estimation, auctions, missing data, and production functions, to name a few chen2016methods. These settings are often considered semiparametric as they involve parameters belonging to an infinite dimensional space, for example, the policy function that characterizes a firm's optimal investment decision given its unobserved productivity and inputs. In addition to these infinite-dimensional components, these models also include a finite-dimensional vector of structural parameters, such as the marginal products of capital and labor, which we are interested in. While our focus is on estimating these structural parameters, it is necessary to estimate the nuisance parameters as well, since they are unknown and play a role in the estimation of the parameters of interest.

Machine Learning tools, e.g., Boosting, Random Forest, Lasso, and Neural Networks, have proven to be suitable algorithms for handling such nonparametric estimations due to their strong prediction capabilities. Consequently, there is great potential to leverage these tools in a wide range of economic models. However, the use of machine learning typically introduces statistical challenges. The strength of these algorithms lies in their specific regularization strategies, which aim to achieve an optimal bias-variance trade-off---crucial for reliable estimation of unknown functions. This, nonetheless, typically leads to biased estimations of the parameters of interest and the impossibility of attaining $\sqrt{n}-$consistency. A practical implication of this issue is poor coverage of the commonly used confidence intervals and tests with sizes exceeding nominal levels.

This paper introduces an algorithm that restores $\sqrt{n}-$consistency when machine learning is used to estimate nuisance parameters in a broad class of models defined by semiparametric CMRs. Thus, it enables the application of modern estimation techniques in various structural economic settings. Specifically, this work presents a general estimation method for Locally Robust (LR)/(Neyman-)Orthogonal/ Debiased Moments chernozhukov2022locally, neyman1959 in these contexts.\footnote{Throughout the paper, we will use these terms interchangeably.} These moments are particularly appealing as estimation based on them is less affected by the regularization involved in learning high-dimensional parameters.

Our work contributes to the existing literature on debiased machine learning by developing a general algorithm for constructing debiased moments in a data-driven way within models characterized by a finite number of CMRs, with possibly different conditioning variables and endogenous regressors (in the sense that the nuisance parameters may be a function depending on variables other than those used for conditioning). To be specific, this paper develops an estimation technique for obtaining Orthogonal Instrumental Variables (OR-IVs) which are the key components to obtain LR moments in models with CMRs, as pointed out by chernozhukovetaldebiased2016, chernozhukov2018double, and, more recently, by arganaraz2025machine. OR-IVs can be interpreted as functions of the conditioning variables, with finite second moments, where the effect of estimating the non-parametric component has been partialed out.

chernozhukovetaldebiased2016, Section 7, first proposed a definition of OR-IVs in models defined by a finite number of semiparametric CMRs, with different conditioning information sets, when the researcher is interested in a finite-dimensional vector of parameters appearing in the residual functions. They focused on a specific OR-IV that minimizes the asymptotic variance of GMM estimators in the class of GMM estimators based on orthogonal instruments. The authors did not discuss how to estimate such optimal OR-IV or other OR-IVs in general. chernozhukov2018double considered the special case with only one conditioning information set under exogeneity (i.e., when the nonparametric component is a function of a sub-vector of the conditioning variable). arganaraz2025machine characterized OR-IVs in the same model of chernozhukovetaldebiased2016, but for the case where the parameter of interest is a general smooth functional of finite-dimensional and nonparametric components, which might not necessarily appear in the residual functions of the model. chernozhukov2018double and arganaraz2025machine worked out estimation of particular OR-IVs in specific applications such as the partially linear model with and without an endogenous treatment. Differently to these works, the goal of the current paper is to introduce a strategy to estimate any OR-IV, in a computationally convenient way, in a general model with CMRs. For simplicity, we construct LR moments for a finite-dimensional vector of parameters appearing in the residual functions; nevertheless, the general principles developed in this paper can be extended to more general functionals.

Building on this strategy, we introduce a Debiased GMM estimator (DGMM) that is locally insensitive to regularization techniques. In our context, the parameter of interest must have a finite efficiency bound, but it can map into the CMRs non-linearly, which is a common situation in structural models in economics.\footnote{Recall that having a finite efficiency bound is equivalent to the (efficient) Fisher information matrix being positive definite, which is a necessary condition for $\sqrt{n}-$ estimability of the parameters of interests in general; see van1991differentiable.} For instance, the object of interest could be interpreted as the parameters of a production function on which the CMRs depend in a non-linear way, as in the semiparametric production function setting considered by ackerberg2014asymptotic. Furthermore, we obtain a rate of consistency of our OR-IV estimators, which is essential to show $\sqrt{n}-$consistency and asymptotic normality of DGMM, as we do in this paper.

Our approach to obtaining debiased moments leverages a characterization of such moments in our general setting for general smooth functionals, previously derived by arganaraz2025machine. They prove that an orthogonal moment is necessarily a sum of the products between the residual functions and the OR-IVs. Since the residual functions are primitives of the economic model, known up to the unknown parameters, the key step in computing an LR moment is the calculation of OR-IVs. This will allow us to construct a moment that is locally insensitive to estimation of the infinite-dimensional parameters, including the OR-IVs, and then proceed with a suitable estimation of the parameters of interest. Unfortunately, OR-IVs are necessarily solutions to functional equations, which might potentially be challenging to solve analytically.

The key idea is based on constructing a linear operator derived from the CMRs.\footnote{See luenberger1997optimization and CARRASCO20075633 for key results on linear operators.} The range of this operator collects all the possible Fréchet derivatives of the initial CMRs with respect to the high-dimensional parameters. OR-IVs must be orthogonal to the range of this operator, under a suitable notion of orthogonality. Our work builds on the observation that this orthogonalization can be accomplished by “residualizing" a known function of the conditional variables, i.e., an IV, from the orthogonal projection onto the range of the operator. This exercise requires obtaining this orthogonal projection, which is challenging as the operator itself is unknown, and thus its orthogonal projection is unknown.

We show that computing such an orthogonal projection can be seen as minimizing a mean squared error, resulting from projecting a vector of known starting IVs onto the linear span of suitable regressors. The part of the starting functions that is not explained by the regressors, i.e., the residual, is shown to be a valid OR-IV. We allow the dimension of the regressors involved in that minimization to be greater than the sample size (and depend on it) by adding an $\ell_1$- norm penalization term. The resulting objective function is then familiar from standard Lasso problems. Nonetheless, this Lasso program is special in that it is based on unknown regressors. This is a by-product of not knowing the linear operator that we mentioned in the previous paragraph. As such, these regressors must be estimated prior to minimizing the corresponding Lasso objective function.

Conveniently, we show that those regressors take the form of conditional expectations (also known as reduced forms), given the conditioning variables of the model. Consequently, they can be estimated by suitable machine learning tools, useful for predicting regression functions. Once these regressors are estimated, the problem reduces to a Lasso problem with generated regressors, which can be implemented straightforwardly by well-known algorithms. We provide theoretical guarantees for our estimation procedure, assuming that some sparsity condition holds, which implies that the orthogonal projection can be well-approximated by a “small" number of terms.

This paper argues that for inference in high-dimensional settings, researchers should use OR-IVs only, as moments based on them are debiased. Hence, our algorithm justifies the use of OR-IVs over other choices of instruments that do not necessarily yield a debiased moment and that are currently employed in applied work. Our approach converts any suitable function of the conditioning variables (e.g., $Z$ or $Z^2$) into an OR-IV that can be used in estimation as all these functions can be “residualized" using the approach that we described above. If multiple such functions are provided, multiple orthogonal moments can be estimated and combined as usual in GMM. A natural approach is to start with the standard choices in applied work (e.g., $Z$ or $Z^2$) and apply our method to transform them into OR-IVs.

We propose estimating the nuisance parameters with machine learning and cross-fitting chernozhukov2018double, and then estimate the parameters of interest with GMM . We develop this estimator in a two-step framework, where the nuisance parameters are estimated in a first step, and the estimation of the finite-dimensional parameter is conducted in a second step, using the first-step estimations. The rationale for employing GMM stems from our theoretical findings, which assure that standard inference remains valid for the finite-dimensional parameter, despite using machine learning tools to estimate the nuisance parameters in a first stage. Specifically, leveraging the asymptotic properties established for the OR-IV estimators, we demonstrate that the estimator we introduce is $\sqrt{n}-$consistent and asymptotically normal. For this, we only require that the nuisance parameters being estimated at a rate faster than $n^{-1/4}$, which can be achieved by a variety of machine learners. When orthogonal moments are not used, standard errors must be correctly computed to account for first-stage estimation newey1994asymptotic. Nonetheless, we show that standard errors based on our debiased moments directly account for such a first-stage estimation by using the usual “sandwich" formula, simplifying their calculation and avoiding using bootstrap methods.

Additionally, we apply our general theory to a semiparametric version of production functions estimation using the proxy variable approach, studied by olleypakes96, levinsohn2003estimating, and ackerberg2015identification, which might be of independent interest.\footnote{Production functions have been at the core of economics since the early 1800's; see the introduction in chambers1988applied for an historical review on production function and productivity.} Production functions and measures of productivity play a central role in several empirical settings in economics. For instance, they have been used to study the effects of trade liberalization, exporting, foreign ownership, competition, importing intermediate goods, investment climate, learning by doing; see ackerberg2007econometric, ackerberg2015identification, and references therein. We apply debiased machine learning in this context using off-the-shelf machine learning tools to conduct inference on the parameters of all the inputs simultaneously. cha2023inference explore the use of the Double Debiased Machine Learning (DML) estimator of chernozhukov2018double in the context of production functions. Since the authors based estimation on the orthogonal score function from robinson1988root, they conduct inference on only one parameter at a time and treat the parameters of the rest of inputs, including capital, as nuisance ones. As ackerberg2015identification pointed out, conducted inference using the Robinson's transformation might not be ideal as it could lead to identification failures from functional form dependence. Moreover, cha2023inference focus on a particular strategy, introduced by the authors, to estimate the nonparametric components. Independently, doraszelski2025productionfunctionestimationinvertibility have derived a Neyman-orthogonal moment for the proxy variable approach by explicitly obtaining the first-step influence function (accounting for regularization bias in the first stage) and imposing restrictions among the conditioning variables, which our more general approach does not require.

We have assessed the finite sample performance of the introduced estimator through several Monte Carlo experiments, in the context of production functions. Despite this setting is not “of high dimension" (as there are only a few inputs), it entails an unknown nonparametric object which is potentially highly nonlinear. It arises from an unknown policy function---a solution to a potentially complex firm's dynamic problem---that links productivity (observed by the firm but unobserved to the econometrician) and a proxy variable, e.g., investment or intermediate inputs, given the other inputs in the production function. We aim to estimate the marginal products of these inputs in the presence of the unknown link function. Estimation is performed using machine learning, and we employ our OR-IV construction to correct for potential regularization bias so that standard inference remains valid. Our results show that DGMM exhibit smaller bias and coverage closer to the nominal level compared to a naive estimator based on the same first stage but not addressing regularization bias, consistent with our theoretical findings.

In our empirical application, we used data from a panel of Chilean firms, which has been extensively studied in the production function literature. We aimed to assess whether using machine learning could lead to different conclusions compared to other estimators that model first stages less flexibly and are commonly used in applied work. The results suggest that the choice of approach for modeling the first stage should be made with caution, as it can significantly influence empirically relevant conclusions, particularly regarding production function parameters. Consequently, inference methods like the one introduced in this paper are potentially appealing because they allow researchers to adopt less parametric specifications in the first stage. This is important for two reasons. First, it provides an additional empirical tool for scholars who wish to avoid imposing restrictive parametric assumptions. Second, it helps those more confident in traditional approaches to evaluate the sensitivity of their results to less flexible specifications.

The rest of the paper is organized as follows. Section (ref) revises the most related works to ours, emphasizing what our contribution is. Section (ref) works out a production functions example to motivate our results and describe the problem we face. Section (ref) defines LR moments in our context and discusses our approach heuristically. Section (ref) is the core of the paper as it introduces our approach to estimate OR-IVs in a general setting. Section (ref) proposes an algorithm that can be followed in practice. Section (ref) obtains a convergence rate for the estimators of the OR-IVs. In Section (ref), we present the estimator of the finite-dimensional parameter in a two-step setting, and Section (ref) provides conditions under which it is $\sqrt{n}-$consistent and asymptotically normal distributed. Section (ref) studies the performance of our estimator through different Monte Carlo experiments. Section (ref) contains our empirical results. Finally, Section (ref) provides final remarks. A Supplementary Appendix collects the proofs of all the theoretical results, states additional formal findings, applies our general asymptotic theory to our example, elaborates on implementation details of the Lasso-type program we introduce and its justification, and provides additional details on the Monte Carlo experiments.

Related Literature

Our work relates to various strands of the literature. First, it contributes to the body of research on CMRs. Our general model encompasses several settings that have been extensively studied in this literature. A seminal contribution in this regard is made by ai2012semiparametric, who establish the semiparametric bound (and propose an efficient estimator) for structural parameters in semiparametric models defined by nested CMRs, and thus involving different conditioning variables. Their work extends previous results by chamberlain1987asymptotic for the purely parametric case with one conditioning variable, as well as by chamberlain1992efficiency and ai2003efficient, who study semiparametric models with a single conditioning variable. It also relates to the findings in chamberlain1992cond and brown1998efficient, who consider parametric cases with nested CMRs; see also ai2007estimation. Furthermore, our model closely relates to the nonparametric setting with general CMRs studied by chen2018overidentification.

As discussed by chen2016methods, nuisance parameters in general CMRs settings have typically been estimated using sieve-based methods chen2007large, with finite-dimensional linear sieve approximations being a popular choice. Estimation of the parameters of interest can be performed simultaneously or in a two-step fashion. While the sieve approach is computationally attractive for both estimation and inference (for example, they already lead to debiased estimators when first steps are efficiently estimated) it can be restrictive. For instance, sieve methods often assume knowledge of the important approximating terms that should be included in estimation, an assumption that might be strong in complex, high-dimensional settings. In contrast, Lasso over a dictionary does not require such an assumption, as it can learn the relevant terms from the data, as discussed by bradic2022minimax.\footnote{A prominent exception is neural networks, which can be seen as a non-linear sieve approach. When neural networks are used to model the first-stage, LR moments and cross-fitting might not be needed, and inference can be alternatively conducted following, e.g., ai2007estimation, ai2012semiparametric, chen2023efficient, and chen2024inference.}

We build upon this literature by allowing for the utilization of a wide array of machine learning tools to estimate nuisance parameters in a first stage, and by automatically constructing LR moments. This provides researchers with additional tools to undertake their analyses using general semiparametric CMRs, invoking more general conditions on the space of functions to be approximated. In this paper, we follow a two-step approach for estimating the parameters of interest. A two-step framework is often more computationally attractive than a one-step estimation ackerberg2014asymptotic.

Second, this paper is connected to the recently developed literature on debiased moments. The general construction and asymptotic theory of these moments have been established by \citet*{chernozhukov2022locally} (CEINR) for high-dimensional settings where machine learners are employed in a first stage. Typically, a debiased moment function consists of two components: the initial identifying moment function and an adjustment term, known as the first-step influence function. This function characterizes the local effect of the first step on the starting (non-debiased) moment and its existence is generally equivalent to the starting moment having a finite semiparametric variance bound. When the functional form of this adjustment term is known, CEINR provide a general method to estimate the components in such functional form.\footnote{Typically, the first-step influence function depends on the data, the finite-dimensional and nonparametric components, and an extra nuisance parameter. CEINR propose a general method to estimate this nuisance object.} In this paper, we specifically examine semiparametric models defined by CMRs that involve varying conditioning variables and we do not need to obtain the functional form of the adjustment term, as it will be entirely estimated from the data, as we explain below. Focusing on CMRs with different conditioning variables allows us to extend the automatic construction of debiased models to new interesting applications such as learning firm-level production functions and productivity measures.\footnote{As CMRs imply unconditional moments, an alternative approach to construct LR moments in our setting, is to first convert conditional restrictions into unconditional ones, and then apply CEINR's approach. This, however, will entail deriving the first-step influence function, which might not be trivial in several structural models in economics. Our method, in contrast, avoids this.}

As we stated in the Introduction, OR-IVs have been introduced by the first version of CEINR, chernozhukovetaldebiased2016, within models with semiparametric CMRs and varying conditioning variables, which leads to debiased moments for a finite-dimensional vector of parameters appearing in the residual functions. chernozhukov2018double consider the special case where there is a common conditioning information set and the nonparametric component is a function of a subset of the IVs of the model. More generally, arganaraz2025machine characterize OR-IVs for general smooth functionals of finite-dimensional and nonparametric unknowns, where the functionals might not index the CMRs. chernozhukov2018double and arganaraz2025machine consider estimation of OR-IVs in specific contexts such as the partially linear model with and without an exogenous treatment, respectively. We complement this line of research by developing an algorithm to estimate any OR-IVs in general contexts with CMRs and providing its asymptotic guarantees. This procedure allows for valid inference on the finite-dimensional components in the residual functions, where there are varying conditioning variables, and potential endogeneity. While we focus on these parameters, our theory can be extended to more general functions, similarly to the algorithms developed by arganaraz2025existence for models with unobserved heterogeneity. However, this is beyond the scope of the current paper and will be explore in future research.

Third, our paper contributes to the literature on the automatic construction of debiased moments, which have typically developed specific algorithms based on that in CEINR. In many applications, the first-step influence function can be expressed as the product of a residual function and the Riesz representer of the initial moment function that identifies the parameter of interest ichimura2022influence. chernozhukov2022automatic propose a Lasso-type minimum distance estimator for such Riesz representer, allowing the nuisance parameter to be estimated using machine learning tools; see also chernozhukov2020adversarial, chernozhukov2021automatic, chernozhukov2022debiaseda, and chernozhukov2022debiasedb for related work. While these papers focus on exogenous cases, where the nuisance parameter depends on conditioning variables, bakhitov2022automatic extends this framework to endogenous settings, proposing a penalized GMM approach for estimating the Riesz representer. Additionally, farrell2021deep offer an automatic method for constructing adjustment terms that produces debiased moments, even when the first-stage parameter is not simply a prediction but any well-defined parameter in an economic model. In these references, the parameter of interest is the average of some smooth function of the first-stage parameter.\footnote{While farrell2021deep focus on parameters that are averages, their first-step correction can potentially be extended to situations where the object of interest solves a general moment condition, see their Remark 3.1.}

Our paper complements these works by working with CMRs and not focusing on situations where the parameter of interest represents an expectation, but on parameters of economic models “with more structure". Specifically, we address the more general case where the model implies conditional moments of arbitrary (but smooth) functions that depend on both finite and infinite-dimensional parameters. Thus, our paper paves the way for applying the automatic construction of debiased moments developed in these previous works to nonlinear and highly complex models with CMRs. This might be especially relevant in semiparametric models prevalent in structural economics such as firm-level production functions olleypakes96, ackerberg2014asymptotic, levinsohn2003estimating, general missing data problems hristache2016semiparametric, hristache2017conditional, graham2011efficiency, auctions kawai2011auction, and simultaneous equations models newey1999nonparametric. Further examples can be found in Section 5 of chen2016methods. In the production function setting, for instance, the parameter of interest are marginal productivities which cannot be written as simple means.

Lastly, our paper develops a robust method for estimation of production functions and productivity at the firm level. Estimation and inference in this context has been studied by extremely popular papers in applied work, including olleypakes96, levinsohn2003estimating, wooldridge2009estimating, ackerberg2015identification, and gandhi2020identification. We contribute to this literature by allowing for the use of off-the-shelf nonparametric estimation routines to handle first-steps in a semiparametric version of this model, while maintaining the ability to perform standard inference on all marginal productivities and related parameters. cha2023inference have recently explored the use of debiased machine learning in this setting. We departure from this paper in two important ways. First, we do not exploit the score in robinson1988root and this allows us to estimate all the parameters of the production function simultaneously without treating the capital marginal product as a nuisance component. In contrast, cha2023inference can only run inference for the labor coefficient with their strategy. More importantly, the Robinson's score might lead to identification failures in this context, as discussed by ackerberg2015identification. Second, we simplify the implementation of machine learning procedures to estimate production functions by invoking several well-known routines already implemented in R or Python, subject to some rate conditions. Differently, cha2023inference estimate the first step using their specific algorithm, but they are able to relax exact sparsity conditions. doraszelski2025productionfunctionestimationinvertibility have independently obtained a Neyman-orthogonal moment in this context by deriving the first-step influence function and imposing restrictions among the conditioning variables. Our approach is more general as we rely on automatic Neyman-orthogonalization, without restricting the conditioning information sets of the model.

Notation: For an arbitrary vector $x \in \mathbb{R}^r$, let $\left|\left|x\right|\right|_1$ and $\left|\left|x\right|\right|$ be the $\ell_1$ and $\ell_2$ norm, respectively. If $A$ is a matrix, $\left|\left|A\right|\right|$ denotes spectral norm. For a measurable function $a$ of an arbitrary random variable $H$, let $\left|\left|a\right|\right|_2 = \sqrt{\mathbb{E}\left[a(H)^2\right]}$. Let $P$ be a probability distribution. Let $L^2$ be the space of functions of $W$ that are square-integrable, when $ W \sim P$, where the precise meaning of $W$ will be established below.\footnote{Technically, we should index $L^2$ by some $\sigma-$finite measure. We avoid this to simplify our notation.} In addition, let $L^2_0$ be a subset of $L^2$ with the additional mean-zero restriction. Similar definitions apply to objects such as $L^2(V)$ and $L^2_0(V)$ for functions of $V$. For a $n \times r$ matrix $A = \left[a_{ik}\right]$, let $\left|\left|A\right|\right|_{\infty} = \underset{i,k}{\operatorname{max}} \left|a_{ik}\right|$, and $\left|\left|A\right|\right|_{\ell_{\infty}} = \underset{i}{\operatorname{max}} \sum^r_{k=1}\left|a_{ik}\right|$. For a set of indexes $S \subseteq \left\{1,\cdots,r\right\}$, let $x_S$ be the modification of $x$ that places zeros in all entries of $x$ whose indexes do not belong to $S$. Moreover, let $S_x$ be a subset of $\left\{1,\cdots,r\right\}$ such that $x_k \neq 0$ for all $k \in S_x$, and $S^c_{x}$ be the complement of $S_x$ in $\left\{1,\cdots,r\right\}$. For a bounded set $Q$, $\# Q$ denotes the cardinality of $Q$. For $a, b \in \mathbb{R}$, $a \vee b = \max\left\{a,b\right\}$. For any arbitrary subset $\mathcal{K}$, let $\overline{\mathcal{K}}$ denote the closure of $\mathcal{K}$ and $\overline{\mathcal{K}}^{\perp}$ be its orthocomplement, when a topology and inner product $\left<\cdot,\cdot\right>$ is defined. Moreover, $\Pi_{\overline{\mathcal{K}}}$ is the orthogonal projection operator onto $\overline{\mathcal{K}}$. Let $J < \infty$ and $Z = \left(Z_1, Z_2, \cdots,Z_J\right)^{'}$, we say that a vector-valued function $f(Z) = \left(f_1(Z_1), \cdots, f_J(Z_J)\right)^{'}$ is in $L^2(Z) \equiv \bigotimes^J_{j=1} L^2(Z_j)$ when each of the elements in a such vector belongs to the corresponding $L^2\left(Z_j\right)$, $j=1,\cdots,J$, and we define $\left|\left|f\right|\right|_{L^2(Z)} := \sqrt{\sum^J_{j=1}\left|\left|f_j \right|\right|^2_{2}}$. Finally, all the CMRs in the sequel are satisfied almost surely (a.s.), but we will not make it explicit often to simplify the exposition.

Production Functions at the Firm Level: A Motivating Example

To illustrate the key ideas of this paper we will consider an interesting setting throughout. We emphasize that this is purely for the sake of clarity in exposition. Our theoretical framework is broad and applicable to models defined by general semiparametric CMRs, as we will show in Section (ref).

Example: Production Functions at the Firm Level. We observe a panel of $n$ firms across $T$ periods, where $i$ and $t$ index firms and periods, respectively. Let $Y_{it}$ be the output of firm $i$ at time $t$, and $X_{it}$ be a vector of inputs, e.g., capital and labor. Output is determined by the following equation:\footnote{Unless otherwise stated, all the variables are expressed in logarithms.}

where $F$ is assumed to be known up to $\theta_{0p}$, where $\theta_{0p} \in \mathbb{R}^{p_1}$, $\omega_{it}$ is firm $i$'s productivity shock (known as anticipated productivity) in period $t$, which is allowed to be correlated with inputs, and $\epsilon_{it}$ is noise in output, which is independent and identically distributed (iid), is assumed to be independent of the current and previous optimal decisions of the firm and anticipated productivities, and with zero mean. Since $\omega_{it}$ is not observed and is correlated with inputs, OLS will provide inconsistent estimates.\footnote{In the case where $F$ is assumed to be linear in inputs.} To address this, we follow the so-called proxy variable approach, started by olleypakes96; see also levinsohn2003estimating and wooldridge2009estimating. We assume that there exists some firm's choice, $I_{it}$, at $t$ that is linked to $\omega_{it}$: $$ I_{it} = I_t\left(\omega_{it}, X_{it}\right). $$ The precise meaning of the variable $I_{it}$, a “proxy", differs across different formulations of the model.\footnote{For instance, olleypakes96 consider in $I_{it}$ the firm's current investment towards future physical capital. In levinsohn2003estimating, $I_{it}$ is the firm's choice of an intermediate input, e.g., electricity or material input.} In addition, let us assume that $I_t$ is strictly monotonic in $\omega_{it}$. No parametric assumptions are imposed on $I_t$. Then, we shall write $$ \omega_{it} = \omega_t\left(I_{it}, X_{it}\right), $$ where $\omega_t$ is also non-parametric. Hence, we were able to express the unobservable productivity in terms of observable inputs. It is immediate that Equation (ref) becomes

Let $$ \eta_{0t}\left(I_{it}, X_{it}\right) = F\left(X_{it},\theta_{0p}\right) + \omega_t\left(I_{it}, X_{it}\right). $$ Then, by the independence and zero mean assumptions,

Note that (ref) does not necessarily exhaust all the information in the model, e.g., we can also condition on $I_{i-t-1}$ and $X_{i,t-1}$. We write this CMR in this way as it is what researchers typically do in applied work, presumably for simplicity; see ackerberg2014asymptotic, Section 4.2. While Equation (ref) identifies $\eta_{0t}$, it is not enough to identify all the parameters of the production function. This is true since $X_{it}$ enters parametrically and non-parametrically in (ref). The last element of the model is the evolution of $\omega_{it}$. Often, this is also treated non-parametrically, as in olleypakes96 and levinsohn2003estimating. We, nonetheless, follow ackerberg2014asymptotic and work with a more “natural" semiparametric model. Let us assume that $\omega_{it}$ follows a First-Order Markov process in the sense that

Equation (ref) is indeed assumed by olleypakes96 and levinsohn2003estimating. What ackerberg2014asymptotic suggests is to parameterize (ref). To keep things simple, let us consider

where $-1 < \theta_{0\omega} < 1$. Then, letting $\Omega_{it}$ be the firm information set at $t$, with $\left(X_{it}, I_{it}\right)$ being a (measurable) function of $\Omega_{it}$, and using the independence assumption, equations (ref), (ref), and (ref), it is not difficult to show that

Suppose that $T=3$, i.e., firms are observed during three periods. Ignoring subscript $i$, and denoting $\eta_{0} = \left(\eta_{01}, \eta_{02}\right)^{\prime}$, the model can be defined by the following CMRs:

Our goal is to learn $\theta_0 = \left(\theta_{0p}^{'},\theta_{0\omega}\right)^{'}$, the parameter of interest, in the presence of an unknown $\eta_0$. To this end, we might exploit (ref)-(ref) to estimate $\eta_{0}$, and then plug the resulting estimator into moments based on (ref)-(ref) for estimation of $\theta_0$. $\square$

In the example above, $\eta_0$ is unknown and ultimately needs to be estimated to estimate $\theta_0$. The traditional literature on production functions estimation estimates the nonparametric component of the model using sieves methods chen2007large, especially linear sieves, as done in olleypakes96, levinsohn2003estimating, and ackerberg2015identification.\footnote{In olleypakes96 and levinsohn2003estimating, the vector of inputs $X_t$ is such that $X_t = \left(X^{\prime}_{1t}, X^{\prime}_{2t}\right)^{\prime}$, and the link function depends on $I_t$ and the subvector $X_{1t}$, so that $\eta_{0t}$ depends on $I_{it}$ and $X_{1it}$ only.} We instead seek to model $\eta_0$ using machine learning tools, which allow us to exploit a wide range of algorithms to estimate $\eta_0$, different from that of the sieves approach. Those tools (e.g., Lasso, random forest, neural networks, and boosting) have proven useful for approximating functions (where the leading case involves conditional expectations) without imposing stringent parametric assumptions. While the sieves approach provides a convenient theoretical and empirical tool, it might be less flexible in modeling unknown nuisance components than machine learning tools. For example, the sieves approach requires knowledge of all the terms necessary to approximate an unknown function, a requirement that might not always be plausible. In contrast, tools like Lasso aim to identify those terms from the data itself---similar notions apply to other machine learning techniques. Therefore, we aim to equip empirical researchers with an additional inference procedure, based on a less parametric modeling for first stages. This is important for two reasons. First, it provides a valuable tool for scholars who prefer to avoid imposing restrictive assumptions. Second, it offers a way to conduct robustness checks for those more confident in traditional (somewhat less flexible) approaches.

Nevertheless, by employing those algorithms, it would be difficult to estimate $\theta_0$ without bias. This is a by-product of the regularization techniques that all those methods impose in estimation. More importantly, this bias typically decays at a rate slower than $\sqrt{n}$ chernozhukov2022locally. Hence, a naive estimator that directly plugs the estimated $\hat{\eta}$ in estimating $\theta_0$, could lead to a biased estimator, say $\hat{\theta}_{Naive}$. Then it would not hold that $\sqrt{n}\left(\hat{\theta}_{Naive}- \theta_0\right )$ is asymptotically normally distributed. Indeed, $\sqrt{n}\left(\hat{\theta}_{Naive}- \theta_0\right )$ would not be $O_p(1)$, invalidating standard inference on $\hat{\theta}_{Naive}$. To illustrate the point, Figure (ref) displays (in orange) the empirical distribution of the standardized $\left(\hat{\theta}_{Naive} - \theta_0\right)$ for one of the parameters of the production function in our example obtained from simulated data on 1,000 firms and 500 Monte Carlo repetitions. This is computed by estimating $\eta_0$ using boosting in the first stage and then estimating $\theta_0$ using GMM.\footnote{For this illustration, we have used the R function gbmt and its default options.} The figure illustrates that this estimator suffers from bias, as its mean is not centered around zero. Furthermore, the shape of the distribution is certainly different from the standard normal distribution (depicted by the red curve), which would be similar to the expected asymptotic distribution if the bias were absent. The situation would likely become worse as the sample size increases.

Our goal is to construct moments that, under regularity conditions, can be used to estimate $\theta_0$, and that will lead to an estimator that remains $\sqrt{n}-$consistent, despite $\eta_{0}$ being estimated flexibly. Such moments, known as LR/orthogonal/debiased moments, play a crucial role in achieving this goal. A key object in this framework will be a function that we denote Orthogonal Instrumental Variables (OR-IVs) depending on the conditioning variables, $\left(I_1,X_1,\Omega_1,I_2,X_2,\Omega_2\right)$ in the example, that satisfies an important property. We will argue that to compute OR-IVs the researcher needs to solve certain functional equations (cf. Equation (ref)). The terms involved in it as well as the plausibility of solving these equations would be extremely linked to the particular setting that one is working on. Our aim is to construct LR moments in a data-driven manner that can be uniformly applied across various settings, e.g., production functions. This is the main contribution of this work. Hence, the researcher does not have to derive explicit expressions for such functions. The algorithm will provide her with a suitable adjustment term that can subsequently be employed in the estimation of $\theta_0$. This approach results in an estimator $\hat{\theta}_{DGMM}$, which is debiased and yields a $\sqrt{n}\left(\hat{\theta}_{DGMM}- \theta_0\right)$ normally distributed. Figure (ref) displays the standardized empirical distribution of such an estimator (in blue). Notably, it closely resembles the expected asymptotic distribution and, crucially, is centered around zero, implying that our construction is less affected by the first-stage bias.

Debiased Moments and Heuristics

We provide a precise notion of an LR moment for $\theta_0$ in our setting. Let $\kappa_0 \in L^2(Z)$ be some vector-valued function, which will be characterized below. Let $\eta_{0} \in \mathbf{B}$, where $\mathbf{B}$ is a possibly infinite dimensional vector space. A debiased moment in our setting is a moment based on a function $\psi: \mathcal{W} \times \Theta \times \mathbf{B} \times L^2(Z) \mapsto \mathbb{R}$ satisfying the following two restrictions:

where $\frac{d}{d\tau}$ denotes derivatives from the right (i.e., from non-negatives values of $\tau$), at $\tau = 0$ . Equation (ref) implies that the Gateaux derivative at $\eta_0$ of a moment based on $\psi$ is zero. Here $b$ represents a possible direction of deviation from $\eta_0$. Intuitively, local perturbations to the nuisance parameter do not affect the moment. This is an appealing property as it would be hard to learn $\eta_0$ exactly, which is particularly true in high-dimensional contexts. This is the reason why estimation of $\theta_0$ based on a moment that is orthogonal is less affected by first-stage bias, typically present in estimating $\eta_0$, and will be a key attribute to establish standard inference results on an estimator of $\theta_0$. As we will show below, we need to introduce the nuisance parameter $\kappa_0$ to achieve (ref). Nonetheless, this would not affect the moment itself as Equation (ref) implies that the expectation of $\psi$ is globally insensitive to deviations from $\kappa_0$.

Now, how can we construct moments that satisfy the key properties (ref)-(ref) in the context of our production function example? The basic idea will be to obtain a function $\kappa_0$ that is orthogonal to a set that contains all the possible derivatives of the CMRs with respect to $\eta$. This corresponds to a more general result of orthogonality for general functionals in models defined by several CMRs, studied by arganaraz2025machine.

Example: Let $W = \left(Y_1,Y_2,Y_3,X_1,X_2,X_3,I_2,I_2\right)$ be the vector of outputs, inputs, and proxy variables. For this case, given some smoothness conditions, we can obtain a debiased moment by means of

where $Z_1 = \left(I_1,X_1\right)$, $Z_2 = \left(I_2,X_2\right)$, $Z = (Z_1,Z_2)$, $\kappa_0 = \left(\kappa_{01}, \kappa_{02}, \kappa_{03}, \kappa_{04}\right)\in L^2\left(Z_1\right) \times L^2\left(Z_1\right) \times L^2\left(Z_2\right) \times L^2\left(Z_2\right)$ is such that

Then, we can verify that (ref) indeed satisfies (ref)-(ref). $\square$

Expression (ref) indicates that we can obtain debiased moments as the sum of the products between the residual functions and the elements in $\kappa_0$. This linear combination is special in the sense that we need to find $\kappa_0 \in L^2(Z)$ such that the moment is invariant to local perturbations to $\eta$ around $\eta = \eta_0$. Therefore, by construction, the resulting moment (ref) is debiased. Since $\kappa_0$ yields an orthogonal moment by properly combining the initial residual functions, chernozhukovetaldebiased2016 refer to these as Orthogonal Instrumental Variables (OR-IVs).

Obtaining a valid $\kappa_0$ entails solving the functional equation (ref). Various approaches can be pursued to achieve this goal. One option is to directly solve (ref) and find a closed from expression for $\kappa_0$, as arganaraz2025machine did for the partially linear model with endogeneity. This strategy, although convenient, may not always be ideal, as its plausibility depends on the specific expressions involved in the functional equations defining $\kappa_0$. Consequently, there is no guarantee that a procedure for solving for $\kappa_0$ in one model can be readily applied in other contexts. For instance, solving Equation (ref) requires finding four terms, and each pair involves functions of the same random variables. In more complex settings, additional terms may emerge. Instead of relying on the tractability of the equations defining the OR-IVs, we propose an approach that can be applied generally. In addition, as it is evident from (ref), these equations might involve unknown quantities, making direct computations potentially infeasible. Our goal is to design a feasible algorithm that addresses unknown quantities intelligently.

Our algorithm converts any suitable vector of instruments $f \in L^2(Z)$, where $Z$ collects all the distinct random elements in the conditioning variables, into a unique $\kappa_0$. In other words, our approach automatically transforms a function of the conditioning variables such that it satisfies (ref). To illustrate the core idea of our procedure, let us again focus on our example. We start with a known function $f = \left(f_1, f_2, f_3, f_4\right) \in L^2\left(Z\right)$, chosen by the researcher. For instance, $f(Z)=\left(K_{1}, K_{1}, I_{2}, I_{2}\right)^{\prime}$---where $K$ denotes capital---provided that these entries have finite second moments. Under some regularity conditions and for suitable vectors $\hat{M}_1$, $\hat{M}_2$, $\hat{M}_3$, and $\hat{M}_4$, our algorithm finds some finite-dimensional vector $\hat{\beta}$ such that an estimator $\hat{\kappa} = \left(\hat{\kappa}_1, \hat{\kappa}_2, \hat{\kappa}_3, \hat{\kappa}_4\right)$ is constructed as follows:

and satisfies (ref), with probability approaching one, where we replace population objects with their sample counterparts. In particular, $\hat{\beta}$ is interpreted as the solution to a Lasso-type problem with generated regressors $\hat{M}_1$-$\hat{M}_4$. Essentially, the vector $\hat{\beta}$ makes $\hat{\kappa}$, constructed as in (ref)-(ref), orthogonal to a linear operator whose range is equal to all the possible Gateux derivatives of a moment based on $\psi$ with respect to $\eta$ at $\eta = \eta_0$. Obtaining such orthogonalization amounts to minimizing a mean squared error. As we allow the dimension of each of the regressors $\hat{M}_1$-$\hat{M}_4$ to be greater than the sample size, we add a $\ell_1-$norm penalization term. Hence, our problem of finding $\hat{\beta}$ can be seen as solving as a Lasso problem. Under an approximate sparsity condition, we show that if this $\hat{\beta}$ is plugged into (ref)-(ref), a valid $\hat{\kappa}$ and thus a debiased moment can be obtained.

Based on this construction, we derive below a convergence rate for our OR-IV estimator, $\hat{\kappa}$. Let $\hat{\eta}$ be a suitable machine learning estimator of $\eta_0 \in \Xi \subseteq \mathbf{B}$.\footnote{While establishing asymptotic results, the distinction between $\Xi$ and $\mathbf{B}$ is important in endogenous settings, as we clarify below.} We provide below conditions under which $$ \sqrt{n}\left|\left|\hat{\eta} - \eta_0\right|\right|_{\Xi}\left|\left|\hat{\kappa} - \kappa_0\right|\right|_{L^2(Z)} \rightarrow 0, $$ where $\left|\left|\cdot \right|\right|_{\Xi}$ is a suitable norm in $\Xi$. Note that the previous display allows for rates slower than $n^{-1/2}$ for $\hat{\eta}$ and $\hat{\kappa}$. Under this key condition, we prove asymptotic normality of a two-step GMM estimator $\hat{\theta}$, implying that standard inference can be conducted straightforwardly. Debiasing plays a pivotal role in this derivation. The following sections develop these ideas in a general framework, derive the exact program that $\hat{\beta}$ solves, and provides the technical conditions required. Readers less interested in these details can skip these sections, except for Section (ref), where we present our estimator.

Computation of the OR-IVs

General Setting

The data $W_i = \left(Y_i, V_i, Z_i\right)$, $i = 1,\cdots,n$, is iid with the same distribution as $W = \left(Y,V,Z\right)$, where $Y$ has support $\mathcal{Y} \subseteq \mathbb{R}^{d_Y}$. Let $\theta \in \Theta \subset \mathbb{R}^{d_\theta}$ denote a finite-dimensional parameter vector. Let $\eta \in \Xi \subseteq \mathbf{B}$ be a vector of real-valued measurable functions of $V$ that might depend on additional unknown parameters that we do not specify, and $\mathbf{B} \subseteq \bigotimes^{d_\eta} L^2\left(V\right)$ is a Hilbert space with inner product $\left<\cdot,\cdot \right>_{\mathbf{B}}$. To be specific, $\eta = \left(\eta_{1}, \cdots, \eta_{d_\eta}\right)$ with $\eta_{s} \equiv \eta_{s}\left(V\right)$.\footnote{With this notation, we also allow for cases where $\eta_s$ depends on a subvector of $V$.} We assume that there is a vector of residual functions $m_j: \mathcal{Y} \times \Theta \times \mathbf{B} \mapsto \mathbb{R}$ such that:

where $\mathbb{E}[\cdot]$ is expectation under the distribution of $Y$ given $Z_j$, $Z$ denotes the union of distinct random elements of the conditioning variables $\left(Z_1,\cdots,Z_J\right)$, and each $m_j$ is known up to the parameters $(\theta_0,\eta_0)$, with finite second moment. To be precise, $m_j$ might depend on $\theta_0$ arbitrarily and the entire function $\eta_0$, not only on its evaluation at a particular realization of $V$. Observe that we are not imposing anything regarding how the conditioning variables relate. These might have all or some elements in common. Note that in the case where $Z_j$ is a constant, for some $j$, we have an unconditional moment. Hereafter, we assume that there exists a unique $\left(\theta_0,\eta_0\right) \in \Theta \times \mathbf{B}$ such that (ref) holds. We say the vector $Z$ gathers the exogenous variables in the model, which might include $V$ or not.

Let $\kappa = \left(\kappa_1, \cdots, \kappa_J\right)$, where $\kappa_j \equiv \kappa_j\left(Z_j\right)$, and $\kappa_j \in L^2(Z_j)$, $1 \leq j \leq J$. Hence, $\kappa \in L^2(Z)$, where $L^2(Z) = \bigotimes^J_{j=1} L^2(Z_j)$. We define $$ h_j\left(Z_j, \theta, \eta\right) := \mathbb{E}\left[\left.m_j\left(Y,\theta,\eta\right) \right|Z_j\right]. $$ Our results will rely heavily on smoothness conditions of the functions $h_j$'s. In particular, we maintain the key assumption throughout:

In our notation above, $b$ denotes a specific direction of deviation and $\tau \in \left[0, \varepsilon\right)$ is the associated length. We next provide conditions in our specific example under which Assumption (ref) holds.

Example: The following lemma states sufficient conditions for the CMRs to be Fréchet differentiable.

Remark that (ref) defines a linear operator $S^{(j)}_{\theta_0,\eta_0}: \mathbf{B} \mapsto L^2(Z_j)$; see luenberger1997optimization and carrasco2000generalization for a theory on linear operators. In addition, let us define $$ S_{\theta_0,\eta_0}b := \left(S^{(1)}_{\theta_0,\eta_0}b, \cdots, S^{(J)}_{\theta_0,\eta_0}b\right). $$ Then, $S_{\theta_0,\eta_0}: \mathbf{B} \mapsto L^2(Z)$ is also a linear operator. We equip $L^2(Z)$ with the inner product $$ \left<f_1,f_2\right>_{L^2(Z)} = \sum^J_{j=1} \mathbb{E}\left[f_{1j}(Z_j)f_{2j}(Z_j)\right], $$ where $f_s = \left(f_{s1},\cdots,f_{sJ}\right)$, $s=1,2$. Therefore, $L^2(Z)$ is a Hilbert space. The range of that operator can be defined as follows: $$ \mathcal{R}\left(S_{\theta_0,\eta_0}\right) := \left\{f \in L^2\left(Z\right): f = S_{\theta_0,\eta_0}b\;\; \text{for some}\;\; b \in \mathbf{B}\right\}. $$ A key object for us is $\overline{\mathcal{R}\left(S_{\theta_0,\eta_0}\right)}^{\perp}$, i.e., the orthogonal complement of the closure of the range of $S_{\theta_0,\eta_0}$ in $L^2(Z)$, which can be defined as $$ \overline{\mathcal{R}\left(S_{\theta_0,\eta_0}\right)}^{\perp} :=\left\{f \in L^2\left(Z\right): \sum^J_{j=1}\mathbb{E}\left[f_j\left(Z_j\right)h_j\left(Z_j\right)\right] = 0, \;\;\;\text{for all}\;\;\; h = \left(h_1,\cdots,h_J\right) \in \overline{\mathcal{R}\left(S_{\theta_0,\eta_0}\right)}\right\}. $$ arganaraz2025machine show that a debiased moment in model (ref) can be constructed as follows:

We point out two important observations. First, $\kappa_0$ might not be unique. In fact, there can potentially exist an infinite number of such $\kappa_0$'s. Despite this, not every IV, i.e., functions in $L^2\left(Z\right)$, belongs to $\overline{\mathcal{R}\left(S_{\theta_0,\eta_0}\right)}^{\perp}$, and thus not every IV may serve as a valid $\kappa_0$. Consequently, choices of instruments commonly made in applied work might not lead to orthogonal moments. For example, $\left(K_{1}, K_{1}, I_{2}, I_{2}\right)^{\prime}$ might not be OR-IVs in our example. Nevertheless, we will demonstrate how to convert suitable functions in $L^2\left(Z\right)$, e.g.,$\left(K_{1}, K_{1}, I_{2}, I_{2}\right)^{\prime}$, into a $\kappa_0$. Therefore, one may begin with the common choices of instrument functions and then apply our transformation to directly obtain an orthogonal moment.

Second, a $\kappa_{0}$ does not necessarily exist. This situation occurs when $\overline{\mathcal{R}\left(S_{\theta_0,\eta_0}\right)}^{\perp} = \left\{0\right\}$. Following the terminology in arganaraz2025existence, when $\overline{\mathcal{R}\left(S_{\theta_0,\eta_0}\right)}^{\perp} = \left\{0\right\}$ we say the model satisfies a local surjectivity property.\footnote{Similar notions have appeared elsewhere for certain panel data models, e.g., in bonhomme2012functional.} A practical implication of this situation is that this would result in a trivial orthogonal moment, i.e., $\psi = 0$. Unfortunately, the plausibility of encountering this situation depends on the particular model at hand.\footnote{Using our approach for estimating OR-IVs, we may check for local surjectivity in a practical way. Let $\hat{\kappa}$ be an estimator of $\kappa_0$. Local surjectivity implies $\hat{\kappa} \approx 0$ for each individual in the sample for a large class of initial functions $f$'s. An alternative (and more rigorous) approach consists of using duality theory and exploiting the fact that $\overline{\mathcal{R}\left(S_{\theta_0,\eta_0}\right)}^{\perp} = \mathcal{N}\left(S^{*}_{\theta_0,\eta_0}\right)$, where $S^{*}_{\theta_0,\eta_0}$ is the adjoint operator of $S_{\theta_0,\eta_0}$ and $\mathcal{N}\left(\cdot\right)$ denotes the null space of an operator luenberger1997optimization. Theoretically checking that $\mathcal{N}\left(S^{*}_{\theta_0,\eta_0}\right) = \left\{0\right\}$ might be easier than studying if $\overline{\mathcal{R}\left(S_{\theta_0,\eta_0}\right)}^{\perp} = \left\{0\right\}$.}

Estimation of OR-IVs

We next explain how to obtain OR-IVs automatically in the general model (ref). As we explained above, the key step to construct an orthogonal moment is to find a function $\kappa_{0} \in L^2\left(Z\right)$ such that it is orthogonal to $\overline{\mathcal{R}\left(S_{\theta_0,\eta_0} \right)}$. This implies that $$ \sum^J_{j=1}\mathbb{E}\left[ \left[S^{(j)}_{\theta_0,\eta_0}b\right](Z_j)\kappa_{0j}\left(Z_j\right)\right] = 0, $$ for all $b \in \mathbf{B}$. We can obtain such $\kappa_0$ by selecting some function $f \in L^2\left(Z\right)$, and then computing $$\kappa_{0} = f - \Pi_{\overline{\mathcal{R}\left(S_{\theta_0,\eta_0}\right)}}f, $$ where recall that $\Pi_{\overline{\mathcal{R}\left(S_{\theta_0,\eta_0}\right)}}$ denotes the orthogonal projection operator onto $\overline{\mathcal{R}\left(S_{\theta_0,\eta_0}\right)}$, and is defined as follows

By the Projection Theorem, (ref) exists and is unique luenberger1997optimization. Hence, for a given $f$, we have a unique $\kappa_0$ a.s. Next, we exploit the following facts. Notice that $\Pi_{\overline{\mathcal{R}\left(S_{\theta,\eta_0}\right)}}f \in \overline{\mathcal{R}\left(S_{\theta_0,\eta_0}\right)}$, and that the range of $S_{\theta_0,\eta_0}S^{*}_{\theta_0,\eta_0}$ is dense in $\overline{\mathcal{R}\left(S_{\theta_0,\eta_0}\right)}$, where $S^{*}_{\theta_0,\eta_0}$ is the adjoint operator of $S_{\theta_0,\eta_0}$. Since $S_{\theta_0,\eta_0}$ is bounded, $S^{*}_{\theta_0,\eta_0}: L^2(Z) \mapsto \mathbf{B}$ exists, is linear, and bounded CARRASCO20075633. Moreover, it satisfies, for all $g \in L^2(Z)$ and $b \in \mathbf{B}$, $$ \left<S_{\theta_0,\eta_0}b, g\right>_{L^2(Z)} = \left<b, S^{*}_{\theta_0,\eta_0}g\right>_{\mathbf{B}}. $$ Let $\Pi_{\overline{\mathcal{R}\left(S_{\theta_0,\eta_0}\right)}}f = f^{*}$, then the previous observations imply that $f^{*} \in \overline{\mathcal{R}\left(S_{\theta_0,\eta_0}S^{*}_{\theta_0,\eta_0}\right)}$. This is potentially useful as $S_{\theta_0,\eta_0}S^{*}_{\theta_0,\eta_0}$ is an operator that maps elements in $L^2(Z)$ into $L^2(Z)$, spaces of exogenous variables in the data. Moreover, $S_{\theta_0,\eta_0}S^{*}_{\theta_0,\eta_0}$ can be characterized by exploiting the CMRs given by the model in a general way, as we show below. In general, we are interested in solving

From a practical standpoint, problem (ref) presents two difficulties. First, we might not able to directly compute $S^{(j)}_{\theta_0,\eta_0}S^{*}_{\theta_0,\eta_0}$ as it might involve unknown objects. Hence, we shall construct an estimator for it.

Second, the problem implied by (ref) requires a search over an infinite-dimensional space, which would be challenging to directly accomplish in practice. To tackle this issue, we exploit a particular regularization strategy. We will look for a solution to (ref) in a set $\mathcal{G}_n$ defined as follows. Let $\gamma\left(Z\right) = \left(\gamma_1\left(Z_1\right)^{'}, \cdots, \gamma_J\left(Z_J\right)^{'}\right)^{'} \in L^2(Z)$ be a vector of basis functions, where $\gamma_j\left(Z_j\right)$ is a $r_j-$dimensional vector of known real-valued functions, with $\mathbb{E}\left[\gamma_{jk}\left(Z_j\right)\right] = 0$ and $\mathbb{E}\left[\gamma^2_{jk}\left(Z_j\right)\right] = 1$, $k=2,\cdots,r_j$.\footnote{We let the first element of $\gamma_j$ be 1.} Let $\beta = \left(\beta_{11}, \cdots, \beta_{1r_1},\cdots,\beta_{1r_J}, \cdots,\beta_{Jr_J}\right)^{'}$ be a $r-$dimensional vector, where $r = \sum^J_{j=1}r_j$. Then, we define $\mathcal{G}_n$ as\footnote{Extending our results to other functional classes such as Reproducing Kernel Hilbert Spaces or Neural Networks is an interesting avenue of future research.} $$ \mathcal{G}_n := \left\{g \in L^2(Z): g_j\left(Z_j\right) = \gamma_j\left(Z_j\right)^{'}\beta_j, \;\;\left|\left|\beta\right|\right|_{\infty} < c \right\}. $$ Remark that since we will look for a solution in $\mathcal{G}_n$ and $\gamma$ is known, we only need to focus on finding an optimal finite-dimensional vector $\beta$. We let $\beta$ be general in the sense that its dimension, $r$, can grow with $n$ and be much larger than it. To this end, we augment the program in (ref) by including an $\ell_1-$norm penalization term on $\beta$. In sum, the need for estimating $S^{(j)}_{\theta_0,\eta_0}S^{*}_{\theta_0,\eta_0}$ and our regularization strategy yields writing a feasible version of (ref) as

where $\gamma(Z)^{'}\beta \equiv \left(\gamma(Z_1)^{'}\beta_1,\cdots,\gamma(Z_J)^{'}\beta_J\right)^{\prime}$, $\hat{S}^{(j)}_{\hat{\theta},\hat{\eta}}\hat{S^{*}}_{\hat{\theta},\hat{\eta}}$ is a suitable estimator of $S^{(j)}_{\theta_0,\eta_0}S^{*}_{\theta_0,\eta_0}$, $\lambda_n \geq 0$ is a regularization term (or tuning parameter) to be determined below, and $\mathbb{E}_n\left[h(Z_i)\right] = \sum^n_{i=1}h(Z_i)/n$. The solution to (ref) is denoted by $\hat{\beta}$.

Let us be more precise about the estimation of $S^{(j)}_{\theta_0,\eta_0}S^{*}_{\theta_0,\eta_0}$. During its estimation, we use cross-fitting, following chernozhukov2018double.\footnote{Cross-fitting has a long tradition in the semiparametric literature; see, e.g., bickel1982adaptive, klaassen1987consistent, vandervaart98, robins2008higher, and Zheng2010AsymptoticTF.} This entails randomly partitioning the sample into $L$ groups, $I_1, \cdots, I_L$, of the same size. Let $I^c_\ell$ be the complement of $I_\ell$. We compute $\hat{\beta}_\ell$, using observations in $I^c_\ell$ only. In other words, we compute, for each $I_{\ell}$, $\hat{\beta}_\ell$ using $\hat{S}^{(j)}_{\hat{\theta}_\ell,\hat{\eta}_\ell}\hat{S^{*}}_{\hat{\theta}_\ell,\hat{\eta}_\ell}$, which is constructed from observations in $I^c_\ell$. Let $n_\ell$ denote the number of observations in $I_\ell$. Hence,

Ultimately, we propose to estimate each component in $\kappa_0$, for each individual $i \in I_\ell$, by means of

The Algorithm

The approach previously developed provide researchers with a general guidance for estimating OR-IVs. We now want to propose a specific algorithm, which can be followed in a wide range of settings that share a common structure in the operator $S_{\theta_0,\eta_0}$. We assume

Assumption (ref) i) implies that there exists a function $\nu_j$ that links the direction $b$ with the operator $S^{(j)}_{\theta_0,\eta_0}$ through the conditional expectation given $Z_j$. In the standard non-parametric IV setting newey2003instrumental, where we have $\mathbb{E}\left[\left.Y - \eta_0\left(V\right)\right|Z_1\right] = 0$, it holds that $\nu_j(Y,\theta_0,\eta_0,b) = -b\left(V\right)$. Hence, i) allows for these derivatives in the first step to depend generally on $\theta_0$ and $\eta_0$. Assumption (ref) ii) requires $\nu_j$ to be linear in the direction $b$ in a specific sense. In the non-parametric IV case, $\tilde{\nu}_j\left(Y,\theta_0,\eta_0\right) = -1$.

As shown in the previous section, Assumption (ref) is not needed, as we could directly work with the generic operators $S^{(j)}_{\theta_0,\eta_0}$'s and its properties. However, in many applications, to find those operators one shall proceed by first finding functions $\nu_j's$ and then showing that their conditional expectations are the corresponding Fréchet derivatives. Additionally, it is more convenient to work with functions $\nu_j$ and $\tilde{\nu}_j$. Our strategy below relies on them to propose a practical algorithm to compute OR-IVs. Notice that if the residual function $m_j\left(Y,\theta_0, \cdot\right): \mathbf{B} \mapsto \mathbb{R}$ is Fréchet differentiable for each fixed $Y = y$ a.s. and its derivative is bounded by some integrable function, Assumption (ref) i) is satisfied, provided that the conditional expectation of such derivative, given $Z_j$, is square-integrable.\footnote{Recall that if $m_j\left(Y,\theta_0, \cdot\right): \mathbf{B} \mapsto \mathbb{R}$ is Fréchet differentiable for each fixed $Y = y$ a.s., then $\frac{d}{d\tau}m_j\left(Y,\theta_0, \eta_0 + \tau b\right) = \nu_j\left(Y,\theta_0,\eta_0,b\right)$ for all $b$, where $\nu_j\left(Y,\theta_0,\eta_0,\cdot\right): \mathbf{B} \mapsto \mathbb{R}$ is linear and continuous. If this expression is bounded by some integrable function, i) is implied by the Dominated Convergence Theorem, provided that $\mathbb{E}\left[\left. \nu_j\left(Y,\theta_0,\eta_0,b\right) \right|Z_j\right] \in L^2(Z_j)$.} Assumption ii) requires that its derivative can be written in the form $ b\left(V\right)^{\prime} \tilde{\nu}_j\left(Y,\theta_0,\eta_0\right)$. Altogether, these conditions simplify our calculations below. We can easily verify that Assumption (ref) holds in our example.

Example: Suppose that the assumptions in Lemma (ref) hold. It is not difficult to obtain that

In addition, $S^{(1)}_{\theta_0,\eta_0}b, S^{(2)}_{\theta_0,\eta_0}b \in L^2(Z_1)$, and $S^{(3)}_{\theta_0,\eta_0}b, S^{(4)}_{\theta_0,\eta_0}b \in L^2(Z_2)$, for any $b$. $\square$

Leveraging Assumption (ref), we can obtain an expression for $S^{*}_{\theta_0,\eta_0}$ that can be applied to various contexts:

We shall now write (ref) using matrix notation. Split $I^c_\ell$ into three mutually exclusive pieces such that $I^c_\ell = A_\ell \cup B_\ell \cup C_\ell$. Possibly, this partition depends on $j$, but we omit this dependence for simplicity. For any estimator, a subscript will indicate the part of $I^c_\ell$ that has been used to compute it. For example, $\hat{\theta}_{A_\ell}$ means that $\theta_0$ has been estimated using observations in $A_\ell$ and so on. Let $\bm{f_{j\ell}}$ be a $n_\ell-$dimensional vector containing each $f_j\left(Z_{ji}\right)$, $i \notin I_{\ell}$. Let $\bm{\hat{M}_{j\ell}}$ be a $n_\ell \times r$ design matrix such that its $(i,l)$-entry is given by

Then, we have $$ \hat{\beta}_\ell = \underset{\beta \in \mathbb{R}^r}{\operatorname{arg\;min}}\; \sum^J_{j=1} \frac{1}{n-n_{\ell}} \left(\bm{f_{j\ell}} - \bm{\hat{M}_{j\ell}}\beta\right)^{'}\left(\bm{f_{j\ell}} - \bm{\hat{M}_{j\ell}}\beta\right) + 2\lambda_n \left|\left|\beta\right|\right|_1, $$ where $\hat{\eta}_{A_\ell}$, $\hat{\eta}_{B_\ell}$, $\hat{\mathbb{E}}_{B_\ell}\left[\left.\cdot\right|V\right]$, and $\hat{\mathbb{E}}_{C_\ell}\left[\left.\cdot\right|Z_j\right]$ are non-parametric estimators, possibly based on some machine learning tool, and $\hat{\theta}_{A_\ell}$ and $\hat{\theta}_{B_\ell}$ are possibly non-LR estimators of $\theta_0$.\footnote{I.e., estimators based on non-LR moments.} Recall that $\tilde{\nu}_j$ is a known function, given estimators of $\theta_0$ and $\eta_0$. Thus, these conditional expectations can be evaluated. Notice that the $\ell_1-$penalization terms allows for $r>n$. As a result, we have transformed our initial problem into a Lasso-type one for which there exist well-known and fast algorithms. Finally, the resulting residuals from such problem is an estimator of an OR-IV. Let $\bm{M_{j}}$ be the population analog of matrix $\bm{\hat{M}_{j\ell}}$. Let $\hat{M}_{j\ell}\left(Z_{ji}\right)$ be a $r-$dimensional vector containing the $i-$ row of $\bm{\hat{M}_{j\ell}}$. A similar definition applies to $M_{j}\left(Z_{ji}\right)$. Our algorithm can be implemented as follows:

Algorithm to estimate OR-IVs:

Step 0: Choose a real-valued function $f \in L^2(Z)$. Choose a basis for each $\gamma_j(Z_j)$, e.g., exponential, Fourier, splines, or power. In addition, specify a low-dimensional dictionary, say $\gamma^{low}(Z)$, which is a sub-vector of $\gamma(Z)$.

Step 1: For each $\ell = 1,\cdots L$, compute (possible) non-LR estimators $\hat{\theta}_{A_\ell}$ and $\hat{\theta}_{B_\ell}$. Moreover, using some machine learning algorithm, compute $\hat{\eta}_{A_\ell}$, $\hat{\eta}_{B_\ell}$, $\hat{\mathbb{E}}_{B_\ell}\left[\left.\cdot\right|V\right]$, and $\hat{\mathbb{E}}_{C_\ell}\left[\left.\cdot\right|Z_j\right]$.

Step 2: Construct design matrix $\bm{\hat{M}_{j\ell}}$ such that its $(i,l)-$entry is (ref).

Step 3: Initialize $\hat{\beta}_\ell$ using $\gamma^{low}(Z)$ such that

Step 4: (While $\hat{\beta}_\ell$ has not converged)

Step 5: Given the optimal $\hat{\beta}_\ell$, compute $\hat{\kappa}_{j\ell}$ as:

In Step 1, we need to specify a low-dimensional dictionary $\gamma^{\text{low}}(Z)$ such that $\sum^J_{j=1}\bm{\hat{M}_{j\ell}^{\prime}\hat{M}_{j\ell}}$ is invertible in Step 3. In our Monte Carlo experiment below, we take the first 5 components of each $\gamma_j(Z_j)$. In our empirical application, we take the first $\tilde{r}_j$ of each $\gamma_j(Z_j)$, where $\tilde{r}_j = r_j/15$.\footnote{How large $\tilde{r}_j$ is depends on the particular setting. Ideally, $\tilde{r}_j$ shall be set as large as possible as long as $\sum^J_{j=1}\bm{\hat{M}_{j\ell}^{\prime}\hat{M}_{j\ell}}$ is non-singular.} Following belloni2012sparse, we include a normalization term, $\hat{D}_{jk\ell}$, in the $\ell_1-$norm above, which is necessary for the good properties of our Lasso estimator. For this same reason, we suggest computing $\lambda_n$ as given in Step 4, as recommended by belloni2012sparse (p. 2380). However, we let $\lambda_n$ be proportional to $n^{-1/4}$ instead of $n^{-1/2}$, as required by our asymptotic theory. Regarding the constants in $\lambda_n$, our Monte Carlo exercises have suggested $c_1 = 1.1$. Several numerical experiments have indicated that a value of at least $c_2 = 0.5/\log(n \vee r)$ works well.\footnote{belloni2012sparse recommend $c_1 = 1.1$ and $c_2 = 0.1/\log(n \vee r)$. We have observed that, in our context, belloni2012sparse's recommendation leads to “too much" regularization. More precisely, their recommendation tends to shrink all coefficients in $\beta$ to zero, which might be explained by the fact that $\lambda_n$ has $\left(n - n_\ell\right)^{1/4}$ in the denominator instead of the larger $\left(n - n_\ell\right)^{1/2}$. To accommodate for this, we let $c_2$ to be larger as well. On strictly theoretical grounds, $c_1$ shall be greater than 1 and $c_2$ needs to be such that it converges to zero and $\log(1/c_2) \leq C \log(n \vee r)$, which is satisfied by our current choices.} To improve numerical stability we follow chernozhukov2022automatic and cap the maximum number of iterations at 10. In addition, we use warm start. This means that in a given iteration, the initial parameter value is equal to the $\hat{\beta}_\ell$ obtained in the previous iteration. In Section (ref) of the Supplementary Appendix we provide a justification for the previous optimization. Notice that Step 4 (b) requires solving for $\hat{\beta}_\ell$. For this, we use an extension of the coordinate descent approach for Lasso Friedman07, friedman2010regularization, fu1998penalized; see Section (ref) of the Supplementary Appendix for details. An important remark is that our theory (and several numerical experiments we conducted) suggest all these choices, as long as the requirements of our asymptotic theory are met, does not play an important role as inference on $\theta_0$ is concerned, which is a result of the OR-IVs being, by construction, in the space of IVs of the model.

We point out that our previous algorithm applies to the most general case. In some applications, it will be implemented with some simplifications, depending on the particular expression of $\tilde{\nu}_j$, and how $V$ and the conditioning variables relate. For example, it might not be necessary to partition $I_\ell$ into three pieces, as it might not be necessary to compute all the estimators in Step 1. This occurs, for example, when $\tilde{\nu}_j$ does not depend on $\theta_0$, or $\eta_0$, or both. Another case is when $\mathbb{E}\left[\left.\tilde{\nu}_j\left(Y,\theta_0,\eta_0\right)\gamma_{jk}(Z_{ji})\right|V\right] = \tilde{\nu}_j\left(Y,\theta_0,\eta_0\right)\gamma_{jk}(Z_{ji})$ as $\left(Y,Z_j\right)$ and the variables that $\eta_0$ depends on are contained in $V$. More simplifications can emerge if $\tilde{\nu}_j$ depends only on $Z_j$. As we will illustrate below, in our example, $\left[\bm{\hat{M}}_{j\ell}\right]_{i\ell}$ has a simpler expression than in (ref).

Example: The user has to provide $f(Z) = \left(f_1(Z_1), f_2(Z_1), f_3(Z_2), f_4(Z_2) \right) \in L^2\left(Z\right)$ and basis $\gamma\left(Z\right) = \left(\gamma_1(Z_1), \gamma_2(Z_1), \gamma_3(Z_2), \gamma_4(Z_2)\right)$. In this example, the regressors have the following expression

As we can see from above, no conditional expectations or estimators of $\eta_0$ appear. Hence, it is not necessary to partition $I^c_{\ell}$. That is, we use all observations in $I^c_{\ell}$ to obtain an estimator of $\theta_{0\omega}$. $\square$

Asymptotic Properties of OR-IVs

This section provides the mean square convergence rate for $\hat{\kappa}$ based on the Lasso estimator introduced above, which is fundamental for deriving the asymptotic properties of $\hat{\theta}$. We define

Then, $\hat{\beta}_\ell$ can equivalently be written as

Interestingly, the previous characterization of the program is similar to the one proposed by chernozhukov2022automatic for automatic estimation of Riesz representers (see Equation (3.7) in this paper), but with different matrices. An important difference is that $\hat{G}_j$, the estimated Gram matrix, depends on regressors $\hat{M}_j$'s that are estimated. This also contrasts with the formulation considered by bakhitov2022automatic. These previous works formulate the estimation problem in terms of known objects analogous to $M_j$'s. On the other hand, notice that we can deal with situations where the parameter of interest is not necessarily an average of some function depending on the vector of nuisance parameters, which is the type of setting considered by chernozhukov2022automatic and bakhitov2022automatic. In what follows, we work with the characterization given in (ref). We start by assuming

Since $\hat{G}_j$ depends on the estimated $\hat{M}_j(Z_j)\hat{M}_j(Z_j)^{'}$, we require a convergence rate for them in a specific sense such that we can assure that $\hat{G}_j$ provides a good approximation to $G_j$. Let $F_0$ be the distribution of the data $W$, then we hypothesize

Assumptions (ref) and (ref) imply that for any $c>0$, $$ \left|\left|\hat{G}_{j\ell} - G_j\right|\right|_{\infty} = O_p\left(\varepsilon^G_{jn}\right),\;\;\;\text{where}\;\;\; \varepsilon^G_{jn} = \max\left\{\bar{\varepsilon}_n, n^c\varepsilon^M_{jn}\right\},\; \bar{\varepsilon}_n = \sqrt{\frac{\log(r)}{n}}. $$ Following similar arguments to the ones in chernozhukov2022automatic and bakhitov2022automatic, who work on situations where regressors are known, we would get $ \left|\left|\hat{G}_{j\ell} - G_j\right|\right|_{\infty} = O_p\left(\bar{\varepsilon}_n\right)$. However, as we deal with unknown regressors, we pay a price in terms of a slower rate for $\hat{G}_{j\ell}$. Hence, the better the rates we can get for $\hat{M}_{j\ell k}(Z_j)$, as required by Assumption (ref), the better the rate for $\hat{G}_{j\ell}$. Remark that $\hat{M}_{j\ell k}(Z_j)$ takes the form of a conditional expectation and as such could be estimated by machine learning tools. This conditional expectation depends on other unknown objects (conditional expectations given $V_i$ and $\eta_0$) that could be estimated with machine learning algorithms as well. Under conditions to be specified below, $\varepsilon^M_{jn}$ will be no faster than the slowest convergence rate among all those other unknown objects in $\hat{M}_{j\ell k}(Z_j)$. The only requirement that we will impose is that such a convergence rate has to be faster than $n^{-1/4}$ but could be slower than $n^{-1/2}$, a reasonable assumption for several machine learners; see for example, neural networks chen1999improved, schmidt2020, farrel21, random forests syrgkanis2020estimation, biau2012analysis, Lasso bickellasso, buneatsybakovwegkamp, and boosting luo2022highdimensional.

In our discussion to follow, let $\varepsilon_n$ be a non-negative sequence that converges to zero no faster than $\varepsilon^G_{jn}$, for any $j$. We impose a sparse approximate condition on the orthogonal projection $f^{*}$, a key assumption for us.

This assumption is arguably standard in the high-dimensional literature; see buneatsybakovwegkamp, belloni2012sparse, belloni2016inference, bradic2022minimax. It controls the squared approximation error from using the linear combination $M_j^{'}\bar{\beta}$ to approximate the orthogonal projection. Note that the rate at which (the square of) this approximation error shrinks is governed by $\varepsilon_n$, which will ultimately influence the rate of OR-IVs, as we show below. Remark that Assumption (ref) does not impose that $f^{*}$ can be written as a linear combination of $s$ terms, i.e., that the orthogonal projection is strictly sparse. Instead, Assumption (ref) only requires the existence of $\bar{\beta}$ with $s$ terms such that the mean of the squared approximation errors across the $J$ elements is bounded by $Cs\varepsilon^2_n$. Moreover, the above assumption does allow the unknown identity of the elements in $M_j$ that give a good approximation, i.e., the researcher does not have to specify which elements are important, a task that will be typically hard to accomplish as we are dealing with highly complex functional objects; see bradic2022minimax. We will see below that a very sparse approximation, with a small number of terms $s$ (and a large $\xi$), will typically lead to faster convergence rates for the OR-IVs. For a more detailed discussion of approximation bias conditions with sparse specifications, see belloni2012sparse. For the remainder of this section let us drop the dependence of random elements on $\ell$ to simplify our notation.

We next impose a sparse eigenvalue condition, following the Lasso literature bickellasso,bellonicherno2013, buhlmann2011statistics:

Notice that the objective function in (ref) depends on a sample counterpart of $F_j$, $\hat{F}_j$, and thus we hypothesize a convergence rate for it.

Finally, following chernozhukov2022automatic, for simplicity of the analysis, we allow the Lasso regularization parameter $\lambda_n$ to shrink slightly slower than $\varepsilon_n$. In the discussion to follow, we argue that choosing a $\lambda_n$ proportional to $n^{-1/4}$ generally suffice. We also restrict the growth of $r$ to be slower than some power of $n$.

Now, we can estate one of the main results of the paper.

Notice that the rate depends on $\xi$, which controls the degree of approximate sparsity in $f^{*}$. The larger $\xi$, i.e., smaller $s$, the faster the rate. Additionally, the rate depends on $\varepsilon_n$, which dominates the rates of $\hat{G}_j$ and $\hat{F}_j$. The better these rates, the faster the rate for $\hat{\kappa}$.

bakhitov2022automatic, in an endogenous setting, finds a slower rate of convergence for his estimator of the Riesz representer than that in chernozhukov2022automatic; see Theorem 1 in bakhitov2022automatic and the author's comments. bakhitov2022automatic proposes a penalized GMM estimator, which can be seen as an estimator from a Lasso extension of the standard GMM program. He employs dictionaries in both $Z$ and $V$, where $V \not \subseteq Z$. Differently, in our case, as the operator $S_{\theta_0,\eta_0}S^{*}_{\theta_0,\eta_0}$ maps elements in $L^2\left(Z\right)$ into $L^2(Z)$, we only use dictionaries in the exogenous variables of the model and the regressors $M_j$ are functions of $Z$, even when $\eta_0$ is a function of variables that are not part of the conditioning ones. This enables us to use similar theoretical devices as in chernozhukov2022automatic, despite we do allow for endogeneity. In other words, our procedure does not pay a price for dealing with situations with endogeneity. We regard this as a theoretical advantage of the way we have formulated our problem.

The rate for the OR-IVs that we obtained in Section (ref) relies crucially on Assumptions (ref) and (ref). We next provide a set of sufficient conditions for these. Additionally, they will conveniently establish a direct relationship between convergence rates of the estimators of all the unknown objects in the estimated regressors $\hat{M}_{j\ell k}$ and the rates in Assumptions (ref) and (ref). Let us define

We next assume

Assumption (ref) i) is a standard high-level condition that will be satisfied by several estimators. The key condition is ii). This assumption requires Lipschitz continuity of $M_{jk}$ in their unknown components. These are $\theta_0$, $\eta_0$, and the conditional expectations given $V$ and $Z_j$. Similar conditions to ii) are not unusual in the semi-parametric literature while dealing with non-linear functions; see, e.g., Assumption 3.9 in ai2003efficient and Assumption 12 in chernozhukov2022automatic. This condition is convenient as we can obtain a connection between rates of convergence for their estimators and the rate in Assumptions (ref) and (ref). When machine learning approaches are employed to compute estimators for $\eta_0$ (when $\eta_0$ is, in turn, a vector of conditional expectations), $\alpha_{0j}\left(V,\theta,\eta\right)$, and $\sigma_{0j}\left(Z_j,\alpha_{j^{\prime}}\left(V,\bar{\theta},\bar{\eta}\right), \theta,\eta\right)$, their rates are known (typically in mean square norm). Notice that, the construction of an estimator of $\theta_0$ that appears in $M_{jk}(Z_j)$ will be typically based on non-orthogonal moments, depending on machining learning estimators of $\eta_0$. Then, under some smoothness conditions (e.g., continuity of the non-orthogonal moments in $\eta$ around $\eta_0$) we can obtain that $\left|\left|\hat{\theta}-\theta_0\right|\right| = O_p\left(\left|\left|\hat{\eta}-\eta_0\right|\right|_{\Xi}\right)$. Hence we will be able to obtain bounds in probability for all the terms above. Assumption (ref) implies then, for any $k$,

Therefore, $\varepsilon^M_{jn} = n^{-d_j}$ in Assumption (ref). By a boundness condition, we can also determine a specific rate for $\hat{F}_{j\ell}$.

Estimation of the Parameter of Interest in a Two-Step Setting

Up to now, we have allowed the functions $m_j$'s to depend on $\eta_0$ and $\theta_0$ in a general fashion. As we have shown, our construction of debiased moments can handle it. Yet, we are interested in the common situation where the researcher works within a two-step setting, in which there are functions $m_j$'s that depend on one component of $\eta_{0}$ only. That is, $m_j\left(y,\theta,\eta\right) \equiv m_j\left(y,\eta_j\right)$ for some $j$. CMRs based on those residual functions can be used to obtain an estimator of $\eta_0$, as it occurs in our example.

We focus on this case as many relevant scenarios in applied work present this feature chen2016methods. Furthermore, a two-step setting might be computationally attractive in practice while estimating $\theta_0$, as pointed out by ackerberg2014asymptotic, as it avoids a simultaneous search over infinite and finite dimensional spaces in settings where $\eta_0$ and $\theta_0$ are jointly estimated by minimizing certain objective function. Our two-step setting can be written as

where $\mathcal{J}_1$ and $\mathcal{J}_1$ are disjoint sets of indexes and $J = |\mathcal{J}_1| + |\mathcal{J}_2|$.

Recall from our previous discussion that for a given instrument $f \in L^2(Z)$, we obtain an OR-IV $\kappa_{0}(Z)$---a vector with $J$ entries. Then, for different choices of $f$'s, say $q$ of them, we can construct $J$ vectors $\bm{\kappa_{0j}(Z_j)}$, of dimension $q$, which collects the $j-$component of the $q$ OR-IVs so constructed. We use the bold notation to emphasize that $\bm{\kappa_{0j}(Z_j)} = \bm{f_j(Z_j)} - \bm{f^{*}_j(Z_j)}$ is a $q-$vector. For the remainder of this paper, let us re-define (ref) such that $$ \psi\left(W, \theta, \eta, \bm{\kappa}\right) = \sum^J_{j=1}m_j\left(Y_i,\theta,\eta\right)\bm{\kappa_j(Z_j)}, $$ where we should notice that $\psi$ is a now a $q-$vector of LR functions. Even though we keep the notation that indicates that each $m_j$ depends on both parameters $\theta$ and $\eta$, $m_j\left(y,\theta,\eta\right) \equiv m_j\left(y,\eta_j\right)$ whenever $j \in \mathcal{J}_1$. Let $\hat{\eta}_\ell$ be an estimator of $\eta_0$, using observations in $I^c_{\ell}$. Let $$ \hat{\psi}\left(\theta\right) = \frac{1}{n} \sum^L_{\ell = 1} \sum_{i \in I_\ell} \psi\left(W_i,\theta,\hat{\eta}_\ell, \hat{\bm{\kappa}}_\ell\right). $$ Our proposed estimator $\hat{\theta}$ is defined as the solution to the GMM program

where $\hat{\Lambda}$ is a positive semi-definite symmetric weighting matrix. A choice that asymptotically minimizes the asymptotic variance is $\hat{\Lambda} = \hat{\Psi}^{-1}$, where $$ \hat{\Psi} = \frac{1}{n} \sum^L_{\ell = 1} \sum_{i \in I_\ell} \hat{\psi}_{i\ell} \hat{\psi}^{\prime}_{i\ell}, \;\;\;\; \hat{\psi}_{i\ell} \equiv \psi\left(W_i,\tilde{\theta}_\ell,\hat{\eta}_\ell, \hat{\bm{\kappa}}_\ell\right), $$ and $\tilde{\theta}_\ell$ is a possibly non-LR estimator of $\theta_0$, based on observations in $I^c_{\ell}$. The estimated $\hat{\psi}_{i\ell}$ directly accounts for estimating $\eta_0$ and $\kappa_0$ in a previous stage. We refer to (ref) as our Debiased GMM Estimator (DGMM). Let us summarize our estimation procedure with the following steps:

Step 1: For each subsample $\ell = 1,\cdots,L$ , compute estimates $\hat{\eta}_\ell$ and $\hat{\kappa}_{\ell}$, using observations not in $I_\ell$.

Step 2: Obtain our estimator $\hat{\theta}$ by means of (ref). The estimator of the asymptotic variance, which accounts for the estimation of the previous objects, takes the “sandwich" form

Figure (ref) provides a diagram that summarizes all the steps and important objects to conduct inference on structural parameters, as proposed by this paper.

Before concluding this section, we would like to remark that $\eta_0$ belongs to the infinite-dimensional space $\Xi$. As such, any suitable machine learning procedure can be used to obtain an estimator $\hat{\eta}$ by exploiting the CMRs that depend on each $\eta_{0j}$. To derive the asymptotic properties of $\hat{\theta}$, some restrictions will need to be imposed on the convergence rate of $\hat{\eta}$. Hence, it is substantial to know the behavior of such a rate asymptotically. A leading case in the semiparametric literature is when $\eta_0$ collects conditional expectations, given some $Z_j^{'}$s (as in our example). In this context, the mean-square convergence rate of $\hat{\eta}$ has been obtained for a variety of machine learners (see references in Section (ref)). We also consider the more general case where $\eta_0$ is a function of variables different from the conditioning ones, allowing for endogeneity. In the endogenous setting, identification of $\eta_0$ is trickier and typically involves an ill-posed problem. In particular, identification of this parameter might require the so-called “completeness condition" newey2003instrumental. The ill-posedness results in slower rates for the mean square norm of estimators of $\eta_0$ than in the exogenous situation. However, if rates are obtained for the so-called projected mean square norm, for example, they typically improved bakhitov2022automatic, gold2020inference, singh2019kernel. The vector $\eta_0$ might include densities as well. In any case, we will assume that a suitable norm $\left|\left|\cdot\right|\right|_{\Xi}$ is defined in $\Xi$, and as we will show below, we will allow the rate for $\hat{\eta}$, to be faster than $n^{-1/4}$ and slower than $n^{-1/2}$, which is arguably plausible for machine learning first steps.

Asymptotic Properties of DGMM

We start by imposing some regularity conditions:\footnote{Recall that $\left|\left|\cdot\right|\right|$ denotes Euclidean norm.}

Assumptions (ref) i) and ii) are mean square convergence conditions for $\hat{\eta}$, while iii) is a convergence condition for $\bm{\hat{\kappa}}$. Sufficient conditions for i) and ii) might involve continuity of $m_j\left(y, \theta_0,\cdot\right)$, for instance. Moreover, note that iii) is implied, e.g., by $\left|m_j\left(Y,\theta_0,\eta_{0j}\right)\right|^2$ being bounded a.s., for all $j$.

We need to assure that the estimators of OR-IVs and $\eta_0$ are well-defined, for this, we impose

Also, we hypothesize a rate of convergence for the estimator $\hat{\eta}$, which is required to be faster than $n^{-1/4}$ only, and the same requirement applies for our estimators of OR-IVs, as the following assumption states:

Lemma (ref) and Theorem (ref) imply that the rate of the OR-IVs is as close as desired to $n^{-2\xi d/(2\xi + 1)}$. Hence, Assumption (ref) is sufficient to assure that $\sqrt{n}\left|\left|\hat{\kappa}- \kappa_0\right|\right|_{L^2(Z)}\left|\left|\hat{\eta} - \eta_0\right|\right|_{\Xi} \rightarrow 0$, which is an essential condition for $\sqrt{n}-$consistency of $\hat{\theta}$.

Unsurprisingly, this depends on the rate of the estimated regressors in our Lasso program. This contrasts with earlier related works where the conditions for $\sqrt{n}$-consistency of $\hat{\theta}$ depends on the rate of consistency of $\hat{\eta}$ only and thus $d = d_\eta$ (see, e.g., Theorem 9 in chernozhukov2022automatic). The reason for this is that the regressors in our Lasso program, captured by $\hat{M}_{j\ell k}$, depends on unknown objects that must be estimated. In contrast, the design matrix in previous works, which have also based estimation on a Lasso-type program, depend on regressors that are known. In particular, they do not depend on $\theta_0$, $\eta_0$, or unknown conditional expectations. As a result, our proofs need to account for the presence of such unknown objects and the possibility that $d \leq d_\eta$ in obtaining a rate for OR-IVs. The practical implication of this result is that to achieve $\sqrt{n}$-consistency, in our context, we might need larger sparseness (i.e., larger $\xi$) than required by these previous papers. Nonetheless, we still allow for $d_\eta$ and $d$ being smaller than $1/2$ simultaneously. It is important to recall that Theorem (ref) requires that $\lambda_n$, the regularization parameter, to shrink slower than $\varepsilon_n$. Since $d > 1/4$, choosing a $\lambda_n$ proportional to $n^{-1/4}$ will generally suffice to achieve this requirement, as we specify in our algorithm of Section (ref).\footnote{Under our assumptions above, we let $\varepsilon_n = n^cn^{-d}$, for any $c>0$. As $c$ is sufficiently close to zero, the rate of $\varepsilon_n$ goes to zero at a rate sufficiently close to $n^{-d}$, which is faster than $n^{-1/4}$. On the other hand, $\lambda_n$ goes to zero at a rate no faster than $n^{-1/4}$.}

Next, let us define the sum of product of residuals $$ \hat{\Delta}_\ell(w) := \sum^J_{j=1}\left(m_j\left(y,\theta_0,\hat{\eta}_{\ell}\right) - m_j\left(y,\theta_0,\eta_{0}\right)\right)\left(\bm{\hat{\kappa}_{j\ell}(z_j)} - \bm{\kappa_{0j}(z_j)}\right). $$ Assumptions (ref), (ref), (ref), and Fréchet differentiability of the CMRs imply

Expression (ref) is a $\sqrt{n}-$convergence result for object $\hat{\Delta}_\ell$, which will be key to derive the asymptotic properties of $\hat{\theta}$. Additionally, let $$ \overline{\psi}\left(\theta,\eta,\bm{\kappa}\right) = \mathbb{E}\left[\psi\left(W,\theta,\eta,\bm{\kappa}\right)\right]. $$

Then, it can be shown that since $\psi$ leads to a debiased moment, there exists a $C >0$ such that $$ \left|\left|\overline{\psi}\left(\theta_0,\eta,\bm{\kappa_0}\right)\right|\right| \leq C \left|\left|\eta - \eta_0 \right| \right|_{\Xi}^2. $$ All the previous conditions yield the most important result of this section:

Interestingly, cross-fitting enables to show (ref) in a simple manner, without the need to impose the so-called Donsker conditions for $\eta_0$, as discussed by chernozhukov2018double and CEINR. Avoiding Donsker conditions is important as it is unknown if machine learners satisfy them. Furthermore, they could be restrictive in high-dimensional contexts.

The next condition is relevant to show consistency of $\hat{\Psi}$, the sample second moments of the debiased functions, which appears in the asymptotic variance of our GMM estimator.

Observe that the previous assumption is trivially satisfied for those $m_j$ depending on $\eta_j$ only, with probability approaching one. Finally, as in any GMM setting, we need conditions for convergence of the Jacobian: $\frac{\partial}{\partial \theta} \hat{\psi}(\bar{\theta}) \overset{p}{\to} \Upsilon = \mathbb{E}\left[\frac{\partial}{\partial \theta} \psi\left(W,\theta_0,\eta_0,\bm{\kappa_0}\right)\right]$ for any $\bar{\theta} \overset{p}{\to} \theta_0$. Let $\bm{\kappa_0(Z)}$ be a $q \times J$ matrix containing all OR-IVs $\bm{\kappa_{0j}(Z_j)}$, $j = 1, \cdots, J$. A similar definition applies to $\bm{\hat{\kappa}_{\ell}(Z)}$.

Given the previous assumptions and findings, the following result shows the asymptotic normality of $\hat{\theta}$:

Theorem (ref) implies that confidence intervals for $\hat{\theta}$ can be obtained straightforwardly in a standard way, using $\hat{\Sigma}$, and the usual quantiles of the standard normal distribution. This holds despite the convergence rates of nuisance parameters being slower than $\sqrt{n}$. Note that Theorem (ref) relies on the consistency of $\hat{\theta}$. We provide sufficient conditions for this in Section (ref) of the Supplementary Appendix.

Monte Carlo

We have conducted several Monte Carlo experiments based on our example, the production function setting, by simulating data from a panel of $n$ firms observed over $T$ periods. Only for this section, let us distinguish between variables in logs and in levels: uppercase letters denote variables in logs (as before), while lowercase letters stand for variables in levels.

Our baseline data generating process (DGP 1) is inspired by that in ackerberg2014asymptotic, Section 4.3. In our experiments, firms are observed over three periods, i.e., $T=3$. We consider a Cobb-Douglas production function in logs with only capital input: $$ Y_{it} = \theta_{01} + \theta_{0k}K_{it} + \omega_{it} + \epsilon_{it}, $$ where $\theta_{01} = 0$ and $\theta_{0k} = 1$. The law of motion of capital is given by $$ k_{it} = \left(1-\delta\right)k_{i,t-1} + \mu_{it}i_{i,t-1}, $$ where $1-\delta = 0.9$, $\mu_{it}$ is a lognormal shock to the capital accumulation process (with the underlying normal distribution having a mean and standard deviation of 1), and $i_{it}$ is the firm's investment decision. The log of investment is determined by

where $\gamma_0 = 0$, $\gamma_1 = -0.7$, and $\gamma_2 = 5$. We consider a large value for $\gamma_2$ to exacerbate the endogeneity bias due to the correlation between inputs and the anticipated productivity. Note that (ref) implies that $\omega_{it} = \omega_{t}\left(I_{it},K_{it}\right)$. We have specified an admittedly ad-hoc process for $I_{it}$ for two reasons. First, this avoids solving potentially complex firms' dynamic programs. Second, this process gives the link between productivity and the proxy variable ($I_{it}$ in this case), which is the source of the non-parametric component of the model. We intentionally model it flexibly and nonlinearly to justify the use of sophisticated machine learning tools.

One essential assumption of the proxy variable approach is the correct specification of the link between productivity and the proxy variable in the sense that $\omega_{it}$ is assumed to be a function of $K_{it}$ and $I_{it}$ only; see the discussion by ackerberg2015identification. In other words, the proxy variable approach rules out the possibility that other unobserved variables or shocks influence $\omega_{it}$. If anticipated productivity is governed by other factors, this is problematic. Controlling only for $I_{it}$ and $K_{it}$ in the estimation of $\eta_{0t}$ could lead to substantial bias in the first stage, resulting in poor estimation of the parameter of interest, both in terms of bias and coverage. As this cannot be tested in practice, if this situation takes place, one should prefer a less sensitive estimator to the first stage, i.e., a debiased estimator, rather than one that is not. This is what we aim to assess with our Monte Carlo simulation. We have considered two additional data generating processes where we add a normal shock to Equation (ref), with mean zero and standard deviations of 0.5 (DGP 2), and 0.7 (DGP 3).\footnote{Note that in this situation both estimators are not necessarily consistent, and thus it is a deeply disadvantageous scenario for both approaches.}

In all cases, productivity is assumed to follow a normal AR(1) process with persistence parameter $\theta_{0\omega} = 0.7$. The variance of the innovation term in this process is specified such that the standard deviation of $\omega_{it}$ is $\sigma_\omega = 0.1$. The unanticipated productivity or measurement error in output is normal and iid over firms and time. To ensure that our DGPs do not depend on the starting values of the variables, we focus on data drawn from the steady-state distribution implied by the model. Specifically, we have simulated the data over one hundred periods and kept the last three.\footnote{The same idea has been followed by ackerberg2014asymptotic and ackerberg2015identification.}

Let us recall that in this case,

where $F$ is Cobb-Douglas. Furthermore, an LR moment is simply

where $Z_t = \left(I_t,K_t\right)$.

We now comment on our estimation strategy. Recall that our algorithm begins with the choice of a vector of functions depending on the conditioning variables in the CMRs (ref)-(ref). Since this example has four CMRs, a initial vector function must have four entries. A given $f(Z)$ will allow us to obtain one debiased moment. In our simulations, we automatically construct four of these moments, and thus we specify four vectors of functions:

Each of these vectors is converted into an OR-IV to be used in estimation. For simplicity, in our experiments, the bases coincide and they are equal to $\gamma_1\left(I_{t}, K_{t}\right)$, which is obtained from the tensor product of univariate exponential bases in $I_{it}$ and $K_{it}$, that is

where $\tilde{\gamma}_s$ denotes the $s-$term of the sequence of exponential bases.\footnote{A one-dimensional exponential basis is of the form $\tilde{\gamma}_k(V) = \exp\left(\alpha_k V\right)$, each with different rate parameter $\alpha_k$. For a discussion of the use of bases with R, we refer the reader to functionaldataR.} Hence, $r_1$, the dimension of $\gamma_1$, is $r_1 = n/25$. Recall that we have to construct the regressors:

To reduce the dimension of the setting, we have imposed a restriction among the $\beta_j$'s. We let $\beta_1 = \beta_j$ for all $j$. As far as our original algorithm is concerned, this restriction can be simply imposed by noticing that since each entry in $\beta_1$ is a common factor across the entries of $\hat{M}_j$, the effective regressors are

and then we simply need to find the solution $\beta \in \mathbb{R}^{r_1}$ using our algorithm above, accordingly accommodated. This is computational appealing as the dimension of our matrices of regressors does not depend on the number of conditioning moments. The low-dimensional dictionary $\gamma^{\text{low}}(Z)$ is based on the first 5 terms of $\gamma_1$. The first stage $\eta_{0t}$ is estimated by Boosting using $L=5$ folds.\footnote{We use the R-function gbmt for boosting. Default options were kept. In particular, weights are equal to 1, the model offset is a vector of zeros, the number of trees is 2000, interaction depth is 3, the minimum number of observations in a node is 10, shrinkage is set at 0.001, bag fraction is 0.5, and the observations in the first half of the database are for training the model while the remainder are used for computing out-of-sample estimates of the loss function. Predictions were conducted using the first 500 iterations of the boosting sequence.} This choice follows a recommendation by chernozhukov2018double (p. C24). We set $\hat{\Lambda}$ as the identity matrix. All results are based on 500 Monte Carlo repetitions.

We remark a last point regarding our numerical computation. Notice that this model has three unknowns: $\theta_{01}$, $\theta_{0k}$, and $\theta_{0\omega}$. While nothing in our theory prevents us from estimating all of them simultaneously, we only focus on $\theta_{0k}$ for two important reasons. First, in the production function literature, researchers are primarily interested in marginal product estimates. Since our goal in the next section is to conduct an empirical exercise, we wanted our Monte Carlo exercise to serve as a preliminary step for applying our algorithm in a more real-world setting. Second, trying to estimate the entire vector of unknown parameters might not be advisable, as it involves a three dimensional nonlinear minimization problem. Our DGPs do not appear to generate strong variation in the instruments (non-OR-IVs and OR-IVs) to ensure a sufficiently convex objective function, potentially leading to very large standard errors, which seems to be potentially serious for $\theta_{0\omega}$, the AR(1) coefficient.\footnote{This issue has been previously noted by ackerberg2015identification (p. 2448).} Conveniently, it is possible to write the GMM objective function in terms of $\theta_{0k}$ only, as the AR(1) process imposes a restriction on the remaining parameters, which must satisfy certain conditions given $\theta_{0k}$; for details see Section (ref). At the same time, this might make our results more reliable.

To assess the performance of DGMM, we have computed a naive plug-in estimator (PI) as our benchmark. This estimator is based on the same first stage as DGMM but uses as instruments natural choices, which are non-OR-IVs, resulting in a non-robust estimator. It is obtained by GMM using the CMRs (ref) and (ref), with instruments: $ K_{t-1}, I_{t-1}, K^2_{t-1}$, and $I^2_{t-1}$. To compute its standard errors, we now have to explicitly take into account the estimation effect from the first stage, unlike DGMM, which automatically correct for this; see newey1994asymptotic.

Table (ref) reports our results. First, let us focus on DGP1, the one that has a correctly specified link between $\omega_{it}$ and $I_{it}$, given $K_{it}$. We observe a clear bias-variance trade-off. Although the bias does not seem to be practically large for the naive estimator, as the sample size increases, it tends to be larger relative to DGMM. Our estimator, across the Monte Carlo repetitions, shows systematically greater standard errors. This is not surprising as it involves the estimation of the OR-IVs, introducing additional sources of uncertainty. However, this greater imprecision does not seem large. In fact, only when $n=250$ and $n=500$, DGMM reports a higher RMSE. In terms of coverage, while our estimator has coverage figures close to the expected nominal level, the naive estimator always undercover, which seems to deteriorate as the sample size increases. All these observations are coherent with our theory.

When we introduce shocks to investment, which increase the bias in the first-stage estimation, we observe that the bias in the plug-in estimator increases relative to DGP1 with no shock, as the first-stage bias becomes greater (captured by larger variability in the investment shock). In contrast, DGMM seems to be less sensitive to the size of the bias in estimating $\eta_{0t}$, which is expected as it uses OR-IVs, resulting in systematically smaller bias compared to the plug-in estimator. Again, although DGMM shows larger standard errors, this price does not seem to be large in this numerical exercise, which can be noticed by the fact that DCRMs can deliver a lower RMSE than the plug-in estimator. Possibly more concerning, the naive approach results in poor coverage, with a convergence level as low as 60$\%$. On the other hand, despite the misspespecification in the first stage, DGMM still provides satisfactory coverage. The implication of this exercise is that, even though DGMM may have higher standard errors, it might be recommendable to use a debiased estimator rather than a non-debiased one, as it is less affected by first-stage bias.

Empirical Application

We now evaluate the performance of our estimator using real data to estimate production functions at the firm level. Our goal is to assess whether the use of machine learning tools makes a significant difference in practice, specifically in delivering different production function estimates, relative to other approaches.

cha2023inference use an adjusted version of the Double Debiased Machine Learning (DML) estimator of chernozhukov2018double within the context of production functions estimation. They exploit the orthogonal score function from robinson1988root. Therefore, the authors conduct inference on only one parameter at a time and treat the parameters of the rest of inputs, including capital, as nuisance ones. Moreover, they focus on a particular strategy to estimate nuisance components. As recognized by the authors (p. 24), their strategy might not lead to identification of the production function parameters as it might suffer from the functional dependence problems discussed by ackerberg2015identification. In contrast, our two-step approach follows ackerberg2014asymptotic and ackerberg2015identification, avoiding functional dependence issues. In addition, we apply debiased machine learning to production functions using off-the-shelf machine learning routines, e.g., random forest, to conduct inference on the parameters of all the inputs in the production function simultaneously. We are not aware of previous works doing this.

Generally speaking, we consider a model similar to our example. We examine a production function in two inputs: $L_{it}$ and $K_{it}$, labor and capital, respectively. The production function $F\left(\cdot,\theta_{0p}\right)$ is Cobb-Douglas and thus Equation (ref) can be written as

As a proxy variable we now use intermediate inputs $E_{it}$, following levinsohn2003estimating.\footnote{The monotonicity assumption, required by the proxy variable approach, is more likely to hold with intermediate inputs than investment. Several firms in our data reports no investment for several periods, which might indicate failure of monotonicity. } Hence, $E_{it} = E_{t}\left(\omega_{it}, L_{it}, K_{it}\right)$ is the variable that allows us to control for unobserved productivity. Note that $E_{it}$ does not appear in (ref), meaning that $Y_{it}$ is “value-added" rather than output. One possible interpretation for this is that output follows a Leontief specification in $E_{it}$; for other possible interpretations see bruno1978analysis and diewert78. As showed by bond2005adjustment for the Cobb-Douglas specification, and gandhi2011identification for a non-parametric production function, under the assumptions of the proxy variable approach, outlined by olleypakes96, levinsohn2003estimating, and ackerberg2015identification, production functions with gross output are not identified unless additional assumptions (e.g., observation on prices that have to vary sufficiently across periods) are imposed. As emphasized by most of the previous authors, the circumstances under which the coefficients in (ref) can be interpreted in a meaningful way are not trivial. Thus, we simply interpret these coefficients as a measure of the “effect" that inputs have on value-added output, not necessarily as marginal productivity as in the usual sense.\footnote{Still, for simplicity, we refer to these coefficients as marginal productivities/products.}

Our starting CMRs are

We follow ackerberg2015identification and let $L_{it}$ be a state variable in the firm's dynamic problem, which can be chosen at $t$, $t-1$, or at some intermediate point (potentially with dynamic implications). Similarly as above, an LR moment can be written as follows

where $\kappa_{0} = \left(\kappa_{012},\kappa_{013}, \cdots,\kappa_{01T}, \kappa_{022}, \kappa_{023} \cdots, \kappa_{02T} \right)$ is a vector of OR-IVs obtained from a vector of functions of starting instruments, to be specified below. These OR-IVs are computed following the algorithm outlined in Section (ref). This setting has $2(T-1)$ CMRs and thus $2(T-1)$ design matrices $\bm{\hat{M}_{j}}$. The $i-$ row of each of these matrices can be written as

where $t = 1,\cdots, T-1$. In our results below we simply construct the same basis for all $t$ (so we can drop the time dependence when referring to the bases) and also we let $\gamma_{1} = \gamma_{2}$, which in our results below is obtained from the tensor product of univariate exponential basis of dimension 5. Specifically, let $\tilde{\gamma}$ be a vector collecting the terms of the sequence of exponential basis, we have

where $\tilde{\gamma}_s$ denotes the $s-$term of the sequence of exponential of bases. Note that $r_1 = 125$. Our starting vectors of functions are, for $2 \leq t \leq T$,

In Section (ref) in the Supplementary Appendix, we give specific conditions for the validity of our inference procedure, as explained in Section (ref), in this setting. This result ensures $\sqrt{n}-$consistency and asymptotic normality of our debiased estimator for the entire vector $\theta_0 = \left(\theta_{01}, \theta_{0l}, \theta_{0k}, \theta_{0\omega}\right)$. However, as we have pointed out in the previous section, in practice researchers are typically interested in inference for $\theta_{0l}$ and $\theta_{0k}$ only. Moreover, estimating all parameters simultaneously might lead to significantly larger standard errors for these key coefficients. Therefore, we only focus on inference for $\theta_{0l}$ and $\theta_{0k}$. We have re-written our problem accordingly, in the same way as we did with our Monte Carlo exercise. That is, we express the GMM objective function such that it depends on $\theta_{0l}$ and $\theta_{0k}$ only. This is advantageous as we only have to proceed by minimizing the GMM objective function in two parameters instead of four, presumably improving the reliability of our results.

We use data from the Instituto Nacional de Estadistica de Chile, which conducted a census of Chilean manufacturing plants with at least ten employees during the period 1979-1987. We are specially interested in this data as it has been extensively employed in the estimation of production functions and productivity at the firm level levinsohn2003estimating, gandhi2020identification, ackerberg2015identification, pavcnik2002trade, alvarez2005exporting.\footnote{We are indebted to David A. Rivers and Salvador Navarro for sharing the data with us and for their assistance with the variables construction details.} The variables in our analysis have been constructed following gandhi2020identification, who use the conventions in greenstreet2007exploiting. In particular, we construct real gross output as deflated revenues. Intermediate inputs are obtained as the sum of the firm's expenditure on electricity, fuels, raw materials, and services. The difference between these two yields value-added output. Capital is formed by the perpetual inventory method. This implies that in each period investment in $t$ is combined with deflated capital from period $t-1$ to produce capital. Labor is aggregated from white-collar and blue-collar workers, with the former weighted by the ratio of the (sector) average white-collar wage to the (sector) average blue-collar wage. Deflators are obtained from bergoeing2003idiosyncratic and all variables are expressed in real 1980 Chilean thousands pesos.\footnote{See Saction D in the Supplementary Appendix of greenstreet2007exploiting for more details; see also liu1991.} Our sample consists of firms from the five largest three-digit ISIC manufacturing industries in Chile: food products (311), textiles (321), apparel (322), wood products (331), and fabricated metal products (381). To maximize statistical power, we pool all firms into a single sample, which has 1,069 establishments.\footnote{This sample size is the one that remains after we drop from the original data those firms with missing values for at least one of the variables in our analysis.} We focus on three subperiods: 1979-1981, 1982-1984, and 1985-1987. These coincide with macroeconomic cycles of the Chilean economy in the data, as identified by levinsohn2003estimating. We will conduct our analysis for each of these time periods separately below.

To reduce the dimensionality of the setting, we have assumed that the coefficients across CMRs are constant as in our Monte Carlo experiments. That is, $\beta_{jk} = \beta_k$ for all $1 \leq j \leq J$ and $k = 1,\cdots, r$, where $r = 125$. The low-dimensional dictionary $\gamma^{\text{low}}(Z)$ is based on the first $r_1/15$ terms of $\gamma_1$. While we could have computed different DGMM estimators using different off-the-shelf machine learning tools to estimate $\eta_{0t}\left(E_t,L_t,E_t\right)$, we focus on first stages computed by random forest, which provided us with the most precise point estimates.\footnote{We have implemented this by the R-function ranger (default options).} We have let $c_1=2$ in our expression for $\lambda_n$, which has also helped us reduce standard errors with this particular data.

We compare the results from DGMM with other two estimators that do not rely on machine learning tools. First, an OLS estimator, obtained by regressing $Y_t$ on $L_t$, $K_t$, and a constant.\footnote{As we discussed, this estimator is inconsistent but gives a simple framework.} Second, the estimator developed by ackerberg2015identification, which is currently implemented by the R function prodestACF.\footnote{See the package at https://github.com/GabrieleRovigatti/prodest/tree/master/prodest.} There are three key differences between this estimator and ours. Firstly, the estimator by ackerberg2015identification (ACF) treats $\eta_{0t}$ non-parametrically and is estimated by a sieve specification chen2007large. In particular, the off-the-shelf routine regresses $Y_t$ on a polynomial of degree $d$ in $E_t$, $L_t$, and $K_t$. This strategy assumes then that all the terms appearing in the polynomial are important in the first stage, which might restrict the unknown complexity of $\eta_{0t}$ and affect the estimates. The default option is $d=2$, a common choice in applied work. We consider this degree and also a more flexible specification with $d=5$. Secondly, while ackerberg2015identification assume that $\omega_t$ follows a first-order Markov process, as we do, they treat it non-parametrically rather than assuming an AR(1) process, and approximate this with sieves. Lastly, standard errors are computed by bootstrap. In our estimation, the optimal GMM weighting matrix is used.

Table (ref) presents the estimation results for the labor and capital parameters in Equation (ref) across the four estimation routines we have implemented.\footnote{These estimation results should be interpreted with caution. Despite their widespread use in applications, identification strategies based on timing and information set assumptions—such as those employed in the proxy variable approach to production functions—may fail and remain not well understood; see kim2019robust and ackerberg2023under for recent discussions.} DGMM, across all time periods, suggests relatively large labor marginal products, in the range of 1.1 and 1.2 (std. err.= 0.38-0.50).\footnote{Note that these are not strictly marginal productivities as the dependent variable is not output but value-added. Yet, they are strongly correlated. We refer to them as marginal productivities to simplify the exposition. These results are intended primarily to illustrate first-stage biases in this dataset.} By contrast, ACF reports lower labor marginal products: between 0.9 and 0.7 (std. err.=0.2–0.3) when the first stage uses $d=2$, and between 0.9 and 0.65 (std. err.=0.2–0.5) when $d=5$. These figures are very similar to those obtained by the inconsistent OLS estimator.

Turning to capital, all methods report systematically lower marginal products. DGMM estimates values around 0.2–0.3, with, for instance, an estimate of 0.32 (std. err.=0.16) in 1982–1984. ACF, in both specifications, yields capital marginal products close to 0.3 in 1979–1981 and 1985–1987. However, in 1982–1984 ACF delivers implausible negative marginal products, likely reflecting severe bias driven by functional-form restrictions. Relative to DGMM, OLS appears downward biased for labor, while the capital coefficient is virtually identical.\footnote{In general, the sign of the bias depends on firms’ ability to adjust inputs after observing productivity.} Taken together, the DGMM results suggest that as firms became more productive, they reduced blue-collar employment while keeping their capital stock essentially unchanged. This aligns with the extensive literature documenting high unemployment in the 1980s and the adverse short-term employment effects of trade liberalization edwards1997trade,riveros86,paus1994economic. In contrast, the ACF results would imply that the transmission bias from ignoring endogeneity is negligible—a conclusion difficult to reconcile with the large structural economic reforms that Chile was undertaken during our period of analysis.

Both DGMM and ACF support constant returns to scale at the 5% level across all periods. For example, in 1985–1987 the sum of labor and capital marginal products is 1.03 under ACF ($d=2$), 1.12 under ACF ($d=5$), and 1.50 under DGMM, with standard errors of similar magnitude. DGMM also suggests that capital–labor intensity ranges from 0.15 (std. err.=0.06) in 1979–1981 to 0.29 (std. err.=0.13–0.17) in 1982–1987. OLS delivers larger ratios, from 0.35 (std. err.=0.02) in 1979–1981 to 0.43 (std. err.=0.03) in 1985–1987, broadly consistent with ACF in some subperiods. Therefore, conclusions differ depending on how we model the first stage.

The data suggests that estimates of labor marginal products are less precise than those of capital: standard errors for labor are typically larger across methods and years. While DGMM requires the additional estimation of OR-IVs, its standard errors are not necessarily larger than those of ACF. In fact, while DGMM typically reports the largest standard errors for $\hat{\theta}_l$, the two methods produce standard errors of comparable magnitude for $\hat{\theta}_k$.

An important observation of this table is that the approach a researcher undertakes for modeling the first stage should be chosen with caution. Failing to flexibly model the first-stage in finite samples might be misleading.\footnote{However, one should not generally observed differences asymptotically when nonparametric approaches are employed.} The point estimates that ACF yields are different from the results using more flexible machine learning tools. In addition, more restrictive specifications might lead to inaccurate results, with negative marginal products reported by ACF. Researchers should therefore evaluate how their results depend on the specific modeling choices for the first stage while maintaining the ability to perform standard inference. We believe that our procedure---which provides a tool to implement more general first-stage estimations---could help them in that regard.

Final Remarks

This paper has developed an approach to estimate OR-IVs, which are the key objects to construct LR/(Neyman-)orthogonal/debiased moments in general semiparametric models defined by a finite number of CMRs, with possible different conditioning variables and endogenous regressors. Our estimation strategy leverages the CMRs implied by the model in a general way, making it applicable to a wide variety of settings. Our approach will hopefully pave the way for the use of machine learning techniques in contexts where the construction of LR has remained unexplored, such as highly non-linear and complex settings with CMRs, prominent in structural economics.

We have demonstrated that OR-IVs can be obtained as the residuals of a Lasso-type program, which can be implemented straightforwardly. We have derived a convergence rate for these functions. Based on this construction, we have introduced an estimator of structural parameters, the DGMM, in a two-step GMM setting. This estimator has been shown to be $\sqrt{n}$-consistent and asymptotically normal, requiring only that nuisance parameters, including OR-IVs, be estimated at a rate faster than $n^{-1/4}$. Our Monte Carlo experimentation has shown that employing such an estimator in high-dimensional contexts offers advantages over naive procedures, particularly in terms of bias and coverage.

Building on our theoretical results and Monte Carlo experiments, we applied our methodology to estimate production functions and productivity at the firm level in our empirical application. A key implication of our analysis is that the choice of approach for estimating first stages can significantly impact empirically relevant conclusions, particularly regarding production function parameters. Consequently, it is appealing to have inference methods, such as the one presented in this paper, that enable researchers to adopt more flexible specifications in the first stage. This is important for two reasons. First, it provides a valuable tool for scholars who prefer to avoid imposing restrictive assumptions. Second, it offers a way to conduct robustness checks for those more confident in traditional (somewhat less flexible) approaches.

We acknowledge that our paper has limitations. We have assumed that our limited number of LR moments are sufficient for identifying $\theta_0$. Moreover, while assuming identification, we have not addressed how efficiency could be improved by selecting other orthogonal moments. These are crucial matters that should not be overlooked. Future work should focus on constructing debiased moments that ensure identification in models defined by a number of CMRs, using modern tools as discussed by muandet2020kernel and zhang2021instrument. Furthermore, as we have seen, several $\kappa_0$'s might exist. This raises the question of whether it is possible to characterize a suitable notion of “optimality" among these OR-IVs and propose an algorithm to select OR-IVs that improve efficiency. From an efficiency perspective, it is well-known that the first best is the optimal IV chamberlain1992efficiency, ai2003efficient. This special IV not only yields an estimator that achieves the efficient semiparametric bound, but it is also a valid OR-IV vandervaart98, newey1990semiparametric. Such a choice, nonetheless, is difficult to implement in general settings. It might be interesting to define a broader criterion for selecting second-best OR-IVs that improve efficiency compared to other OR-IVs while remaining computationally feasible. Additionally, a more general theory for the estimation of OR-IVs can be derived. We have performed our construction exclusively for the space of sparse functions. It might be promising to develop a general framework for different functional spaces $\mathcal{G}_n$, including Reproducing Kernel Hilbert Spaces and Neural Networks. Finally, our algorithm focuses on debiasing structural parameters $\theta_0$. It might be important to extend this approach to more general parameters that depend on both $\theta_0$ and $\eta_0$, that do not necessarily appear in the CMRs. Hopefully, these ideas will be addressed in subsequent works.

\phantomsection