Automatic Debiased Estimation with Machine Learning-Generated Regressors

Automatic Locally Robust GMM with Machine-Learning-Generated Regressors

Introduction

Many parameters of interest depend on predicted or generated regressors. Leading examples include structural parameters in models with endogenous variables estimated by control functions stock1989nonparametric,stock1991nonparametric,blundell2004endogeneity,imbens2009identification, average partial effects in sample selection models ahn1993semiparametric,das2003nonparametric,newey2009two, propensity score matching heckman1998matching,abadie2006large, and marginal treatment effects heckman2005structural. More recently, machine learning (ML) is routinely used to generate regressors for imputing missing covariates fongtyler, dimension reduction sorzano2014survey, learned proxies, confounders, and treatments knox2022testing, and feature engineering with unstructured data such as text, images, or audio feder2022causal, among many others.\footnote{knox2022testing estimate that about two thirds of recent computational work in political science uses predictions of unobserved concepts as regressors in their analyses.} In all of these settings, the parameter of interest depends on a regressor that is itself estimated in a preliminary step, often by flexible or high-dimensional methods.

Despite the prevalence of such problems in modern empirical work, there is currently no general inference framework that remains valid when regressors are generated by flexible or high-dimensional ML methods and then used again in downstream estimation. A common practice is to treat the generated regressors as if they were known, and to apply standard Generalized Method of Moments (GMM) or double/debiased machine-learning (DML) methods as if one were in a two-step setting. This practice typically yields invalid inference: the influence functions and asymptotic variances of such plug-in estimators have complicated analytic forms hahn2013asymptotic, and ignoring first-step estimation (i.e., the estimation of the generated regressor) generally leads to distorted standard errors and large regularization or model-selection bias in the final estimates.\footnote{For other biases induced by the repeated use of ML-generated data, see shumailov2024ai.} These difficulties are exacerbated in ML settings, where preliminary estimators are high-dimensional, nonparametric, and often non-Donsker chernozhukov2018double.

This paper develops automatic locally robust/debiased estimation and inference for structural parameters in three-step models with ML-generated regressors, generalizing the two-step setting of chernozhukov2022locally. A new idea is downstream local robustness: valid inference must neutralize not only the direct impact of estimating generated regressors, but also the indirect downstream effects that arise because these regressors are themselves inputs to later nuisance functions. Indirect effects are annihilated by making the moment robust to the second step. This paper shows how to automatically achieve downstream local robustness.

A simple toy example may help to fix ideas. Let $V(g)$ denote a ML-generated regressor produced by a first-step $g$ with true value $g_0$, and let the second-step nuisance $h(g)$ denote the optimal linear predictor of $Y$ on $V(g)$, with slope coefficient $\beta_g$, i.e., $h(g)(v)=\beta_g v$. The direct effect of the first step is related to the mapping $g\mapsto \beta_0 V(g)$, where $\beta_0=\beta_{g_0}$, while the indirect effect operates through the second step via the mapping $g\mapsto \beta_g v$. The indirect effect is more complex than the direct effect, as can be seen from $\beta_g=\mathbb{E}[YV(g)]/\mathbb{E}[V(g)^2]$. Let the parameter of interest be a functional of $(g,h)$, say $\theta_0=\theta(g_0,h_0)$ with $h_0=h(g_0)$. Downstream local robustness means that, by the functional chain rule, if $\partial\theta/\partial h (g_0,h_0)=0$, then \[ \left. \frac{\partial}{\partial g} \theta(g_0,h(g)) \right|_{g=g_0} =\frac{\partial\theta}{\partial h}(g_0,h_0)\cdot \frac{\partial h}{\partial g}(g_0) =0, \] so orthogonality with respect to the second step removes the indirect effect of $g$. The direct effect $\partial\theta/\partial g(g_0,h_0)$ still remains. This paper provides an automatic construction of functionals that delivers zero derivatives for both direct and indirect effects.

Automatic debiased estimators with generated regressors are useful for two main reasons. First, debiased estimators deliver downstream local robustness and correct for the large regularization and model-selection biases that arise when ML-generated regressors are plugged into subsequent stages. In our simulations, with a moderately large sample ($n=1000$), they reduce the bias of the DML estimator by up to 95%. Second, in three-step procedures, the analytic form of influence functions and asymptotic variances becomes complex and hard to derive (cf. hahn2013asymptotic). Our estimators and tests are automatic in the sense that these objects are estimated directly from data and identifying moments, without requiring analytic derivations or bootstrap approximations whose theoretical justification is delicate in the presence of ML-generated regressors.

A key feature of these problems is that the use of generated regressors induces a natural three-step structure. We therefore generalize the existing debiasing literature from a two-step to our three-step framework. In the first step, some regressors are predicted (for example, via imputation, ML-estimated propensity scores, or control functions with high-dimensional covariates). In the second step, a nuisance function is constructed as a (potentially high-dimensional) least-squares projection using the generated regressors and possibly other covariates. In the third step, the parameter of interest is identified by a GMM criterion involving the first two steps and the data. Existing debiasing methods could be applied by treating either the generated regressor or the second-step nuisance as known, effectively reducing the problem to two steps, but this generally leads to invalid inference. Additionally, the three-step structure induces a constrained, non-product parameter space in which the second-step nuisance depends on first-step generated regressors, thereby invalidating standard local-robustness arguments that rely on product-space perturbations in the two-step literature (see Remark (ref)).

We now summarize our main contributions.

First, we develop a general three-step locally robust GMM framework for models with generated regressors. We fully and separately account for the first and second steps and characterize their contributions to the parameter's influence function, including an indirect effect of the first step that operates through the generated regressors as conditioning variables in the second step. We show that when the second-step effect is zero, this indirect effect is also zero, extending a remark in hahn2013asymptotic to a more general class of three-step procedures that include leading ML methods. This establishes the broader orthogonality principle of downstream local robustness: by the functional chain rule, moment functions constructed to be orthogonal to the second step eliminate the indirect effects of generated regressors (see Proposition (ref)).

Second, we provide automatic estimation of influence functions and asymptotic variances for models with generated regressors. Under a linearization assumption newey1994asymptotic,ichimura2022influence, we show how the Riesz representers in the first- and second-step influence functions can be identified and estimated separately without knowing their analytic form. This is achieved via cross-fitted auxiliary regressions that remain valid for generic non-Donsker ML methods in the first and second steps; see, e.g., chernozhukov2021automatic,chernozhukov2022locally,chernozhukov2023automatic,chernozhukov2022automatic. Automatic estimation is particularly well motivated for generated regressors, where the Riesz representers typically have complex forms hahn2013asymptotic,mammen2016semiparametric,escanciano2014uniform. Together with the first contribution, we establish feasible standard errors and valid asymptotically normal inference for debiased estimators with ML-generated regressors. Relative to the DML literature, the presence of generated regressors makes the asymptotic analysis---and, in particular, the control of higher-order terms in the asymptotic expansions---more delicate, and we address this issue.

Third, we propose novel automatic three-step debiased estimators for leading applications such as high-dimensional propensity score (Hd-PS) regression adjustment, treatment effects with learned confounders (autoencoders), nonparametric Average Treatment Effect (ATE) estimation on a boosted propensity score, and the nonparametric Counterfactual Average Structural Function (CASF). In these settings, the generated regressors arise, for example, from a control-function approach using Lasso, Random Forest, or Deep Learning; from Logit-Lasso or Boosting Hd-PS; or from learned confounders via autoencoders bengio2013representation. The nonparametric ATE estimator with a Hd-PS generalizes heckman1998matching, hahn2013asymptotic, and mammen2016semiparametric to a ML setup with debiasing and automatic inference, reducing regularization bias from both first and second steps. The application to the CASF with a control-function approach appears to be novel even in low dimensions, and it is related to the literature on domain adaptation, transfer learning, and covariate shift. Relative to that literature, we allow for endogeneity and a flexible non-separable structural model, which is important in applications where counterfactuals involve endogenous variables such as prices.

Our work builds on two strands of the literature. The first is the classical literature on semiparametric estimation with generated regressors ichimura1991semiparametric,ahn1993semiparametric,heckman1998matching,newey1999nonparametric,li2002semiparametric,rothe2009semiparametric,imbens2009identification,escanciano2010testing,song2012smoothness. In an important work, hahn2013asymptotic derive the influence function of three-step estimators that are averages of evaluation functionals of nonparametric regressions with generated regressors. We build on these influence-function calculations by considering a more general class of first, second, and third steps, including high-dimensional regressions (e.g., Logit-Lasso) and targets that may depend on the entire second step (not only evaluation functionals, as in, the CASF example). For estimation, mammen2012nonparametric,mammen2016semiparametric and escanciano2014uniform study the asymptotic properties of (non--locally robust) estimators using empirical process methods. These existing results are formulated for nonparametric first and second steps in Donsker classes and are generally not applicable to ML estimators, which often fall outside Donsker classes chernozhukov2018double. We contribute to this literature by providing automatic debiased GMM estimators that explicitly account for ML-generated regressors and reduce regularization and model-selection biases, and by proving their asymptotic properties accounting for ML-generated regressors.

The second strand is the literature on locally robust/debiased estimators chernozhukov2018double,chernozhukov2022locally,chernozhukov2022automatic. With the exception of sasaki2021estimation, the DML literature prior to our work has not considered or accounted for generated regressors in inference. Our results complement sasaki2021estimation by providing a general three-step framework and automatic estimation of adjustment terms for a broad class of models with generated regressors, including empirically relevant settings such as the partially linear model with ML-generated regressors. Relative to the Automatic DML literature, we innovate by (i) working in a three-step setting where the second step depends on the generated regressor and the product-space structure of chernozhukov2022locally fails; (ii) exploiting novel partial and downstream local robustness results that allow separate identification and automatic estimation of individual Riesz representers; and (iii) accounting for generated regressors in the estimation of Riesz representers and the bounds for higher-order terms in functional derivatives with respect to the high-dimensional generated regressors.

The rest of the paper is organized as follows. Section (ref) introduces the setting and examples. Section (ref) describes the debiased moment functions and defines the debiased GMM estimator in the presence of ML-generated regressors, illustrating its performance in two Monte Carlo experiments. Section (ref) gives the separate identification and automatic estimation of the Riesz representers. Debiased automatic estimators for the examples are presented in Section (ref). The asymptotic theory is developed in Section (ref). Section (ref) concludes. Appendix (ref) summarizes the estimation steps. Appendix (ref) provides a further application to a ML implementation of the nonparametric ATE estimator of heckman1998matching. Appendix (ref) contains details about the Monte Carlo simulations. Appendix (ref) discusses regularity conditions, and Appendix (ref) gathers the proofs of the main results.

Setting and examples

Three-step setting

We observe data $W=(Y,D,Z)$ from a cumulative distribution function (cdf) $F_{0}$. We describe our three-step setting as follows:

First step. There is a first-step nuisance function $g_{0}(Z)$ satisfying the moment restrictions

where $\epsilon (W,g_{0})$ is a generalized error depending on the data $W$ and the nuisance $g_{0}\in \Delta _{1}$, where $\Delta _{1}$ is a linear and closed subspace of $L_{2}(Z)$. Henceforth, for a generic random variable $U$, we denote by $L_{2}(U)$ the Hilbert space of square-integrable functions of $U$, i.e., $g\in L_{2}(U)$ iff $\mathbb{E}[g^2(U)]<\infty$.

This setting covers a wide variety of semiparametric and nonparametric first steps. For example, when $\epsilon (W,g_{0})=D-g_{0}(Z)$ and $\Delta _{1}=L_{2}(Z)$, we have $g_{0}(Z)=\mathbb{E} [D|Z]$, as in hahn2013asymptotic. However, if $\operatorname{dim}(Z)$ is high, fully nonparametric first steps may not be feasible to implement. We could then consider a high-dimensional additive regression model with the same error but with $\Delta _{1}=\sum_{j =1}^{\operatorname{dim}(Z)}\Delta _{1,j}$, where $\Delta _{1,j}$ is a subset of $L_{2}(Z_{j})$ for the $j$-th component of $Z$ wainwright2019high. When $\Delta _{1}$ is the mean-square limit of linear combinations $\sum_{k=1}^{K}\beta _{0k}c_{k}(Z)$ for $K\in\mathbb{N}$, a sequence of real numbers $(\beta _{0k})_{k=1}^{\infty }$, a dictionary $(c_{k})_{k=1}^{\infty }$ of functions in $L_{2}(Z)$, and $\epsilon (W,g_{0})=D-\Lambda (g_{0}(Z))$ for the logistic cdf $\Lambda$, this setting covers high-dimensional logistic regression (Logit-Lasso), which is commonly used for propensity-score and classification modeling in high dimensions. These ML-generated regressors complement the fully nonparametric mean-regression first steps in hahn2013asymptotic. For numerous other examples of $\epsilon (W,g_{0})$, including quantile regression, see Section 3 of ichimura2022influence. For general parametric first steps, see Remark (ref); and for other first steps not covered by our setting, see Remark (ref).

The first-step nuisance $g_{0}$ in ((ref)) is used to construct the population generated regressors

where $\varphi$ is a known function of observed variables $(D,Z)$ and the unknown function $g_{0}.$ Note the simplified notation $V\equiv V(g_{0}).$ A high-dimensional extension of propensity score matching in heckman1998matching has $V=\Lambda (g_{0}(Z))$; some dimension-reduction methods have $\varphi (D,Z,g_{0})=g_{0}(Z)$, as in hahn2013asymptotic, or some components of $g_{0}$ (as with autoencoders); imputation for conditionally missing-at-random regressors has $\varphi (D,Z,g_{0})=Z_{1}D+(1-Z_{1})g_{0}(Z_{2})$, where $Z_{1}$ is a “not missing" indicator for the covariate $D$ and $g_{0}(Z_{2})=\mathbb{E} [D|Z_{1}=1,Z_{2}]$ for observed covariates $Z_{2}$; and control-function methods often lead to $\varphi (D,Z,g_{0})=D-g_{0}(Z)$, for an endogenous variable $D$ and exogenous variables $Z.$ Our setting covers these and other generated regressors.

Second step. Let $S$ and $X$ denote some components (or all) of $(Y,D)$ and $(D,Z)$, respectively. The second step links $S$ with $X$ and the generated regressor $V$ through the moment restrictions

where $\Delta _{2}(g_{0})$ is a linear, closed subspace of $L_{2}(X,V)$ (note $\Delta _{2}$ depends on $g_{0}$ because $V$ depends on $g_{0}$). When $S$ (and hence $h_{0}$) has dimension $\operatorname{dim}(S)>1$, we understand ((ref)) as being applied to each component of $S.$ This dependence of the parameter space $\Delta _{2}(g_{0})$ on the first step $g_{0}$ is a point of departure from the existing debiasing literature in, e.g., chernozhukov2022locally. hahn2013asymptotic and mammen2016semiparametric consider cases where the second step $h_{0}$ is a nonparametric regression of $Y$ on $(X,V)$, corresponding to $\Delta _{2}(g_{0})=L_{2}(X,V)$. In contrast, we also allow $\Delta _{2}(g_{0})$ to be a strict subset of $L_{2}(X,V)$ (e.g., with sparse or sieve restrictions).

Third step. Let $\Theta \subseteq \mathbb{R}^{p}$ denote the parameter space where the parameter of interest lies. Consider the moment function $m\colon \mathbb{R}^{\operatorname{dim}(W)}\times L_{2}(Z)\times L_{2}(X,V)^{\operatorname{dim}(S)}\times \Theta \rightarrow \mathbb{R}^{q}$, $q\geq p.$ The parameter of interest $\theta _{0}$ is identified in a third step by a GMM moment condition

Here we assume that $\theta _{0}$ is identified by these moments, i.e., that $\theta _{0}$ is the unique solution to $\mathbb{E}[m(W,g_{0},h_{0},\theta )]=0$ over $\theta \in \Theta $.

Examples

The following examples are used to illustrate the main results of this paper.

Debiased estimation with ML-generated regressors

A fundamental property that allows us to develop debiased estimators is Neyman orthogonality, also referred to as local robustness neyman1959optimal,chernozhukov2018double,chernozhukov2022locally. In our three-step setting, Neyman-orthogonal moments are obtained by augmenting the original identifying moments with influence-function (IF) corrections associated with the first ($g_0$) and second ($h_0$) steps. The second-step IF accounts for the effect of estimating the second-step nuisance $h_0$ and corresponds to the classical correction in newey1994asymptotic.

A key difference relative to standard two-step problems is that, here, the first-step nuisance $g_0$ enters the moment condition in two ways: directly through $m(W,g_0,h_0,\theta)$ and indirectly through the fact that estimation of $h_0$ depends on the generated regressor $V=\varphi(D,Z,g_0)$. Thus, estimation error in $g_0$ affects the target parameter through a direct (or evaluation) effect (in the toy example, $\beta_0 v$ evaluated at $v=V(g)$) and an indirect (or conditioning) effect that operates through $h_0$ (in the toy example, $g\mapsto h(g)(v)=\beta_g v$); see Figure (ref) and Section (ref). This indirect effect is absent in standard two-step locally robust problems but is unavoidable whenever the conditioning variable in a regression is itself ML-generated.

We show that the debiased moment function takes the generic form

where $\alpha_0 \equiv (\alpha_{01}, \alpha_{02})$ are the Riesz representers associated with the first and second steps, respectively. In the case of multiple moment conditions ($q>1$), each component of $m$ is debiased separately.

The function $\phi_{1}$ in (ref) is the first-step IF and captures the effect of the generated regressors on the identifying moments. It is generally nonzero, so inference that ignores generated regressors is typically invalid.\footnote{One instance where $\phi _{1}=0$ and inference that does not account for generated regressors is valid is when the sample size used to construct the generated regressors is asymptotically larger than the sample size used to estimate the main parameter (see Remark (ref) for a formal statement).} The explicit analytic expression for $\alpha_{01}$ is typically complicated (see equation (ref) in Appendix (ref)), but we construct automatic estimators that do not require this expression. The second-step IF $\phi _{2}$ is of the usual form newey1994asymptotic, but automatic estimation of the corresponding Riesz representer $\alpha_{02}$ must be generalized to allow for generated regressors as inputs.

A central insight of this paper is that the indirect contribution of $g_0$ operates entirely through the second-step. By the functional chain rule, this implies that, when the moment is orthogonal with respect to $h$ (so $\alpha_{02}=0$), the indirect effect of the generated regressor is zero and the first-step IF $\phi_1$ simplifies. We refer to this property as downstream local robustness. It generalizes an observation in hahn2013asymptotic to a general three-step ML framework and to a much broader class of problems beyond generated regressors (cf.\ Proposition (ref)).

The debiased estimator

Automatic debiased estimation with generated regressors is based on the moment condition in equation (ref), where the Riesz representers $\alpha_{01}$ and $\alpha_{02}$ are estimated automatically (see Section (ref)). We construct sample analogues using cross-fitting, as in chernozhukov2018double: the sample is split into $L$ folds $I_\ell$, and for each fold we evaluate $\psi(W_i,g_0,h_0,\alpha_0,\theta)$ only on observations $i\in I_\ell$ that were not used to estimate $(g_0,h_0,\alpha_0)$. Formally, we partition $(W_i)_{i=1}^n$ into $L$ groups $I_\ell$, for $\ell=1,\dots,L$. For each group, we have estimators $\hat{g}_\ell$, $\hat{h}_\ell$, and $\hat{\alpha}_\ell=(\hat{\alpha}_{1\ell},\hat{\alpha}_{2\ell})$ based only on observations outside $I_\ell$.

The debiased sample moment function is

with

for $\hat{V}_{i\ell }\equiv \varphi (D_{i},Z_{i},\hat{g}_{\ell })$. When there is more than one moment condition, each component of $m$ is debiased by its own Riesz representers, so as many $\hat\alpha_{\ell}$'s must be estimated as there are moment conditions.

The three-step debiased GMM estimator is then defined as

where $\hat{\Upsilon}$ is a positive semi-definite weighting matrix of dimension $q\times q$. Under regularity conditions (see Section (ref)), $\hat{\theta}$ is asymptotically normal with the usual GMM asymptotic variance.

Monte Carlo simulations

We give an overview of two Monte Carlo studies: estimation of Hd-PS regression adjustment in the partially linear model and estimation of the CASF with a control-function approach. We evaluate the finite-sample performance of several estimation procedures. First, the plug-in estimator that uses the original moment condition. For inference based on the plug-in estimator, we consider both accounting and not accounting for estimation effects in the asymptotic variance. Second, the natural application of the DML procedure of chernozhukov2018double,chernozhukov2022locally, which corrects only for the second step in parameter and asymptotic-variance estimation. Third, our proposed three-step debiased (3SD) estimator with an asymptotic-variance estimator (see Section (ref)). A detailed description of the setups, estimation procedures, and results is provided in Appendix (ref).

Hd-PS regression adjustment in the partially linear model

The available data are $(Y, D, Z)$, with $Z \equiv (Z_j)_{j=1}^{10}$. The outcome and treatment equations are:

The error terms $\varepsilon$ and $\nu$ are independent, with $\varepsilon \sim N(0, 1)$. The distribution of $\nu$ varies with the specification: it can be logistic or standard normal. The regressors $Z$ are independent of each other and are uniformly distributed on $[-1,1]$. The regressors are also independent of $\nu$. On the other hand, the first two regressors $(Z_1, Z_2)$ and $\varepsilon$ are correlated, rendering the treatment $D$ endogenous. The constant $C$ is chosen so that the propensity score is supported on $[0.01, 0.99]$. Here, $\theta_0 = 1$.

Results regarding the mean bias are similar across specifications. For a small sample ($n=100$), the DML estimator (the same as the plug-in estimator here) is heavily biased, while the three-step debiased estimator performs well (see Table (ref) in Appendix (ref)). When $n=100$, debiasing removes around 85% of the bias present in the DML estimator. As the sample size increases, the bias of the DML estimator becomes smaller. Nevertheless, the bias of the DML estimator remains orders of magnitude larger than that of the three-step debiased estimator. Figure (ref) displays histograms of both estimators for the logistic-$\nu$ specification. We see that, when $n=1000$, the DML estimator is still biased. The distribution of the three-step debiased estimator is centered around $\theta_0 = 1$.

Table (ref) in Appendix (ref) shows that the coverage of the three-step debiased estimator is close to the nominal 95%, even for $n=100$. On the other hand, the DML asymptotic-variance estimator tends to overestimate the true asymptotic variance. This leads to coverage rates that exceed the nominal level, except for the $n=100$ case, where the bias dominates. In addition, the plug-in estimator for $\theta_0$ (which equals the DML estimator) shows poor performance even when using the correct asymptotic variance for inference. Its coverage is below 90% even when $n=1000$.

CASF with a control-function approach

The available data are $(Y, D, Z)$, with $Z \equiv (Z_j)_{j=1}^6$. The variables $D$ and $Y$ are generated by:

The error terms $U$ and $V$ are correlated, with $U, V \sim N(0,1)$, so $D$ is endogenous. The regressors $Z$ are standard normal and are independent of each other and of the errors $(U,V)$. In this case, $X = (Z_1, \dots, Z_5, D)$. We estimate the CASF for the following counterfactual distribution $F^*$: (i) the distribution of $(Z_1,\dots, Z_5)$ remains unchanged and (ii) $D$ is normal with mean $1$ (instead of $0$) and the same variance as in the DGP. Therefore, the true parameter is $\theta_0=2$.

The plug-in estimator is severely biased across all sample sizes (see Table (ref) in Appendix (ref) and Figure (ref) below). However, the comparison between the DML and the three-step debiased estimator differs from that in the previous example. For small samples ($n=100$), both estimators have similar bias. As the sample size increases, the bias of the three-step debiased estimator decreases, while the bias of the DML estimator remains sizable. This confirms the presence of an asymptotic bias in the DML estimator.

Results regarding coverage also differ from those of the previous example (see Table (ref) in Appendix (ref)). The three-step debiased estimator shows good coverage, close to the nominal 95% level when $n=500$ or $n=1000$. In this case, the DML asymptotic-variance estimator underestimates the true asymptotic variance. Thus, its coverage is well below the nominal 95% level across all sample sizes. In estimating the CASF, the plug-in estimator performs poorly due to the large asymptotic bias. Even when using the correct asymptotic variance for inference, its coverage ranges from 59.4% when $n=100$ to 24.5% when $n=1000$.

Summarizing, not accounting for the generated regressor leads to large biases in finite samples. In contrast, our three-step debiased procedure substantially reduces bias and delivers robust and valid inference. The following sections show how the Riesz representers needed to build the debiased moment function are identified and estimated. These sections are more technical than the previous ones; thus, an applied reader may wish to jump directly to Section (ref), where additional details about the examples are gathered.

First- and second-step influence functions

This section provides a detailed construction of orthogonal moment functions in our three-step setting with generated regressors. We begin by introducing additional concepts and notation. Let $F$ denote a possible cdf for a data observation $W$. We denote by $g(F)$ the probability limit of an estimator $\hat{g}_\ell$ of the first step when the true distribution of $W$ is $F$, i.e., under general misspecification newey1994asymptotic. That is, $F$ is unrestricted except for regularity conditions such as existence of $g(F)$ and finiteness of the expectation of certain functions of the data. For example, if $\hat{g}_\ell(z)$ is a nonparametric estimator of $\mathbb{E}[D\mid Z=z]$, then $g(F)(z)=\mathbb{E} _{F}[D\mid Z=z]$ is the conditional expectation function when $F$ is the true distribution of $W$. We denote expectation under $F$ by $\mathbb{E}_{F}$, which is well defined under the regularity condition that $\mathbb{E}_{F}[|D|]$ is finite. We assume that $g(F)$ is identified as the solution in $g$ to

Our notation is consistent with $g(F_{0})=g_{0}$ being the probability limit of $\hat{g}$ when $F_{0}$ is the cdf of $W$.

To study the effect of the second step, suppose again that $W$ is distributed according to $F$, but the first-step nuisance is independently fixed to $g\in \Delta _{1}$. Let $h(F,g)$ be the solution in $h\in \Delta _{2}(g)$ to

where $V(g)\equiv \varphi (D,Z,g)$ and $\Delta _{2}(g)$ is a linear and closed subspace of $L_{2}(X,V(g))$ for each $g\in \Delta _{1}$. The solution of the above equation is a function of $(x,v)$, written as $h(F,g)(x,v)$. In the toy example, $h(F,g)(x,v)=\beta_g(F)v$, where $\beta_g(F)=\mathbb{E}_{F}[YV(g)]/\mathbb{E}_{F}[V(g)^2]$. We use the short notation $h_{0}(x,v)\equiv h(F_{0},g_{0})(x,v)$. Thus, henceforth, a subscript $0$ in $h$ means that the conditioning variable is the true generated regressor $V\equiv V(g_{0})$; for example, $h_{0}(x,v)=\mathbb{E}[Y\mid X=x,V=v]$ when $\Delta_{2}(g)=L_{2}(X,V(g))$. We may think of the mapping $h(F,g)$ as the probability limit of an estimator of $h_{0}$ under the following conditions: (i) the true distribution of $W$ is $F$ and (ii) the estimator is constructed with the first-step nuisance fixed at $g\in \Delta_{1}$. A feasible estimator $\hat{h}_\ell$ of $h_{0}$ will, however, rely on the estimator $\hat{g}_\ell$ with probability limit $g(F)$. Therefore, we assume that the probability limit of $\hat{h}_\ell$ under general misspecification is $h(F,g(F))$.

Let $H$ be some alternative distribution that is unrestricted except for regularity conditions, and define $F_{\tau }\equiv (1-\tau )F_{0}+\tau H$ for $\tau \in [0,1]$. We assume that $H$ is chosen so that $g(F_{\tau })$ and $h(F_{\tau },g(F_{\tau }))$ exist for sufficiently small $\tau$, and that other regularity conditions are satisfied. The effect of both first- and second-step estimation on the moment condition is measured by the derivative with respect to $\tau$ at $\tau=0$ of $\bar{m}(g(F_\tau), h(F_\tau, g(F_\tau)))$, with

We study these effects separately. By the chain rule,

where, henceforth, $d/d\tau$ denotes the right derivative with respect to $\tau$, evaluated at $\tau=0$. In the display above, the first derivative on the right-hand side (RHS) accounts for the first step. As in hahn2013asymptotic, the first step affects the moment condition in two ways (see Figure (ref)). We have a direct impact on $\bar{m}$, quantified by the derivative of $\bar{m}(g(F_{\tau }),h_{0})$. This direct impact includes the effect of evaluating $h$ at the generated regressor. We also have an indirect effect on the moment that arises because $g$ affects estimation of $h_{0}$ in the second step (through conditioning), quantified by the derivative of $\bar{m}(g_{0},h(F_{0},g_{\tau }))$. Both effects (direct and indirect) are considered in (ref). The derivative in (ref) accounts for the effect of the second step. This effect is independent of the first step and therefore treats $g_{0}$ as known.

To debias the moment conditions, we compute separate IFs for each estimation step. That is, we seek functions $\phi _{1}(w,g,\alpha _{1})$ and $\phi _{2}(w,g,h,\alpha _{2})$ such that, all $H$ defining a regular path $F_{\tau }\equiv (1-\tau )F_{0}+\tau H$,

Additionally, we require the IFs to have zero mean and finite variance. Note that $(\alpha_{01}, \alpha_{02})$ are the Riesz representers of the above derivatives, which are evaluated at $(g_0, h_0, \theta_0)$, and thus may depend on $(g_0, h_0, \theta_0)$. For functions $\phi_1$ and $\phi_2$ satisfying the above conditions, the moments

are orthogonal/locally robust/debiased. The next section shows that, under some conditions, the IFs have the form displayed in equation (ref). It also illustrates the separate automatic estimation of each Riesz representer $\alpha_{01}$ and $\alpha_{02}$.

It is worth highlighting that when the second-step effect is zero, the first-step indirect effect is zero by the chain rule. This is a special case of a more general result that applies beyond generated regressors.

Automatic estimation of the Riesz representers

The orthogonal moments require a consistent estimator $\hat{\alpha }_\ell$ of the Riesz representers $\alpha _{0}\equiv (\alpha _{01},\alpha _{02})$ . When the shape of $\alpha _{0}$ is known, one can plug-in nonparametric estimators of the unknown components of $\alpha _{0}$ to form $\hat{\alpha}_\ell$ . In the generated regressors setup, however, the nuisance parameters (especially $ \alpha _{01}$) have a complex analytical shape (see the result in equation (ref) in Appendix (ref)). Therefore, the plug-in estimators may be cumbersome to compute in practice.

To ease exposition and without loss of generality, in this section, we consider that there is a single moment condition ($p=q=1$). Recall that in the multi-dimensional case one must estimate Riesz representers $\alpha _{0}$ for each moment condition.

Separate identification of Riesz representers

We provide separate orthogonality conditions that will serve as a basis for the identification and automatic estimation of the Riesz representers $\alpha _{01}$ and $ \alpha _{02}$. Define the following moment functions: $\psi _{1}(w,g,\alpha _{1},\theta )\equiv m(w,g,h(F_{0},g),\theta )+\phi _{1}(w,g,\alpha _{1})$ for the first step, and $\psi _{2}(w,h,\alpha _{2},\theta )\equiv m(w,g_{0},h(F,g_0),\theta )+\phi _{2}(w,g_0,h(F,g_0),\alpha _{2})$ for the second step. Since, individually, the spaces $\Delta _{1}$ and $\Delta _{2}(g_{0})$ are linear, an application of Theorem 3 in chernozhukov2022locally to each step leads to

where $\delta_{1}$ represents a possible direction of deviation of $g(F)$ from $g_{0}$ and $\delta _{2}$ represents a possible deviation of $h(F,g_{0})$ from $h_{0}$. The innovation relative to chernozhukov2022locally is that we can compute the IFs $\phi_1$ and $\phi_2$ by separately studying $\psi _{1}$ and $\psi _{2}$, respectively. This means we can separately identify $\alpha _{01}$ and $\alpha _{02}$ from ((ref)) and ((ref)), even though $\psi _{1}$ and $\psi _{2}$ are not LR moment functions ($\psi _{1}$ and $\psi _{2}$ are not LR to $h_0$ and $g_0$, respectively).

Likewise, rather than joint identification from the analytical derivatives of the original identifying moments as in chernozhukov2022locally, which are not be available with generated regressors, we propose an approach that uses the linearization and orthogonality of $\psi _{1} $ and $\psi _{2}$ with respect to $g$ and $h$, respectively, to construct separate estimators of $\alpha _{01}$ and $\alpha _{02}$. This approach does not require knowing the shape of $\alpha _{0}$. It is \textquotedblleft automatic" in only requiring the orthogonal moment functions and data for the construction of $\hat{\alpha}_\ell$. Moreover, an automatic estimator can be constructed separately for each step.

The key ingredients for our approach are (i) the shape of the IFs and (ii) a consistent estimator of the linearization of the moment condition with respect to each parameter ---$g$ for the first step and $h$ for the second. Section (ref) provides the formal development. For a detailed construction of the automatic estimators, we refer to Section (ref).

First- and second-step linearization

We start with the linearization of the second-step effect because this will show up in the first-step linearization. The linearization of the second step with a known first step is a well-established result in the literature newey1994asymptotic, and it will follow immediately if $\bar{m}(g_{0},h)$ is linear in $h$.

Before introducing the result, we note that throughout this section, for $F_\tau\equiv (1-\tau)F_0+\tau H$, we consider that $ \tau\mapsto h_\tau\equiv h(F_\tau,g_0)$ and $\tau \mapsto g_\tau \equiv g(F_\tau)$ denote differentiable paths in $L_2(X,V)$ and $L_2(Z)$, respectively; i.e., $0\mapsto h_0$ and $ dh_\tau/d\tau$ exists (equivalently for $g_\tau$). When an assumption is stated for $h_0$ or $h_\tau$, it is understood that it applies to each of its components.

We assume that $\bar{m}$ can be linearized with respect to the second step parameter:

The same assumption has been considered in newey1994asymptotic. A necessary and sufficient condition for the linearity and continuity part is the existence of $r_{02} \in L_2(X,V)$ such that $\mathbb{E}[D_{02}(W,h)]=\mathbb{E}[r_{02}(X,V)h(X,V)]$ for all $h\in L_2(X,V)$. We can then get the shape of the second step IF:

An important observation is that if $r_{02}$ is zero, then $\alpha _{02}$ (and hence $\phi_{2}$) is also zero. We also note that $\bar{m}$ is linearized at $(g_{0},h_{0},\theta_0)$, so $D_{02}$, $r_{02}$, and $\alpha _{02}$ may also depend on $(g_{0},h_{0},\theta_0)$. This is omitted for notational simplicity, but it is of course accounted for in the theory of this paper, and it will become relevant to construct feasible automatic estimators (see Section (ref)).

We now move to the more complicated linearization of the first-step effect. Note that if the chain rule can be applied:

The first derivative in the RHS can be easily analyzed if we linearize $\bar{ m}(g,h_{0})$ in $g$:

The term $D_{dir}$ is responsible for the direct effect of the first step (the evaluation effect). Again, for simplicity of notation, we drop the dependence of $D_{dir}$ on $(g_0,h_0, \theta_0)$, though our theory accounts for this dependence.

To study the indirect effect, $d\bar{m}(g_0,h(F_0, g(F_\tau)))/d\tau$, we generalize the key Lemma 1 in hahn2013asymptotic to allow for ML second steps as in equation (ref). The lemma is stated for one-dimensional $S$. For higher dimensions, it must be applied component-wise.

The condition that the function $\delta_2$ belongs to every set $\Delta(g)$ for $g$ close to $g_0$ is related to “regularity" of $\Delta_2(g)$. If the functions in the sets $\Delta_2(g)$, with $g \in \Delta_{1}$, have the same shape, one would expect that many $\delta_2$'s satisfy the condition in the above lemma. The condition allows to take derivatives in equation (ref) along the path $(F_0, g_\tau)$.

To linearize the first step, we ask $\alpha_{02}$ to satisfy the condition for $\delta_{2}$ in Lemma (ref). We also impose some additional assumptions on the paths $\tau \mapsto h(F_0, g_\tau)$. This allows us to express $d\bar{m}(g_0,h(F_0, g(F_\tau)))/d\tau$ as an inner product.

As we have emphasized, this assumption is related to \textquotedblleft regularity" in the shape of the functions in $\Delta _{2}(g)$. It is needed to deal with a non-linear parameter space for $(g, h)$. In both the nonparametric case $\Delta _{2}(g)=L_{2}(X,V(g))$ and the partially linear case $\Delta _{2}(g)=\{\beta ^{\prime }x+\kappa (v)\colon \beta \in \mathbb{R }^{p},\kappa \in L_{2}(V(g))\}$ the assumption translates into square-integrability conditions (see Appendix (ref) for a detailed discussion). What Assumption (ref) rules out, for example, it is to specify a partially linear model for some $g$ and a nonparametric regression for others.

Once we can apply Lemma (ref), the remaining step is to linearize the terms $h_0(X,\varphi(D,Z,g(F_\tau)))$ and $\alpha_{02}(X, \varphi(D,Z,g(F_\tau)))$. To achieve this, we require $h_0$, $\alpha_0$, and $\varphi$ to be differentiable in an appropriate sense:

The Hadamard derivative of $\varphi $ is a linear and continuous map $ D_{\varphi }\colon L_{2}(Z)\rightarrow L_{2}(D,Z)$ such that

To illustrate, if $\varphi (d,z,g)=g(z)$ (first step prediction) or $\varphi (d,z,g)=d-g(z)$ (first step residual), then $D_{\varphi }g=g$ or $D_{\varphi }g=-g$, respectively.

To identify the first-step Riesz representer $\alpha_{01}$ while allowing for general residuals $\epsilon(W,g)$, we need the following assumption ichimura2022influence.

The usual first-step error $\epsilon(W,g)=D-g(Z)$ has $r_e(Z)=-1$, satisfying the above assumption. The Logit-Lasso error $\epsilon(W, g)=D-\Lambda(g_0(Z))$ has $r_e(Z) = -\Lambda(g_0(Z))(1-\Lambda(g_0(Z))$ and satisfies the above assumption if the propensity score is bounded away from zero and one.

The next theorem gives the shape of the first-step IF:

The shape of the first step Riesz representer $\alpha _{01}$ has a rather complex form. Indeed, the linearization with respect to the first step effect is also complex (c.f. equation (ref)). The first term corresponds to the linearization of the direct effect of $g$. It is given by $D_{dir}$, the linearization of $d\bar{m} (g_{\tau },h_{0})/\tau $. The second term corresponds to the indirect effect. Consistent estimation of the second term generally requires estimators for (i) $g_{0}$, (ii) $h_{0}$, (iii) $\partial h_{0}/\partial v$, (iv) $\alpha _{02}$, and (v) $\partial \alpha _{02}/\partial v$. Section (ref) provides the details on how to estimate $\mathbb{E}[D_{01}(W,g)]$. We also note that some simplifications and variations on the expression for $\mathbb{E}[D_{01}(W,g)]$ and for $\phi_{1}$ occur under different scenarios.

Cross-fitted automatic estimators

The debiased sample moment functions are estimated using cross-fitting, where we partition the sample $(W_i)_{i=1}^n$ into $L$ groups $I_\ell$, for $\ell = 1, \dots, L$. Estimation of the debiased moment function $\psi$ for an observation $i \in I_\ell$ requires estimators of the Riesz representers $(\hat{\alpha}_{1\ell},\hat{\alpha}_{2\ell})$ based only on observations not in $I_\ell$. This section is devoted to the construction of automatic estimators satisfying this property. Through the section, we consider that the researcher has at her disposal first and second step estimators, $\hat{g}_{\ell\ell^{\prime }}$ and $\hat{h}_{\ell\ell^{\prime }}$, and a preliminary estimator $\tilde{ \theta}_{\ell\ell^{\prime }}$, that use only observations not in $I_\ell \cup I_{\ell^{\prime }}$; and estimators $(\hat{g}_{\ell\ell^{\prime }\ell^{\prime \prime }}$, $\hat{h}_{\ell\ell^{\prime }\ell^{\prime \prime }}, \tilde{\theta}_{\ell\ell^{\prime }\ell^{\prime \prime }})$ that use only observations not in $I_\ell\cup I_{\ell^{\prime }} \cup I_{\ell^{\prime \prime }}$. Depending on the application, some of the preliminary estimators may not be needed (see, e.g., the debiased estimator in the partially linear model).

Our approach to automatically estimate the Riesz representers relies on the orthogonality conditions discussed in Section (ref). We can combine the orthogonality conditions with the linearization results in Section (ref) to obtain sets of moment conditions for the estimation of the Riesz representers. In particular, a combination of equation (ref) and Theorem (ref) gives

where $D_{01}$ is the linearization of the identifying moment function $\bar{m}$ with respect to the first step and $r_e$ gives the linearization of the generalized error function $\epsilon(w,g)$ (cf. Assumption (ref)). Varying $\delta_1$, the above equation provides a set of moment conditions that identify $\alpha_{01}$. Likewise, identification of the second-step Riesz representer $\alpha _{02}$ follows from equation (ref) and Proposition (ref):

where $D_{02}$ is the linearization of the identifying moment function $\bar{m}$ with respect to the second step. The shape of the linearizations $D_{01}$, $D_{02}$, and $r_{e}$ may vary with the problem (see Section (ref) for some examples).

Equations (ref) and (ref) form the basis for automatic estimation of the Riesz representers. These require finding consistent estimators of the linearizations of the identifying moment functions and the generalized error. In this section, we will write $D_{02}(w,h|g_{0},h_{0},\theta_0)$ to make explicit that the linearization with respect to $h$ may depend on $(h_{0},g_{0},\theta_0)$. For the linearization of the effect of first-step estimation, we will write $D_{01}(w,g|g_{0},h_{0},\alpha _{02},\theta_0 )$, to emphasize that it may also depend on the second-step Riesz representer. $D_{01}$ generally also depends on the derivatives $\partial h_{0}/\partial v$ and $ \partial \alpha _{02}/\partial v$. We do not make this explicit, but we will address the issue in this section. We also write $r_e(z|g_0)$ to express that the generalized error is linearized at $g_0$.

Automatic estimation of the second-step Riesz Representer

We start with the automatic estimator for $\alpha_{02}$. We assume that there is a dictionary $(b_{j})_{j=1}^{\infty } $ whose closed linear span is $\Delta _{2}(g_{0})$. That is, any function in $\Delta _{2}(g_{0})$ can be approximated, in the $L_{2}$ sense, by a linear combination of the atoms. Thus, $\alpha _{02}$ can be approximated by $\mathbf{b}_{J}^{\prime }\boldsymbol{\rho }_{0J}$, where $\mathbf{b} _{J}=(b_{1},...,b_{J})^{\prime }$ and $\boldsymbol{\rho }_{0J}=(\rho _{01},...,\rho _{0J})^{\prime }$. We can now plug in $\mathbf{b}_{J}^{\prime }\boldsymbol{\rho }_{0J}$ into equation (ref) for $ \delta _{2}=b_{j}$, $j=1,...,J$. This gives the following $J$ moment conditions:

where $D_{02}(w,\mathbf{b}_{J})\equiv (D_{02}(w,b_{1}),...,D_{02}(w,b_{J}))^{\prime }$.

The above moment conditions can be used to construct an OLS-like estimator of $\boldsymbol{\rho }_{0J}$. Note, however, that in high-dimensional settings $\mathbb{E}[\mathbf{b}_{J}(X,V)\mathbf{b}_{J}(X,V)^{\prime }]$ may be near singular. Therefore, we use the regularized estimator solving

where $\lVert \boldsymbol{\rho }_{J}\rVert _{q}\equiv (\sum_{j=1}^{J}|\rho _{j}|^{q})^{1/q}$ for $q\geq 1$ and $\lambda \geq 0$ is a tuning parameter. For $q=1$, the above is the Lasso objective function, while $q=2$ corresponds to Ridge Regression. Additionally, we could consider elastic-net-type penalties, where $\lambda (\xi \lVert \boldsymbol{\rho }_{J}\rVert _{2}^{2}+(1-\xi )\lVert \boldsymbol{\rho }_{J}\rVert _{1})$, for $\xi \in \lbrack 0,1]$, replaces the $L_q$ penalization.

For a given $\ell\in{1,...,L}$, the automatic estimator $\hat{\alpha}_{2\ell}$ is based on the sample version of the objective function in equation (ref). We estimate $\mathbb{E}[D_{02}(W,\mathbf{b}_{J})]$ by

where $n_{\ell }$ is the number of observations in $I_{\ell }.$ In turn, $ \mathbb{E}[\mathbf{b}_{J}(X,V)\mathbf{b}_{J}(X,V)^{\prime }]$ is estimated by

With this, we can build an automatic estimator of the second-step Riesz representer that only uses observations not in $I_{\ell }$. It is given by

The tuning parameter $\lambda $ can be chosen by cross-validation.

Automatic estimation of the first-step Riesz representer

We also assume that there is a dictionary $(c_{k})_{k=1}^{\infty }$ that spans $ \Delta _{1}$. This means that $ \alpha _{01}$ can be approximated by $\mathbf{c}_{K}^{\prime }\boldsymbol{ \beta }_{0K}$, where $\mathbf{c}_{K}=(c_{1},...,c_{K})^{\prime }$ and $ \boldsymbol{\beta }_{0K}=(\beta _{01},...,\beta _{0K})^{\prime }$. We can now plug in $\mathbf{c}_{K}^{\prime }\boldsymbol{\beta }_{0K}$ into equation (ref) for $\delta _{1}=c_{k}$, $k=1,...,K$. This gives the following $K$ moment conditions:

where $D_{01}(w,\mathbf{c}_{K})\equiv (D_{01}(w,c_{1}),...,D_{01}(w,c_{K}))^{\prime }$. Recall that $r_e$ gives the derivative of the generalized error $\epsilon$ (see Assumption (ref)).

We use these conditions as a basis to construct the objective function to estimate $\boldsymbol{\beta }_{0K}$:

where the tuning parameter $\lambda $ may be different from that of the second step. The automatic estimator for the first-step Riesz representer is built with the sample version of the above equation. The estimator is given by

with

and

Note that the estimator for the first-step linearization $\hat{D}_{1\ell}$ may require estimators of the second-step Riesz representer $\hat{\alpha}_{2\ell\ell'}$ that do not include observations in $I_\ell \cup I_{\ell'}$. These estimators can be obtained using the methodology of the previous section. To construct $\hat{\alpha}_{2\ell \ell ^{\prime }}=\mathbf{b}_{J}^{\prime }\widehat{\boldsymbol{\rho }}_{J\ell \ell ^{\prime }}$, we let $\widehat{ \boldsymbol{\rho }}_{J\ell \ell ^{\prime }}$ solve the optimization problem in equation (ref), with $\hat{D}_{2\ell}$ and $\hat{ B}_{\ell}$ replaced by

respectively.

Furthermore, $D_{01}$ may depend on the derivatives $\partial h_{0}/\partial v$ and $ \partial \alpha _{02}/\partial v$ (see equation (ref)). We thus need to provide consistent estimators of these derivatives to build $\hat{D}_{1\ell }$. It is straightforward to construct an estimator $\partial\hat{\alpha}_{2\ell\ell^{\prime }}/\partial v$ of the derivative of $\alpha_{02}$ based on the cross-fitted Lasso estimator. Since we have already estimated $\hat{\alpha} _{2\ell\ell}= \mathbf{b}_J^{\prime }\widehat{\boldsymbol{\rho}} _{J\ell\ell^{\prime }}$, if each $b_j$ is differentiable w.r.t. $v$, we have that $\partial\hat{\alpha}_{2\ell\ell^{\prime }}/\partial v\equiv (\partial \mathbf{b}_J/\partial v)^{\prime }\widehat{\boldsymbol{\rho}}_{J\ell\ell^{\prime }}$.

Estimation of $\partial h_{0}/\partial v$ may be trickier. It will depend on the shape of the estimator $\hat{h}_{\ell \ell }$. Note that, since $ h_{0}\in \Delta _{2}(g_{0})$, we may use the dictionary $(b_{j})_{j=1}^{ \infty }$ to approximate the parameter. In this case, $\hat{h}_{\ell \ell }$ will be a Lasso or Ridge Regression estimator and we can estimate the derivative of $h_{0}$ as we have estimated the derivative of $\alpha _{02}$. Moreover, if estimating $h_{0}$ involves a nonparametric regression problem with a low-dimensional covariate, we can often take $\hat{h}_{\ell \ell }$ as a Kernel or a Local Linear Regression estimator, as in heckman1998matching. Then, the derivatives of $h_{0}$ can be estimated by finding the analytical expression of the derivatives of the kernel function.

For a general ML estimator $\hat{h}_{\ell\ell^{\prime }}$ (e.g., Random Forest), we propose a numerical derivative approach to estimate $\partial h_0/\partial v$. Let $t_n$ be a tuning parameter depending on the sample size with $t_n \downarrow 0 $. We propose to estimate $\partial h_0(x,v)/\partial v$ by

This approach has been used and justified theoretically in bravo2020two in a two-step setting. Note that, usually, we need to compute the derivative evaluated at $(X_i, \varphi(D_i,Z_i, \hat{g}_{\ell\ell^{\prime }}))$. Alternative ML estimators that achieve optimal rates for partial derivatives are discussed in dai2016optimal.

Examples: Three-step debiased estimators

Hd-PS regression adjustment

The debiased three-step estimator is

where $\hat{V}_{i\ell} = \Lambda(\hat{g}_\ell(Z_i))$ and $\Lambda$ is the Logistic function. Cross-fitted estimators $(\hat{g}_\ell, \hat{h}_\ell)$ are discussed in Example (ref) (p. (ref)). In this section, we detail the estimation of the Riesz representer $\alpha_{01}$. As discussed in Section (ref), we propose to estimate the Riesz representer by $\hat{\alpha}_{1\ell }(z)=\mathbf{c}_{K}(z)^{\prime }\widehat{\boldsymbol{\beta } }_{K\ell }$ with $\widehat{\boldsymbol{\beta }}_{K\ell }$ solving (ref). We show how to construct $\hat{C}_\ell$ and $\hat{D}_{1\ell}$.

The term $\hat{C}_\ell$ depends on the linearization of the first-step generalized error $\epsilon(w, g) = d - \Lambda(g(z))$. Since

Therefore,

To find $\hat{D}_{1\ell}$, note that, since $\alpha_{02}=0$, we have that $D_{01} = D_{dir}$ and

where $\dot{h}_0 \equiv dh_0/dv$ and $\dot{\Lambda} \equiv d\Lambda/du = \Lambda \cdot (1 - \Lambda)$. From this representation and $\varepsilon=Y - h_0(V) -\theta_0(D-V)$, it follows that the effect of the first step is zero if $\mathbb{E}[\varepsilon|D,Z] = 0$, which we do not assume as it imposes strong restrictions on heterogeneity. To estimate $\alpha_{01}$, the linearization $D_{dir}$ is projected onto $\Delta_1 \subseteq L_2(Z)$. Therefore, as $\mathbb{E}[D|Z] = V$, the expression for $D_{dir}$ simplifies, and we consider the following estimator for the linearization of the first step:

Note that, in this case, the linearization $D_{01}$ does not depend on $\theta_0$ and $\alpha_{02}$. Hence, no additional estimators $\tilde{\theta}_{\ell\ell'}$ and $\hat{\alpha}_{2\ell\ell'}$ are needed.

For the Hd-PS regression adjustment estimator, the estimation of the first-step Riesz representer can be reframed as a weighted Lasso regression, which is defined in equation (ref).

Partially linear model with generated regressors: Autoencoders

The partially linear model is a workhorse for debiased machine learning methods, see ahrens2025introduction and references therein. Here we propose a debiased estimator for the partially linear model that is robust to ML-generated regressors. We first consider a general $V=\varphi(D,Z,g_0)$, where $g_0$ is identified by (ref). Then, we illustrate the framework with learned confounders via autoencoders.

Suppose $\dim(D)=p$. For the partially linear model with generated regressors, introduce the second-step nuisances $h_{0S}(v)=\mathbb{E}[S\mid V=v]$ for $S=Y$ or $S=D_j$, for $j=1,\dots,p$. Let $h_{0D} \equiv (h_{0D_1}, \dots, h_{0D_p})'$. The DML-type identifying moment is

where $h_0\equiv(h_{0Y},h_{0D})$. We assume $\kappa_0\in\Delta_2(g_0)$ so this moment identifies $\theta_0$ for the relevant second-step space. In this example, $\alpha_{02}=0$, but $\alpha_{01}$ is generally nonzero (cf.\ Section (ref)), so standard DML inference that ignores generated regressors is not generally valid.

We therefore use a debiased three-step estimator that is robust to the first step. The estimator solves

with $\hat{\psi}(\theta) = n^{-1} \sum_{\ell=1}^L\sum_{i\in I_\ell} \hat{\psi}_{i\ell}(\theta)$, weighting matrix $\hat\Upsilon$, and debiased moments

where $\widehat{\boldsymbol{\alpha}}_{1\ell}(z)$ is the $p$-vector of automatic first-step Riesz-representer estimators (one per moment component). Each component is constructed as in Section (ref): $\hat\alpha_{1j\ell}(z)=\mathbf c_K(z)'\widehat{\boldsymbol\beta}_{Kj\ell}$, with $\widehat{\boldsymbol\beta}_{Kj\ell}$ solving (ref). The construction of $\hat C_\ell$ depends on $\epsilon(w,g)$. For example, in a control-function setup with $V=D-g_0(Z)$ one has $\epsilon(w,g)=d-g(z)$ and $r_e(z)=-1$.

Construction of $\hat D_{1j\ell}$, for each $j=1,\dots, p$, parallels Section (ref). Since $\alpha_{02}=0$, the indirect effect is zero and only the direct effect remains. The linearization of the $j$-th moment condition is

where $\dot{h}_{0S} \equiv dh_{0S}/dv$. To build $D_{01j}(W_i, \mathbf{c}_K|\hat{g}_{\ell\ell'}, \hat{h}_{\ell\ell'}, \tilde{\theta}_{\ell\ell'})$, we replace these terms in equation (ref): (i) $\varepsilon$ by $\hat\varepsilon_{i\ell\ell'}=Y_i-\hat h_{\ell\ell',Y}(\hat V_{i\ell\ell'})-\tilde\theta_{\ell\ell'}'(D_i-\hat h_{\ell\ell',D}(\hat V_{i\ell\ell'}))$, (ii) $\dot\varepsilon$ by $\dot\varepsilon_{i\ell\ell'}=-\dot h_{\ell\ell',Y}(\hat V_{i\ell\ell'})+\tilde\theta_{\ell\ell'}'\dot h_{\ell\ell',D}(\hat V_{i\ell\ell'})$, (iii) $\dot{h}_{0S}$ by the corresponding cross-fitted estimators $\dot h_{\ell\ell',S}$, and (iv) $D_\varphi\mathbf{c}_K$ by a cross-fitted estimator $\hat D_{\varphi i\ell\ell'}\mathbf{c}_K$ of the linearization of $\varphi$ w.r.t. $g$; e.g., for $\varphi(d,z,g)=\Lambda(g(z))$, $\hat D_{\varphi i\ell\ell'}\mathbf{c}_K=\Lambda(\hat g_{\ell\ell'}(Z_i))(1-\Lambda(\hat g_{\ell\ell'}(Z_i)))\mathbf{c}_K(Z_i)$, while for $\varphi(d,z,g)=d-g(z)$, $\hat D_{\varphi i\ell\ell'}=-\mathbf{c}_K(Z_i)$.

In general, $\alpha_{01}$ is nonzero. It is the orthogonal projection onto $\Delta_1$ of $D_\varphi^*r_{dir}$, where $D_\varphi^*$ is the adjoint of $D_\varphi$. Even if $\mathbb{E}[\varepsilon\mid D,Z]=0$, $\alpha_{01}$ typically remains nonzero, so inference that ignores generated regressors is invalid. For $V=g_0(Z)$, $\Delta_1=L_2(Z)$ and $\Delta_2=L_2(V)$, these influence function calculations are covered by hahn2013asymptotic.

\noindentLearned confounders via autoencoders.

An autoencoder consists of an encoder $e_0(\cdot)$, a low-dimensional representation $V=e_0(Z)$ (our generated regressor), and a decoder $d_0(\cdot)$, identified by

where $\mathcal{L}$ is a loss and $\mathcal E$ and $\mathcal D$ are function classes (see Figure (ref)). For concreteness, we take $\mathcal{L}(Z,f)=|Z-f|^2$ and feed-forward neural networks indexed by $\zeta\in\mathbb{R}^{\dim(\zeta)}$, so that $e_0(Z)=e_{\zeta_0}(Z)$ and $d_0(V)=d_{\zeta_0}(V)$ for some minimizer $\zeta_0$. A key feature of autoencoders is the bottleneck $\dim(V)\ll\dim(Z)$, which yields nonlinear dimension reduction bengio2013representation.

To write this example in our setting, define $g_0=(d_0,e_0)$ and $V=\varphi(D,Z,g_0)=e_0(Z)$, and assume w.l.o.g.\ $\dim(V)=1$. The generalized error is $\epsilon(W,g)=Z-d(e(Z))$, $g=(d,e)$, and the first-step identifying condition is

where $\Delta_1$ is the linear span generated by the columns of the Jacobian

Define also $\mathbf e_K \equiv \left.\partial e_\zeta/\partial\zeta\right|_{\zeta=\zeta_0}$, which shows up in the linearization of the moment condition w.r.t. $g$ (the direct effect). The implementation follows the generic construction with $\mathbf c_K(Z)=\dot f_0(Z)$, $K=\dim(\zeta)$, and objective

with $D_{01j}$ given by equation (ref) with $D_\varphi\mathbf{e}_K = \mathbf{e}_K$.

Since $\zeta_0$ is unknown, we use $\hat\zeta_{\ell\ell'}$, which is estimated without observations in $I_\ell\cup I_{\ell'}$, and compute Jacobians by backpropagation:

Then,

with $\hat V_{i\ell\ell'}=e_{\hat\zeta_{\ell\ell'}}(Z_i)$. Solving the sample analog of (ref) yields $\widehat{\boldsymbol\beta}_K$ and $\widehat{\boldsymbol\alpha}_{1\ell}(Z)=\hat{\mathbf c}_{K,\ell\ell'}(Z)'\widehat{\boldsymbol\beta}_K$. The three-step debiased estimator uses the moment function in (ref), with $\epsilon(W_i, \hat{g}_\ell) = Z_i - f_{\hat{\zeta}_\ell}(Z_i)$ and $\hat{\zeta}_\ell$ estimated without observations in $I_\ell$.

CASF in a non-separable model

The three-step debiased estimator of the CASF is given in equation (ref). We provide the ingredients to build the estimators $\hat{\alpha}_{1\ell}$ and $\hat{\alpha}_{2\ell}$. Recall that the moment function defining the CASF is

This moment is already linear in $h$ and hence

for each atom $b_{j}$ in the dictionary. We follow the same strategy as before and approximate $D_{02}$ by Monte Carlo integration. Let $(X_{s}^{\ast })_{s=1}^{S}$ be a sample drawn from $F^{\ast }$. To construct the objective function to estimate $\widehat{\boldsymbol{\rho }} _{J\ell }$, for an observation $i\in I_{\ell ^{\prime }}$, we set

for each $j=1,\dots ,J$. Here we emphasize that the linearization does not depend on $h_0$ and $\theta_0$, it only depends on $g_0$. With this we construct $\hat{\alpha}_{2\ell }= \mathbf{b}_{J}^{\prime }\widehat{\boldsymbol{\rho }}_{J\ell }$ following ((ref)).

It is straightforward to show that the linearization of the moment condition w.r.t. $g$ is $D_{dir}(w,g)=r_{dir}(w)g(z)$, with

We can now plug in the expression for $D_{dir}$ into equation (ref), where the linearization of the first step effect is defined. Recall that $D_{\varphi }g=-g$. Then, for the CASF, equation (ref) becomes

The linearization depends on $h_{0}$ and $\alpha _{02}$, and the derivative of $h_{0}$ w.r.t. $v$. It also depends on $g_{0}$, as $v\equiv d-g_{0}(z)$. However, it does not depend on $\theta_0$. Note that $\mathbb{E}[\partial\alpha_{02}/\partial v \cdot (Y - h_0)]=0$ by the control-function assumption.

We approximate $r_{dir}(W_i)$, with $i\in I_{\ell ^{\prime }}$, by

where, ${\partial \hat{h}_{\ell \ell ^{\prime }}}/{\partial v} =(\partial \mathbf{b}_{J}/\partial v)^{\prime }\widehat{\boldsymbol{\eta}}_{\ell \ell ^{\prime }}$. The parameters $\boldsymbol{\hat{\eta}}_{\ell \ell ^{\prime }}$ are Lasso cross-fitted slope estimates for the second step $ h_{0}$. To estimate $D_{1\ell }$, it remains to show how to estimate $\alpha_{02} \cdot \partial h_0/\partial v$ for an observation $i\in I_{\ell ^{\prime }}$. Being $ V_{i\ell \ell ^{\prime }}\equiv D_{i}-\hat{g}_{\ell \ell ^{\prime }}(Z_{i})$, we can estimate it by

Therefore, we have that, for $i\in I_{\ell ^{\prime }}$,

for each $k=1,...,K$. Finally, note that $\epsilon(w,g)=d-g(z)$ and, hence, $r_e(z) = -1$. These results can then be used to construct the objective function to estimate $\widehat{\boldsymbol{\beta }}_{K\ell }$ in ((ref)) and then $\hat{\alpha}_{1\ell }=\mathbf{c}_{K}^{\prime }\widehat{\boldsymbol{\beta }}_{K\ell }$.

Asymptotic theory

General results

This section gives general conditions for asymptotic normality of the automatic debiased GMM and conditions for consistent estimation of its asymptotic variance. The conditions are based on the mean-square consistency, small interaction of estimation biases, and locally robust conditions. These asymptotic results generalize chernozhukov2022locally to our three-step setting with generated regressors. Estimation rates for the Riesz representers $(\alpha_{01}, \alpha_{02})$ require (i) that the dictionaries approximate well the Riesz representers and (ii) being able to estimate the linear approximations of $\bar{m}(g,h)$ given by $ \mathbb{E}[D_{01}(W,g)]$ and $\mathbb{E}[D_{02}(W,h)]$ at a certain rate chernozhukov2022automatic.

In the presence of generated regressors, the theory needs to account for the fact that the estimator of the correction term (and probably that of the moment condition) evaluates the estimators $\hat{h}_\ell$ and $\hat{\alpha} _{2\ell}$ in the generated regressor $\hat{V}_{i\ell} \equiv \varphi(D_i, Z_i, \hat{g}_\ell)$ (c.f., equation (ref)). We modify the expansion of $\hat{\psi}_{i\ell}(\theta_0)-\psi(W_i, g_0,h_0,\alpha_{0},\theta_0)$ given by chernozhukov2022locally to deal with this fact. After a first order expansion, which forms the basis of the local robustness property, remainders implying the generated regressor are of a particularly complex form. In the case of downstream local robustness, when $\alpha_{02} = 0$, the remainder simplifies. In any other cases, we rely on smoothness conditions on the dictionaries and $g \mapsto \varphi(D,Z, g)$ to bound the remainder (c.f. Assumption (ref)).

We begin with assumptions on the dictionaries. The first assumption formally states that the dictionaries $(b_j)_{j=1}^\infty$ and $(c_k)_{k=1}^\infty$ span $\Delta_{2}(g_0)$ and $\Delta_1$, respectively.\footnote{ In this section, for a measurable function $f$, $\lVert f \rVert_2 \equiv \sqrt{\mathbb{E}[f(W)^2]}$ denotes its $L_2$-norm. Also, for a $m\times n$ matrix $A=(A_{i,j})_{i=1,j=1}^{m,n}$, $\lVert A\rVert_\infty \equiv \max_{i,j} |A_{ij}|$.}

We also assume bounded dictionaries newey1997convergence:

The assumption translates into consistency of $\hat{B}_\ell$ and $\hat{C} _\ell$. Also, on top of the following assumption, it will guarantee that the Riesz representers are bounded:

This assumption keeps the $L_1$-norm of the coefficient of the Lasso penalized regression under control. The result is relevant to estimate the asymptotic variance chernozhukov2022automatic. We also note that the absolute summability of the coefficients imposes a sparsity condition on the relevant terms to approximate $\alpha_{01}$ and $\alpha_{02}$ chernozhukov2022automatic.

We require the following estimation rates:

This assumption imposes standard rate conditions on the estimators of the nuisance parameters and on the linearization of the moment condition. For general results on rates with generated regressors see mammen2012nonparametric; for Lasso rates, see bickel2009simultaneous,bunea2007sparsity,zhang2008sparsity, and references therein; for $L_2$-rates with boosting with high-dimensional regressors see kueck2023estimation; for deep neural networks with a ReLU activation function, see farrell2021deep. Under some regularity conditions on the linearizations chernozhukov2022automatic, Assumption (ref) can be derived from the rate conditions on the estimators of the nuisance parameters.

We also ask for the following rates for the Lasso penalty and the number of terms in the dictionaries:

This assumption asks for the Lasso penalty to go to zero slightly slower than $n^{-r}$, where $r$ is the rate from Assumption (ref). For instance, a rate of $\log(n)/n^{r}$ is allowed. Moreover, it requires polynomial rates in the growth of the number of terms in the dictionaries.

The above are general conditions imposed on the dictionaries and the tuning parameters for the Lasso penalized regression. The specific problem at hand only appears in two instances. First, Assumption (ref) requires that the dictionaries approximate well the correction-term Riesz representers (living in $\Delta_1$ and $\Delta_2(g_0)$, respectively). Second, Assumption (ref) requires (i) mean-square rates for the estimators of $g_0$ and $h_0$ and (ii) to be able to estimate the linearizations at the same rate. As discussed before, these conditions provide rates of estimators of the Riesz representers $\alpha_{01}$ and $\alpha_{02}$ chernozhukov2022automatic. For instance, the convergence rate of $\hat{\alpha}_{1\ell}$ will be fast enough to guarantee that the interaction term satisfies $\lVert \hat{\alpha} _{1\ell}-\alpha_{01}\rVert_2 \cdot \lVert \hat{g}_\ell-g_0\rVert_2 = o_p(n^{-1/2})$ chernozhukov2022locally.

We now provide assumptions on the moment condition. The first is a mean-square consistency condition similar to Assumption 1 in chernozhukov2022locally:

Assumption (ref) is necessary for regular estimation of $\theta_0$. Assumptions (ref) and (ref) are mean-square consistency conditions for the moment condition. Boundedness of the conditional errors (Assumption (ref)) translates into mean-square consistency conditions for the debiasing term $\phi$. We repeat here that the boundedness of $\alpha_{01} $ and $\alpha_{02}$ is implied by Assumptions (ref) and (ref).

We also need to strengthen Assumption (ref) to control the remainder for linealizing the generalized error:

The following assumption is standard in the literature, see, e.g., newey1994asymptotic. It imposes a quadratic remainder bound for the first-order linearization of $\bar m(g,h)$ and therefore strengthens Assumptions (ref) and (ref). Consider the linearization \[ \bar{\psi}(g,h) \equiv \mathbb{E}\!\left[ m(W,g,h,\theta_0)-m(W,g_0,h_0,\theta_0) -D_{dir}(W,g-g_0) -D_2(W,h-h_0) \right]. \] Note that the linearization treats both $g$ and $h$ as “independent" nuisance, i.e., it does not account for the fact that $g$ affects estimation of $h$. We assume the following:

To account for the generated regressors, we introduce the following assumption. Its goal is to guarantee that the remainder of the chain rule in our Lemma (ref), which accounts for the indirect effect, is quadratic. To state the assumption, we introduce the mapping $\nu(h, \alpha_2)\equiv[\partial\phi_2/\partial v](w, g_0, h, \alpha_2) = \partial/\partial v\{\alpha_2(x,v) \cdot (s - h(x,v)\}$ (c.f. Lemma (ref)). Define $\hat{\nu}_\ell=\nu(\hat{h}_\ell, \hat{\alpha}_{2\ell})$ and $\nu_0=\nu(h_0, \alpha_{02})$.

First, under downstream local robustness, the indirect effect is zero, and the above conditions are not needed. Regarding these conditions, Assumption (ref) asks for a quadratic remainder in the linearization of the generated regressor. Assumptions (ref)-(ref) strengthen Assumption (ref). Assumption (ref) also strengthens Assumption (ref) (see Appendix (ref) for a general discussion of this assumption). Assumption (ref) is a product-rate condition to handle higher-order components from the generated regressors. When estimation of $h_0$ is conducted by (penalized) regression onto the dictionary $(b_j)_{j=1}^\infty$, this assumption may be understood as smoothness conditions on the dictionary. This is the case of the CASF example. We provide a detailed discussion of Assumption (ref) in Section (ref), where we verify it for the CASF, and in Appendix (ref).

Finally, the GMM procedure requires consistent estimation of the Jacobian of the moment condition. Being the following assumption specific to the GMM procedure, it is stated for the case with an arbitrary number of parameters and moment conditions.

Assumptions (ref), (ref), (ref), and (ref) are stated for a single moment condition. In the presence of more than one condition, they must be understood to hold componentwise. Assumption (ref), since it refers to a GMM-specific situation, is already formulated in the general case. The remaining assumptions do not depend on the dimension of the moment condition (they depend, on the other hand, on the dimension of $Y$ and $D$).

Let $\Xi \equiv (M'\Upsilon M)^{-1}M'\Upsilon' \Psi \Upsilon M (M'\Upsilon M)^{-1}$ be the usual asymptotic variance of the GMM estimator based on the debiased moment functions, where

Define the plug-in estimator $\hat{\Xi}\equiv (\hat{M}^{\prime }\hat{ \Upsilon}\hat{M})^{-1}\hat{M}^{\prime }\hat{\Upsilon}^{\prime }\hat{\Psi} \hat{\Upsilon}\hat{M}(\hat{M}^{\prime }\hat{\Upsilon}\hat{M})^{-1}$, where $\hat{M}$ and $\hat{\Psi}$ are given by the corresponding cross-fitted sample analogs

The following theorem ensures asymptotic normality of $\sqrt{n}(\hat{\theta}-\theta _{0})$:

Regularity conditions for some examples

Hd-PS regression adjustment

Here we verify Assumptions (ref), (ref), (ref), (ref), and (ref). These are the assumptions that explicitly depend on the identifying moment condition and the first and second estimation steps. We also provide sufficient conditions for Assumptions (ref) and (ref) required for linearizing the moment condition. The moment condition that identifies $\theta_0$ in the partially linear model is:

were recall that $\Lambda$ stands for the logistic cdf and $V=\Lambda(g_0(z))$.

The following assumption gives the result:

Assumption (ref) is the usual overlap assumption. Assumption (ref) bounds $\mathbb{E}[Y^2|D,Z]$ (note that $Y$ may still be supported on $\mathbb{R}$). If regressors $Z$ have compact support, continuity of $E[Y^2|D=d,Z=z]$ would be sufficient for Assumption (ref). Assumption (ref) imposes smoothness conditions on the atoms. Note that if $h_0(v) = \mathbb{E}[Y|V=v]$ is smooth enough, the econometrician can always choose a dictionary with smooth atoms to estimate it.

We show that the assumptions for Theorem (ref) holds in the Hd-PS setting:

CASF in a non-separable model

Here we verify Assumptions (ref), (ref), (ref), (ref), and (ref) for the CASF. To achieve this, we require some regularity on the distribution of $(X,V)$, where $V = D - g_0(Z)$, on the second step $h_0(x,v) \equiv \mathbb{E}[Y|X=x, V=v]$, and on the dictionary that is used for second-step estimation.

Assumption (ref) is guarantees regular identification of the CASF. Note that, in the case of the CASF, the linearization of the second step can be written as $\mathbb{E}[D_2(w, h)] = \mathbb{E}[r_2(X,V)h(X, V)]$, with $r_2 = f^* f_{v} / f_{xv}$. This assumption ensures finite variance of $r_2$ and ask for additional smoothness conditions. Assumption (ref) is the usual bounded conditional variance assumption. Assumption (ref) also imposes smoothness conditions on $h_0$, while Assumption (ref) requires the dictionary used to estimate $h_0$ to satisfy the same smoothness conditions.

Assumption (ref) requires product-rate conditions involving the estimation error of the derivatives of the second-step nuisance functions. These conditions hold under standard sparse high-dimensional assumptions when $\hat h_\ell$ and $\hat\alpha_{2\ell}$ are estimated by Lasso on the dictionary $\mathbf b_J(x,v)$ introduced in Section (ref). In particular, if $h_0$ and $\alpha_{02}$ admit sparse approximations on $\mathbf b_J$, if the Gram matrices $\mathbb{E}[\mathbf b_J(X,V)\mathbf b_J(X,V)']$ and $\mathbb{E}[(\partial \mathbf b_J(X,V)/\partial v) (\partial \mathbf b_J(X,V)/\partial v)']$ have eigenvalues bounded away from zero and infinity, and if the Lasso estimators achieve the usual $L_1$ coefficient rates, then the derivative estimation errors satisfy $\|\partial\hat h_\ell/\partial v-\partial h_0/\partial v\|_2 =O_p(s_h\sqrt{\log J/n})$ and $\|\partial\hat\alpha_{2\ell}/\partial v-\partial \alpha_{02}/\partial v\|_2 =O_p(s_\alpha\sqrt{\log J/n})$ bickel2009simultaneous. Hence, Assumption (ref) holds whenever these rates multiplied by the first-step rate $\|\hat g_\ell-g_0\|_2$ are $o_p(n^{-1/2})$. A detailed verification is given in Appendix (ref).

We can then show that the assumptions for Theorem (ref) hold for the moment condition defining the CASF.

Conclusion

We propose Automatic Locally Robust estimators for structural parameters in the presence of ML-generated regressors. We show that the debiasing correction term can be decomposed into terms accounting for the first-step and the second-step estimation. Each of the first- and second-step IFs depends on an additional Riesz representers, which can be automatically estimated (i.e., estimated without finding their analytic shape).

We apply our results to construct Automatic Locally Robust estimators for causal treatment effects and the CASF under different modelling assumptions (partially linear and nonparametric models) and different generated regressors (Hd-PS, autoencoders, control function, etc). The analytic shape of the Riesz representers in these cases is particularly complex. For the partially linear model, our automatic debiased estimator overcomes the large biases of the state-of-the-art method, the DML, which does not account for the generated regressors. For the CASF parameter, the moment condition depends on the whole shape of the second-step nuisance parameter (not only its pointwise value), making existing results on generated regressors not applicable even in low-dimensional scenarios. Therefore, automatic estimation is particularly well suited for these problems. We have shown that commonly used plug-in or DML methods lead to highly biased inferences with ML-generated regressors. Three-step debiased estimators correct the bias and deliver much more accurate inference in a complex setting with ML-generated regressors.

\addcontentsline{toc}{section}{References} \makeatletter \makeatother

\cleardoublepage