Debiased Machine Learning for Unobserved Heterogeneity: High-Dimensional Panels and Measurement Error Models

Debiased Machine Learning for Unobserved Heterogeneity: High-Dimensional Panels and Measurement Error Models

Introduction

One of the most important insights of empirical economics is the pervasiveness of heterogeneity in observed economic behavior heckman2001micro. In particular, Unobserved Heterogeneity (UH) plays a central role in bridging economic modeling with empirical data, with significant implications for policy Browning_Carro_2007. Structural analysis often focuses on functionals of the nonparametric UH, which are key summary objects of interest in many empirical applications. These structural parameters may include average treatment effects, average marginal effects, variances, counterfactual outcomes, linear projections of UH on covariates, welfare measures, bounds, and other downstream functionals. Many of these parameters can be expressed as functionals of the unknown distribution of UH and other parameters of the model.

This paper introduces a systematic approach to conduct inferences for functionals of the UH based on Neyman-orthogonal moments neyman1959.\footnote{Throughout the paper, we use indistinguishably the terms Neyman-orthogonal, orthogonal, debiased, and Locally Robust (LR) moments.} Neyman-orthogonal moments for a target functional are robust to local deviations of the nuisance parameters of the target at the truth. The application of locally robust inference to contexts with nonparametric UH is naturally important for (at least) four reasons. First, the distribution of UH is typically not nonparametrically identified, e.g., in short panel data models with discrete outcomes. This implies that nuisance parameters of target functionals are often not identified in these settings, so obtaining moments that are locally robust to them is important. Second, even when the distribution of UH is point-identified and estimated, plug-in estimators may suffer from large regularization biases---so much so that they render inference on the target functional invalid. Third, nonparametric estimators of UH may have very slow (e.g., logarithmic) rates of convergence, invalidating critical assumptions in semiparametric estimation theory.\footnote{Plug-in estimators with these properties are available from a large literature on nonparametric estimation of the distribution of UH in arellano2012identifying, bonhomme2010generalized, carroll1988optimal, heckman1984method, horowitz1996semiparametric, kato2021robust, kiefer1956consistency, koenker2014convex, li1998nonparametric, and szekely2000identifiability, among many others.} Orthogonal moments are free from the distribution of the UH if an additional support condition holds, which overcomes the problem of the regularization bias and slow rate of convergence of nonparametric estimators of the distribution of the UH. Fourth, efficient moments are locally robust moments. It is thus important to provide inference methods in models with UH that are asymptotically valid (provable), robust to the UH and additional nuisance parameters, and informative (relevant) about the target parameter. We show that our Debiased Machine Learning (DML) for models with UH (DML-UH in short) is a systematic approach to achieving these goals.

As the main result of this paper, we provide necessary and sufficient conditions for the existence and relevance of orthogonal moments for general smooth functionals in models with nonparametric UH. The precise sense of smoothness of the functional is defined in Section (ref). Relevant Neyman-orthogonal moments, a concept we introduce in this paper, are informative about the target functional even when its Fisher information matrix is singular but not zero (e.g., in some partially identified cases). We thus characterize all relevant Neyman-orthogonal moments in models with nonparametric UH. Our necessary and sufficient conditions are constructive, providing new tools for identification, estimation, and inference on functionals of UH in a broad range of settings.

Based on the necessary and sufficient conditions for orthogonality, we propose a three-step algorithm for their construction. Specifically, we show that relevant orthogonal moments must solve a system of two functional equations. The first is an integral equation that is free of the distribution of UH when certain additional support conditions hold. This means that locally robust moments are globally robust to the distribution of UH. In some cases, this first equation admits a closed-form solution; more generally, it can be solved numerically or via functional approximation. Such support conditions and similar equations to our first one have notably appeared in the functional differencing literature, see the important work by bonhomme2012functional, and more recently honore2024moment. Yet, ours is the first paper to underscore the role played by support conditions to obtain global robustness. We show that similar arguments invoked by such literature apply more broadly, e.g., to semiparametric/nonparametric conditional likelihoods given the UH. Thus, the first equation characterizing orthogonal moments generalizes estimating equations appearing in the functional differencing literature.

More importantly, we show that this first equation alone does not give an orthogonal moment for general functionals. For example, existing moments for average marginal effects (e.g., chamberlain1992efficiency, aguirregabiria2024identification) that solve the first equation are not orthogonal in general. This lack of robustness may be particularly problematic when additional high-dimensional nuisance parameters are present, as in the high-dimensional random coefficient model we consider in this paper. Orthogonal moments must also satisfy a second functional equation to ensure orthogonality with respect to other nuisance parameters different from the UH. We develop a three-step algorithm that finds a solution to both equations simultaneously, and hence delivers an orthogonal moment.

Once a relevant orthogonal moment has been constructed using our method, we introduce a score-type test based on it. We derive its asymptotic properties and show how it can be used for valid inference. Importantly, our characterization and inference results allow for multivariate partially identified target parameters and the use of machine learning methods to estimate high-dimensional nuisance components, incorporating cross-fitting to mitigate overfitting bias, as in chernozhukov2018double.\footnote{While our characterization of local robustness allows for non-identified high-dimensional nuisance parameters, the implementation of our inference procedures requires the orthogonal moment to be identified. This situation is analogous to that considered in santos2011instrumental, escanciano2021optimal, and, more recently, in bennett2023source. As in these references, regularization of the nuisance parameters helps in identifying one element of the nuisance parameter's identified set.} We also present a number of novel results on identification and LR estimation.

We illustrate our three-step construction and inferences with three distinct applications involving nonparametric UH. These specific applications build on our general characterization of relevant orthogonal moments. The first example application develops DML inferences in a random coefficient panel data model with high-dimensional exogenous controls and a finite number of periods. Such specification may be well-motivated for several reasons, such as the availability of large data sets that contain a myriad of controls, concerns about strong functional form assumptions, or to make causal statements more credible.\footnote{In Section (ref) of the Supplementary Appendix, we build on borusyak2024revisiting and argue that our procedure complements the literature on event-study designs by allowing for flexible modeling of parallel trends assumptions. We are not aware of inference procedures---using LR moments---that allow for high-dimensional regressors in this context, which might be of independent interest.} We build on the important work by chamberlain1992efficiency and arellano2012identifying to generalize their results to a setting with high-dimensional covariates and robust inferences. We construct in Section (ref) novel orthogonal moments and develop DML-UH procedures for estimating linear combinations of common parameters, average marginal effects, and variances of the UH. We are currently developing an R package, dmlpanel.R, that implements all the procedures described in this section.\footnote{A preliminary version of this function can be obtained from https://github.com/argafacu/dmlpanel.}

A second application is to a canonical measurement error model--the Kotlarski model with a factor loading. The identification of this model has been recently investigated in lewbel2022kotlarski, who also provides a rich historical perspective. We develop in Section (ref) the first debiased inferences for non-central moments of the common factor in this setting. This model has been used, for example, in studies of income dynamics bonhomme2010generalized, horowitz1996semiparametric, the formation of cognitive and noncognitive skills in childhood cunha2010estimating, studies on the impact of grade retention fruehwirth2016timing, the intergenerational mobility literature nybom2016heterogeneous, haider2006life, panel data evdokimov2010identification,lewbel2022kotlarski, and the triangular two-equation system without instruments lewbel2024identification. Inference on this model is notoriously difficult, see kato2021robust, as the conditional likelihood of the data given the common factor contains two high-dimensional nuisance parameters (the idiosyncratic errors' distributions). Our orthogonal moments are free from these nuisance parameters.

The third example application that we work on belongs to the literature on the quality of institutions, such as hospitals, doctors, nurses, teachers, schools, neighborhoods, and firms kane2008estimating, fletcher2014estimating, yakusheva2014nurse, chetty2014measuring, angrist2016interpreting, angrist2017leveraging, finkelstein2017adjusting, doyle2019evaluating, kline2022systemic, chetty2018impacts, abowd1999high. These studies typically employ linear random coefficient models with multidimensional UH, drawing on “value-added” estimation methods originally developed in the teacher value-added literature koedel2015value, and focus on functionals of UH using plug-in procedures; see bonhomme2024estimating for discussion of this framework.

We consider the standard teacher value-added model with a nonparametric distribution of “value-added” or random effects, as in, e.g., gilraine2020new. We make two main contributions within this framework in Section (ref). First, we use the necessity of our conditions for orthogonality to show the nonexistence of orthogonal moments for prominent policy parameters in the teacher value-added literature. This result has important policy implications, as it shows the non-robustness of certain counterfactual measures that inform teacher hiring and retention decisions. Second, we develop robust inference methods for smooth functionals of “value-added”, leveraging Neyman-orthogonality. Notably, our approach to constructing orthogonal moments does not require estimating the distribution of UH, and therefore, it provides a robust alternative to commonly employed regularized plug-in methods.

There is nothing fundamentally special about the previous settings. Our results can be applied to other contexts involving some form of UH or latent variables, where Neyman-orthogonality principles for general smooth functionals (and possibly under more general model assumptions) are not yet known. Examples include, among many others, panel data models with censored, truncated, or missing data honore1992trimmed; nonparametric panel data with interactive fixed effects freyberger2018non; non-classical measurement error models with instrumental variables hu2008instrumental; mixed logit models rafi2024nonparametric; production function models with firm specific technologies duffy2004capital, MairesseGriliches1990, dobbelaere2013panel, structural models of firm investment and capital formation cooper2006nature; duration models heckman1984identifiability, heckman1984method, horowitz1999semiparametric, including multiple-spell duration models with individual heterogeneity in different parameters van2001duration, alvarez2016decomposing; auctions with nonparametric UH krasnokutskaya2011identification; linear and nonlinear random coefficient models, lewbel2017unobserved, masten2018random; individual-specific income profile/earning dynamic models lillard1979components, guvenen2009empirical, arellano2017earnings; and learning models bunting2024heterogeneity .

The most common approach to nonparametric inference in these applications is to use a nonparametric regularized maximum likelihood (sieve) approach followed by a plug-in principle. However, this approach has several important limitations relative to the methods proposed here. First, it requires fully nonparametric identification, which may entail strong assumptions. Second, it can be computationally intensive, as it needs to solve a high-dimensional nonlinear optimization problem. Finally, it can be very sensitive to the choice of the sieve basis, since estimating these nuisance parameters is an ill-posed inverse problem. Nevertheless, sieve methods could be combined with our orthogonal moments to provide DML inferences with desirable properties. We illustrate this point in Monte Carlo simulations, where we show the large efficiency gains from using DML-UH with efficient moments (which debias a plug-in moment). This application emphasizes the practical importance of relevant orthogonal moments.

An additional important setting where our results may be useful is the analysis of heterogeneous treatment effects under flexible specifications. In our empirical application, we study heterogeneity in the effect of smoking during pregnancy on child birth weight, an outcome that has received considerable attention in economics due to its strong correlation with later-life outcomes such as health, education, and labor market performance hack1995long, corman1995effects, corman1998effect, currie1999impact. Our approach extends abrevaya2006estimating by allowing for heterogeneous slopes, and builds on arellano2012identifying by incorporating high-dimensional controls and constructing novel moments that are locally robust to the covariance structure of the idiosyncratic errors. Our empirical findings are broadly consistent with those in arellano2012identifying and show, in particular, that their conclusions remain robust to flexible specifications of the control variables. Yet, our results suggest slightly more variability of the treatment effects of smoking across mothers than previously determined.

The rest of the paper is organized as follows. Section (ref) reviews the related literature. Section (ref) introduces the examples, target functionals, and a heuristic summary of the main results. Section (ref) establishes conditions for the existence and relevance of orthogonal moments in models with nonparametric UH. Section (ref) proposes a three-step method for constructing relevant debiased moments. Section (ref) introduces a test statistic for DML inference, and Section (ref) applies the general results to our examples. Section (ref) presents Monte Carlo simulations. Section (ref) reports an empirical application to the effect of smoking on birth weight. Finally, Section (ref) concludes. The Supplementary Appendix contains proofs, asymptotic results, regularity conditions, and additional discussion.

Literature Review

This paper contributes to two different strands of the literature---the literature on DML and the literature on models with UH or latent variables. We explain our contributions relative to these two literatures.

While the literature on DML has been growing rapidly in the last years, the work on orthogonal moments for models with UH is relatively scarce, and more so for general functionals. The most recent literature on DML has shown that inferences based on LR moments have several advantages over plug-in approaches. Most notably, LR moments are useful to reduce model selection and regularization biases when Machine Learning (ML) estimators are used as first steps chernozhukov2018double, they lead to optimal inferences lee2024locallyregularefficienttests, and they provide valid confidence intervals for structural parameters under more general conditions than plug-in methods chernozhukov2022locally.

In this paper, we contribute to the DML literature by studying functionals of the UH, a setting for which existing DML methods remain underdeveloped. In the context of linear panel data models with additive and scalar UH, our work complements recent DML inference methods proposed by clarke2025double, klosin2022estimating, and semenova2023inference. Like some of these papers, we build on the estimation framework of belloni2016inference. However, we depart from this literature in the specific application to high-dimensional panel data by allowing for multivariate UH and considering general functionals of UH (e.g., average marginal effects and variances).

Moreover, our approach is conceptually distinct from the existing DML literature, which requires nuisance and target parameters to be identified from a starting GMM moment condition, see, e.g., chernozhukov2022locally. Instead, we adopt a score operator approach to characterize Neyman-orthogonal moments, which complements the identification results of escanciano2022semiparametric and the operator-based identification and orthogonality results for causal models in navjeevan2023identification. In particular, our paper is the first to characterize all (relevant) Neyman-orthogonal moments in models with UH.

Our paper also relates to the important literature on robust inference in the presence of UH in panel data settings with a finite number of time periods or repeated measures, and an unrestricted distribution of individual effects given covariates (i.e., a fixed effects approach).\footnote{See the differencing methods proposed in, e.g., andersen1970asymptotic, bonhomme2012functional, chamberlain1984panel, honore1992trimmed, honoré2023momentconditionsdynamicpanel, and mundlak1978pooling. There is also an extensive literature on debiasing for large $T$ panel data, see arellano2007understanding for a survey, where $T$ denotes the number of periods. Our DML inferences are valid for fixed $T$.} In the linear random coefficient model with strictly exogenous covariates, chamberlain1992efficiency introduced a GMM estimator for the common parameters and average marginal effects that does not rely on the unknown distribution of UH. graham2012identification further examined the irregular identification and estimation of these parameters.\footnote{In this paper, we focus on regular identification for target parameters, in the general sense that the information matrix for them is non-zero but possibly of reduced rank; the irregular case—characterized here by an identically zero information matrix—is left for future research.}

Building on chamberlain1992efficiency, arellano2012identifying used covariance and additional independence restrictions on the idiosyncratic errors to identify variances or even the full distribution of UH. They proposed an estimator for the distribution of UH and reported plug-in estimates for functionals such as quantiles in their empirical application.

In this paper, we show that the moments commonly used to estimate average marginal effects, variances, and other functionals of the UH are generally not orthogonal. We derive novel locally robust (Neyman-orthogonal) moments for these parameters—when they exist—and provide valid asymptotic inference procedures based on them.

This paper also relates to the important literature on functional differencing (bonhomme2012functional). Under parametric assumptions on the conditional distribution of the dependent variables given the exogenous covariates and the UH—a setting we refer henceforth to as conditionally parametric—bonhomme2012functional developed a systematic approach, known as functional differencing, to construct moment conditions for common parameters that are free from the distribution of individual effects. This method has seen several successful applications, including recent work by honore2024moment and dobronyi2024identificationdynamicpanellogit in the context of dynamic logit models, and bonhommeDano2024functional in network models.

We show that when our method is applied to: (i) the common parameters as the functional of interest; (ii) conditionally parametric models; and (iii) under additional support conditions, it yields functional differencing as a special case. Thus, ours is the first paper to link the literature of functional differencing with the DML literature. This link is useful, as it enables the application of the full suite of theoretical and computational tools developed in the DML literature to functional differencing. In particular, this link provides a formal framework for the theoretical justification of cross-fitted efficient implementations of functional differencing. We illustrate the substantial efficiency gains of employing cross-fitted efficient functional differencing in a Monte Carlo simulation study.

As bonhommeDano2024functional succinctly write: “The scope of functional differencing is not limited to panel data, and the approach can be applied to any setting with a parametric conditional distribution involving latent variables”. Our paper shows that the seminal contribution by bonhomme2012functional on functional differencing has much broader applicability, beyond conditionally parametric models. However, the moments it produces must be modified to achieve orthogonality for general functionals. This is our main contribution relative to the literature of functional differencing. We fully characterize these additional restrictions that functional differencing moments must satisfy through score operators and their adjoints.\footnote{The first version of this paper, arganaraz2023existence, has characterization results on general models, beyond models with UH.}

Recently, bonhomme2024neyman has extended the concept of Neyman (first-order) orthogonality to higher-order orthogonal moments within the conditionally parametric panel data and large $T$ setting. In contrast, our analysis is restricted to first-order orthogonal moments, but applies to more general models and fixed $T$, including high-dimensional panel data models and the Kotlarski model. In this sense, our results are complementary to those of bonhomme2024neyman.

More broadly, our paper contributes to the recent literature on how to systematically construct Neyman-orthogonal moments. Our approach complements the influence function and GMM-based approach of chernozhukov2022locally by considering a likelihood approach for models with nonparametric UH. The construction of orthogonal moments is the fundamental step in the implementation of popular DML inferences, see chernozhukov2018double. We propose novel orthogonal moments and DML inferences in panel data models with high-dimensional covariates, the Kotlarski model and the teacher value-added models. These applications demonstrate the utility of the general characterization of orthogonal moments in this paper.

Examples, Functionals and Heuristics

A High-Dimensional Linear Random Coefficient Model

Consider the following linear Random Coefficient model

where $Y_{i}=(Y_{i1},...,Y_{iT})$ is a $T\times1$ vector of dependent variables, for a finite and fixed $T,$ $X_{i}=(W_{i},V_{i})$ are strictly exogenous covariates, for $T\times p$ and $T\times q$ matrices $W_{i}\ $and $V_{i},$ respectively, $\alpha_{i}$ is a $q\times1$ vector of individual heterogeneity (UH), and $\varepsilon_{i}$ are idiosyncratic (unobserved) errors.\footnote{Other names for $\alpha_{i}$ include fixed effects, individual-specific effects, incidental parameters, etc.} The data observation is $Z_{i}=(Y_{i},X_{i})$ for a large number of units $i=1,...,n.$ We assume that ${(Z_{i},\alpha_{i}):i=1,...,n}$ are independent and identically distributed (iid), with $Z_{i}$ following the distribution $\mathbb{P}_0$. To simplify notation, we drop the reference to $i$ when is convenient. This model is parametrized by $\lambda_{0}=\left( \theta_{0},\eta_{0}\right)$ where $\theta_{0}=(\beta_{0},f_{\left. \varepsilon\right\vert \alpha,X})\in\mathbb{R}^{p}\times\mathcal{F} _{\mathcal{\varepsilon}}$ and $\eta_{0}=f_{\left. \alpha\right\vert X},$ and where henceforth, for generic random vectors $U_{1}$ and $U_{2},$ $f_{\left. U_{1}\right\vert U_{2}}$ and $f_{U_{1}}$ denote, respectively, the conditional density of $U_{1}$ given $U_{2}$ and the marginal density of $U_{1}$. The set $\mathcal{F}_\varepsilon$ of error distributions might be large, i.e., we are not necessarily imposing parametric assumptions on $f_{\left. \varepsilon\right\vert \alpha,X}$. We will be concerned with situations where $p$, the dimension of $W$, might be large, specifically, larger than the sample size. This is consistent with modern applications of panel data that may contain a myriad of time-varying characteristics or with flexible specifications of a low-dimensional set of controls that can be well approximated by a dictionary given by the columns of $W$.\footnote{We note that our results can be extended to nonlinear (in covariates and common parameters) models, as in Equation (4.1) of chamberlain1992efficiency, but with high-dimensional regressors; however, this extension is beyond the scope of the current paper.}

In the more classical setting with $p$ fixed, chamberlain1992efficiency proposed GMM estimators for $\beta_{0}$ and $\mathbb{E} \left[ \alpha\right] $ when the following strict exogeneity condition on $f_{\left. \varepsilon\right\vert \alpha,X}\ $holds,

He also considered the key assumption that $\mathbb{P}_{0}\left( rank(V^{\prime}V)=q\right) =1$ and $q<T$. This is a type of non-surjectivity condition which plays a key role for the existence of orthogonal moments, see the first version of this paper arganaraz2023existence for further discussion on existence. Under additional covariance restrictions (independence conditions of UH) of the idiosyncratic errors $\varepsilon$ given the exogenous regressors, arellano2012identifying showed identification of the variance of $\alpha$ (and its distribution, respectively). All these results consider a fixed $p$.

We present in Section (ref) novel orthogonal moments, DML inferences, and estimators for linear combinations of $\beta_{0}$, average marginal effects, and quadratic forms of the UH (e.g., variances) in this model with high-dimensional regressors, i.e., when $p$ is potentially larger than the sample size.

Nonlinear Panel Data Models and Functional Differencing

A strand of the literature on panel data has considered a model for $f_{Y|X}$ given by\footnote{Henceforth, for the sake of exposition, we drop the region of integration, as it will be clear from the context.}

where $f_{Y|\alpha,X}\left( y|\alpha,x;\theta_{0}\right)$ is the conditional density of $Y$ given $(\alpha,X),$ which is assumed to be known up the finite-dimensional parameter $\theta_{0}$. This setting includes static and dynamic models, as $X$ may include exogeneous variables and initial conditions. We refer to this class of models as conditionally parametric models. It is in this setting where the recent and important literature on functional differencing has been developed, bonhomme2011,bonhomme2012functional,honore2024moment,dobronyi2021identificationdynamicpanellogit,davezies2021identification,bonhomme2024neyman.

We show that functional differencing is a special case of our approach for the common parameters as the functional of interest, when $\theta_{0}$ is a finite-dimensional parameter. Our characterization of orthogonal moments applies to more general models and other functionals, such as average marginal effects. We note that the moments resulting from functional differencing are generally not locally robust for general functionals (e.g., for average marginal effects). We thus provide novel orthogonal moments for those. We also develop DML-UH inferences allowing for partial identification of the multivariate target and nuisance parameters, and provide a theoretical framework for efficient implementations of functional differencing. We illustrate in a small Monte Carlo exercise the efficiency gains of optimal implementations of functional differencing and the importance of relevant orthogonal moments in practice.

Kotlarski Model with a Factor Loading

This example illustrates that we allow for nonparametric specifications of $f_{\left. Y\right\vert \alpha}$. Specifically, suppose the observed data $Z=(Y_{1},Y_{2})$ satisfies

where $\alpha$ is a common unobserved factor, independent of the bivariate error $\varepsilon=(\varepsilon_{1},\varepsilon_{2})$, which also has independent components with zero mean and densities $f_{\varepsilon_1}$ and $f_{\varepsilon_2}$, respectively. The scalar parameter $\beta_{0}$ is the unknown factor loading, with $\beta_{0}\neq0$. The density of the data is given by ((ref)) with $X$ empty, $Z=Y=(Y_{1},Y_{2})$, and \[ f_{Y/\alpha}(y_{1},y_{2};\theta_{0})=f_{\varepsilon_{1}}(y_{1}-\alpha)f_{\varepsilon_{2}}(y_{2}-\beta_{0}\alpha). \] Thus, here $\theta_{0}=(\beta_{0},f_{\varepsilon_{1}},f_{\varepsilon_{2}})$ is infinite-dimensional. This model is called the Kotlarski model with a factor loading in lewbel2022kotlarski, and it generalizes the classical Kotlarski model first used in econometrics by li1998nonparametric. See rao1983nonparametric for a description of the basic Kotlarski's lemma and evdokimov2012some for extensions. The identification of the model above has been studied in detail by lewbel2022kotlarski.

The Kotlarski model has several prominent applications in economics. It has been used to study income dynamics, where $Y_1$ and $Y_2$ denote individual earnings at two different periods and $\alpha$ represents a persistent component of earnings bonhomme2010generalized, horowitz1996semiparametric. cunha2010estimating analyze the formation of cognitive and noncognitive skills in childhood, where $Y_1$ and $Y_2$ are different measurements, e.g., test scores and $\alpha$ is unobserved skill.\footnote{cunha2010estimating do not necessarily assume independence between $\varepsilon_1$ and $\varepsilon_2$.} fruehwirth2016timing focus on the treatment effect of grade retention (or grade repetition), where $Y_1$ is an observed measure that is involved in the retention decision and $Y_2$ is some outcome of interest, e.g., test score, and $\alpha$ is an unobservable common “cause" that determines the effect of grade retention and whether a student is retained. Other interesting applications have appeared to measure the importance of information sets while making optimal decisions navarro2017identifying and in the intergenerational mobility literature nybom2016heterogeneous, haider2006life. lewbel2022kotlarski discusses further applications of the model (e.g., to panel data), and lewbel2024identification has extended these applications to a triangular two-equation system without instruments.\footnote{Note that some of the aforementioned applications control for strictly exogenous regressors that appear in (ref)-(ref). For simplicity, we focus on the case where these regressors are absent from the analysis. Our results could be extended to one-step or two-step methods that account for exogenous covariates.}

The Kotlarski model is representative of a large class of semiparametric models with a similar mathematical structure for which our results apply. To illustrate, consider the general class of non-classical measurement error models put forward by hu2008instrumental. In these models $\alpha$ is a latent (unobserved) regressor, with a contaminated counterpart $Y_2$ and a dependent variable $Y_1$, so $Y=(Y_1,Y_2)$. The vector $X$ is a set of instrumental variables. Under their conditions, the model is a semiparametric mixture model as above, but with \[ f_{Y/\alpha}(y_{1},y_{2};\theta_{0})=f_{Y_1/\alpha}(y_{1};\beta_{0},\alpha)f_{Y_2/\alpha}(y_{2};\alpha), \] where $\theta_{0}=(\beta_{0},f_{Y_2/\alpha})$ is again infinite-dimensional. hu2008instrumental proposed a joint sieve method for $\theta_{0}$ and $\eta_{0}$. Other examples with a similar mathematical structure include nonparametric panel data with interactive fixed effects in freyberger2018non and learning models in bunting2024heterogeneity, among many others.

In this application, we obtain debiased inferences for the $k$-order moment of the common factor and the parameter $\beta_{0}$ that do not require estimating the distribution of UH and other high-dimensional nuisance parameters. Thus, our methods in this application overcome the aforementioned limitations of plug-in functionals.

Teacher Value-Added

In this model, $ij$ indexes student-teacher pairs (teacher $i$, student $j$), $t=1,\dots, T$ denotes periods, $Y_{ijt}$ captures $j$-student's test scores of a particular teacher $i$ at time $t$, which follows the model $$ Y_{ijt}=X_{ijt}^{\prime}\beta_{0}+\alpha_{i}+\varepsilon_{ijt}, $$ where $X_{ijt}$ is a $q-$dimensional vector of observed characteristics of the student and the teacher. The literature assumes that the entries of $\varepsilon$ are iid normal (across $i$ and $T$) with unknown and common variance and that the linear projection of $Y_{ijt}$ on $X_{ijt}$ identifies $\beta_{0}$. Recently, gilraine2020new has proposed a nonparametric maximum likelihood estimation method for the distribution of UH from the equation that averages over $j$ and $t$, denoted by $$ Y_i = \alpha_i + \varepsilon_i. $$ We propose LR moments for smooth (analytic) functionals of the teacher value-added, which are free from its distribution, and also show that LR moments do not exist for some empirically relevant policy parameters, see Section (ref). These results abstract from the estimation of $\beta_{0}$ and the selection of controls $X_{ijt}$, as the results from the high-dimensional random coefficient model apply here to provide robust inferences to the selection of controls.

Functionals of Interest

Researchers working with the previous models might be interested in various parameters, which typically depend on the UH, i.e., functionals of UH. We will denote a generic parameter by $\psi_0 \equiv \psi(\lambda_0)$, with $\lambda_0 = (\theta_0,\eta_0)$. Some examples of these functionals include:

Most of the previous quantities are special cases of the following generic functional

where $r$ is known up to $\theta_0$. Our results apply to more general functionals of UH than in ((ref)), but we focus for the most part on this generic functional for simplicity of exposition.

A Heuristic Summary of the Main Results

The main result of this paper can be heuristically stated as follows. For a general functional $\psi_{0}=\mathbb{E}\left[r\left(\alpha,X,\theta_0\right)\right]$, and for a model with nonparametric UH, the set of relevant Neyman-orthogonal moments for $\psi_{0}$ is characterized as the set of non-zero moment functions $g(Z,\rho_{0},\psi_{0})$ with mean zero and finite variance, depending on additional nuisance parameters $\rho_{0}$ (which may include $\theta_0$ as well as other parameters different from $\psi_0$), satisfying the following two equations:

and

The first equation ((ref)) depends on the distribution of UH only through its support. This equation makes the moment $\mathbb{E}\left[ g(Z,\psi_{0},\rho_{0})\right]$ locally robust to the distribution of UH. In fact, under some commonly employed additional support conditions, the solutions to this equation will be globally robust to the distribution of UH, i.e., $\rho_{0}$ will not contain $\eta_0$.

The second equation is a heuristic and commonly used locally robust condition for the orthogonal moment with respect to the additional nuisance parameters. We provide a formal (and corrected) version of ((ref)) in terms of pathwise derivatives, or equivalently, in terms of a certain score operator and the derivative of the functional $\psi(\lambda_0)$.

The main result shows that these two equations are necessary and sufficient for a relevant orthogonal moment, up to scale, with the necessary part being the most difficult part to establish. We then propose a three-step procedure to construct a moment function $g$ that satisfies these two equations. We illustrate the wide applicability of our three-step procedure in our examples and their finite sample performance through a Monte Carlo and an empirical application.

The next section provides a formal framework and regularity conditions to prove the main result. Applied researchers can skip this and the next sections and move directly to the examples in Section (ref), where the general results are applied.

Existence and Relevance of Debiased Moments

This section formally introduces the model, regularity conditions, and the main theoretical result of the paper: the necessary and sufficient conditions for relevant orthogonal moments in models with nonparametric UH. We introduce some notation that will be used throughout the paper.

Notation: Let $f_{\lambda}$ denote the density of $\mathbb{P}_{\lambda}$ with respect to (wrt) a $\sigma$-finite measure $\mu$, and let $f_{0}\equiv f_{\lambda_0}$ be the density pertaining to $\mathbb{P}_{0}$. Let $L_{2}\equiv L_{2}(f_{0})$ denote the Hilbert space of $\mathbb{P}_{0}-$square integrable measurable functions with inner product $\langle g_{1},g_{2}\rangle=\mathbb{E}\left[ g_{1}\left( Z\right) g_{2}\left( Z\right) \right] $ and norm $\left\Vert g\right\Vert_{2} ^{2}=\mathbb{E}\left[ g^{2}\left( Z\right) \right] $. The set $L^{0}_{2}\equiv L^{0}_{2}(f_{0})$ is the subspace of zero-mean functions in $L_2,$ i.e., $g$ with $\mathbb{E}\left[ g\left( Z\right) \right] =0$ and $\mathbb{E}\left[ g^{2}\left( Z\right) \right] <\infty.$ More generally, $L_{2}(f)$ and $L_{2}^{0}(f)$ are defined analogously for any density $f.$ For a vector $\mathbf{v}$ in $\mathbb{R}^{p}$, we denote $\left|\mathbf{v}\right|_1$ and $\|\mathbf{v}\|$ the $l_1$ and $l_2$ norms, respectively. If $A$ is a matrix, let $\left|\left|A\right|\right|$ be its spectral norm, $vec(A)$ its vectorization, let $A^{\prime}$ be the transpose of $A$, and let $A^{\dagger}$ denote the generalized Moore-Penrose inverse of $A$. Henceforth, $I_{r}$ denotes the identity matrix of order $r\ $($r\times r$ identity matrix). For any linear operator $S$, let $\mathcal{R}(S)$ and $\mathcal{N}\left(S\right)$ denote its range and null space, respectively. For any linear space $V$, let $\overline{V}$ be its closure and $\Pi_{\overline{V}}$ be the orthogonal projection operator onto $\overline{V}$ (the topology will be clear from the context). For any subset $B$ of a generic Hilbert space $\mathcal{H}$, let $B^{\perp}$ be the orthogonal complement of $B$ in $\mathcal{H}$. Henceforth, $C$ is a generic positive constant that can change from expression to expression.

Model, Regularity and Scores

The model we consider has a density with respect to $\mu=\mu_{Y}\times v_{X},$ where $\mu_{Y}$ is a $\sigma$-finite measure and $v_{X}$ is the probability measure of $X$, given by

where $f_{Y|\alpha,X}\left( y|\alpha,x;\theta_{0}\right)$ is the conditional density of $Y$ given $(\alpha,X)$, which is known up to the parameter $\theta_0$, and $\eta_{0}(\alpha|x)$ is the conditional density of the UH given the exogenous regressors $X$, which is unknown and unrestricted (i.e., we consider a fixed effects approach). The model is indexed by the parameter $\lambda_{0}=\left(\theta_{0},\eta_{0}\right) \in \Lambda\equiv\Theta\times\Xi$, where $\Theta$ and $\Xi$ are subsets of possibly infinite-dimensional spaces. For simplicity of presentation, we consider the UH to be continuously distributed, but this can be trivially relaxed. Additionally, we note that the dimension of the UH is unrestricted for the moment.

The observed data is an iid sample $Z_{i}=(Y_{i},X_{i})$, $i=1,...,n,$ from $\mathbb{P}_{0}\in\mathcal{P}=\{d\mathbb{P}_\lambda/d\mu =f_{\lambda} \text{ in } (\ref{mixt2}):\lambda=(\theta ,\eta)\in\Lambda\}$. The parameters that generated the data are denoted by $\lambda_{0}=\left( \theta_{0},\eta _{0}\right) \in\Lambda,$ i.e., $\mathbb{P}_{0}=\mathbb{P}_{\theta_{0},\eta _{0}}.$ However, the model may be partially identified in the sense that the equation $\mathbb{P}_{0}=\mathbb{P}_{\lambda}$ may have more than one solution in $ \lambda \in \Lambda.$

Our first assumption is that the model in ((ref)) is regular, as defined in Section (ref) of the Supplementary Appendix. Regularity includes the concept of Differentiability in Quadratic Mean (DQM), which is a standard smoothness assumption in the literature, see newey1990semiparametric.

To explain what regularity entails, we need to introduce some notation that will also be useful for subsequent developments. To model deviations from $\theta_0$ and $\eta_0$, we introduce the Hilbert space $\mathbf{H}=\mathcal{H}_{\theta}\times \mathcal{H}_{\eta} $, which has the inner product $\langle(\delta_{1},b_{1} ),(\delta_{2},b_{2})\rangle_{\mathbf{H}}:=\langle\delta_{1},\delta_{2} \rangle_{\mathcal{H}_{\theta}}+\langle b_{1},b_{2}\rangle_{\mathcal{H}_{\eta} }$, and where $\mathcal{H}_{\theta}$ and $\mathcal{H}_{\eta}\ $are Hilbert spaces endowed with the inner products $\langle\cdot,\cdot\rangle _{\mathcal{H}_{\theta}}$ and $\langle\cdot,\cdot\rangle_{\mathcal{H}_{\eta}},$ respectively$.$ Since we are considering a fixed effects approach, the space of deviations from $\eta_0$ is unrestricted and given by $\mathcal{H}_{\eta}=L_2(\eta_0\times f_X)$, where $f_X$ is the density of $X$ (with respect to some dominating measure). The inner product associated to $\mathcal{H}_{\eta}$ is then simply $\langle b_{1},b_{2}\rangle_{\mathcal{H}_{\eta} }=\mathbb{E}\left[ b_{1}(\alpha,X)b_{2}(\alpha,X)\right]$.

We will then consider paths $\tau\in\lbrack0,\epsilon)\rightarrow\lambda_{\tau}=(\theta_{\tau},\eta_{\tau})\in \Lambda$ such that $d\log\eta_{\tau}(\alpha)/d\tau= b\in\mathcal{H}_{\eta}$ and $d\theta_{\tau}/d\tau= \delta\in\mathcal{H}_{\theta}$, where henceforth $d/d\tau$ is the derivative from the right (i.e., for nonnegative values of $\tau$) at $\tau=0$. To achieve this, it will often suffice to consider $\theta_{\tau}=\theta_0+\tau \delta$ and $\eta_{\tau}=(1+\tau b)\eta_{0}$. Note that $b$ has the interpretation of the (conditional) score for $\eta_0$. If $b$ is bounded with $\int b(\alpha,x) \eta_0(\alpha|x) d\alpha = 0$ a.s., then $\eta_{\tau}$ will be a density for a small enough $\tau$.

Under regularity, scores are well-defined in a quadratic mean sense, but for simplicity of notation and exposition, in the main text, we will express them pointwise, as in \[ S_{\theta} \delta = \frac{d}{d\tau} \log f_{\theta_{\tau}, \eta_{0}}(z), \] \[ S_{\eta} b = \frac{d}{d\tau} \log f_{\theta_{0}, \eta_{\tau}}(z). \] Since these derivatives may not exist for all deviations $h=(\delta,b)$ in $\mathbf{H}=\mathcal{H}_{\theta}\times \mathcal{H}_{\eta}$, we consider a subset $\Delta(\lambda _{0})=B(\theta_{0})\times B(\eta_{0})\subseteq\mathbf{H,}$ with the property that $B(\eta_{0})$ is dense in $\mathcal{H}_{\eta}\cap L^{0}_2(\eta_0)$ (consistent with the fixed effects approach). Thus, an implication of regularity and the chain rule is that the score of the model corresponding to the path $\lambda_{\tau}$ is

Moreover, from a direct calculation with the likelihood of the model, in Section (ref) of the Supplementary Appendix, we show that the scores for $\theta$ and $\eta$, have the expressions given, respectively, by \[ S_{\theta} \delta (z)= \mathbb{E}\left[ \left. \ell_\theta(Y\vert\alpha,X,\theta_0 )[\delta]\right\vert Z=z\right], \] \[ S_{\eta} b (z)= \mathbb{E}\left[ \left. b(\alpha,X )\right\vert Z=z\right], \] where $\ell_\theta(Y\vert\alpha,X,\theta_0 )[\delta]$ is the conditional score operator for $\theta$ at $\theta_0$, defined formally in Section (ref) of the Supplementary Appendix. These formulas for the scores, and in particular that for $S_\eta$, will have fundamental consequences for the characterization of orthogonal moments. The score operators $S_{\theta}$ and $S_{\eta}$ transform \textquotedblleft scores\textquotedblright\ of the parameter spaces in $\Delta(\lambda _{0})\subseteq\mathbf{H,}$ into scores of the model in $L_{2}$.\footnote{See, e.g., begun1983information for an introduction to score operators.}

To illustrate the calculation of $S_{\theta}\delta$ for a simple case, consider conditionally parametric models where $\theta$ is finite-dimensional (say $p$-dimensional). In this case, $\mathcal{H}_{\theta}$ is $\mathbb{R}^p$, the inner product $\langle\cdot,\cdot\rangle _{\mathcal{H}_{\theta}}$ is the standard Euclidean inner product, and $S_{\theta}\delta(z)=s_\theta(z)'\delta$, where $s_\theta(z)=\partial \log f_{\theta_0,\eta_0}/\partial\theta$ is the classical score for $\theta$ at $\theta_0$. Here, $\ell_\theta(y\vert\alpha,x,\theta_0 )[\delta]=(\partial \log f_{Y|\alpha,X}\left( y|\alpha ,x;\theta_{0}\right)/\partial\theta)'\delta$ is the conditional score in the “direction” $\delta\in\mathbb{R}^p$. Further discussion on regularity is given in Section (ref) of the Supplementary Appendix.

Relevant orthogonal moments

Inference on $\psi(\lambda_{0})$ is based on a non-zero vector of moment functions $g\in L_2$ such that, for all $\lambda_0\in\Lambda$,

where $Z$ is distributed as $Z_{i}\ $and $\mathbb{E}$ denotes expectation under $\mathbb{P}_{0}$, i.e., under $\lambda_0$.

Our goal is to construct relevant orthogonal moments for inference on $\psi(\lambda_{0})$. To explain what relevant orthogonal moments for $\psi(\lambda_{0})$ are, let again $\lambda_{\tau}$ denote a path in $\Lambda,$ i.e., a mapping $\tau \in\lbrack0,\epsilon)\rightarrow\lambda_{\tau}\in\Lambda,$ with $\epsilon>0,$ such that $\lambda_{\tau}$ equals $\lambda_{0}$ at $\tau=0$ and where other regularity conditions given in Section (ref) of the Supplementary Appendix are satisfied.

In addition to the regularity of the model, we also require regularity conditions on the moments to be able to interchange derivatives and expectations. This will be a maintained assumption throughout for all the moments we consider. Let $\mathbb{E}_{\tau}$ denote the expectation under the path $\mathbb{P}_{\lambda_{\tau}}\equiv\mathbb{P}_{\tau}.$

Then, we give a formal definition of an orthogonal moment function for $\psi(\lambda)$ at $\psi_{0}$.

We will show below that the lack of sensitivity of the moment $\mathbb{E}\left[ g(Z,\lambda_{\tau})\right] $ for all restricted paths fixing the parameter of interest, i.e., with $\psi (\lambda_{\tau})=\psi_{0},$ is the fundamental notion to extend DML to general functionals. The DML literature has considered moments with a partition of the parameters as $\lambda_0=(\psi_0,\gamma_0)$, for some $\gamma_0$, but such partition is not possible for most functionals $\psi(\lambda_0)$.

We now define relevant orthogonal moments.

Note that the concept of relevance is in a weak sense, as the non-zero derivative has to hold for some path for which $\psi(\lambda_{\tau})\neq\psi_{0}$ and not for all such paths. Relevant orthogonal moments give rise to score tests with a non-trivial local power function. Further discussion on local power analysis and relevance is given in Section (ref) of the Supplementary Appendix.

Smooth Functionals

After discussing the regularity of the model and the moments, this section introduces smoothness assumptions on the functionals considered. For simplicity of exposition, we focus on a scalar functional

Extensions to vector functionals are trivial.\footnote{The scalar case is without loss of generality. The case $\psi(\lambda_{0})\in\mathbb{R}^{d_{\psi}}$, with $d_{\psi}>1$, can be analyzed by considering one-dimensional projections $C'\psi(\lambda_{0})$, for a selected basis of $C\in\mathbb{R}^{d_{\psi}}$. In fact, our results also hold for $d_{\psi}=\infty$, but we focus on a finite-dimensional target parameter for simplicity of exposition.} Also, extensions to nonlinear functionals of UH are considered in Section (ref), and are permitted in our theory. Examples of empirically relevant nonlinear functionals are quantiles of the distribution of the UH.

This assumption is weaker than differentiability of the functional in the sense of van1991differentiable, since in addition to Assumption (ref), differentiability requires a special representation for $r$ and $r_\theta$; see Remark (ref) for further discussion.

To illustrate this assumption, in the case of a finite-dimensional $\theta$, it requires the existence of the standard derivative $r_{\theta}=\partial \mathbb{E}\left[ r(\alpha,X,\theta_{0})\right]/\partial \theta$.

Main result: existence of relevant orthogonal moments

Under the regularity of the model (Assumption (ref)), the adjoint operator of $S_{\theta}$ is well-defined: there exists a continuous linear operator $S_{\theta}^{\ast} :L_2\rightarrow\overline{B(\theta_{0})}$ satisfying for all $\delta\in\overline{B(\theta_{0})}\subset{\mathcal{H}_{\theta}}$ and all $g\in L_2,$

To illustrate, when $\theta$ is finite-dimensional, since $S_{\theta}\delta(z)=s_\theta(z)'\delta$, with $s_\theta(z)=\partial \log f_{\theta_0,\eta_0}/\partial\theta$, we have that $\mathbb{E}\left[\delta's_\theta(Z) g(Z)\right] =\delta'\mathbb{E}\left[s_\theta(Z) g(Z)\right]$, and hence $S_{\theta}^{\ast}g=\mathbb{E}\left[s_\theta(Z) g(Z)\right]$. Therefore, $S^{*}_\theta$ is continuous when $\mathbb{E}\left[\|s_\theta(Z)\|^2 \right]<\infty$, which is a standard assumption in the literature. Furthermore, by the information equality, for a regular $g\in L^{0}_2$,

For future reference, the adjoint operator for $\eta$ is given by

We are ready to give the main theoretical result of this paper. Define the $\eta$-influence function

and recall $r_{\theta}$ has been defined in Assumption (ref).

We give a corollary of this theorem, which may be very useful for identification and for a simpler verification of the orthogonality conditions.

To summarize the main result of this section, Theorem (ref) characterizes all Neyman-orthogonal moments for general smooth functionals (ref) in models with nonparametric UH.\footnote{Extensions to semiparametric UH are also possible, see the first version of this paper in arganaraz2023existence.} The existing DML literature had only limited results for such functionals.

Construction of Relevant LR\ moments

Theorem (ref) is constructive, as it provides us with a method to construct debiased moments. In this section, we describe a three-step procedure that directly computes an LR moment function without estimating UH, as we next show. For simplicity of exposition, in this section we focus on the regular scalar case where $\psi(\lambda_0)$ has positive information.\footnote{This gives a strong sense of relevance, where the corresponding test has non-trivial power in all directions.}

Suppose that we can find a $g_0 \in L_2$ such that

where $\rho_{01}=\rho_{1}(\lambda_0)$ may contain $\theta_0$ and additional nuisance parameters.

Moreover, suppose that we can find a $g_1 \in L_2$ such that

Finally, let $\Gamma_0$ be such that, for all $h\in \Delta(\lambda_0)$ and $\rho_{\tau1}=\rho_{1}(\lambda_\tau)$,

In all our examples, we have found that it suffices to take $\Gamma_{0}(X)$ as a constant $\Gamma_{0}$. Therefore, for simplicity of exposition, we henceforth consider this case. Following the proof of Corollary (ref), condition ((ref)) is equivalent to

For example, for a finite-dimensional $\rho_{01}=\theta_0$ equation ((ref)) reads

This equation has a unique solution $\Gamma_0$ when, for example, $g_1$ satisfies a full column rank condition (i.e., $\theta_0$ is locally identified). We only need the existence of $\Gamma_0$ (a minimum norm solution $\Gamma_0$ can be selected, for example).

In the general case, we can obtain an orthogonal moment for $\psi(\lambda_{0})=\mathbb{E}\left[ r(\alpha,X,\theta_{0})\right]$ as

where $\rho_0=(\rho_{01},\rho_{02})$ and $\rho_{02}=\Gamma_{0}$.

To see that (ref) leads to an orthogonal moment, we just need to verify that it satisfies (ref)-(ref) in Theorem (ref), or (ref)-(ref) in Corollary (ref) with $m_0(Z,\rho_{0})=g_{0}(Z,\rho_{01})-\Gamma _{0}g_{1}(Z,\rho_{01}).$ Indeed,

and by the definition of $\Gamma_0$,

where we have used that $S^{*}_\theta \psi_0=0$ by the zero mean property of scores. Hence, (ref)-(ref) hold, and (ref) is an orthogonal moment for $\psi(\lambda_0)$.

We summarize this in the following three-step algorithm:

Algorithm for obtaining an orthogonal function for $\psi(\lambda_{0})=\mathbb{E}\left[ r(\alpha ,X,\theta_{0})\right]$:

A orthogonal moment function for $\psi(\lambda_{0})$ can be then constructed as

The moments resulting from Step 1, however, will not be LR\ for $\psi(\lambda_{0})$ in general, as they will be affected by the estimation of $\theta_0$, which motivates our three-step procedure. For finite-dimensional $\theta_0$, honore2024moment have obtained analytical solutions for Step 2 in dynamic discrete choice panel data models, and bonhomme2012functional has proposed numerical approximations for solutions in Step 2, and estimators for $\theta_{0}$ when this parameter is identified. Step 3 involves Jacobians of the moments, whose estimation is discussed in bonhomme2012functional, p. 1366. For high-dimensional $\theta_0$, we show how these steps are solved in our applications to the high-dimensional random coefficient model, Kotlarski model, and the teacher value-added example.

DML Inferences

In this section, we propose a score-type test based on orthogonal moments for the hypothesis $H_{0}:\psi(\lambda_{0})=\psi_{0}$ vs $H_{0}:\psi(\lambda_{0})\neq\psi_{0}$, where we some abuse of notation, in this section $\psi_{0}$ denotes a specified value under the null (e.g., $\psi_{0}=0$). The proposed test is in the spirit of the Anderson-Rubin test, see anderson1949estimation. For ease of our exposition, we collect all the involved nuisance parameters, i.e., parameters different from $\psi_0$, in the vector $\rho_0 = \left(\rho_{01},\rho_{02}\right)$, for some suitable partition. The nuisance parameter $\rho_0$ may contain $\eta_0$ or be part of it. We write the debiased moment as follows:

such that $\phi_0\left(Z,\rho_{01},\psi_0\right), \phi_1\left(Z,\rho_{0},\psi_0\right) \in L^0_2$. Decomposition (ref) is useful to develop our asymptotic theory in the Appendix. The correction term could be $\phi_1\left(Z,\rho_{0},\psi_0\right)=-\Gamma_{0}g_{1}(Z,\rho_{01})$ when $\phi_0\left(Z,\rho_{01},\psi_0\right) = g_{0}(Z,\rho_{01}) - \psi_0$. In this case, $\rho_{02}=\Gamma_{0}$, which satisfies the key global robustness property used in the asymptotic theory that $\mathbb{E}\left[\phi_1\left(Z,\rho_{01},\rho_{2},\psi_0\right)\right]=0$ for all $\rho_{2}$. Our setting here generalizes chernozhukov2022locally to a partial identification framework.\footnote{Recently, bennett2023source has studied locally robust inference on linear functionals in nonparametric instrumental variables models.}

As in the current DML literature, an important element in our construction of orthogonal moments is cross-fitting. It has been commonly used in the semiparametric literature, see, e.g., bickel1982adaptive, klaassen1987consistent, and schick1986asymptotically. Cross-fitting allows us to deal with high-dimensional estimators and avoid an “own observation” bias; see the discussion by chernozhukov2018double. It is implemented as follows. Given the data $\left\{ Z_{i}\right\} _{i=1} ^{n}$, we randomly partition the observation indices $(i=1,...,n)$ into $L$ groups $I_{\ell},$ $(\ell=1,...,L)$ of equal size. Let $\hat{\rho}_{\ell}$ denote an estimator of $\rho_{0}$ using all observations not in $I_{\ell}$; see below for the specific choices of $\hat{\rho}_{\ell}$ in our examples. Thus, we here assume that there exists an estimator of the components in $\rho_0$ to implement our DML inferences.\footnote{Recall that $\rho_0$ might not necessarily contain the entire vector of nuisance parameters $(\theta_0,\eta_0)$.} If $\rho_{0}$ is low-dimensional there is no need to use cross-fitting, and we can set $\hat{\rho}_{\ell}=\hat{\rho}$ based on all the observations. The estimated orthogonal moment is then \[ \hat{g}_{i\ell}(\psi_{0})\equiv g(Z_{i},\hat{\rho}_{\ell},\psi_{0}). \] For a given (estimated) $k-$dimensional LR\ moment $\hat{g},$ a debiased cross-fitted sample moment functions is

Based on our estimated LR moments, (ref) we propose a test to conduct robust inference on $\psi(\lambda_0)$. Let the cross-fitted weighting matrix be \[ \check{W}_{n}=\frac{1}{n}\sum_{\ell=1}^{L}\sum_{i\in I_{\ell}}\hat {g}_{i\ell}(\psi_{0})\hat{g}_{i\ell}^{\prime}(\psi_{0}), \] which under regularity conditions consistently estimates $W_{g}=\mathbb{E} \left[ g(Z,\rho_{0},\psi_0)g(Z,\rho_{0},\psi_0)^{^{\prime}}\right] .$

We allow for the possibility that $W_{g}$ is singular, as when there is partial identification. Define $0<k_{g}:=rank(W_{g})\leq k$. To construct a consistent estimator for $W_{g}^{\dagger},$ we follow a common nonparametric strategy in the literature, see andrews1987asymptotic, or more recently lee2024locallyregularefficienttests.\footnote{When information of the rank $k_{g}$ is available, more efficient and simpler to implement estimators of $W_{g}^{\dagger}$ might be available. Here, we focus on the most difficult case where $k_{g}$ is unknown, but of course, our results can also be applied to the simpler case where $k_{g}$ is known, with obvious modifications. Similar remarks apply to the set of eigenvalues and eigenvectors of $W_g$.} Let $\check{W}_{n}=\check{U}_{n}\check{\Lambda}_{n}\check{U} _{n}^{^{\prime}}$ be the eigendecomposition of $\check{W}_{n}$,\footnote{Using either $\hat{g}_{i\ell}\left(\psi_0\right)$ or $\hat{g}_{i\ell}\left(\hat{\psi}_\ell\right)$.} with $\check{\lambda}_{n,i}$ being the $i-th$ element of the diagonal matrix $\check{\Lambda}_{n}\ $of eigenvalues, where we drop the dependence on $\psi_{0}$ for simplicity of notation, and define

where $\check{\Lambda}_{n}(\nu_{n})$ is a diagonal matrix with $i-$th element equal to $\check{\lambda}_{n,i}$ if $\check{\lambda}_{n,i} \geq\nu_{n}$ and zero otherwise, and $\nu_{n}$ is non-negative user-specific sequence that $\nu_{n}\longrightarrow0$ at a suitable rate; see Section (ref) of the Appendix. We derive specific $\nu_n$ for some of our examples below.

Our feasible cross-fitted debiased score test statistic is

The proposed test for $H_{0}:\psi(\lambda_{0})=\psi_{0}$ vs $H_{0}:\psi(\lambda_{0})\neq\psi_{0}$ rejects the null hypothesis at $\zeta$ nominal level if $\hat{C}_{n}(\psi_{0})>c_{\zeta,\hat{k}_{g}}$, where $\hat{k}_{g}=rank(\hat{W})$ and $c_{\zeta,d}$ is the $1-\zeta$ quantile of a chi-square distribution with $d$ degrees of freedom. To define confidence regions for $\psi(\lambda_{0})$, let the parameter space be $\Psi=\{\psi(\lambda ):\lambda\in\Lambda\}$ and the asymptotic confidence region \[ CR_{\zeta}=\{\psi\in\Psi:\hat{C}_{n}(\psi)\leq c_{\zeta,\hat{k}_{g}(\psi )}\}, \] where we emphasize the dependence of the rank $\hat{k}_{g}$ on the parameter $\psi.$ Note that when the conditions of Remark (ref) hold, $\hat{k}_{g}$ and $\hat{W}$ do not depend on $\psi$. In Section (ref) of the Appendix, we give primitive conditions for the asymptotic validity of the proposed tests as well as a detailed local power analysis of it.

Application to the Examples

High-Dimensional Random Coefficient Model

In this section, we apply our theory to the High-Dimensional Random Coefficient Model.

Common Parameters

We first construct novel LR moments for $\psi(\lambda_0)=C_{1}^{\prime}\beta_{0}$ for a known $p\times p_{1}$ matrix $C_{1}$ of full column rank. In an important paper, chamberlain1992efficiency proposed a GMM estimator for $\beta_{0}$ that is locally robust for a fixed $p$. He assumed $\mathbb{P}_{0}\left(rank(V^{\prime}V)=q\right) =1$, and introduced the generalized within-group and between operators, respectively,

Then, by the exogeneity condition ((ref)) and $QV=0$ a.s., it follows that

which leads to the LR moment function $g(Z,\beta_{0})=W^{\prime}Q(Y-W\beta _{0})$ for $\psi(\lambda_{0})=\beta_{0}$ when $p$ is fixed.

In high-dimensional settings where $p$ is large, possibly much larger than the sample size, Chamberlain's GMM estimator requires regularization, resulting in the large biases common of high-dimensional estimators. Moreover, in many applied settings, researchers are not interested in the whole vector $\beta_{0}$ but rather in some linear combination (e.g., a particular coefficient) $\psi(\lambda_{0})=C_{1}^{\prime}\beta_{0}$. The moment derived in chamberlain1992efficiency is not LR for $\psi(\lambda_{0})=C_{1}^{\prime}\beta_{0}$, which is particularly problematic for cases where $p>>nT$. To see this, define the matrices $H_{C_{1}}:=C_{1}(C_{1}^{\prime}C_{1})^{-1}$ and $Q_{C_{1}}:=I-C_{1}(C_{1}^{\prime}C_{1})^{-1}C^{\prime}_{1}$, and note the identity $\beta_0=H_{C_{1}}C_{1}^{\prime}\beta_0+Q_{C_{1}}\beta_0$. In Chamberlain's moment, now $Q_{C_{1}}\beta_0$ becomes a high-dimensional nuisance parameter, and the moment $g(Z,\beta_{0})$ in ((ref)) is sensitive to it.

The form of the LR moment for $\psi(\lambda_{0})=C_{1}^{\prime}\beta_{0}$ and the assumptions needed for it follow from our general necessary and sufficient conditions for the existence of a relevant orthogonal moment. Here, we omit the technical details and give the LR moment directly. A more general definition of the generalized within-group operator is $Q=I_{T}-VV^{\dagger}$. Define also the $p\times p$ matrix $M:=\mathbb{E}\left[ W^{\prime}QW\right]$ and the $p_{1}\times p$ matrix $\rho_{02} := C^{\prime}_1M^{\dagger}$. We need the following assumption.

Assumption (ref) iii) requires that the variation in $W$ not collinear with $V$ (variation in $QW$) is rich enough to identify the component $C'_1\beta_0$. This assumption implies the order conditions $rank(M)\geq p_{1}$ and $rank(V)<T$ a.s. Under Assumption (ref) iii), $\rho_{02}$ is well defined. Before providing a formal result, it is worth noticing what the functions $g_0$, $g_1$, and $\Gamma_0$ look like in this situation. We have $r\left(\alpha,X,\beta_0\right) = C^{\prime}_1\beta_0 = \psi_0$. Next, if we let $$ g_0(Z,\rho_{0},\psi_0) = \rho_{02}W^{\prime}Q\left(Y-WH_{C_1}\psi_0 - WQ_{C_1}\beta_0\right) + \psi_0, $$ where $\rho_0=(\beta_0,\rho_{02})$. By (ref), we obtain $$ \mathbb{E}\left[\left.g_0(Z,\rho_{0},\psi_0) \right|\alpha,X\right] = C^{\prime}_1\beta_0, $$ as required by Step 1 of our algorithm above. Interestingly, direct calculation yields $$ \frac{d}{d\tau}\mathbb{E}\left[g_0(Z,\beta_\tau)\right] = 0. $$ Thus, $g_0$ is already orthogonal to $\beta_0$, which implies that the resulting moment is LR to, e.g., bias selection/bias regularization in $\beta_0$. Hence, we shall let $g_1 = 0$ and $\Gamma_0 = 0$, using our notation above. For these, Step 2 and Step 3 are therefore satisfied. We are now ready to give an identification result and a relevant LR moment for $\psi(\lambda_{0})=C_{1}^{\prime}\beta_{0}$.

The following algorithm describes the steps to estimate the LR moment. Define $Q_i:=I_{T}-V_iV_i^{\dagger}$.

Algorithm for obtaining an LR moment and inferences for $\psi(\lambda_{0})=C_{1}^{\prime}\beta_{0}$:

Our penalized least squares estimator $\hat{\beta}_\ell$ extends the Lasso estimator of belloni2016inference of high-dimensional (linear) panel models to account for multidimensional, possibly non-additive, UH.

Inference is then implemented following Section 5. This LR for $\psi(\lambda_{0})=C_{1}^{\prime}\beta_{0}$ is LR with respect to all nuisance parameters, including $\rho_{02}$ and $Q_{C_{1}}\beta_{0}$. When $p_{1}=1$, under Assumption (ref) the asymptotic variance $W_g$ is positive. We recommend implementing the following cross-fitted variance estimator \[ \hat{W}_{n}=\frac{1}{n}\sum_{\ell=1}^{L}\sum_{i\in I_{\ell}} ({g}_{i\ell}(\hat{\psi}_{\ell})-\bar{g}_{\ell}(\hat{\psi}_{\ell}))^2, \] where $\hat{\psi}_{\ell}=\hat{\rho}_{2\ell}\hat{R}_{\beta,\ell}$, with \[ \hat{R}_{\beta,\ell}=\frac{1}{n-n_{\ell}} \sum_{l\neq\ell}^{L}\sum_{i\in I_{l}}W'_iQ_iY_i,\hspace{1cm} \] \[ \bar{g}_{\ell}(\hat{\psi}_{\ell})=\frac{1}{n}\sum_{\ell=1}^{L}\sum_{i\in I_{\ell}}\hat {g}_{i\ell}(\hat{\psi}_{\ell}). \]

A debiased estimator for $\psi_0$ is \[ \hat{\psi}=\frac{1}{n}\sum_{\ell=1}^{L}\hat{\rho}_{2\ell}\hat{R}_{\beta,\ell}. \] This estimator is asymptotically normal under the conditions of Proposition (ref), as it follows from its proof in Section (ref) of the Supplementary Appendix.

Average Marginal Effects

chamberlain1992efficiency also derived the following moment function $g(Z,\beta_{0},\psi _{0})=H(Y-W\beta_{0})-\psi_{0}$ for the average marginal effects $\psi_{0}=\mathbb{E}\left[ \alpha\right]$ when $p$ is fixed. This moment is motivated from the exogeneity condition, since $HV=I_{q}$ a.s., and hence \[ \mathbb{E}\left[ \left. H(Y-W\beta_{0})\right\vert \alpha,X\right] =\alpha,\text{ a.s.} \] Taking expectations both sides, we conclude $\mathbb{E}\left[ H(Y-W\beta _{0})-\psi_{0}\right] =0.$ However, this moment is not LR for $\psi _{0}=\mathbb{E}\left[ \alpha\right]$, and it is in general, quite sensitive to the estimation of $\beta_{0}$ when $W$ and $V$ are correlated (which is often the case in applications). This sensitivity is particularly problematic for inference on $\psi_{0}$ when $p$ is large, which will be the case in many applications.

We now derive novel LR moments for $\psi(\lambda_{0} )=\mathbb{E}\left[ C_{2}^{\prime}\alpha\right],$ for a known $q\times q_{2}$ matrix $C_{2}$ of full column rank, when possibly $p>>nT$. Define the generalized between operator $H:=V^{\dagger}$, and also define the $q\times p$ matrix $S_1:=\mathbb{E}\left[HW\right]$, and the $p_{2}\times p$ matrix $\Gamma_{0}:=-C_{2}^{\prime}S_1M^{\dagger}$.

This assumption implies the order conditions $rank(M)\geq q_{2}$ and $rank(V) > q_2$. Assumption (ref) iii) requires that the variation in $W$ not colinear with $V$ (variation in $QW$) is rich enough to identify the component $C'_2S_1\beta_0$ (as when $C'_1=C'_2S_1$ in the previous section. This assumption is needed for Steps 2 and 3 of the three-step algorithm. Assumption (ref) iv) is equivalent to the column space of $C_2$ beging contained in that of $V'$ (which is weaker than $rank(V)=q$ a.s.). This is needed for the identification of $\psi(\lambda_0)=\mathbb{E}\left[ C_{2}^{\prime}\alpha\right]$. Note that identification of $\beta_{0}$ or $\mathbb{E}\left[ \alpha\right]$ is not required for $\mathbb{E}\left[ C_{2}^{\prime}\alpha\right]$ to be identified.

We describe an algorithm for the derivation of an LR moment for $\psi(\lambda_{0})=\mathbb{E}\left[ C_{2}^{\prime}\alpha\right]$ as well as for LR inferences based on it. Define $H_i:=V_i^{\dagger}$.

Algorithm for obtaining an LR moment and inferences for $\psi(\lambda_{0} )=\mathbb{E}\left[ C_{2}^{\prime}\alpha\right]$:

Inference is then implemented following Section 5. This moment for $\psi(\lambda_{0})=\mathbb{E}\left[ C_{2}^{\prime}\alpha\right]$ is LR with respect to all nuisance parameters, including $\Gamma_{0}$ and $\beta_{0}$. We recommend implementing the following cross-fitted variance estimator \[ \hat{W}_{n}=\frac{1}{n}\sum_{\ell=1}^{L}\sum_{i\in I_{\ell}} ({g}_{i\ell}(\hat{\psi}_{\ell})-\bar{g}_{\ell}(\hat{\psi}_{\ell}))^2, \] where $\hat{\psi}_{\ell}=\hat{R}_{\alpha,\ell}$ for $\psi(\lambda_{0} )=\mathbb{E}\left[ C_{2}^{\prime}\alpha\right]$, with \[ \hat{R}_{\alpha,\ell}=\frac{1}{n-n_{\ell}} \sum_{l\neq\ell}^{L}\sum_{i\in I_{l}} (C_{2}^{\prime}H_i-\hat{\Gamma}_{\ell}W'_iQ_i)(Y_i-W_i\hat{\beta}_{\ell}), \] \[ \bar{g}_{\ell}(\hat{\psi}_{\ell})=\frac{1}{n}\sum_{\ell=1}^{L}\sum_{i\in I_{\ell}}\hat {g}_{i\ell}(\hat{\psi}_{\ell}). \]

An LR estimator for $\psi_0 = \mathbb{E}\left[C_2^{\prime} \alpha\right]$ can be constructed as $\hat{\psi}=\frac{1}{L}\sum^L_{\ell=1}\hat{R}_{\alpha,\ell}$. This estimator is asymptotically normal, as it follows from the proof of Proposition (ref).

Second Moments

We aim to construct LR moments for quadratic functionals of the UH $\psi_0 = \mathbb{E}\left[\alpha^{\prime}\Omega \alpha\right]$, where $\Omega$ is a $q \times q$ known matrix. We consider the setting of arellano2012identifying. These authors derive a system GMM for common parameters, variances of idiosyncratic errors, and variances of UH by employing the following additional generalized within-group and between operators, respectively,

where $\otimes$ denotes the Kronecker product. arellano2012identifying consider the following covariance restrictions on the conditional variance $\Sigma_i:=\mathbb{V}\left[ \left. \varepsilon_i\right\vert X_i\right]$,

where $vec$ is the vectorization operation, $S_2$ is a known selection matrix and $\omega_{0i}$ is an $m$ dimensional vector of parameters, possibly depending on $X_i$.

Define the parameter $\zeta_0:=(\beta_{0},\omega_0)$ and the generalized errors

Also, to simplify the notation, denote $u(\beta)=Y-W\beta$. We apply the general algorithm for constructing an LR for $\psi_0$ as follows. First, we show that the moment function

satisfies the following conditional moment restriction \[ \mathbb{E}\left[ \left. g_0(Z,\beta_0,\omega_0)\right\vert \alpha,X\right] =\alpha^{\prime}\Omega \alpha,\text{ a.s.} \]

Thus, taking expectations both sides, we conclude $\mathbb{E}\left[ g_0(Z,\beta_0,\omega_0)-\psi_{0}\right] =0.$ However, this moment is not LR for $\psi_{0}$ because it is in general quite sensitive to the estimation of $\beta_{0}$ and the covariance parameters $\omega_0$. To orthogonalize this moment, we require the following assumptions. Recall the definition $M=\mathbb{E}\left[ W^{\prime}QW\right]$, let $\dot{\omega}_0$ be the Jacobian of $\omega_0$ and define the matrices $A:=vec(\Omega)'\mathbb{E}\left[ \mathbb{H}S_2\dot{\omega}_0\right]$, $B:=\mathbb{E}\left[ \mathbb{Q}S_2\right]$. Define also $\Gamma_{0\omega}:=AB^{\dagger}$, and $L:=-\mathbb{E}\left[ (vec(\Omega)'\mathbb{H}-\Gamma_{0\omega}\mathbb{Q})\{(W\otimes u(\beta_0))+(u(\beta_0)\otimes W)\}\right]$.

To implement this inference procedure, construct estimators: $$ \hat{A}_\ell = \frac{1}{n - n_\ell} \sum^L_{l \neq \ell} \sum_{i \in I_\ell} vec\left(\Omega\right)^{\prime}\mathbb{H}_iS_2\dot{\omega}_i,\;\;\;\; \hat{B}_\ell = \frac{1}{n - n_\ell} \sum^L_{l \neq \ell} \sum_{i \in I_\ell} S^{\prime}_2\mathbb{Q}_i,\;\;\; \hat{\Gamma}_{\omega \ell} = \hat{B}^{\dagger}_\ell \hat{A}_\ell, $$ where $\hat{B}^{\dagger}_\ell$ is a suitable estimator of the Moore-Penrose inverse of $B$ (e.g., spectral hard thresholding, as in Section (ref)).\footnote{In our empirical emplication below, we set the threshold to be $(\log(T^2)/n)^{1/2}$, which might be a conservative choice. Undertaking an optimal choice for the estimator of $B^{\dagger}$ is beyond the scope of this paper; see wainwright2019.} Moreover, let\footnote{Possibly, results can be improved by using double cross-fitted estimators $\hat{\Gamma}_{\omega \ell, \ell^{\prime}}$.} $$ \hat{L}_\ell = - \frac{1}{n - n_\ell} \sum^L_{l \neq \ell} \sum_{i \in I_\ell} (vec(\Omega)'\mathbb{H}_i-\hat{\Gamma}_{\omega \ell}\mathbb{Q}_i)\{(W_i\otimes u(\hat{\beta}_\ell))+(u(\hat{\beta}_\ell)\otimes W_i)\},\;\;\; \hat{\Gamma}_{\beta \ell} = \hat{L}_\ell\hat{M}^{\dagger}_\ell. $$ In addition, let

where $\hat{\omega}_\ell$ is an estimator of $\omega_0$, e.g., as in arellano2012identifying, using observations not in $I_\ell$. A feasible LR moment is then

An algorithm to implement LR inferences for $\psi_0 = \mathbb{E}\left[\alpha^{\prime}\Omega \alpha\right]$ follows in an analogous way as for average marginal effects. It is also possible to show the validity of an inference procedure based on (ref). We omit the details for the sake of space.

Kotlarski model with factor loading

We provide in this section a recursive construction of LR moments for $\psi _{0k}:=\mathbb{E}\left[ \alpha^{k}\right] ,$ $k\in\mathbb{N}$, in the Kotlarski model with a factor loading given by

where $\alpha$ is a common unobserved factor, and $\varepsilon = (\varepsilon_1, \varepsilon_2)$ has independent nonobservable components with zero mean and densities $f_{\varepsilon_1}$ and $f_{\varepsilon_2}$, respectively. The variables $\alpha$ and $\varepsilon$ are independent. Define for $k\in\mathbb{N},$ $\overrightarrow{\psi}_{0k}=(\psi_{01},...,\psi_{0k})^{\prime}\ $and the moments \[ a_{0k}(Z,\beta_{0},\overrightarrow{\psi}_{0k})=\beta_{0}\psi_{0k}-Y_{2} Y_{1}^{k-1}+\sum_{j=1}^{k-1}\left( \binom{k-1}{j-1} \beta_{0}Y_{1}^{k-j}-\binom{k-1}{j}Y_{2}Y_{1}^{k-j-1}\right) \psi_{0j}. \] The following lemma provides orthogonal moments for the multivariate functional $\overrightarrow{\psi }_{0k}$ if $\beta_{0}$ is known. Define \[ c_{k,j}:=\left\{

\right. . \]

The moment functions $a_{0k}$ are not LR when $\beta_0$ is unknown and the parameter of interest is just the $k-$th moment. We provide LR moments in a recursive way as follows. Define, for $k\in\mathbb{N},$

and

where the sum from 1 to zero is zero, \[ b_{k}=\psi_{0k}+(k-1)!\sum_{j=1}^{k-1}\frac{\mathbb{E}\left[ Y_{1}^{k-j}\right] }{j!(k-1-j)!}\psi_{0j}. \] and the constants $\gamma_{k,j}$ are defined recursively as, for $j=k-1,k-2,...,1$ \[ \gamma_{k,j}=-\frac{c_{k,j}+\sum_{h=j+1}^{k-1}\gamma_{k,h}c_{h,j}}{\beta_{0}}, \] and \[ \gamma_{k,0}=\frac{b_{k}+\sum_{h=1}^{k-1}\gamma_{k,h}b_{h}}{\mathbb{E}\left[ Y_{1}\right] }. \] For example,

Note that

Then, we have the following proposition. Define $\rho_{0k}:=(\beta_{0},\overrightarrow{\psi}_{0k-1},\Gamma_0)$ and write $g_{k}(Z,\beta_{0},\overrightarrow{\psi}_{0k},\Gamma_0)=g_{k}(Z,\rho_{0k},\psi_{0k})$.

The moment function $g_{k}$ satisfies the conditions of Theorem (ref) with $c=-\beta_{0}\neq0$, and hence, it is a relevant orthogonal moment. This construction avoids the nonparametric estimation of the distribution of latent variables, which are known to have very slow (logarithmic) rates of convergence, see, e.g., horowitz1996semiparametric.

Teacher Value-Added

In this application we follow gilraine2020new and consider the model for estimated teacher fixed effects, given by $Y=\alpha+\varepsilon$, $\varepsilon=\theta_{0}u$, $u\sim N(0,1),$ $\alpha$ independent of $u,$ with $\alpha\sim\eta_{0}$ for an unknown Lebesgue density $\eta _{0}$. This model with auxiliary moment conditions for $\theta_{0}$ is the model used in gilraine2020new. In its most general version without auxiliary moments for $\theta_{0}$, this model has been extensively studied in the literature, see, e.g., carroll1988optimal, fan1991optimal, and ishwaran1999information. In particular, it is well known that, without auxiliary moments and further assumptions, the model is not identified, see kiefer1956consistency. In any case, we are not aware of any investigation on LR moments in this model, which is the goal of our paper.

To that end, a key condition for the existence of a relevant orthogonal moment for the functional $\psi _{0}=\mathbb{E}\left[r(\alpha ,\theta _{0})\right] $ is the existence of a $g_{0}\in L_{2},$ such that, a.s.

where $\theta _{0}$ is the unknown error's standard deviation. In its integral form, this equation is

where, henceforth, we use the short notation $\phi _{\sigma }(x)=\phi \left( x/\sigma \right) /\sigma$ and recall $\phi\left( \cdot\right)$ denotes the standard normal density. We consider two types of functionals in this application: (i) analytic functionals, as defined below; and (ii) CDFs, quantiles and policy functionals.

Analytic Functionals

Analytic functionals correspond to moments $r(\alpha ,\theta _{0})$ that are an entire function of $\alpha$, that is, a real analytic function in $\alpha$ with a Maclaurin series expansion $r(\alpha )=\sum_{j=0}^{\infty }r_{j}\alpha^{j}$, where we drop the dependence on $\theta _{0}$ in the coefficients $r_j$ for simplicity of notation. This class includes moments and other smooth functionals of teacher value-added. To implement our approach, we first solve the equation for moments of $\alpha$. To that end, let us introduce the Hermite polynomials $\{H_{j}(\cdot)\}_{j=0}^{\infty }.$\footnote{ The $j$-th Hermite polynomial is defined by $H_{j}(x)=(-1)^{j}\exp (-x^{2}/2)\partial ^{j}\exp (-x^{2}/2)/\partial x^{j}.$} Then, the convolution identity for Hermite polynomials, see, e.g., davis2024general, yields

and thus, $a_{0k}(Y,\theta _{0})=\theta_0^{k}H_{k}(Y/\theta _{0})$ is a solution of Step 1 in the direct approach for the $k-$th moment. Exploiting the Maclauring series, we conjecture that

is a solution to ((ref)) in the general analytic case, provided this solution is in $L_2$. Our next result gives primitive conditions for this and gives an LR moment for the functional. Recall, that after solving for $g_0$, the LR moment will be given by \[ g(Y,\theta_{0},\psi_{0})=g_{0}(Y,\theta_{0})-\psi_{0}-\Gamma_{0}g_1(Y,\theta_{0}), \] where $\Gamma_{0}$ solves the equation \[ \mathbb{E}\left[ \frac{\partial g_{0}(Y,\theta_{0},\psi_{0})}{\partial\theta }\right]=\Gamma_{0}\mathbb{E}\left[ \frac{\partial g_{1}(Y,\theta_{0} )}{\partial\theta}\right], \] and $g_{1}$ is the auxiliary estimating equation for $\theta_0$ with finite variance and $\mathbb{E}\left[ \frac{\partial g_{1}(Y,\theta_{0} )}{\partial\theta}\right]\neq0$.\footnote{Note that we can let $\rho_{01} = \theta_0$ and $\rho_{02} = \Gamma_0$.} We also assume that the derivative $\mathbb{E}\left[ \frac{\partial g_{0}(Y,\theta_{0},\psi_{0})}{\partial\theta }\right]$ is well defined (this derivative can be obtained from well-known expressions of derivatives of Hermite polynomials and the derivatives of the coefficients of $r$). The following assumption is also required

for some $R$ such that $0<R<1$.

CDFs, quantiles and policy functionals

The literature of teacher value-added has long been interested in policy parameters such as the marginal gain in test scores resulting from replacing the bottom $\varphi$% (typically, 5%) of teachers. More generally, CDFs and quantiles are common parameters employed in the literature of models with nonparametric UH. In this section, we show that, as an important application of the necessity of our conditions, relevant orthogonal moments do not exist for these important parameters in teacher value-added models. This raises the question of the suitability of these policy parameters if the goal is to derive robust policy implications.

The following proposition shows that if $r$ satisfies equation ((ref)), then it has to be an analytic function of $\alpha$. An implication of this is that CDFs, quantiles, and policy parameters with non-analytic Riesz representers with respect to $\eta$ cannot satisfy equation ((ref)), and thus do not have a relevant orthogonal moment.

This proposition also motivates the approach of the previous section.

Let us consider some examples, and their corresponding Riesz representer $r$. This Riesz representer is defined in general as

For the CDF at $\bar{\alpha}$, the corresponding Riesz representer is $r(\alpha)=\bm{1}(\alpha \leq \bar{\alpha})$. For $\varphi$-quantile it is well-known that $r(\alpha)=\bm{1}(\alpha \leq F^{-1}_\alpha(\varphi))/\eta_{0}(F^{-1}_\alpha(\varphi))$. For the policy parameter $\psi(\lambda_0) = -\mathbb{E}\left[\alpha \bm{1}\left\{ \alpha < F^{-1}_\alpha(\varphi)\right\}\right]$, by the chain rule and the fundamental theorem of calculus, $r(\alpha)$ $=(F^{-1}_\alpha(\varphi)-\alpha)\bm{1}(\alpha \leq F^{-1}_\alpha(\varphi))$. In all these cases, the corresponding $r$ is not analytic. Therefore, by our results, there does not exist an LR moment for these functionals.

Monte Carlo

The aim of this section is to illustrate that constructing relevant orthogonal moments is important. In particular, we present a numerical exercise that shows that those moments might lead to substantial gains in terms of efficiency, although it might involve estimating UH.

Consider a simpler version of the Kotlarski model of Section (ref). Specifically,

where $\alpha$ is independent of $\varepsilon= (\varepsilon_1, \varepsilon_2)$, which has independent components, and $\beta_0 $ is unknown (with a true value of one). In this section, for a better comparison with existing procedures, we consider the case where the marginal densities $f_{\varepsilon_1}$ and $f_{\varepsilon_2}$ are known. Particularly, their distributions are standard normals. However, $\eta_0 = f_{\alpha}$ is unknown. The UH $\alpha$ is a discrete approximation to a standard normal with unknown support points; see below for the specification of $\eta_0$. We are interested in conducting LR inference on $\psi(\lambda_0) = \beta_0$, which is a special functional, e.g., it does not depend on the UH. This setting, with continuous $\alpha$, has been used by bonhomme2012functional (cf. Section 6) to evaluate the numerical performance of functional differencing, which we will exploit to have a benchmark for comparison. Functional differencing constructs moments orthogonal to UH (under a support condition), but does not necessarily deliver relevant moments.\footnote{This is a point acknowledged by Bonhomme; see Section 4.3 in bonhomme2012functional.} Our construction algorithms allow us to purposely compute informative moments. Indeed, in this exercise, we can construct the most efficient of them all, as we show below.

Notice that in this case, $$ f_{\beta_0,\eta_0}(y_1,y_2) = \int f_{z|\alpha}(z|\alpha;\beta_0)\eta_0 (\alpha)d\nu(\alpha), $$ with $$ f_{z|\alpha}(z|\alpha;\beta_0) = \phi\left(y_1 - \alpha\right)\times \phi\left(y_2 -\beta_0\alpha\right), $$ where $\phi$ is the standard normal pdf and $\nu$ is the counting measure.

Since our target functional is $\psi(\lambda_{0}) = \beta_0 :=\theta_0$ and finite-dimensional, we can apply Remark (ref). We aim to find a function $g$ such that\footnote{Observe that our second equation, $\mathbb{E}\left[s_\theta(Z,\lambda_0)g\left(Z,\lambda_0\right)\right] = c$, with $c \in \mathbb{R}$, holds in this case, provided that $s_\theta$ is in $L_2$.}

Notice that we are writing $g$ as a function of $\lambda_0$, instead of $\theta_0$ only. Our algorithm, explained in Section (ref), suggests that is not necessary. The reason why we want $g$ to depend on $\lambda_0$ will be come clear below. To solve for $g$ in (ref), we implement a projection approach, where $g=\Pi_{\mathcal{N}(S_{\eta_{0}}^{*})}m$, $$ S_{\eta}^{*}g = \mathbb{E}\left[ \left. g(Z, \lambda_0)\right\vert \alpha\right], $$ and $m$ is a given moment function, $m \in L^0_2$. A convenient candidate is $$ m\left(Y_1,Y_2,\lambda_0\right) = \frac{\partial}{\partial \beta} \log f_{\beta_0,\eta_0}(Y_1,Y_2), $$ i.e., the score of the model with respect to $\beta$. This is an appealing choice since the resulting relevant LR $g$ will correspond to the efficient (semiparametric) score of $\beta_0$, which could derive in efficiency gains when inference is based on our test statistic; see vandervaart98, Section 25.4, and lee2024locallyregularefficienttests. In Section (ref) of the Supplementary Appendix, we provide details regarding how to approximate $g=\Pi_{\mathcal{N}(S_{\eta}^{*})}m$ in practice.

In our simulations, the distribution of $\alpha$ is specified as follows

where $\Phi(.)$ is the standard normal distribution. Notice that (ref) indicates that the true density of individual effects is discrete but it closely behaves as a standard normal; Section (ref) of the Supplementary Appendix discusses regularized estimation of (ref). All results are based on $1,000$ Monte Carlo repetitions.

We are interested in comparing the performance of our test statistic and the naive plug-in one, which is based on $m(Z)$ without correcting for the first stage bias. To have a precise notion of the impact of the bias for inference, both tests are based on the same estimator of the asymptotic variance $\hat{W}$. Additionally, we consider a Wald-test statistic using the functional differencing approach of bonhomme2012functional. Specifically, we start with the function $\tilde{h}(y) = \varphi_2\left(y - \varrho\right)$ where $\varphi_2$ is the standard bivariate normal p.d.f. and $\varrho$ belongs to one of three increasing sets containing 9, 25, and 49 points, respectively.\footnote{These sets are $\left\{(0,0), (0,1), (0,-1),\cdots,(-1,-1)\right\}$, $\left\{(0,0), (0,1), (0,-1), (0,2), (0,-2)\cdots,(-2,-2)\right\}$, $\left\{(0,0), (0,1), (0,-1),(0,2), (0,-2), (0,3), (0,-3)\cdots,(-3,-3)\right\}$. } We can then obtain 9, 25, and 49 different functions $g_{FD}(y,\beta_0)$ by functional differencing. In fact, these are the same functions that bonhomme2012functional considers in his simulations.\footnote{bonhomme2012functional's code is publicly available, which we have closely followed to reproduce the results here.} Then, we estimate $\beta_0$ using GMM based on $g_{FD}(y,\beta_0)$, and compute Wald-type statistics $W_{9}$, $W_{25}$ and $W_{49}$, where $W_j = (\hat{\beta}-\beta_0)^2/\widehat{AsyVar}(\hat{\beta})$, $j=9,25,49$, respectively.

Table (ref) compares rejection rates under the null hypothesis $H_0: \beta_0 =1$ of all the test statistics considered and different sample sizes. The most important point from this table is the importance of controlling for regularizarion/first stage bias to conduct inference on parameters of interest. In particular, we observe that failing to construct LR moments results in test statistics with rejection rates sufficiently above the nominal values to be a serious concern in applied settings. In contrast, our construction of debiased moments yields a test statistic that presents controlled sizes figures, even though we based inference in regularized estimators $\hat{\eta}_\ell$, as our theory suggests. Functional differencing presents good performance in terms of size, uniformly in the number of moments constructed by such an approach. This is unsurprising since the Wald statistics are not affected by $\hat{\eta}_\ell$.

To evaluate power figures we compute (size-adjusted) rejection rates under several local alternatives. To be specific, we consider equally spaced points for true $\beta_0$'s such that $\beta_0 = \bar{\beta} + \delta$, with $\delta \in \left\{-0.5,-0.475,\cdots,0.425,0.450\right\}$ and $\bar{\beta}=1$. In this case, we only show the performance of our LR test statistics and the Wald-type ones. All them present larger rejection rates as the departure from the null becomes more important, although when a small number of moments is considered and for $n=250$, this behavior is not necessarily monotonous. More importantly, our test statistic is deemed more efficient as it presents power figures larger than the Wald-statistics. This is the gain that our relevant LR moment delivers in this case. It is not surprising since our test, in this case, is based on an estimator of the efficient score for $\beta_0$, which is the optimal LR moment. Overall, our numerical experimentation leads us to conclude that the finite sample performance of our test, both in terms of size and power, is in line with our theoretical results, and shows the advantage of correcting for the first-stage bias and using relevant LR moments.

Application: The Effect of Smoking on Birth Weight

We consider the following model:

where $i$ and $t$ index mothers and children, respectively, $y_{it}$ measures the child $t$'s birth weight of mother $i$, $s_{it}$ is equal to 1 if mother $i$ was smoking when she was pregnant with child $t$ and 0 otherwise, and $w_{it}$ is a vector collecting other determinants of weight at birth, including the sex of the child and the age of the mother. Notice that we can write (ref) in our linear random coefficient model notation as $$ Y_i = W_{i}\beta_0 + V_{i}\alpha_i + \varepsilon_i, \;\;\; i=1,\cdots,n, $$ where $Y_i = \left(y_{i1},\cdots,y_{iT}\right)^{\prime}$, $W_{i} = \left(w_{i1}, \cdots,w_{iT}\right)^{\prime}$, $V_i = \left(\imath_{T},\left(s_{i1},\cdots,s_{iT}\right)^{\prime}\right)$, where $\imath_T$ is a vector of ones with dimension $T=3$, and $\alpha_i = \left(\alpha_{i1}, \alpha_{i2}\right)^{\prime}$ varies across mothers ($q=2$). The intercept $\alpha_{i1}$ measures mother-specific unobserved characteristics (which could be partly genetic) that affect the weight of the child at birth whereas $\alpha_{i2}$ is the marginal effect of smoking.

abrevaya2006estimating, assuming homogeneous $\alpha_{2i}$, finds large and significant negative effects of smoking on birth weight. This coincides with the results in arellano2012identifying (AB, hereafter), who allow for UH across mothers. In both cases, the authors work with a restrictive specification for the additional controls in (ref). Differently, taking advantage of our theory for high-dimensional $W_i$, we will allow for a richer specification in covariates, addressing potential bias due to functional form assumptions on how controls enter in the model. Nevertheless, our results below generally confirm the findings obtained by the aforementioned works, suggesting a positive result: the existing empirical conclusions on the treatment effect of smoking on child birth weight are robust to the inclusion of a flexible model in covariates.

Some comments regarding the data are in order. We use the same sample of mothers from abrevaya2006estimating employed by AB. This data comes from the Natality Data Sets for the U.S. for the years 1990 and 1998. Abrebaya implements a method to match mothers to children as there is no unique identifiers in the original data. We use the “matched panel $\#$3", which is the subsample less likely to suffer from large matching error. The restriction $rank(V_i) = q < T$ makes us focus on mothers who had at least three children during 1989-1998. In our data, all mothers have had three children during the period of analysis. In addition, we need the smoking status $s_{it}$ to be time-varying. Hence, we keep only mothers for whom $s_{it}$ changes at least once. Our final sample cointains 1,445 mothers.

In the original specification of AB, $w_{it}$ includes the following variables: the sex of the child, the age of the mother at time of birth, dummy variables indicating number of conducted prenatal visits, and the value of an index---the “Kessner" index, which measures the quality of prenatal care abrevaya2006estimating. In this empirical exercise, we want to evaluate the robustness of this specification while still allowing for heterogeneous treatment. Specifically, we control for the aforementioned covariates along with a (cubic) B-spline basis expansion of five terms in age.\footnote{Our approach allows us to go beyond this specification. We, for instance, also considered a larger number of terms in the basis construction and interaction terms with the rest of the covariates in $w_{it}$; yet, some point estimates became reasonably implausible.} Our algorithm above will select the relevant terms in $w_{it}$.

We first conduct inference on some common parameters that appear in $\beta_0$, following the algorithm outlined in Section (ref). The results are shown in the first panel of Table (ref). Columns “Est. (AB)" and “CI$_{95}$ (AB)" display point estimates and $95\%$ confidence intervals by AB, respectively. Columns “Est. (AE)" and “CI$_{95}$ (AE)" show our own results. Our strategy, despite allowing for more flexibility, generally confirms the main conclusion obtained by AB: the sex of the child seems to be the only important control. In addition, both slope estimates for this variable are virtually equal: Boys tend to have, on average, 130 grams larger birth weights than girls. Although the rest of the regressors show the expected signs, they do not seem to be statistically significant.

Turning to heterogeneous treatment effects, we have estimated the average UH in the population of mothers with three children. Our results can be seen in the rows “Mean $\alpha_{1i}$" and “Mean $\alpha_{2i}$", indicating the average of the random intercepts and slopes, respectively. For this, we have followed the algorithm in Section (ref). While our results virtually coincide with those in AB for the average treatment effect of smoking, our intercept coefficient differs. In particular, we estimate an average effect of -161.65 grams of smoking during pregnancy. This effect is statistically different from zero, and our $95\%$ confidence intervals are almost the same as the one implied by the results obtained by AB. Our intercept coefficient is $22\%$ larger and equal to 3,404 grams (and statistically different from zero). In addition, our confidence intervals do not fully overlap with those in AB and suggest larger estimates for the random intercept.

Finally, we have estimated second moments for both UH parameters, following the guidelines in Section (ref). From these calculations and our previous estimates, we have computed implied variance estimates, which appear in the rows “Var $\alpha_{1i}$" and “Var $\alpha_{2i}$". We first assume iid homoskedastic errors.\footnote{In this case, $S_2 = vec(I_3)$ in (ref) and $\omega_{0i} = \sigma^2_0$, where $\sigma^2_{0}$ is the variance of the idiosyncratic errors.} Then, we allow for a more flexible variance-covariance matrix for the errors, following the same specification that Arellano and Bonhomme considered. In this case, the variance of errors for the $jth$ child is $a+b(j-1)$, where $a$ and $b$ are scalars.\footnote{Notice that the selection matrix $S_2$ must also change.} We refer to this model as the non-stationary errors case. In general, our results suggest significant dispersion,\footnote{Except for Var $\alpha_{1i}$ when iid errors are considered.} which is typically larger than that determined by AB. The most noticeable difference is observed for the variance (or standard deviation) of $\alpha_{2i}$ in the non-stationary errors case---our point estimate of its standard deviation is 18% larger. Interestingly, while the estimates by AB indicate a reduction in both variances, our method implies an increase in the dispersion of the random intercept and slope.\footnote{This responds to the fact that the two approaches are based on different moments to estimate the variances.} All in all, similarly to AB, our results suggest the presence of UH in the treatment effects of smoking on birth weights for this data; however, we find evidence of slightly more variability across mothers than previously determined.

Concluding Remarks

In the recent decade, there has been an increased interest in orthogonal moments and DML inferences, due to their convenient properties in contexts with high-dimensional first-step estimators. Yet, their extensions to models with UH have remained limited. In this context, several important parameters are functionals of common components and UH. Typically, the distribution of UH is not identified, and if regularized plug-in approaches are used, this might result in slow rates for the target functionals and large biases that invalidate inferences. All of this justifies the importance of applying Neyman-orthogonality to these settings, which is the main goal of this paper.

We have shed light on DML for models with UH and general parameters of interest by providing necessary and sufficient conditions for the existence and relevance of Neyman-orthogonal moments. Crucially, our conditions are constructive for identification, estimation, and inference. We have obtained novel identification results and DML inferences for high-dimensional random coefficient models, the Kotlarski model with a factor loading, and teacher value-added models. We have introduced an algorithm to construct relevant LR moments for smooth functionals, which might not necessarily require estimating UH, as illustrated with our examples. We hope these ideas will be applied in different contexts. Our paper opens the door for the application of valid asymptotic inferences in models with UH and other parameters estimated by machine learning methods---including high-dimensional regressors---and efficient inferences, as illustrated by our Monte Carlo and empirical exercises.

There are important avenues for future research on this subject. On an applied note, we strongly believe our results are useful to a whole host of applications, which were not covered by us in this paper. These might include applications to other random coefficient models, production functions with heterogeneous technologies, duration models, individual-specific income profile models, and heterogeneous treatment effects. It would be interesting to apply our theory to these contexts to obtain reliable conclusions in relevant empirical settings, allowing for flexible specifications and UH. Furthermore, on a theoretical side, we have shown that some important functionals, such as CDFs, counterfactual measures, and quantiles of UH, fail to have relevant orthogonal moments, as defined in this paper. In the future, we plan to explore the extensions of our results to such parameters, further extending the applicability of our results.

\phantomsection