The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
156,087 characters
Automatic Locally Robust GMM with Machine-Learning-Generated Regressors
\title{Automatic Locally Robust GMM with Machine-Learning-Generated Regressors\thanks{ First version: 25 Jan 2023 (arXiv:2301.10643v1). Research supported by MICIN/AEI/10.13039/501100011033, grant
CEX2021-001181-M, Comunidad de Madrid, grants EPUC3M11 (V PRICIT) and
H2019/HUM-5891, and grant PID2021-127794NB-I00 (MCI/AEI/FEDER, UE), Programa Primas y Problemas de la Fundaci\'on BBVA 2023. }}
\author{Juan Carlos Escanciano \\
\textit{Universidad Carlos III de Madrid} \and Telmo P\'{e}rez-Izquierdo
\\
\textit{University of the Basque Country} }
\date{\today}
\maketitle
\begin{abstract}
Machine-learning (ML) methods now routinely generate regressors used in subsequent econometric analyses—for example, estimated propensity scores, control-function residuals, imputed covariates, learned proxies, or low-dimensional embeddings of high-dimensional data. As these ML-generated regressors become ubiquitous, the lack of general inference methods for models that use them has become a critical limitation. Standard plug-in and Double ML procedures ignore how generated regressors enter later stages, leading to large biases and invalid inference. We develop a three-step locally robust GMM framework for inference with ML generated regressors. A key new insight is \textit{downstream local robustness}: by a functional chain rule, moment functions that are constructed to be orthogonal to the second step eliminate the complicated indirect (conditioning) effects from the ML-generated regressors. We show how to implement this automatically by estimating the associated Riesz representers through cross-fitted auxiliary regressions, allowing for generic non-Donsker ML in both early steps. In leading treatment-effect and counterfactual settings, simulations demonstrate severe bias in existing methods and reductions of 85–95\% using our procedures.
\vspace{1cm}
\noindent Keywords: Locally robust; Machine learning; Generated regressors; GMM; Orthogonal scores; High-dimensional estimation; Causal inference.
\vspace{0.5cm}
\noindent JEL Classification: C13; C14; C18; C21.
\end{abstract}
\newpage
\section{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 \citep[see, e.g.,][]{stock1989nonparametric,stock1991nonparametric,blundell2004endogeneity,imbens2009identification}, average partial effects in sample selection models \citep{ahn1993semiparametric,das2003nonparametric,newey2009two}, propensity score matching \citep{heckman1998matching,abadie2006large}, and marginal treatment effects \citep{heckman2005structural}. More recently, machine learning (ML) is routinely used to generate regressors for imputing missing covariates \citep{fongtyler}, dimension reduction \citep{sorzano2014survey}, learned proxies, confounders, and treatments \citep{knox2022testing}, and feature engineering with unstructured data such as text, images, or audio \citep{feder2022causal}, among many others.\footnote{\cite{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 \citep{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 \cite{shumailov2024ai}.} These difficulties are exacerbated in ML settings, where preliminary estimators are high-dimensional, nonparametric, and often non-Donsker \citep{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 \citet{chernozhukov2022locally}. A new idea is \textit{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. \citealp{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{twostep}).
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 \citet{hahn2013asymptotic} to a more general class of three-step procedures that include leading ML methods. This establishes the broader orthogonality principle of \textit{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{prop:zero_indirect}).
Second, we provide automatic estimation of influence functions and asymptotic variances for models with generated regressors. Under a linearization assumption \citep[as in, e.g.,][]{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., \citet{chernozhukov2021automatic,chernozhukov2022locally,chernozhukov2023automatic,chernozhukov2022automatic}. Automatic estimation is particularly well motivated for generated regressors, where the Riesz representers typically have complex forms \citep[see, for instance,][]{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 \citep{bengio2013representation}. The nonparametric ATE estimator with a Hd-PS generalizes \cite{heckman1998matching}, \cite{hahn2013asymptotic}, and \cite{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 \citep[see, among many others,][]{ichimura1991semiparametric,ahn1993semiparametric,heckman1998matching,newey1999nonparametric,li2002semiparametric,rothe2009semiparametric,imbens2009identification,escanciano2010testing,song2012smoothness}. In an important work, \citet{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, \citet{mammen2012nonparametric,mammen2016semiparametric} and \citet{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 \citep[see][]{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 \citep[e.g.,][]{chernozhukov2018double,chernozhukov2022locally,chernozhukov2022automatic}. With the exception of \citet{sasaki2021estimation}, the DML literature prior to our work has not considered or accounted for generated regressors in inference. Our results complement \citet{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 \citet{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{Setting} introduces the setting and examples. Section~\ref{sec:debiased_est} 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{sec:Automatic} gives the separate identification and automatic estimation of the Riesz representers. Debiased automatic estimators for the examples are presented in Section~\ref{sec:Examples}. The asymptotic theory is developed in Section~\ref{sec:asymptotic}. Section~\ref{sec:conclusion} concludes. Appendix~\ref{sec:algorithm} summarizes the estimation steps. Appendix~\ref{sec:AdditionalExamples} provides a further application to a ML implementation of the nonparametric ATE estimator of \cite{heckman1998matching}. Appendix~\ref{sec:MC_details} contains details about the Monte Carlo simulations. Appendix~\ref{sec:app_inclusion} discusses regularity conditions, and Appendix~\ref{sec:Proofs} gathers the proofs of the main results.
\section{Setting and examples}
\label{Setting}
\subsection{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: \bigskip
\textbf{First step}. There is a first-step nuisance function $g_{0}(Z)$ satisfying the moment restrictions
\begin{equation}
\mathbb{E}[\delta _{1}(Z)\epsilon (W,g_{0})]=0\text{ for all }\delta _{1}\in
\Delta _{1}, \label{orth1}
\end{equation}
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 \cite{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$ \citep[see Chapter 7 in][]{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 \cite{hahn2013asymptotic}. For numerous other examples of $\epsilon (W,g_{0})$, including
quantile regression, see Section 3 of \citet{ichimura2022influence}. For
general parametric first steps, see Remark \ref{Parfirststep}; and for other first steps not covered by our setting, see Remark \ref{OtherFS}.
The first-step nuisance $g_{0}$ in (\ref{orth1}) is used to construct the population generated
regressors
\begin{equation*}
V\equiv \varphi (D,Z,g_{0}),
\end{equation*}
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 \cite{heckman1998matching} has $V=\Lambda (g_{0}(Z))$; some dimension-reduction methods have $\varphi (D,Z,g_{0})=g_{0}(Z)$, as in \cite{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.\bigskip
\textbf{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
\begin{equation}
\mathbb{E}[\delta _{2}(X,V)(S-h_{0}(X,V))]=0\text{ for all }\delta _{2}\in
\Delta _{2}(g_{0}), \label{orth2}
\end{equation}
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{orth2}) 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., \cite{chernozhukov2022locally}. \citet{hahn2013asymptotic} and \citet{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).
\bigskip
\textbf{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
\begin{equation*}
\mathbb{E}[m(W,g_{0},h_{0},\theta _{0})]=0.
\end{equation*}
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 $.
\subsection{Examples}
The following examples are used to illustrate the main results of this paper.
\begin{ex}[Partially linear model with ML-generated regressors]
\label{ex:PLM_ML}
We first consider a general partially linear model with a generated regressor.
We observe $W\equiv (Y, D, Z)$, where $Y$ is an outcome, $D$ is a (possibly
vector-valued) treatment or regressor of interest, and $Z$ is a (potentially)
high-dimensional covariate. The starting point is the partially linear model
\begin{equation*}
Y = \theta_0' D + \kappa_0(V) + \varepsilon
\quad \text{with } \mathbb{E}[\varepsilon\mid D,V] = 0,
\end{equation*}
where $\theta_0$ is the parameter of interest, $V \equiv \varphi(D,Z,g_0)$
is a generated regressor constructed from a first-step nuisance $g_0$ and the
covariates $(D,Z)$, and $\kappa_0(\cdot)$ is an unknown nuisance function.
\medskip
\noindent\emph{First step and generated regressor.}
The first step $g_0$ solves \eqref{orth1} for some generalized error
$\epsilon(W,g_0)$ and a linear, closed subspace $\Delta_1\subseteq L_2(Z)$. This encompasses, for example:
\begin{itemize}
\item \emph{Hd-PS:}
$\epsilon(W,g_0)=D-\Lambda(g_0(Z))$, with $\Delta_1$ the mean-square limit
of sparse linear combinations $\sum_{k=1}^K \gamma_{0k} c_k(Z)$, and
$V = \Lambda(g_0(Z))$ (for treatment effects, sample selection, etc.).
\item \emph{Dimension reduction / learned confounders (e.g., autoencoders):}
$\epsilon(W,g_0)=Z-d_0(e_0(Z))$, $g_0=(d_0,e_0)$, with $\Delta_1$ defined in Section~\ref{sec:ex_PLMauto}, and
$V = e_0(Z)$ equal to the encoder.
\item \emph{Control-function residuals:}
$\epsilon(W,g_0)=D-g_0(Z)$, with $\Delta_1=L_2(Z)$, and $V = D-g_0(Z)$ equal to the residual of the first-stage regression of $D$ on $Z$.
\end{itemize}
\medskip
\noindent\emph{Second and third steps.}
In the partially linear model, we set $S=(Y,D)$ and
$X=\emptyset$ in \eqref{orth2}, and define
\begin{equation} \label{eq:PLM_2step}
h_{0Y}(v)\equiv \mathbb{E}[Y\mid V=v],
\qquad
h_{0D}(v)\equiv \mathbb{E}[D\mid V=v],
\end{equation}
so that $h_0\equiv (h_{0Y},h_{0D})$,
$\Delta_2(g_0)=L_2(V)^{1+\dim(D)}$. Following the downstream local robustness principle, the identifying moment is \citep[cf.][]{robinson1988root}
\begin{equation}
\label{eq:PLM_general_moment}
m(W,g_{0},h_{0},\theta _{0})
=
\big( Y-h_{0Y}(V)-\theta
_{0}' ( D-h_{0D}(V))\big) \cdot (D - h_{0D}(V)).
\end{equation}
When $h_{0Y}$ and $h_{0D}$ are estimated by cross-fitted ML, the estimator based on \eqref{eq:PLM_general_moment} is the DML estimator of \cite{chernozhukov2018double}.
We now indicate two leading special cases.
\medskip
\noindent\textbf{(a) Hd-PS regression adjustment.}
Let $D$ be a binary treatment and $Y$ satisfy the potential-outcome model
$Y=Y_0 + D(Y_1-Y_0)$, where $Y_0$ and $Y_1$ are the potential outcomes under
control and treatment, respectively. Under strong ignorability,
$(Y_0, Y_1)\perp D\mid Z$, the propensity score $\mathbb{E}[D\mid Z]$ is a
balancing score \citep{rosenbaum1983central}, and classical matching and reweighting methods can be based on it
\citep{heckman1998matching,hirano2003efficient}.
With high-dimensional $Z$, we estimate the propensity score by Logit--Lasso (Hd-PS), where $\Delta_1$ is the mean-square limit of sparse linear combinations of a dictionary $(c_k)_{k=1}^\infty$ in $L_2(Z)$, and the generated regressor is $V=\varphi(D,Z,g_0)=\Lambda(g_0(Z))$. For expositional clarity, we consider that the generated regressor recovers the propensity score: $V=\Lambda(g_0(Z))=\mathbb{E}[D|Z]$. This simplifies the second step (now $h_0=h_{0Y}$, since $h_{0D}(V)=V$ is known), while it accommodates a rich set of controls through the high-dimensional first step. The estimand $\theta_0$ has a transparent causal (weighted-variance) interpretation because $V$ is a balancing score. In the general case of $V\neq\mathbb{E}[D|Z]$, our three-step inference procedure applies, though interpreting the target as causal requires additional assumptions.
\begin{comment}
\emph{Identification and overlap-weighted interpretation.}
Let $\tau(v)\equiv \mathbb{E}[Y_1-Y_0\mid \pi_0(Z)=v]$ denote the conditional average
treatment effect given the propensity score. The partially linear coefficient
$\theta_0$ can be written as an overlap-weighted average of $\tau(V)$:
\begin{equation*}
\theta_0
=\frac{\mathbb{E}\!\left[Var(D\mid V)\,\tau(V)\right]}{\mathbb{E}\!\left[Var(D\mid V)\right]}
=\frac{\mathbb{E}\!\left[V(1-V)\,\tau(V)\right]}{\mathbb{E}\!\left[V(1-V)\right]}.
\end{equation*}
The weighting kernel $V(1-V)$ emphasizes regions with good overlap and
automatically downweights propensity-score extremes, providing a stability
rationale for regression adjustment relative to estimators that rely on
$1/V$ or $1/(1-V)$ weighting.
\end{comment}
This regression-adjustment formulation is closely related to the estimating
equations in \citet{robins1992estimating}, which exploit the ``exposure residual''
$D-V$. In particular,
orthogonalized moment conditions can be based on products of
$(Y-\theta D)$ and $(D-V)$, yielding robustness to first-step estimation error.
Our three-step locally robust GMM framework constructs such orthogonal moments
for $\theta_0$ when $V$ is estimated by ML and $\kappa_0$ is
flexible, delivering valid inference in this widely used workflow. The debiased GMM estimator can be easily implemented as a (cross-fitted) Ordinary Least Squares (OLS) estimator, see equation~\eqref{DATE} and below.
A nonparametric version of the Hd-PS regression adjustment example is provided in Section~\ref{ATEPSM} of Appendix~\ref{sec:AdditionalExamples}. This example generalizes \cite{heckman1998matching}, \cite{hahn2013asymptotic}, and \cite{mammen2016semiparametric} to a machine-learning propensity score estimator and locally robust estimation and inference.
\medskip
\noindent\textbf{(b) Learned confounders via autoencoders and other embeddings.}
\smallskip
Researchers controlling for high-dimensional unstructured data (images, text, audio, or video) often employ a low-dimensional learned representation. Let
$e_0(Z)$ be such an embedding (e.g., the encoder from an autoencoder), let $g_0=(d_0,e_0)$, and define
\[
V = \varphi(D,Z,g_0) = e_0(Z).
\]
The second and third steps remain as in
\eqref{eq:PLM_2step}--\eqref{eq:PLM_general_moment}. We construct locally
robust estimators that account for the estimation of these embeddings in Section~\ref{sec:ex_PLMauto}.
There is a growing literature using deep latent-variable models and learned representations for causal inference. \citet{louizos2017causal} use variational autoencoders to learn latent confounders from proxy variables, while \citet{klaassen2024doublemldeep} and \citet{schulte2025adjustment} study treatment-effect estimation with multimodal or non-tabular data based on pre-trained or jointly trained neural-network representations within a DML framework. Surveys such as \citet{scholkopf2021toward} review causal representation learning more broadly. Unlike this literature, we treat the learned embedding $e_0(Z)$ as an ML-generated regressor and derive three-step locally robust GMM estimators that explicitly account for its estimation in downstream inference. Further details about this example, including the construction of locally robust estimators, are provided in Section~\ref{sec:ex_PLMauto}.
\end{ex}
\begin{comment}
\begin{ex}[Nonparametric ATE]
\label{sec:ATE_non}
Among existing causal methods that use machine learning to estimate the propensity score, twang (gradient boosting) is widely used in empirical work, particularly in combination with doubly robust augmented inverse probability weighted (AIPW) estimation; see \cite{leite2024machine}. However, the high sensitivity of AIPW to the estimation of the propensity score has been well documented in the literature; see, in particular, the recent analysis by \cite{yadlowsky2022explaining} based on high-dimensional asymptotics. Here, we propose an alternative nonparametric estimator of the ATE based on nonparametric regression on a gradient-boosting propensity score estimator. That is, the first step learns $V=g_0(Z)$ from a gradient-boosting nonparametric estimation of $D$ onto $Z$, corresponding to $\epsilon(W,g)=D-g(Z)$. The second step fits $h_0(d,v) = \mathbb{E}[Y\mid D=d, V=v]$, so $\Delta_2(g_0) = L_2(D,V)$.
We also investigate the case where the second step is misspecified, which does not fall under the influence function calculations of \cite{hahn2013asymptotic} and is quite relevant for applied work. Our proposal here relates to the $L_2$-boosting treatment effect estimation results of \cite{kueck2023estimation} and to the literature on balancing weights, see, e.g., \cite{athey2018approximate}, with the key difference that, in our setting, the estimated propensity score is an ML-generated regressor.
\end{ex}
\end{comment}
\begin{ex}[CASF with a control-function approach in a non-separable model]
\label{ex:CF}
We observe $W=(Y,D,Z)$ satisfying the model $Y=H(X,U)$, for an
unknown function $H$ and unobserved error term $U$. The main feature of this model is that $D$, a
component of $X$, may be an endogenous regressor. We assume that the
endogenous regressor satisfies $D=g_{0}(Z)+V$, with $U$ and $V$ being
unobserved correlated error terms. The function $g_{0}$ can be identified
by a conditional mean restriction, as in equation~\eqref{orth1} with $\epsilon(W, g_0) = D - g_0(Z)$. We assume a
control-function approach: $U\mid D,Z\sim U\mid X,V\sim U\mid V$, where $\sim$ denotes
equality in distribution. Thus, the generated regressor is a first-step residual:
\begin{equation*}
V\equiv \varphi (X,Z,g_{0}) = D-g_{0}(Z).
\end{equation*}
As in \citet{blundell2003endogeneity}, the control-function
assumption implies
\begin{align*}
\mathbb{E}[Y\mid X=x,V=v]& =\mathbb{E}[H(X,U)\mid X=x,V=v] \\
& =\mathbb{E}[H(x,U)\mid V=v]\equiv h_{0}(x,v).
\end{align*}
This defines the second step, which satisfies \eqref{orth2} with $\Delta
_{2}(g_{0})=L_{2}(X,V)$.
The control-function assumption allows us to identify the Average Structural
Function (ASF) at a point $x\in \mathbb{R}^{\operatorname{dim}(X)}$:
\begin{equation*}
\operatorname{ASF}_{0}(x)\equiv \mathbb{E}[H(x,U)]=\mathbb{E}[\mathbb{E}[H(x,U)\mid V]]=
\mathbb{E}[h_{0}(x,V)].
\end{equation*}
Some well-known conditions on the support of the random vectors are needed
for the above equation to hold
\citep[see][]{blundell2004endogeneity,imbens2009identification}.
In this setup, a parameter of interest is the CASF, given by
\begin{equation*}
\theta _{0}=\int \operatorname{ASF}_{0}(x^{\ast })dF^{\ast }(x^{\ast }),
\end{equation*}
for a counterfactual distribution $F^{\ast }$. When $F^{\ast }$ is implied
by a certain policy, the CASF may be used to measure the effect of the
policy
\citep[see][]{stock1989nonparametric,stock1991nonparametric,blundell2004endogeneity}.
By Fubini's Theorem, the CASF can be written as
\begin{equation*}
\theta _{0}=\int \mathbb{E}[h_{0}(x^{\ast },D-g_{0}(Z))]dF^{\ast }(x^{\ast
})=\mathbb{E}\left[ \int h_{0}(x^{\ast },D-g_{0}(Z))dF^{\ast }(x^{\ast })
\right] .
\end{equation*}
Hence, the moment function that identifies the CASF is:
\begin{equation}
\label{eq:CASF_moment}
m(w,g_0,h_0,\theta_0 )=\int h_0(x^{\ast },d-g_0(z))dF^{\ast }(x^{\ast })-\theta_0 .
\end{equation}
We propose in (\ref{eq:CASF_D}) a debiased estimator for the CASF, allowing for and accounting for ML first and second steps. The CASF estimator here generalizes the automatic debiased estimator of trained regression averages under covariate shifts, as in \citet{chernozhukov2023automatic}, to the empirically important case in which the covariate shift arises through an endogenous regressor.
A remarkable feature of the CASF example is that, to evaluate the moment condition at a point $w=(y,d,z)$, one needs the entire second-step nuisance function $h_0$. Thus, even in the low-dimensional case, it is not encompassed by the setup of \citet{hahn2013asymptotic,hahn2019three} (nor are Examples~\ref{ex:PLM_ML}(a)--(b)).
\end{ex}
\section{Debiased estimation with ML-generated regressors}
\label{sec:debiased_est}
A fundamental property that allows us to develop debiased estimators is Neyman orthogonality, also referred to as local robustness \citep[see][]{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 \citet{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 \emph{direct} (or \emph{evaluation}) effect (in the toy example, $\beta_0 v$ evaluated at $v=V(g)$) and an \emph{indirect} (or \emph{conditioning}) effect that operates through $h_0$ (in the toy example, $g\mapsto h(g)(v)=\beta_g v$); see Figure~\ref{fig:chainrule} and Section~\ref{sec:first_second_IF}. 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
\begin{equation}\label{orth}
\begin{aligned}
\psi (w,g_{0},h_{0},\alpha _{0},\theta ) &= m(w,g_{0},h_{0},\theta ) \\
&\quad+ \underbrace{\alpha _{01}(z)\cdot
\epsilon (w,g_{0})}_{=\phi _{1}(w,g_{0},\alpha _{01})} \\
&\quad+ \underbrace{
\alpha _{02}(x,\varphi (d,z,g_{0}))^{\prime }[s-h_{0}(x,\varphi (d,z,g_{0}))]
}_{=\phi _{2}(w,g_0,h_{0},\alpha _{02})},
\end{aligned}
\end{equation}
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 \eqref{orth} 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{Parfirststep} for a formal statement).} The explicit analytic expression for $\alpha_{01}$ is typically complicated (see equation~\eqref{eq:alpha1_defi} in Appendix~\ref{sec:Proofs}), but we construct automatic estimators that do not require this expression. The second-step IF $\phi _{2}$ is of the usual form \citep{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 \emph{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 \emph{downstream local robustness}. It generalizes an observation in \citet{hahn2013asymptotic} to a general three-step ML framework and to a much broader class of problems beyond generated regressors (cf.\ Proposition~\ref{prop:zero_indirect}).
\subsection{The debiased estimator}
Automatic debiased estimation with generated regressors is based on the moment condition in equation~\eqref{orth}, where the Riesz representers $\alpha_{01}$ and $\alpha_{02}$ are estimated automatically (see Section~\ref{sec:est_auto}). We construct sample analogues using cross-fitting, as in \citet{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
\begin{equation*}
\hat{\psi}(\theta )\equiv \frac{1}{n}\sum_{\ell =1}^{L}\sum_{i\in I_{\ell }}
\hat{\psi}_{i\ell }(\theta ),
\end{equation*}
with
\begin{equation}
\hat{\psi}_{i\ell }(\theta )\equiv m(W_{i},\hat{g}_{\ell },\hat{h}_{\ell
},\theta )+\hat{\alpha}_{1\ell }(Z_{i})\cdot \epsilon(W_{i},\hat{g}_{\ell })+
\hat{\alpha}_{2\ell }(X_{i},\hat{V}_{i\ell })^{\prime } (S_{i}-\hat{h}
_{\ell }(X_{i},\hat{V}_{i\ell })), \label{eq:dmon_i}
\end{equation}
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
\begin{equation} \label{dgmm}
\hat{\theta}=\operatornamewithlimits{argmin}_{\theta\in\Theta}\hat{\psi}
(\theta)^{\prime}\hat{\Upsilon}\hat{\psi}(\theta),
\end{equation}
where $\hat{\Upsilon}$ is a positive semi-definite weighting matrix of
dimension $q\times q$. Under regularity conditions (see
Section~\ref{sec:asymptotic}), $\hat{\theta}$ is asymptotically normal with the usual GMM asymptotic variance.
\begin{ex}[continues=ex:PLM_ML]
\label{ex:ATE_PLM1}
We illustrate the construction of a three-step debiased and cross-fitted ML estimator for the partially linear model in Example~\ref{ex:PLM_ML}. For expositional clarity, we consider $\operatorname{dim}(D)=1$. Suppose we have a generated regressor $\hat{V}_{i\ell} = \varphi(D_i, Z_i, \hat{g}_\ell)$. Within each fold, the second step estimates $h_{0Y}(v)=\mathbb{E}[Y\mid V=v]$ and $h_{0D}(v)=\mathbb{E}[D\mid V=v]$ by regressing $Y_i$ and $D_i$ on $\hat V_{i\ell}$ using a dictionary $\mathbf{b}_{J}(v)=(b_{1}(v),\ldots,b_{J}(v))^{\prime }$; for instance, one may take $b_j(v)=v^{j-1}$ and conduct an $\ell_1$-penalized least squares. This yields $\hat{h}_{\ell,Y}$ and $\hat{h}_{\ell,D}$ and defines $\hat h_\ell=(\hat h_{\ell,Y},\hat h_{\ell,D})$.
The partially linear moment in \eqref{eq:PLM_general_moment} is orthogonal with respect to $h$, implying $\alpha_{02}=0$. Therefore, a natural application of the DML estimator of \citet{chernozhukov2018double} yields the closed-form expression
\begin{equation}
\label{eq:ATE_PLM_DML}
\hat{\theta}_{DML} = \frac{\sum_{\ell =1}^{L}\sum_{i\in I_{\ell }}
\left( Y_i - \hat{h}_{\ell,Y}(\hat{V}_{i\ell})\right)\left( D_{i}-\hat{h}_{\ell,D}(\hat{V}_{i\ell})\right)}{\sum_{\ell
=1}^{L}\sum_{i\in I_{\ell }}\left( D_{i}-\hat{h}_{\ell,D}(\hat{V}_{i\ell})\right) ^{2}}.
\end{equation}
This estimator is locally robust with respect to the second step but does not account for the effect of learning the generated regressor $\hat V_{i\ell}$.
Using the general debiased moment \eqref{orth}, the three-step debiased estimator takes the form
\begin{equation}
\label{DATE}
\hat{\theta}=\hat{\theta}_{DML} + \frac{\sum_{\ell =1}^{L}\sum_{i\in I_{\ell }}
\hat{\alpha}_{1\ell }(Z_{i})\cdot \epsilon(W_{i},\hat{g}_{\ell})}{\sum_{\ell
=1}^{L}\sum_{i\in I_{\ell }}\left( D_{i}-\hat{h}_{\ell,D}(\hat{V}_{i\ell})\right) ^{2}}.
\end{equation}
Since the partially linear moment is orthogonal to $h$, downstream local robustness implies that the indirect contribution of the generated regressor through the second step vanishes. Hence, the first-step Riesz representer $\alpha_{01}$ also simplifies considerably. We construct an automatic cross-fitted estimator $\hat\alpha_{1\ell}$ in Section~\ref{sec:auto_1step}, with a special case provided below in equation~\eqref{eq:alpha1ell_riesz_lasso}.
\medskip
\noindent\textbf{Hd-PS regression adjustment.}
Let $\mathbf{c}_K(z) \equiv (c_1(z), \dots, c_K(z))'$ be a dictionary with $K$ atoms (in the high-dimensional case, this can simply collect the regressors $z_j$). The first step estimates the propensity score $\mathbb{E}[D\mid Z]$ via Lasso--Logit: $\hat{g}_\ell(z) =\mathbf{c}_K(z)'\widehat{\boldsymbol{\gamma}}_{K\ell}$, where
\begin{equation*}
\widehat{\boldsymbol{\gamma}}_{K\ell}
=
\operatornamewithlimits{argmin}_{\boldsymbol{\gamma }_{K}\in \mathbb{R}^{K}}\left\{ -\sum_{i\notin I_{\ell }} \left[ D_{i}\log\Lambda(\mathbf{c}_{K}(Z_{i})^{\prime }\boldsymbol{\gamma }_{K})+(1-D_{i})\log(1-\Lambda(\mathbf{c}_{K}(Z_{i})^{\prime }\boldsymbol{\gamma }_{K})) \right]+\lambda \lVert \boldsymbol{\gamma }_{K}\rVert
_{1}\right\},
\end{equation*}
$\lVert \cdot\rVert _{1}$ is the $\ell_1$ norm, and $\lambda$ is a penalization parameter. The score of this problem leads to the orthogonality condition in equation~\eqref{orth1} with $\epsilon(W,g_0)=D-\Lambda(g_0(Z))$ and $\Delta_1$ the mean-square limit of sparse linear combinations of $(c_k)_{k=1}^\infty$. The generated regressor is the estimated propensity score
\[
\hat{V}_{i\ell}=\Lambda(\hat{g}_{\ell }(Z_{i})).
\]
In the Hd-PS regression adjustment, estimation of the second step simplifies to $\hat{h}_{\ell,D}(\hat{V}_{i\ell})=\hat{V}_{i\ell}$, since $\mathbb{E}[D\mid V]=V$. We provide a simple weighted Lasso estimator $\hat{\alpha}_{1\ell }$ as follows:
\begin{equation}\label{eq:alpha1ell_riesz_lasso}
\hat{\alpha}_{1\ell}(z)\equiv\mathbf{c}_K(z)'\widehat{\boldsymbol{\beta}}_{K\ell},\quad
\widehat{\boldsymbol{\beta}}_{K\ell}\in\operatornamewithlimits{argmin}_{\boldsymbol{\beta}_K\in\mathbb{R}^K}\left\{\sum_{\ell'\neq\ell}\sum_{i\in I_{\ell'}}\omega_{i\ell\ell'}\bigl(\mathcal{E}_{i\ell\ell'}-\mathbf{c}_K(Z_i)'\boldsymbol{\beta}_K\bigr)^2+\lambda\|\boldsymbol{\beta}_K\|_1\right\}.
\end{equation}
where the ``dependent'' variable is $\mathcal{E}_{i\ell\ell'}=-[Y_{i}-\hat{h}_{\ell \ell
^{\prime }}(\hat{V}_{i\ell \ell ^{\prime }})]$ and the weights are $\omega_{i\ell\ell'} \equiv \hat{V}_{i\ell \ell^{\prime }} (1-\hat{V}_{i\ell \ell^{\prime }})$. Here, $\hat{h}_{\ell \ell^{\prime }}$ and $\hat{V}_{i\ell \ell^{\prime }}$ are estimators that use only observations not in $I_\ell\cup I_{\ell^{\prime }}$. The arguments leading to this construction are detailed in Section~\ref{sec:ex_ATE_PLM}. In general, $\alpha_{01}\neq 0$, and hence, without our correction, inference is generally invalid. The DML estimator $\hat{\theta}_{DML}$ must be debiased to obtain an estimator that is locally robust to the generated propensity score. This estimator $\hat{\theta}$ can be easily implemented as a (cross-fitted) OLS of $Y_i - \hat{h}_{\ell,Y}(\hat{V}_{i\ell})+\hat{\alpha}_{1\ell }(Z_{i})$ on $D_{i}-\hat{V}_{i\ell}$.
\end{ex}
\begin{ex}[continues=ex:CF]
The moment condition defining the CASF is not orthogonal to the second step, so a debiasing term for each step is needed. We illustrate how to build a three-step debiased estimator of the CASF, starting from a plug-in estimator and comparing it to the natural extension of the DML estimator.
The first step recovers the control function $\hat{V}_{i\ell} = D_i - \hat{g}_\ell(Z_i)$, with $\hat{g}_\ell(z) =\mathbf{c}_K(z)'\widehat{\boldsymbol{\gamma}}_{K\ell}$ being a Lasso fit of $D_i$ on a dictionary $\mathbf{c}_K(Z_i)$. The second step estimates $h_0(x,v) = \mathbb{E}[Y\mid X=x, V=v]$ via Lasso. For a dictionary $\mathbf{b}_J(x, v)$ with $J$ atoms, we get $\hat{h}_\ell(x, v) = \mathbf{b}_J(x, v)'\widehat{\boldsymbol{\eta}}_J$ by $L_1$-penalized least squares of $Y_i$ on the dictionary evaluated at $X_i$ and the generated control function $\hat{V}_{i\ell}$.
To estimate the CASF according to equation~\eqref{eq:CASF_moment}, we compute the integral by Monte Carlo integration, since the counterfactual distribution $F^*$ is fixed by the researcher. Let $(X_s^*)_{s=1}^S$ be a sample from $F^*$, independent of the original sample ($S\gg n$). The cross-fitted plug-in estimator for the CASF is:
\begin{equation} \label{eq:CASF_plugin}
\hat{\theta}_{PI}=\frac{1}{nS}\sum_{\ell =1}^{L}\sum_{i\in I_{\ell
}}\sum_{s=1}^{S}\hat{h}_{\ell }(X_{s}^{\ast },\hat{V}_{i\ell }).
\end{equation}
A DML estimator that accounts for estimation of $h_0$ in the second step, but not for the generated control function, it is given by
\begin{equation} \label{eq:CASF_DML}
\hat{\theta}_{DML}= \hat{\theta}_{PI} + \frac1n \sum_{\ell =1}^{L}\sum_{i\in I_{\ell}} \hat{\alpha}_{2\ell}(X_i,\hat{V}_{i\ell}) \cdot (Y_i - \hat{h}_\ell(X_i, \hat{V}_{i\ell})).
\end{equation}
The three-step debiased estimator accounts for the generated regressor by extending the moment condition:
\begin{equation} \label{eq:CASF_D}
\hat{\theta}= \hat{\theta}_{DML} + \frac1n \sum_{\ell =1}^{L}\sum_{i\in I_{\ell}} \hat{\alpha}_{1\ell}(Z_i) \cdot (D_i - \hat{g}_\ell(Z_i)).
\end{equation}
Automatic estimation of the Riesz representers $\alpha_{01}$ and $\alpha_{02}$ is detailed in Section~\ref{sec:CASF_details}.
\end{ex}
\subsection{Monte Carlo simulations}
\label{sec:MC}
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 \citet{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{sec:asymptotic}). A detailed description of the setups, estimation procedures, and results is provided in Appendix~\ref{sec:MC_details}.
\subsubsection{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:
\begin{align*}
Y &= D + Z_1 + Z_2 + \varepsilon, \\
D &= \mathbf{1}\left( C \nu \leq Z_1 + Z_2 + Z_3 + Z_4 + Z_5 + Z_6 \right).
\end{align*}
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{tab:MC_ATE} in Appendix~\ref{sec:MC_details}). 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{fig:histogram_ATE} 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$.
\begin{figure}[h!]
\centering
\begin{subfigure}{.33\textwidth}
\centering
\includegraphics[width=\textwidth]{Plots/samedgp_logit_lin_unif_dep_subset_n100_p10.png}
\caption{$n=100$.}
\end{subfigure}\hfill
\begin{subfigure}{.33\textwidth}
\centering
\includegraphics[width=\textwidth]{Plots/samedgp_logit_lin_unif_dep_subset_n500_p10.png}
\caption{$n=500$.}
\end{subfigure}\hfill
\begin{subfigure}{.33\textwidth}
\includegraphics[width=\textwidth]{Plots/samedgp_logit_lin_unif_dep_subset_n1000_p10.png}
\caption{$n=1000$.}
\end{subfigure}
\caption{Histograms of the estimators for $\theta_0$ in a partially linear model framework (logistic $\nu$) with Hd-PS first-step. Number of replications is $2000$. \textit{DML} = Double/Debiased Machine Learning estimator, \textit{3SD} = Three-Step Debiased estimator. Note that the DML estimator equals the plug-in estimator.}
\label{fig:histogram_ATE}
\end{figure}
Table~\ref{tab:MC_ATE} in Appendix~\ref{sec:MC_details} 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$.
\subsubsection{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:
\begin{align*}
Y = \sum_{k=1}^5 Z_k + 2D + U \text{ and } D = \sum_{k=1}^6 Z_k + V.
\end{align*}
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{tab:MC_CASF} in Appendix~\ref{sec:MC_details} and Figure~\ref{fig:histogram_CASF} 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.
\begin{figure}[h!]
\centering
\begin{subfigure}{.33\textwidth}
\centering
\includegraphics[width=\textwidth]{Plots/L_MODEL_glin_hlin_n100.png}
\caption{$n=100$.}
\end{subfigure}\hfill
\begin{subfigure}{.33\textwidth}
\centering
\includegraphics[width=\textwidth]{Plots/L_MODEL_glin_hlin_n500.png}
\caption{$n=500$.}
\end{subfigure}\hfill
\begin{subfigure}{.33\textwidth}
\includegraphics[width=\textwidth]{Plots/L_MODEL_glin_hlin_n1000.png}
\caption{$n=1000$.}
\end{subfigure}
\caption{Histograms of the CASF estimators with a Control-Function Approach. Number of replications is $2000$. \textit{PI} = Plug-in estimator, \textit{DML} = Double/Debiased Machine Learning estimator, and \textit{3SD} = Three-Step Debiased estimator.}
\label{fig:histogram_CASF}
\end{figure}
Results regarding coverage also differ from those of the previous example (see Table~\ref{tab:MC_CASF} in Appendix~\ref{sec:MC_details}). 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{sec:Examples}, where additional details about the examples are gathered.
\subsection{First- and second-step influence functions}
\label{sec:first_second_IF}
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
\citep[see][]{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
\begin{equation*}
\mathbb{E}_{F}[\delta _{1}(Z)\epsilon (W,g)]=0\text{ for all }\delta _{1}\in
\Delta _{1}.
\end{equation*}
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
\begin{equation}
\mathbb{E}_{F}\left[ \delta _{2}(X,V(g))\{S-h(X,V(g))\}\right] =0\text{ for
all }\delta _{2}\in \Delta _{2}(g), \label{eq:ortho_genFg}
\end{equation}
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))$.
\begin{ex}[continues=ex:PLM_ML]
\label{ex:PLM_ortho}
In the partially linear model with generated regressors, $h$ has two components,
$h=(h_{Y},h_{D})$. For $S$ equal to $Y$ or $D$, denote $h_{S}(F,g)(v)=\mathbb{E}
_{F}[S\mid V(g)=v]$. The first step $g$ enters each second step $h_{S}$ in
two ways: (i) indirectly, through the conditioning variable $V(g)$, and (ii)
directly, when we evaluate $v$ at $V(g)$. Following our notation,
$h_{S}(F_{0},g_{0})(V(g_{0}))$ simplifies to $h_{0S}(V)$, and
$h_{0}=(h_{0Y},h_{0D})$.
\end{ex}
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 \emph{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
\begin{equation*}
\bar{m}(g,h) \equiv \mathbb{E}[m(W,g,h,\theta_0)].
\end{equation*}
We study these effects separately. By the chain rule,
\begin{align}
\frac{d}{d\tau }\bar{m}(g(F_{\tau }),h(F_{\tau },g(F_{\tau })))& =
\frac{d}{d\tau }\bar{m}(g(F_{\tau }),h(F_{0},g(F_{\tau })))
\label{eq:FSD} \\
& \quad +\frac{d}{d\tau }\bar{m}(g_{0},h(F_{\tau },g_{0})), \label{eq:SSD}
\end{align}
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 \citet{hahn2013asymptotic}, the first
step affects the moment condition in two ways (see Figure~\ref{fig:chainrule}). We have a \textit{direct impact} on $\bar{m}$, quantified by the derivative of $\bar{m}(g(F_{\tau }),h_{0})$. This direct impact includes the \textit{effect of evaluating $h$} at the generated regressor. We also have an
\textit{indirect effect} on the moment that arises because $g$ affects
\emph{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 \eqref{eq:FSD}. The derivative in
\eqref{eq:SSD} accounts for the effect of the second step. This effect is
independent of the first step and therefore treats $g_{0}$ as known.
\begin{figure}[h!]
\centering
\begin{tikzpicture}[->,>=stealth',auto,node distance=2.6cm,
thick]
\node (par) {$\tau$};
\node (dist) [align=center, right of=par] {$F_\tau$};
\node (g) [align=center,right of=dist] {$g(F_\tau)$};
\node (h) [align=center, right of=g] {$h(F_\tau, g(F_\tau))$};
\node (exp) [align=center, below of=h] {$\bar{m}(g(F_\tau), h(F_\tau, g(F_\tau)))$};
\draw [->] (par) -- (dist);
\draw [->] (dist) -- (g);
\draw [->] (g) -- (h) node [midway, fill=white, above=0.25em] {(I)};
\draw [->] (h) -- (exp) node [midway, fill=white, right = 0.25em] {(2S)};
\draw [->] (g) -- (exp) node [midway, fill=white, above right =0.01em and 0.01em] {(D)};
\draw [->] (dist) to [out=90,in=90] (h);
\end{tikzpicture}
\caption{The effect of a deviation $F_\protect\tau$ on the moment condition.
(2S) represents the second-step effect. (D) represents the direct effect of
the first step. The path (I)-(2S) represents the indirect estimation effect
of the first step.}
\label{fig:chainrule}
\end{figure}
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$,
\begin{align}
\frac{d}{d\tau }\bar{m}(g(F_{\tau }),h(F_{0},g(F_{\tau })))& =\int
\phi _{1}(w,g_{0},\alpha _{01})\,dH(w) \text{ and} \label{eq:IF1} \\
\frac{d}{d\tau }\bar{m}(g_{0},h(F_{\tau },g_{0}))& =\int \phi
_{2}(w,g_0, h_{0},\alpha _{02})\,dH(w). \label{eq:IF2}
\end{align}
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
\begin{equation*}
\psi (w,g,h,\alpha ,\theta )\equiv m(w,g,h,\theta )+\phi _{1}(w,g,\alpha
_{1}) + \phi _{2}(w,g,h,\alpha _{2})
\end{equation*}
are orthogonal/locally robust/debiased. The next section shows that, under some conditions, the IFs have the form displayed in equation~\eqref{orth}. 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.
\begin{prop}[Downstream local robustness]
\label{prop:zero_indirect}
Assume that $\bar m$ is Hadamard differentiable in the second step $h$ at $(g_0,h_0)$, with derivative $D_{02}$, and that $h(F, g)$ is Hadamard differentiable in $g$ at $(F_0, g_0)$. If the moment is locally robust with respect to the second step, i.e., $D_{02} = 0$,
then, for any regular path $\tau \mapsto F_\tau$ through $F_0$,
\begin{equation*}
\frac{d}{d\tau}\bar m\big(g(F_0),h(F_0,g(F_\tau))\big)
= 0.
\end{equation*}
\end{prop}
\begin{rem}[Comparison with two-step approaches]
\label{twostep}
The dependence of the second step on the first step makes the results in \cite{chernozhukov2022locally} not applicable in our setting. In particular, we require that for each $g\in \Delta _{1}$, $h(X,V(g))\in \Delta _{2}(g)$, so the parameter space of the second step depends on the first step. In a two-step setting with multiple parameters, one would need $(g,h)$ to live in a linear product space \citep[see Theorem~3 in][]{chernozhukov2022locally}. The dependence of $\Delta_2(g)$ on $g$ breaks this structure. Our results rely on alternative assumptions (see Assumption~\ref{ass:inclusion}) to handle this parameter space.
\end{rem}
\section{Automatic estimation of the Riesz representers}
\label{sec:Automatic}
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~
\eqref{eq:alpha1_defi} in Appendix \ref{sec:Proofs}). 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.
\subsection{Separate identification of Riesz representers}
\label{sec:ident_alphas}
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 \cite{chernozhukov2022locally} to each step leads to
\begin{align}
\frac{d}{d\tau }\mathbb{E}[\psi _{1}(W,g_{0}+\tau \delta _{1},\alpha
_{01},\theta_0 )]& =0\text{ for all }\delta _{1}\in \Delta _{1}\text{ and}
\label{eq:auto_momentF} \\
\frac{d}{d\tau }\mathbb{E}[\psi _{2}(W,h_{0}+\tau \delta _{2},\alpha
_{02},\theta_0 )]& =0\text{ for all }\delta _{2}\in \Delta _{2}(g_{0}),
\label{eq:auto_momentS}
\end{align}
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 \cite{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{eq:auto_momentF}) and (\ref{eq:auto_momentS}), 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 \cite{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
{sec:Linearization} provides the formal development. For a detailed construction of the automatic estimators, we refer to Section~
\ref{sec:est_auto}.
\subsection{First- and second-step linearization}
\label{sec: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 \citep[see,
e.g.,][Equation~4.1]{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:
\begin{ass}
\label{ass:Diff_m_h} There exists a function $D_{02}(w,h)$ such that
\begin{equation*}
\frac{d}{d\tau}\bar{m}(g_0,h_\tau)=\frac{d}{d\tau}\mathbb{E}[D_{02}(W,h_\tau)].
\end{equation*}
Moreover, $h\mapsto \mathbb{E}[D_{02}(W,h)]$
is linear and continuous in $L_2(X,V)$.
\end{ass}
The same assumption has been considered in \cite{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:
\begin{prop}
{\label{prop:second_step}} Under Assumption~\ref{ass:Diff_m_h}, there exists an $\alpha _{02}\in \Delta
_{2}(g_{0})$, given by the orthogonal projection of $r_{02}$ onto $\Delta_{2}(g_{0})$, such that the function
\begin{equation*}
\phi _{2}(w,g_0,h_{0},\alpha _{02})=\alpha _{02}(x,\varphi
(d,z,g_{0}))^{\prime }\{s-h_{0}(x,\varphi (d,z,g_{0}))\},
\end{equation*}
satisfies equation~\eqref{eq:IF1} and is thus the second-step IF.
\end{prop}
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{sec:est_auto}).
We now move to the more complicated linearization of the first-step effect. Note that if the chain rule can be applied:
\begin{equation}\label{eq:1st_step_decomposition}
\begin{aligned}
\frac{d}{d\tau} \bar{m}(g(F_\tau),h(F_0, g(F_\tau))) &= \frac{d}{d\tau
} \bar{m}(g(F_\tau),h_0) \\
&+\frac{d}{d\tau} \bar{m}(g_0,h(F_0, g(F_\tau))).
\end{aligned}
\end{equation}
The first derivative in the RHS can be easily analyzed if we linearize $\bar{
m}(g,h_{0})$ in $g$:
\begin{ass}
\label{ass:Diff_m_g} There exists a function $D_{dir}(w,g)$ such that
\begin{equation*}
\frac{d\bar{m}(g_\tau,h_0)}{d\tau}=\frac{d\mathbb{E}
[D_{dir}(W,g_\tau)]}{d\tau}.
\end{equation*}
Moreover, $g\mapsto \mathbb{E}
[D_{dir}(W,g)]$ is linear and continuous in $L_{2}(Z)$.
\end{ass}
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 \cite{hahn2013asymptotic} to allow for ML second
steps as in equation~\eqref{orth2}. The lemma is stated for one-dimensional $S$. For higher dimensions, it must be applied component-wise.
\begin{lma}
{\label{lma:generalized_HR}} Assume that the chain rule can be applied along
the path $\tau \mapsto g_{\tau }$. Then, for every $\delta _{2}\in L_{2}(X,V)
$ satisfying that there exists an $\varepsilon >0$ such that $\delta _{2}\in
\cap _{\tau <\varepsilon }\Delta _{2}(g_{\tau })$:
\begin{equation*}
\frac{d}{d\tau }\mathbb{E}[\delta _{2}(X,V)\cdot h(F_{0},g_{\tau })(X,V)]= \frac{d}{d\tau }\mathbb{E}[\delta_{2}(X,V(g_{\tau})) \cdot (S-h_{0}(X,V(g_{\tau})))]
\end{equation*}
\end{lma}
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~\eqref{eq:ortho_genFg}
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{lma:generalized_HR}. 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.
\begin{ass}
\label{ass:inclusion} For every path $\tau \mapsto g_\tau$ there exits an $
\varepsilon>0$ such that \mbox{} \\[-20pt]
\begin{enumerate}[label=\textbf{\alph*.},ref=\ref{ass:inclusion}.\alph*]
\item \label{ass:inclusion:alpha} $\alpha_{02}\in \cap_{\tau<\varepsilon}\Delta_2(g_\tau)$, and
\item \label{ass:inclusion:h} $h(F_0, g_\tau)\in \Delta_2(g_0)$ for all $0\leq\tau<\varepsilon$.
\end{enumerate}
\end{ass}
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{sec:app_inclusion} for a detailed discussion). What
Assumption~\ref{ass:inclusion} 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{lma:generalized_HR}, 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:
\begin{ass}
\label{ass:diff_h_alpha_phi} $h_0(x,v)$ and $\alpha_{02}(x,v)$ are almost surely differentiable w.r.t. $v$ with square-integrable derivatives. Moreover, the mapping $g\mapsto \varphi(d,z,g)$, from $L_2(Z)$ to $L_2(D,Z)$
, is Hadamard differentiable at $g_0$, with derivative $D_\varphi$.
\end{ass}
The Hadamard derivative of $\varphi $ is a linear and continuous map $
D_{\varphi }\colon L_{2}(Z)\rightarrow L_{2}(D,Z)$ such that
\begin{equation*}
\frac{d}{d\tau }\varphi (d,z,g_{\tau })=\frac{d}{d\tau }D_{\varphi }g_{\tau
}.
\end{equation*}
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 \citep[cf.][Ass.~2 and the discussion below]{ichimura2022influence}.
\begin{ass}
\label{rho}
The mapping $g \mapsto \mathbb{E}[\epsilon(W, g)]$ is Hadamard differentiable at $g_0$. The Riesz representer of the derivative ($r_e$) satisfies $r_{e}(z)<0$ and is bounded and bounded away from zero.
\end{ass}
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:
\begin{thm}
{\label{thm:first_step}} Under Assumptions \ref{ass:Diff_m_h}-\ref{rho}:
\begin{itemize}
\item[\textsc{\textbf{(Lin)}}] The function
\begin{equation}
{\label{eq:D1_defi}}D_{01}(w,g)\equiv D_{dir}(w,g)+\frac{
\partial }{\partial v}\left[ \alpha_{02}(x,v)' (s-h_{0}(x,v))
\right] \cdot D_{\varphi }g,
\end{equation}
where the derivative is evaluated at $v=\varphi (d,z,g_{0})$, satisfies
\begin{equation*}
\frac{d}{d\tau }\bar{m}(g(F_{\tau }),h(F_{0},g(F_{\tau })))=\frac{d}{
d\tau }\mathbb{E}[D_{01}(W,g(F_{\tau }))].
\end{equation*}
\item[\textsc{\textbf{(IF)}}] There exists an $\alpha _{01}\in \Delta _{1}$, given by equation~\eqref{eq:alpha1_defi} in Appendix \ref{sec:Proofs}, such that the function
\begin{equation*}
\phi _{1}(w,g_{0},\alpha _{01})=\alpha _{01}(z)\cdot \epsilon
(w,g_{0}),
\end{equation*}
satisfies equation~\eqref{eq:IF2} and is thus the first-step IF.
\end{itemize}
\end{thm}
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~\eqref{eq:D1_defi}).
The first term corresponds to the linearization of the \textit{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 \textit{
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{sec:auto_1step} 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.
\begin{rem}[Relation to Hahn and Ridder (2013)]\label{rem:HR}
Theorem~5 in \cite{hahn2013asymptotic} studies a three-step generated-regressor problem with
$V=g_0(Z)$ (i.e.\ $\varphi(d,z,g)=g(z)$), $\Delta_1=L_2(Z)$, $\Delta_2=L_2(X,V)$, and third-step moment
\[
m(w,g,h,\theta)=\eta\!\left(w,\,h(x,g(z))\right)-\theta,
\]
so that $(g,h)$ enters through the scalar $h(x,g(z))$.
In our notation, their second-step derivative corresponds to $D_{02}$. If $\eta$ is differentiable in its second argument,
\[
D_{02}(w,\delta)
=
\frac{\partial \eta}{\partial y}\!\left(w,h_0(x,v)\right)\delta(x,v),
\qquad v=g_0(z),
\]
with Riesz representer
\[
\alpha_{02}(x,v)
=
\mathbb{E}\!\left[\left.\frac{\partial \eta}{\partial y}\!\left(W,h_0(X,V)\right)\right|X=x,V=v\right].
\]
Let $r_{02}(w)=\partial\eta/\partial y\bigl(w,h_0(x,v)\bigr)$, so that
$\alpha_{02}(x,v)=\mathbb{E}[r_{02}(W)\mid X=x,V=v]$.
In this evaluation-functional setting the first-step derivative is separable,
\[
D_{01}(w,g)=r_{01}(w)\,g(z),
\]
where
\[
r_{01}(w)
=
\bigl(y-h_0(x,v)\bigr)\,\partial_v\alpha_{02}(x,v)
+
\bigl(r_{02}(w)-\alpha_{02}(x,v)\bigr)\,\partial_v h_0(x,v).
\]
Since $\Delta_1=L_2(Z)$, the first-step Riesz representer is
\[
\alpha_{01}(z)=\mathbb{E}[r_{01}(W)\mid Z=z].
\]
As noted by \citet[Remark~3]{hahn2013asymptotic}, if $\alpha_{02}=0$ then
$r_{01}(w)=r_{02}(w)\,\partial_v h_0(x,v)$, illustrating a special case of downstream local robustness.
Relative to \cite{hahn2013asymptotic}, our results:
(i) allow general generated regressors $V=\varphi(D,Z,g_0)$ for Hadamard differentiable $\varphi$;
(ii) allow general Hadamard differentiable functionals $\bar m(g,h)=\mathbb{E}[m(W,g,h,\theta_0)]$, not only evaluation functionals; and
(iii) allow more general spaces $\Delta_1$ and $\Delta_2(g)$ (e.g.\ sparse or sieve structures).
Beyond influence functions, we provide an automatic Riesz-representer implementation suitable for high-dimensional/non-Donsker ML estimators and establish asymptotic normality of the resulting debiased estimators.
\end{rem}
\begin{rem}[Other First Steps]
\label{OtherFS}
There are examples of first steps that are not included in (\ref{orth1}), such as some parametric estimators, functions
identified by orthogonality conditions with instruments, where $g$ depends on
other variables different from $Z$, or the control function approach of
\citet{imbens2009identification}, among others. Nevertheless, much of our analysis is still useful for these other cases. In particular, the expression for $D_{01}(w,g)$ remains the same, and our results can be readily extended to other first steps by characterizing the corresponding first step IF $\phi_{1}$ such that
\begin{equation*}
\frac{d}{d\tau }\mathbb{E}[D_{01}(W,g_{\tau })]=\int \phi _{1}(w,g_{0},\alpha _{01})dH(w).
\end{equation*}
We illustrate the application of this equation with parametric first steps in the next remark.
\end{rem}
\begin{rem}[General Parametric First Steps]
\label{Parfirststep} Suppose we replace our definition of $g_{0}$ and $\Delta_{1}$ in (\ref{orth1}) by a generic parametric fit $g_{0}(z)=G(z,\zeta_{0})$, where $G$ has a known functional form and $\zeta_{0}$ is an unknown finite-dimensional parameter in a parameter space $B\subseteq \mathbb{R}^{\operatorname{dim}(\zeta)}$. We allow for $\zeta_{0}$ to be identified by parametric or semiparametric restrictions. For example, this setting includes semiparametric estimators (e.g., single-index models with $g_{0}(z)=z'\zeta_{0}$). Let $\hat{\zeta}$ be a regular estimator for $\zeta_{0} $ satisfying
\begin{equation*}
\sqrt{n}(\hat{\zeta}-\zeta_{0})=\frac{1}{\sqrt{n}}\sum_{i=1}^{n}\psi_{\zeta}
(D_{i},Z_{i},\xi_0)+o_{P}(1),
\end{equation*}
where the IF $\psi_\zeta$ has zero mean and finite variance, and $\xi_0$ contains $\zeta_{0}$ and may contain additional nuisance parameters. We assume that the pathwise derivative of $\mathbb{E}[\psi_\zeta
(D_{i},Z_{i},\xi)]$ with respect to these additional nuisance parameters at $\xi_{0}$ is zero. Then, by standard arguments in regular estimation, see \cite{newey1994asymptotic}, and the previous remark, all our results apply in the parametric first step case with the adjustment term $\phi
_{1}=-\alpha_{01}'\psi_\zeta$, where
\begin{equation*}
\alpha_{01}=\mathbb{E}[D_{01}(W,\dot{G})],
\end{equation*}
and $\dot{G}(Z,\zeta_{0})=\partial G/\partial \zeta (Z,\zeta_{0})$. In particular, if $\psi_\zeta=0$, i.e., if $\sqrt{n}(\hat{\zeta}-\zeta_{0})=o_{P}(1),$ and other mild conditions are satisfied (to apply a Delta method), then there is no estimation effect from the generated regressors. This is typically the case when $\hat{\zeta}$ is constructed from a large sample (with a large sample size relative to $n$).
\end{rem}
\begin{rem}[Simplifications]
\label{DRcase} If the original identifying moments are such that $\alpha
_{02}= 0,$ then the first step linearization simplifies to:
\begin{equation*}
D_{01}(w,g) = D_{dir}(w,g).
\end{equation*}
Another simplification occurs under the index restriction $\mathbb{E}[S|D,Z]=\mathbb{E}
[S|X,V]$, which implies:
\begin{equation*}
D_{01}(w,g) = D_{dir}(w,g)-\frac{\partial h_{0}}{\partial v
}(x,v)\alpha _{02}(x,v)D_{\varphi }g.
\end{equation*}
In both cases, the corresponding $\alpha _{01}\in \Delta _{1}$, given by (\ref{eq:alpha1_defi}) in Appendix \ref{sec:Proofs}, simplifies accordingly.
\end{rem}
\subsection{Cross-fitted automatic estimators}
\label{sec:est_auto}
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{sec:ident_alphas}. We can combine the orthogonality conditions with the linearization results in Section~\ref{sec:Linearization} to obtain sets of moment conditions for the estimation of the Riesz representers. In particular, a combination of equation~\eqref{eq:auto_momentF} and Theorem~\ref{thm:first_step} gives
\begin{equation}
\mathbb{E}[D_{01}(W,\delta _{1})]=\mathbb{E}[-\alpha
_{01}(Z)r_{e}(Z)\delta _{1}(Z)],\text{ for each }\delta _{1}\in \Delta _{1},
\label{eq:Lin_projection1}
\end{equation}
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{rho}). 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~\eqref{eq:auto_momentS} and Proposition~\ref{prop:second_step}:
\begin{equation}
\mathbb{E}[D_{02}(W,\delta _{2})]=\mathbb{E}[\alpha
_{02}(X,V)\delta _{2}(D,Z)]=0,\text{ for each }\delta _{2}\in \Delta
_{2}(g_{0}), \label{eq:Lin_projection2}
\end{equation}
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{sec:Examples} for some examples).
Equations \eqref{eq:Lin_projection1} and \eqref{eq:Lin_projection2} 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$.
\subsubsection{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~\eqref{eq:Lin_projection2} for $
\delta _{2}=b_{j}$, $j=1,...,J$. This gives the following $J$ moment
conditions:
\begin{equation*}
\mathbb{E}[D_{02}(W,\mathbf{b}_{J})]=\mathbb{E}[\mathbf{b}_{J}(X,V)\mathbf{b}
_{J}(X,V)^{\prime }]\boldsymbol{\rho }_{0J},
\end{equation*}
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
\begin{equation}
\min_{\boldsymbol{\rho }_{J}\in \mathbb{R}^{J}}\left\{ -2\mathbb{E}[D_{02}(W,
\mathbf{b}_{J})^{\prime }]\boldsymbol{\rho }_{J}+\boldsymbol{\rho }
_{J}^{\prime }\mathbb{E}[\mathbf{b}_{J}(X,V)\mathbf{b}_{J}(X,V)^{\prime }]
\boldsymbol{\rho }_{J}+\lambda \lVert \boldsymbol{\rho }_{J}\rVert
_{q}^{q}\right\}, \label{eq:Auto2_objective}
\end{equation}
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~\eqref{eq:Auto2_objective}. We estimate $\mathbb{E}[D_{02}(W,\mathbf{b}_{J})]$ by
\begin{equation*}
\hat{D}_{2\ell }\equiv \frac{1}{n-n_{\ell }}\sum_{\ell ^{\prime }\neq \ell
}\sum_{i\in I_{\ell ^{\prime }}}D_{02}(W_{i},\mathbf{b}_{J}|\hat{g}_{\ell
\ell ^{\prime }},\hat{h}_{\ell \ell ^{\prime }},\tilde{\theta}_{\ell \ell
^{\prime }}),
\end{equation*}
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
\begin{equation*}
\hat{B}_{\ell }\equiv \frac{1}{n-n_{\ell }}\sum_{\ell ^{\prime }\neq \ell
}\sum_{i\in I_{\ell ^{\prime }}}\mathbf{b}_{J}(X_{i},\varphi (D_{i},Z_{i},
\hat{g}_{\ell \ell ^{\prime }}))\mathbf{b}_{J}(X_{i},\varphi (D_{i},Z_{i},
\hat{g}_{\ell \ell ^{\prime }}))^{\prime }.
\end{equation*}
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
\begin{equation}\label{eq:min_rho}
\hat{\alpha}_{2\ell}=\mathbf{b}_J'\widehat{\boldsymbol{\rho}}_{J\ell},\quad
\widehat{\boldsymbol{\rho}}_{J\ell}=\operatornamewithlimits{argmin}_{\boldsymbol{\rho}_J\in\mathbb{R}^J}
\left\{-2\hat{D}_{2\ell}'\boldsymbol{\rho}_J+\boldsymbol{\rho}_J'\hat{B}_\ell\boldsymbol{\rho}_J
+\lambda\|\boldsymbol{\rho}_J\|_q^q\right\}.
\end{equation}
The tuning parameter $\lambda $ can be chosen by cross-validation.
\subsubsection{Automatic estimation of the first-step Riesz representer}
\label{sec:auto_1step}
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~
\eqref{eq:Lin_projection1} for $\delta _{1}=c_{k}$, $k=1,...,K$. This gives
the following $K$ moment conditions:
\begin{equation*}
\mathbb{E}[D_{01}(W,\mathbf{c}_{K})]=-\mathbb{E}[r_e(Z)
\mathbf{c}_{K}(Z)\mathbf{c}_{K}(Z)^{\prime }]\boldsymbol{\beta }_{0K},
\end{equation*}
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{rho}).
We use these conditions as a basis to construct
the objective function to estimate $\boldsymbol{\beta }_{0K}$:
\begin{equation*}
\min_{\boldsymbol{\beta }_{K}\in \mathbb{R}^{K}}\left\{ -2\mathbb{E}[D_{01}(W,
\mathbf{c}_{K})^{\prime }]\boldsymbol{\beta }_{K} + \boldsymbol{\beta }
_{K}^{\prime }\mathbb{E}[-r_e(Z)\mathbf{c}_{K}(Z)\mathbf{c}
_{K}(Z)^{\prime }]\boldsymbol{\beta }_{K}+\lambda \lVert \boldsymbol{\beta }
_{K}\rVert _{q}^{q}\right\} ,
\end{equation*}
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
\begin{equation}\label{eq:min_beta}
\hat{\alpha}_{1\ell}=\mathbf{c}_K'\widehat{\boldsymbol{\beta}}_{K\ell},\quad
\widehat{\boldsymbol{\beta}}_{K\ell}=\operatornamewithlimits{argmin}_{\boldsymbol{\beta}_K\in\mathbb{R}^K}
\left\{-2\hat{D}_{1\ell}'\boldsymbol{\beta}_K+\boldsymbol{\beta}_K'\hat{C}_\ell\boldsymbol{\beta}_K
+\lambda\|\boldsymbol{\beta}_K\|_q^q\right\}.
\end{equation}
with
\begin{equation*}
\hat{C}_{\ell }\equiv \frac{-1}{n-n_{\ell }}\sum_{\ell ^{\prime }\neq \ell
}\sum_{i\in I_{\ell ^{\prime }}}r_e(Z_{i}|\hat{g}_{\ell \ell ^{\prime
}})\mathbf{c}_{K}(Z_{i})\mathbf{c}_{K}(Z_{i})^{\prime },
\end{equation*}
and
\begin{equation*}
\hat{D}_{1\ell }\equiv \frac{1}{n-n_{\ell }}\sum_{\ell ^{\prime }\neq \ell
}\sum_{i\in I_{\ell ^{\prime }}}D_{01}(W_{i},\mathbf{c}_{K}|\hat{g}_{\ell
\ell ^{\prime }},\hat{h}_{\ell \ell ^{\prime }},\hat{\alpha}_{2\ell \ell
^{\prime }},\tilde{\theta}_{\ell \ell ^{\prime }}).
\end{equation*}
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~\eqref{eq:min_rho}, with $\hat{D}_{2\ell}$ and $\hat{
B}_{\ell}$ replaced by
\begin{align*}
\hat{D}_{2 \ell\ell^{\prime }} &\equiv \frac{1}{n-n_{\ell}-n_{\ell^{\prime }}
} \sum_{\ell^{\prime\prime}\notin \{\ell,\ell^{\prime }\}}\sum_{i\in
I_{\ell^{\prime\prime}}} D_{02}(W_i, \mathbf{b}_J | \hat{g}_{\ell\ell^{\prime
}\ell^{\prime \prime }}, \hat{h}_{\ell\ell^{\prime }\ell^{\prime \prime }},
\tilde{\theta}_{\ell\ell^{\prime }\ell^{\prime \prime }}) \text{ and } \\
\hat{B}_{\ell\ell^{\prime }} &\equiv \frac{1}{n-n_{\ell}-n_{\ell^{\prime }}}
\sum_{\ell^{\prime\prime}\notin \{\ell,\ell^{\prime }\}}\sum_{i\in
I_{\ell^{\prime\prime}}} \mathbf{b}_J(X_i,\varphi(D_i,Z_i,\hat{g}
_{\ell\ell^{\prime }\ell^{\prime \prime }}))\mathbf{b}_J(X_i,\varphi(D_i,Z_i,
\hat{g}_{\ell\ell^{\prime }\ell^{\prime \prime }}))^{\prime },
\end{align*}
respectively.
Furthermore, $D_{01}$ may depend on the derivatives $\partial h_{0}/\partial v$ and $
\partial \alpha _{02}/\partial v$ (see equation~\eqref{eq:D1_defi}). 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 \cite{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
\begin{equation}
\frac{\partial \hat{h}_{\ell\ell^{\prime }}}{\partial v}(x,v)\equiv \frac{
\hat{h}_{\ell\ell^{\prime }}(x,v+t_n)-\hat{h}_{\ell\ell^{\prime }}(x,v)}{t_n}
. \label{eq:deriv_h_est}
\end{equation}
This approach has been used and justified theoretically in \cite{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 \cite{dai2016optimal}.
\section{Examples: Three-step debiased estimators}
\label{sec:Examples}
\subsection{Hd-PS regression adjustment}
\label{sec:ex_ATE_PLM}
The debiased three-step estimator is
\begin{equation*}
\hat{\theta} = \frac{\sum_{\ell =1}^{L}\sum_{i\in I_{\ell }} \left[
\left( Y_i - \hat{h}_\ell(\hat{V}_{i\ell})\right)\left( D_{i}-\hat{V}_{i\ell}\right) + \hat{\alpha}_{1\ell }(Z_{i})\cdot \left(D_{i}-\hat{V}_{i\ell}\right) \right]}{\sum_{\ell
=1}^{L}\sum_{i\in I_{\ell }}\left( D_{i}-\hat{V}_{i\ell}\right) ^{2}},
\end{equation*}
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{ex:ATE_PLM1} (p.~\pageref{ex:ATE_PLM1}). In this section, we detail the estimation of the Riesz representer $\alpha_{01}$. As discussed in Section~\ref{sec:est_auto}, 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 \eqref{eq:min_beta}. 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
\begin{equation*}
\frac{d}{d\tau} \mathbb{E}[\epsilon(W, g_\tau)] = \mathbb{E}\left[-\Lambda(g_0(Z))[1-\Lambda(g_0(Z))] \frac{d}{d\tau} g_\tau\right] \text{ and } r_e(Z)= -\Lambda(g_0(Z))[1-\Lambda(g_0(Z))].
\end{equation*}
Therefore,
\begin{equation*}
\hat{C}_\ell = \frac{1}{n-n_{\ell }}\sum_{\ell ^{\prime }\neq \ell
}\sum_{i\in I_{\ell ^{\prime }}}\hat{V}_{i\ell\ell'}(1-\hat{V}_{i\ell\ell'})\mathbf{c}_{K}(Z_{i})\mathbf{c}_{K}(Z_{i})^{\prime }, \text{ where } \hat{V}_{i\ell\ell'} \equiv \Lambda(\hat{g}_{\ell\ell'}(Z_i)).
\end{equation*}
To find $\hat{D}_{1\ell}$, note that, since $\alpha_{02}=0$, we have that $D_{01} = D_{dir}$ and
\begin{align*}
\frac{\partial}{\partial\tau}\mathbb{E}[m(W, g_\tau, h_0, \theta_0)]=\frac{\partial}{\partial\tau}\mathbb{E}\left[ \{Y - h_0(\Lambda(g_\tau(Z))) -\theta_0 (D - \Lambda(g_\tau(Z)))\} \cdot (D - \Lambda(g_\tau(Z)))\right] \\
= \mathbb{E} \left[ \left\{ -\left(Y - h_0(V) -\theta_0(D-V)\right) + \left( -\dot{h}_0(V) + \theta_0 \right) \cdot (D - V) \right\} \dot{\Lambda}(g_0(Z)) \frac{\partial g_\tau}{\partial\tau} \right]
\\ = \mathbb{E} \left[ -\left(Y - h_0(V) \right) \dot{\Lambda}(g_0(Z)) \frac{\partial g_\tau}{\partial\tau} \right],
\end{align*}
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:
\begin{equation*}
D_{01}(W_{i},\mathbf{c}_{K}|\hat{g}_{\ell \ell ^{\prime }},\hat{h}_{\ell \ell
^{\prime }})=-\left[Y_{i}-\hat{h}_{\ell \ell
^{\prime }}(\hat{V}_{i\ell \ell ^{\prime }}) \right] \hat{V}_{i\ell \ell
^{\prime }} \left( 1- \hat{V}_{i\ell \ell
^{\prime }} \right) \mathbf{c}_{K}(Z_i).
\end{equation*}
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~\eqref{eq:alpha1ell_riesz_lasso}.
\subsection{Partially linear model with generated regressors: Autoencoders}
\label{sec:ex_PLMauto}
The partially linear model is a workhorse for debiased machine learning methods, see \cite{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 \eqref{orth1}. 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
\begin{equation*}
m(W,g_{0},h_{0},\theta _{0})=( Y-h_{0Y}(V)-\theta
_{0}' ( D-h_{0D}(V))) \cdot (D - h_{0D}(V)),
\end{equation*}
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{sec:ex_ATE_PLM}), 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
\begin{equation*}
\hat{\theta}=\operatornamewithlimits{argmin}_{\theta\in\Theta}\hat{\psi}
(\theta)^{\prime}\hat{\Upsilon}\hat{\psi}(\theta),
\end{equation*}
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
\begin{equation}\label{eq:PLM_moments}
\hat{\psi}_{i\ell}(\theta)=
\left(Y_i-\hat h_{\ell, Y}(\hat V_{i\ell})-\theta'(D_i-\hat h_{\ell, D}(\hat V_{i\ell}))\right)
(D_i-\hat h_{\ell, D}(\hat V_{i\ell}))
+\widehat{\boldsymbol{\alpha}}_{1\ell}(Z_i)\,\epsilon(W_i,\hat g_\ell),
\end{equation}
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{sec:auto_1step}:
$\hat\alpha_{1j\ell}(z)=\mathbf c_K(z)'\widehat{\boldsymbol\beta}_{Kj\ell}$, with
$\widehat{\boldsymbol\beta}_{Kj\ell}$ solving \eqref{eq:min_beta}. 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{sec:ex_ATE_PLM}. Since $\alpha_{02}=0$, the indirect effect is zero and only the direct effect remains. The linearization of the $j$-th moment condition is
\begin{equation} \label{eq:PLM_linearization}
\begin{aligned}
D_{01j}(W_i, g) = r_{dir}(W_i)D_\varphi g, \text{ with } r_{dir}(W_i) = -\dot h_{0D_j}(V_i) \cdot\varepsilon+\dot\varepsilon \cdot [D_{ji}-h_{0D_j}(V_i)], \\
\varepsilon=Y_i-h_{0Y}(V_i)-\theta_0'(D_i-h_{0D}(V_i)) \text{ and } \dot\varepsilon \equiv \frac{\partial\varepsilon}{\partial v} = \dot{h}_{0Y}(V_i) + \theta_0'\dot{h}_{0D}(V_i),
\end{aligned}
\end{equation}
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~\eqref{eq:PLM_linearization}: (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 \cite{hahn2013asymptotic}.
\medskip
\noindent\textbf{Learned 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
\begin{equation*}
(e_{0},d_{0})=\operatornamewithlimits{argmin}_{e\in \mathcal{E},d\in \mathcal{D}}\mathbb{E}[\mathcal{L}(Z,d(e(Z)))],
\end{equation*}
where $\mathcal{L}$ is a loss and $\mathcal E$ and $\mathcal D$ are function classes (see Figure~\ref{fig:autoencoders}). 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 \citep{bengio2013representation}.
\begin{figure}[htbp]
\centering
\includegraphics[width=0.65\linewidth]{Plots/figure_autoencoders.pdf}
\caption{Autoencoder architecture.}
\label{fig:autoencoders}
\end{figure}
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
\begin{equation*}
\mathbb{E}[\delta_1(Z)\epsilon(W,g_0)]=0\quad\text{for all }\delta_1\in\Delta_1,
\end{equation*}
where $\Delta_1$ is the linear span generated by the columns of the Jacobian
\begin{equation*}
\dot f_0(Z)=\left.\frac{\partial f_\zeta}{\partial\zeta}(Z)\right|_{\zeta=\zeta_0},
\text{ with } f_\zeta(Z)=d_\zeta(e_\zeta(Z)).
\end{equation*}
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
\begin{equation}\label{AutoObj}
\min_{\boldsymbol{\beta }_{K}\in \mathbb{R}^{K}}\left\{ -2\mathbb{E}[D_{01j}(w,\mathbf{e}_{K})']\boldsymbol{\beta }_{K}+\boldsymbol{\beta }
_{K}^{\prime }\mathbb{E}[\mathbf{c}_{K}(Z)\mathbf{c}_{K}^{\prime }(Z)]
\boldsymbol{\beta }_{K}+\lambda \lVert \boldsymbol{\beta }_{K}\rVert
_{q}^{q}\right\},
\end{equation}
with $D_{01j}$ given by equation~\eqref{eq:PLM_linearization} 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:
\begin{equation*}
\hat{\mathbf c}_{K,\ell\ell'}=\left.\frac{\partial f_\zeta}{\partial\zeta}\right|_{\zeta=\hat\zeta_{\ell\ell'}}'
\text{ and }
\hat{\mathbf e}_{K,\ell\ell'}=\left.\frac{\partial e_\zeta}{\partial\zeta}\right|_{\zeta=\hat\zeta_{\ell\ell'}}.
\end{equation*}
Then,
\begin{align*}
\hat{C}_{\ell }&=\frac{1}{n-n_{\ell }}\sum_{\ell ^{\prime }\neq \ell
}\sum_{i\in I_{\ell ^{\prime }}}\hat{\mathbf c}_{K,\ell\ell'}(Z_i)\hat{\mathbf c}_{K,\ell\ell'}(Z_i)^{\prime} \text{ and} \\
\hat{D}_{1j\ell }&=\frac{1}{n-n_{\ell }}\sum_{\ell ^{\prime }\neq \ell
}\sum_{i\in I_{\ell ^{\prime }}}D_{01j}(W_{i},\hat{\mathbf e}_{K,\ell\ell'}\mid\hat g_{\ell\ell'},\hat h_{\ell\ell'},\tilde\theta_{\ell\ell'}),
\end{align*}
with $\hat V_{i\ell\ell'}=e_{\hat\zeta_{\ell\ell'}}(Z_i)$. Solving the sample analog of \eqref{AutoObj} 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 \eqref{eq:PLM_moments}, 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$.
\subsection{CASF in a non-separable model}
\label{sec:CASF_details}
The three-step debiased estimator of the CASF is given in equation~\eqref{eq:CASF_D}. 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
\begin{equation*}
m(w,g,h,\theta )=\int h(x^{\ast },\varphi (d,z,g))dF^{\ast }(x^{\ast
})-\theta .
\end{equation*}
This moment is already linear in $h$ and hence
\begin{equation*}
D_{02}(w,b_{j})=\int b_{j}(x^{\ast }, \varphi(d,z,g_{0}))dF^{\ast }(x^{\ast }),
\end{equation*}
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
\begin{equation*}
D_{02}(W_{i},b_{j}|\hat{g}_{\ell \ell ^{\prime }})=\frac{1}{S}
\sum_{s=1}^{S}b_{j}(X_{s}^{\ast },D_{i}-\hat{g}_{\ell \ell ^{\prime
}}(Z_{i})),
\end{equation*}
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
{eq:min_rho}).
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
\begin{equation*}
r_{dir}(w)= -\int \frac{\partial h_{0}}{\partial v}(x^{\ast
},d-g_{0}(z))dF^{\ast }(x^{\ast }).
\end{equation*}
We can now plug in the expression for $D_{dir}$ into equation~\eqref{eq:D1_defi}, where the linearization of the first step effect is defined. Recall that $D_{\varphi }g=-g$. Then, for the CASF, equation~\eqref{eq:D1_defi} becomes
\begin{equation*}
D_{01}(w,g)\equiv \left\{ r_{dir}(w)+\frac{\partial h_{0}}{\partial v}
(x,v)\alpha _{02}(x,v)\right\} g(z).
\end{equation*}
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
\begin{equation*}
-\frac{1}{S}\sum_{s=1}^{S}\frac{\partial \hat{h}_{\ell \ell ^{\prime }}}{
\partial v}(X_{s}^{\ast },D_{i}-\hat{g}_{\ell \ell ^{\prime }}(Z_{i})),
\end{equation*}
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
\begin{equation*}
\mathbf{b}_{J}(X_{i},\hat{V}_{i\ell \ell ^{\prime }})^{\prime }\widehat{
\boldsymbol{\rho }}_{J\ell \ell ^{\prime }}\cdot \frac{\partial \hat{h}
_{\ell \ell ^{\prime }}}{\partial v}(X_{i},\hat{V}_{i\ell \ell ^{\prime }}).
\end{equation*}
Therefore, we have that, for $i\in I_{\ell ^{\prime }}$,
\begin{align*}
D_{01}(W_{i},c_{k}|\hat{g}_{\ell \ell ^{\prime }},\hat{h}_{\ell \ell ^{\prime
}},\hat{\alpha}_{2\ell \ell ^{\prime }})=c_{k}(Z_{i})\cdot & \left\{ -\frac{1}{S}\sum_{s=1}^{S}\frac{\partial
\hat{h}_{\ell \ell ^{\prime }}}{\partial v}(X_{s}^{\ast },\hat{V}_{i\ell \ell
^{\prime }}) +\mathbf{b}_{J}(X_{i},\hat{V}_{i\ell \ell ^{\prime }})^{\prime }
\widehat{\boldsymbol{\rho }}_{J\ell \ell ^{\prime }}\cdot \frac{\partial \hat{h}
_{\ell \ell ^{\prime }}}{\partial v}(X_{i},\hat{V}_{i\ell \ell ^{\prime
}})\right\} ,
\end{align*}
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
{eq:min_beta}) and then $\hat{\alpha}_{1\ell }=\mathbf{c}_{K}^{\prime }\widehat{\boldsymbol{\beta }}_{K\ell }$.
\begin{comment}
\begin{remark}[Profiling]
Estimation of the correction term with a profiled-out $h_0$ is based on the initial estimators $\tilde{\theta}_\ell$. For each $\ell$, the $\hat{h}_\ell(\cdot, \theta)$ is estimated for $\theta=\tilde{\theta}_\ell$. Additionally, we note that the estimators $\hat{h}_{\ell\ell'}(\cdot,\theta)$ and $\hat{h}_{\ell\ell'\ell''}(\cdot,\theta)$ (required for automatic estimation), that do not use observations in $I_\ell\cup I_{\ell'}$ or not in $I_\ell\cup I_{\ell'} \cup I_{\ell''}$, respectively, are estimated for initial estimators $\tilde{\theta}_{\ell\ell'}$ and $\tilde{\theta}_{\ell\ell'\ell''}$ satisfying those same properties.
To sum up, in the presence of profiling, the debiased moment functions in equation~\eqref{eq:dmon_i} are estimated by:
\begin{equation*}
\hat{\psi}_{i\ell}(\theta) \equiv m(W_{i},\hat{g}_{\ell},\hat{h}_\ell(\cdot, \theta), \theta)+\phi(W_{i},\hat{g}_{\ell}, \hat{h}_\ell, \hat{\alpha}_{\ell}(\cdot,\tilde{\theta}_\ell),\tilde{\theta}_{\ell}).
\end{equation*}
Also, the Jacobian of $m$ w.r.t. $\theta$ must be extended:
\begin{equation*}
M \equiv \mathbb{E}\left[\frac{\partial m}{\partial\theta}(W,g_0,h_0(\cdot, \theta),\theta) \Big\rvert_{\theta=\theta_0}\right],
\end{equation*}
which may be estimated by
\begin{equation*}
\hat{M} \equiv \frac{1}{n}\sum_{\ell=1}^{L}\sum_{i\in I_{\ell}} \frac{\partial m}{\partial\theta}(W_i, \hat{g}_\ell,\hat{h}_\ell(\cdot, \tilde{\theta}_\ell), \tilde{\theta}_\ell).
\end{equation*}
\end{remark}
\end{comment}
\begin{comment}
\begin{remark}{\label{rem:initial_est}}
Following \cite{chernozhukov2022locally}, we can base the the initial estimator $\tilde{\theta}_{\ell}$ on the original moment conditions (without the bias correction). Note, however, that for each group $I_\ell$, the initial estimator must only use observations that are not in $I_\ell$. We therefore conduct standard GMM, but only with the observations not in $I_\ell$.
Furthermore, for each group $I_\ell$, the first and second step estimators must not use observations in $I_\ell$. Therefore, we construct estimators $\hat{g}_{\ell\ell'}$ and $\hat{h}_{\ell\ell'}$ that use only observations not in $I_\ell \cup I_{\ell'}$. Then, we can define the moment condition
\begin{equation*}
\hat{m}_{\ell}(\theta)=\frac{1}{n-n_{\ell}}\sum_{\ell^{\prime}\neq\ell}\sum_{i\in
I_{\ell^{\prime}}}m(W_{i},\hat{g}_{\ell\ell'}, \hat{h}_{\ell\ell'},\theta),
\end{equation*}
where $n_\ell$ is the number of observations in $I_\ell$. Consider a weighting matrix $\hat{\Upsilon}_{\ell}$ that does not use observations in $I_\ell$. The initial estimator for group $I_\ell$ is defined by the GMM condition:
\[
\tilde{\theta}_{\ell}=\operatornamewithlimits{argmin}_{\theta\in\Theta}\hat{m}_{\ell}(\theta
)^{\prime}\hat{\Upsilon}_{\ell}\hat{m}_{\ell}(\theta).
\]
\end{remark}
\end{comment}
\begin{comment}
\subsection{Semiparametric CASF}
A practical semiparametric specification for $h_{0}$ when $X$ is $p$
-dimensional and $p$ is high is a PLM, where
\begin{equation*}
h_{0}(x,v)=x^{\prime }\beta _{0}+\kappa _{0}(v),
\end{equation*}
with $\beta _{0}$ and $\kappa _{0}$ unknown finite and infinite-dimensional
parameters, respectively. In this specification, $X$ contains an intercept,
and hence, we can assume that $\kappa _{0}$ belongs to the subspace of zero
mean functions in $L_{2}(V)$, denoted as $L_{2}^{0}(V)$. This setting
corresponds to the semiparametric orthogonality conditions \eqref{orth2}
where
\begin{equation*}
\Delta _{2}(g_{0})=\{\delta (x,v)=x^{\prime }\beta +\kappa (v)\colon \beta
\in \mathbb{R}^{p},\kappa \in L_{2}^{0}(V)\}\subseteq L_{2}(X,V).
\end{equation*}
This specification generalizes the classical linear structural control
function approach to a semiparametric setting. In this PLM, the CASF is
given by $\theta _{0}=\beta _{0}^{\prime }\mathbb{E}^{\ast }[X]$, where $
\mathbb{E}^{\ast }[X]$ denotes the mean of $X$ under the counterfactual
distribution $F^{\ast }$.
To estimate $h_{0}$ in this semiparametric specification use a dictionary $
\mathbf{b}_{J}=(b_{1},...,b_{J})^{\prime }$ such that the first $p$
components give $X$, i.e. $X=(b_{1},...,b_{p})^{\prime }$, and $\mathbf{b}
_{p+1,J}(v)=(b_{p+1}(v),...,b_{J}(v))^{\prime }$ are functions of $v$. Then,
let $\hat{h}_{\ell }(X_{i},\hat{V}_{i\ell })=X_{i}^{\prime }\hat{\boldsymbol{
\beta }}_{J\ell }+\mathbf{b}_{p+1,J}^{\prime }(\hat{V}_{i\ell })\hat{
\boldsymbol{k}}_{J\ell }$ denote the fitted value of a lasso regression of $
Y_{i}$ on $\mathbf{b}_{J}(X_{i},\hat{V}_{i\ell })\ $based on all
observations not in $I_{\ell }$. The debiased semiparametric CASF estimator
is
\begin{equation*}
\hat{\theta}=\hat{\boldsymbol{\beta }}_{J\ell }^{\prime }\bar{X}^{\ast }+
\frac{1}{n}\sum_{\ell =1}^{L}\sum_{i\in I_{\ell }}\phi (W_{i},\hat{g}_{\ell
},\hat{h}_{\ell },\hat{\alpha}_{\ell }),
\end{equation*}
where
\begin{equation*}
\bar{X}^{\ast }=\frac{1}{S}\sum_{s=1}^{S}X_{s}^{\ast },
\end{equation*}
and $\phi $ is given in (\ref{adju}) with some simplified expressions for $
\hat{\alpha}_{1\ell }(Z_{i})$ and $\hat{\alpha}_{2\ell }(X_{i},\hat{V}
_{i\ell }).$
For semiparametric specifications of the second step, such as the PLM, the
nuisance parameter $\alpha _{02}$ is the orthogonal projection of $r_{2}$
onto $\Delta _{2}(g_{0})$. These objects are generally hard to estimate and
may cause the plug-in estimator for $\alpha _{02}$ to behave poorly. We
advocate automatic estimation as a potential
solution to this issue.
We illustrate how some simplifications for estimating the linearizations may
occur in semiparametric settings with the PLM and the CASF. The zero mean
restriction of the nonparametric component $\kappa (v)$ in the partial
linear specification implies that, for $\delta _{2}(x,v)=x^{\prime }\beta
+\kappa (v)\in \Delta_{2}(g_{0}),$
\begin{equation*}
\mathbb{E}[D_{02}(W,\delta _{2})]=\beta ^{\prime }\mathbb{E}^{\ast }[X].
\end{equation*}
The linear approximations necessary for the automatic estimation of $\alpha
_{02}$ are known in this example, with
\begin{align}
\mathbb{E}[D_{02}(W,X)]& =\mathbb{E}^{\ast }[X] \\
\mathbb{E}[D_{02}(W,b_{j})]& =0,\text{ for }j=p+1,\ldots ,J.
\end{align}
Therefore, $D_{02}(W_{i},b_{j}|\hat{g}_{\ell \ell ^{\prime }},\hat{h}_{\ell
\ell ^{\prime }},\tilde{\theta}_{\ell \ell ^{\prime }})=\bar{X}^{\ast }$ for
$j=1,...,p,$ and $D_{02}(W_{i},b_{j}|\hat{g}_{\ell \ell ^{\prime }},\hat{h}
_{\ell \ell ^{\prime }},\tilde{\theta}_{\ell \ell ^{\prime }})=0$ for $
j=p+1,\ldots ,J.$ With this we construct $\hat{\alpha}_{2\ell }=\mathbf{b}
_{J}^{\prime }\hat{\boldsymbol{\rho }}_{J\ell }$ following (\ref{eq:min_rho}
).
Likewise, the expression for the linearization w.r.t. the first step
simplifies to
\begin{equation*}
\mathbb{E}[D_{01}(W,c_{k})]=\mathbb{E}[(\alpha _{02}(X,V)-1)\dot{\kappa}
_{0}(V)c_{k}(Z)],
\end{equation*}
where $\dot{\kappa}_{0}(v)=\partial h_{0}/\partial v$ can be approximated by
$(\partial \mathbf{b}_{J}/\partial v)^{\prime }\boldsymbol{\eta }_{0J}$ when
$h_{0}$ is approximated by $\mathbf{b}_{J}^{\prime }\boldsymbol{\eta }
_{0J}\equiv X_{i}^{\prime }\mathbf{\beta }_{J\ell }+\mathbf{b}
_{p+1,J}^{\prime }(\hat{V}_{i\ell })\kappa _{J\ell }.$ Therefore, to
estimate $D_{1\ell }$ according to equation~\eqref{eq:D1_estimation}, we
have that, for $i\in I_{\ell ^{\prime }}$,
\begin{equation}
D_{01}(W_{i},c_{k}|\hat{g}_{\ell \ell ^{\prime }},\hat{h}_{\ell \ell ^{\prime
}},\hat{\alpha}_{2\ell \ell ^{\prime }},\tilde{\theta}_{\ell \ell ^{\prime
}})=\left( \mathbf{b}_{J}(X_{i},\hat{V}_{i\ell \ell ^{\prime }})^{\prime }
\hat{\boldsymbol{\rho }}_{J\ell \ell ^{\prime }}-1\right) \cdot \frac{
\partial \hat{h}_{\ell \ell ^{\prime }}}{\partial v}(X_{i},\hat{V}_{i\ell
\ell ^{\prime }})\cdot c_{k}(Z_{i}),
\end{equation}
for each $k=1,...,K$. This can then be used to construct the objective
function to estimate $\hat{\boldsymbol{\beta }}_{K\ell }$.
\end{comment}
\section{Asymptotic theory}
\label{sec:asymptotic}
\subsection{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 \cite{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
\citep[see][]{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~\eqref{eq:dmon_i}). We modify the
expansion of $\hat{\psi}_{i\ell}(\theta_0)-\psi(W_i,
g_0,h_0,\alpha_{0},\theta_0)$ given by \cite{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{ass:gen_reg}).
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}|$.}
\begin{ass} \label{ass:span} \mbox{} \\[-20pt]
\begin{enumerate}[label=\textbf{\alph*.},ref=\ref{ass:span}.\alph*]
\item \label{ass:span:b} For every $j$, $b_j\in \Delta_{2}(g_0)$. Also, $\forall \delta_{2}\in \Delta_{2}(g_0)$ and for every $\varepsilon>0$, there exist $J$ and $\boldsymbol{\rho}_J$ such that $\lVert \delta_2 - \mathbf{b}_J'\boldsymbol{\rho}_J \rVert_2 < \varepsilon$.
\item \label{ass:span:c} For every $k$, $c_k\in \Delta_{1}$. Also, $\forall \delta_{1}\in \Delta_{1}$ and for every $\varepsilon>0$, there exist $K$ and $\boldsymbol{\beta}_K$ such that $\lVert \delta_1 - \mathbf{c}_K'\boldsymbol{\beta}_K \rVert_2 < \varepsilon$.
\end{enumerate}
\end{ass}
We also assume bounded dictionaries
\citep[see, for
instance,][]{newey1997convergence}:
\begin{ass} \label{ass:bounded_dicts}
$\sup_{j\in\mathbb{N}} |b_j(X,V)|<\infty$ and $\sup_{k\in\mathbb{N}} |c_k(Z)|<\infty$.
\end{ass}
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:
\begin{ass}
\label{ass:abs_sumability} For the real-valued sequences $(\rho_{0j})_{j=1}^\infty$ and $(\beta_{0k})_{k=1}^\infty$ such that $\alpha_{02}(x,v)=\sum_{j=1}^{\infty} \rho_{0j} b_j(x,v)$ and $\alpha_{01}(z)=\sum_{k=1}^{\infty} \beta_{0k} c_k(z)$:
\begin{enumerate}[label=\textbf{\alph*.},ref=\ref{ass:abs_sumability}.\alph*]
\item \label{ass:abs_sumability:main} $\sum_{j=1}^{\infty} |\rho_{0j}|<\infty$ and $\sum_{k=1}^{\infty} |\beta_{0k}|<\infty$.
\item \label{ass:abs_sumability:largest_terms} For a $C>0$, the atoms $b_j$ and $c_k$ corresponding to the largest $C\sqrt{n}$ values of $\rho_{0j}$ and $\beta_{0k}$ are included in $\mathbf{b}_J$ and $\mathbf{c}_K$.
\end{enumerate}
\end{ass}
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 \citep[see][]{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}$
\citep[see][p.~985]{chernozhukov2022automatic}.
We require the following estimation rates:
\begin{ass} \label{ass:est_convergence} There is $1/3 < r < 1/2$ such that
\begin{enumerate}[label=\textbf{\alph*.},ref=\ref{ass:est_convergence}.\alph*]
\item \label{ass:est_convergence:estimators} $\lVert \hat{g}_\ell-g_0 \rVert_2 = O_p(n^{-r})$ and $\lVert \hat{h}_\ell-h_0 \rVert_2 = O_p(n^{-r})$.
\item \label{ass:est_convergence:linearizations} $\lVert\hat{D}_{1\ell}-\mathbb{E}[D_{01}(W, \mathbf{c}_K)]\rVert_\infty=O_p(n^{-r})$ and $\lVert\hat{D}_{2\ell}-\mathbb{E}[D_{02}(W, \mathbf{b}_J)]\rVert_\infty=O_p(n^{-r})$.
\end{enumerate}
\end{ass}
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 \cite{mammen2012nonparametric}; for Lasso rates, see \cite{bickel2009simultaneous,bunea2007sparsity,zhang2008sparsity}, and references
therein; for $L_2$-rates with boosting with high-dimensional regressors see \cite{kueck2023estimation}; for deep neural networks with a ReLU activation function, see \cite{farrell2021deep}. Under some regularity conditions on the linearizations
\citep[see][Ass.~12]{chernozhukov2022automatic}, Assumption~\ref{ass:est_convergence:linearizations} 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:
\begin{ass}
\label{ass:rates} \mbox{} \\[-20pt]
\begin{enumerate}[label=\textbf{\alph*.},ref=\ref{ass:rates}.\alph*]
\item \label{ass:rates:penalty} The Lasso penalty term $\lambda=\lambda(n)$ for estimation of $(\alpha_{01},\alpha_{02})$ satisfies: $n^{-r}=o(\lambda)$ and $\lambda=o(n^{c-r})$ for every $c>0$.
\item \label{ass:rates:numterms} The number of terms in the dictionaries satisfy $J, K = O(n^\kappa)$ for a constant $\kappa>0$.
\end{enumerate}
\end{ass}
This assumption asks for the Lasso penalty to go to zero slightly slower
than $n^{-r}$, where $r$ is the rate from Assumption \ref{ass:est_convergence}. 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{ass:span} requires
that the dictionaries approximate well the correction-term Riesz representers (living in $\Delta_1$ and $\Delta_2(g_0)$, respectively). Second,
Assumption~\ref{ass:est_convergence} 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}$ \citep[see][]{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})$ \citep[c.f. Assumption~2 in][]{chernozhukov2022locally}.
We now provide assumptions on the moment condition. The first is a
mean-square consistency condition similar to Assumption~1 in \cite{chernozhukov2022locally}:
\begin{ass} \label{ass:sq_convergence} \mbox{} \\[-20pt]
\begin{enumerate}[label=\textbf{\alph*.},ref=\ref{ass:sq_convergence}.\alph*]
\item \label{ass:sq_convergence:m_bound} $\mathbb{E}[|m(W, g_0,h_0,\theta_0)|^2]<\infty$.
\item \label{ass:sq_convergence:m_conv} $\int |m(w,\hat{g}_\ell,\hat{h}_\ell,\theta_0)-m(w,g_0,h_0,\theta_0)|^2 dF_0(w)\xrightarrow{P} 0$.
\item \label{ass:sq_convergence:m_conv_theta} $\int |m(w,\hat{g}_\ell,\hat{h}_\ell,\tilde{\theta}_\ell)-m(w,\hat{g}_\ell,\hat{h}_\ell,\theta_0)|^2 dF_0(w)\xrightarrow{P} 0$.
\item \label{ass:sq_convergence:var_bound} $\mathbb{E}[(S-h_0(X,V))^2|D,Z]$ and $\mathbb{E}[\epsilon(W, g_0)^2|Z]$ are bounded almost surely.
\end{enumerate}
\end{ass}
Assumption~\ref{ass:sq_convergence:m_bound} is necessary for regular
estimation of $\theta_0$. Assumptions~\ref{ass:sq_convergence:m_conv} and
\ref{ass:sq_convergence:m_conv_theta} are mean-square consistency conditions
for the moment condition. Boundedness of the conditional errors
(Assumption~\ref{ass:sq_convergence:var_bound}) 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
{ass:bounded_dicts} and \ref{ass:abs_sumability:main}.
We also need to strengthen Assumption~\ref{rho} to control the remainder for linealizing the generalized error:
\begin{ass}
\label{ass:general_erro}
The mapping $g \mapsto \mathbb{E}[\epsilon(W, g)]$ is Frechet differentiable at $g_0$, with derivative $D_e$. Moreover, the Riesz representer of the derivative ($r_{e}$) satisfies $r_{e}(z)<0$ and is bounded and bounded away from zero. The remainder from the linearization is quadratic: there exists $\varepsilon > 0$ and $C> 0$ such that, for $\lVertg-g_0\rVert_2 < \varepsilon$,
\begin{equation*}
\lvert \mathbb{E}[\epsilon(W, g) - \epsilon(W, g_0) - D_e(g-g_0)] \rvert \leq C \lVertg-g_0\rVert_2^2
\end{equation*}
\end{ass}
The following assumption is standard in the literature, see, e.g.,
\cite{newey1994asymptotic}. It imposes a quadratic remainder bound for the
first-order linearization of $\bar m(g,h)$ and therefore strengthens
Assumptions~\ref{ass:Diff_m_h} and \ref{ass:Diff_m_g}. 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:
\begin{ass}\label{ass:linear_approx}
For each $\ell=1,\dots,L$, one of the following conditions holds:
\begin{enumerate}[label=(\roman*)]
\item For a $C>0$, with probability tending to one, $|\bar{\psi}(\hat g_\ell,\hat h_\ell)|
\le
C(
\|\hat g_\ell-g_0\|_2^2
+
\|\hat h_\ell-h_0\|_2^2)$ or
\item $\sqrt{n}\,\bar{\psi}(\hat g_\ell,\hat h_\ell) \xrightarrow{P} 0$.
\end{enumerate}
\end{ass}
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{lma:generalized_HR}, 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{lma:generalized_HR}). Define $\hat{\nu}_\ell=\nu(\hat{h}_\ell, \hat{\alpha}_{2\ell})$ and $\nu_0=\nu(h_0, \alpha_{02})$.
\begin{ass}
\label{ass:gen_reg} \mbox{}
Either $\alpha_{02}(x, v) = 0$ or
\begin{enumerate}[label=\textbf{\alph*.},ref=\ref{ass:gen_reg}.\alph*]
\item \label{ass:gen_reg:phi_Fderiv} The mapping $g \mapsto \varphi(d, z, g)$, from $\Delta_1$ to $L_2(D,Z)$, is twice Frechet differentiable with continuous second derivative.
\item \label{ass:gen_reg:twice_deriv} $h_0$ and $\alpha_{02}$ are almost surely twice continuously differentiable with respect to $v$.
\item \label{ass:gen_reg:bounded_deriv} $\partial h_{0}/\partial v$ and $\partial \alpha_{02}/\partial v$ are almost surely bounded.
\item \label{ass:gen_reg:bounded_alpha2} $\hat{h}_\ell \in \Delta_2(g_0)$, almost surely. Moreover, $\hat{h}_\ell$ is almost surely twice continuously differentiable with respect to $v$.
\item \label{ass:gen_reg:continuous_deriv} $\lVert\hat{\nu}_\ell-\nu_0\rVert_2 \lVert \hat{g}_\ell-g_0 \rVert_2 = o_p(n^{-1/2}).$
\end{enumerate}
\end{ass}
First, under downstream local robustness, the indirect effect is zero, and the above conditions are not needed. Regarding these conditions, Assumption~\ref{ass:gen_reg:phi_Fderiv} asks for a quadratic remainder in the linearization of the generated regressor. Assumptions~\ref{ass:gen_reg:phi_Fderiv}-\ref{ass:gen_reg:bounded_deriv} strengthen Assumption~\ref{ass:diff_h_alpha_phi}. Assumption~\ref{ass:gen_reg:bounded_alpha2} also strengthens Assumption~\ref{ass:inclusion} (see Appendix~\ref{sec:app_inclusion} for a general discussion of this assumption). Assumption~\ref{ass:gen_reg:continuous_deriv} 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{ass:gen_reg:continuous_deriv} in Section~\ref{sec:CASF_asymp}, where we verify it for the CASF, and in Appendix~\ref{sec:lasso_derivatives}.
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.
\begin{ass} \label{ass:jacobian}
There exists a neighborhood $\mathcal{N}$ of $\theta_0$ such that, for small $\lVert g - g_0 \rVert_2$ and $\lVert h - h_0 \rVert_2$:
\begin{enumerate}[label=\textbf{\alph*.}, ref=\ref{ass:jacobian}.\alph*, series=JAC]
\item \label{ass:jacobian:diff} $m(W, g, h, \theta)$ is almost surely differentiable in $\mathcal{N}$.
\item \label{ass:jacobian:bound}There exists a $C>0$ and a function $d(W, g,h)$, with $\mathbb{E}[d(W,g,h)]<C$, such that for $\theta\in\mathcal{N}$
\begin{equation*}
\left\lVert \frac{\partial m}{\partial\theta}(W, g,h, \theta) - \frac{\partial m}{\partial\theta}(W, g,h, \theta_0)\right\lVert_\infty \leq d(W, g,h) \lVert \theta - \theta_0 \rVert_\infty^{1/C} \text{ almost surely}.
\end{equation*}
\end{enumerate}
Moreover, we assume that:
\begin{enumerate}[resume*=JAC]
\item \label{ass:jacobian:existence} The expectation of the Jacobian, $\partial m/\partial\theta$, exists.
\item \label{ass:jacobian:convergence} It holds that
\begin{equation*}
\int\left\lVert\frac{\partial m}{\partial\theta}(w, \hat{g}_\ell,\hat{h}_\ell, \theta_0) - \frac{\partial m}{\partial\theta}(w, g_0,h_0, \theta_0)\right\rVert_\infty dF_0(w)\xrightarrow{P} 0.
\end{equation*}
\end{enumerate}
\end{ass}
Assumptions~\ref{ass:abs_sumability}, \ref{ass:est_convergence}, \ref
{ass:sq_convergence}, and \ref{ass:linear_approx} are stated for a single
moment condition. In the presence of more than one condition, they must be understood to hold componentwise. Assumption~\ref{ass:jacobian}, 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
\begin{equation*}
M \equiv \mathbb{E}\left[\frac{\partial m}{\partial\theta}
(W,g_0,h_0,\theta_{0})\right] \text{ and }
\Psi \equiv \mathbb{E}[\psi(W,g_0,h_0,\alpha_{0},\theta_{0})\psi(W,g_0,h_0,
\alpha_{0},\theta_{0})^{\prime }].
\end{equation*}
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
\begin{equation*}
\hat{M} \equiv \frac{1}{n}\sum_{\ell =1}^{L}\sum_{i\in I_{\ell }}\frac{
\partial m}{\partial \theta }(W_{i},\hat{g}_{\ell },\hat{h}_{\ell },\tilde{
\theta}_{\ell })\text{ and }
\hat{\Psi} \equiv \frac{1}{n}\sum_{\ell =1}^{L}\sum_{i\in I_{\ell }}\hat{
\psi}_{i\ell }(\tilde{\theta}_{\ell })\hat{\psi}_{i\ell }(\tilde{\theta}
_{\ell })^{\prime }.
\end{equation*}
The following theorem ensures asymptotic normality of $\sqrt{n}(\hat{\theta}-\theta _{0})$:
\begin{thm} \label{thm:asymptotics}
Consider that Assumptions \ref{ass:span}-\ref{ass:jacobian} are satisfied, $\hat{\Upsilon}\xrightarrow{P}\Upsilon$, and $M'\Upsilon M$ is non-singular. Then, the three-step debiased GMM estimator in equation~\eqref{dgmm} satisfies
\begin{equation*}
\sqrt{n}(\hat{\theta}-\theta_0) \xrightarrow{D} N(0, \Xi).
\end{equation*}
Moreover, the plug-in estimator for the asymptotic variance is consistent: $\hat{\Xi} \xrightarrow{P} \Xi$.
\end{thm}
\subsection{Regularity conditions for some examples}
\subsubsection{Hd-PS regression adjustment}
\label{sec:asymp_hdps}
Here we verify Assumptions~\ref{ass:sq_convergence}, \ref{ass:general_erro}, \ref{ass:linear_approx}, \ref{ass:gen_reg}, and \ref{ass:jacobian}. 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{ass:Diff_m_h} and \ref{ass:Diff_m_g} required for linearizing the moment condition. The moment condition that identifies $\theta_0$ in the partially linear model is:
\begin{equation*}
m(w, g, h, \theta) = \left[ y - h(\Lambda(g(z))) - \theta (d - \Lambda(g(z))) \right] \cdot \left[ d - \Lambda(g(z)) \right],
\end{equation*}
were recall that $\Lambda$ stands for the logistic cdf and $V=\Lambda(g_0(z))$.
The following assumption gives the result:
\begin{ass} \label{ass:regularity_ATE} \mbox{} \\[-20pt]
\begin{enumerate}[label=\textbf{\alph*.},ref=\ref{ass:regularity_ATE}.\alph*]
\item \label{ass:regularity_ATE_compact} The propensity score $V=\Lambda(g_0(Z))$ is bounded away from $0$ and $1$.
\item \label{ass:regularity_ATE_expectation} $\mathbb{E}[Y^2 |D, Z]$ is bounded almost surely.
\item \label{ass:regularity_ATE_dictionary} The atoms in the dictionary $\mathbf{b}_J$ are continuously differentiable.
\end{enumerate}
\end{ass}
Assumption~\ref{ass:regularity_ATE_compact} is the usual overlap assumption. Assumption~\ref{ass:regularity_ATE_expectation} 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{ass:regularity_ATE_expectation}. Assumption~\ref{ass:regularity_ATE_dictionary} 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{thm:asymptotics} holds in the Hd-PS setting:
\begin{prop} \label{prop:ATE_conditions}
Suppose that the convergence conditions in Assumption~ \ref{ass:est_convergence:estimators} hold and assume the existence of a consistent preliminary estimator $\tilde{\theta}_\ell \xrightarrow{P} \theta_0$. Then, Assumption~\ref{ass:regularity_ATE} guarantees that Assumptions~\ref{ass:sq_convergence}, \ref{ass:general_erro}, \ref{ass:linear_approx}, \ref{ass:gen_reg}, and \ref{ass:jacobian} are satisfied for the moment condition identifying $\theta_0$ in the Hd-PS regression adjustment within the partially linear model.
\end{prop}
\subsubsection{CASF in a non-separable model}
\label{sec:CASF_asymp}
Here we verify Assumptions~\ref{ass:sq_convergence}, \ref{ass:general_erro}, \ref{ass:linear_approx}, \ref{ass:gen_reg}, and \ref{ass:jacobian} 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.
\begin{ass} \label{ass:CASF_conditions} \mbox{} \\[-20pt]
\begin{enumerate}[label=\textbf{\alph*.},ref=\ref{ass:CASF_conditions}.\alph*]
\item \label{ass:CASF_conditions_alpha2} $(X, V)$ has joint density $f_{xv}$ w.r.t. an absolutely continuous measure $\mu$, with $\mu(x,v) = \mu_x(x) \times \mu_v(v)$. Also, $F^*$ has density $f^*$ w.r.t. $\mu_x$. With $f_v$ being the marginal density of $V$, it holds that $f^* f_{v} / f_{xv}$ is almost surely bounded and twice continuously differentiable with bounded first derivative w.r.t. $v$.
\item \label{ass:CASF_conditions_expectation} $\mathbb{E}[Y^2 |D, Z]$ and $\mathbb{E}[D^2|Z]$ are bounded almost surely.
\item \label{ass:CASF_conditions_h} $h_0$ is twice differentiable w.r.t. $v$, with $\partial h_0/\partial v$ and $\partial^2 h_0/\partial v^2$ bounded almost surely.
\item \label{ass:CASF_conditions_dictionary} The atoms in the dictionary $\mathbf{b}_J$ are twice differentiable w.r.t. $v$ with $\partial\mathbf{b}_J/\partial v$ and $\partial^2\mathbf{b}_J/\partial v^2$ bounded almost surely.
\item \label{ass:CASF_conditions_deriv} $\lVert\partial\hat{h}_\ell/\partial v - \partial h_0 / \partial v\rVert_2 \lVert\hat{g}_\ell - g_0\rVert_2=o_p(n^{-1/2})$ and $\lVert\partial\hat{\alpha}_{2\ell}/\partial v - \partial \alpha_{02} /\partial v\rVert_2 \lVert\hat{g}_\ell - g_0\rVert_2=o_p(n^{-1/2})$.
\end{enumerate}
\end{ass}
Assumption~\ref{ass:CASF_conditions_alpha2} 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{ass:CASF_conditions_expectation} is the usual bounded conditional variance assumption. Assumption~\ref{ass:CASF_conditions_h} also imposes smoothness conditions on $h_0$, while Assumption~\ref{ass:CASF_conditions_dictionary} requires the dictionary used to estimate $h_0$ to satisfy the same smoothness conditions.
Assumption~\ref{ass:CASF_conditions_deriv} 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{sec:CASF_details}.
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})$ \citep{bickel2009simultaneous}.
Hence, Assumption~\ref{ass:CASF_conditions_deriv} 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{sec:lasso_derivatives}.
We can then show that the assumptions for
Theorem~\ref{thm:asymptotics} hold for the moment condition defining the CASF.
\begin{prop} \label{prop:CASF_conditions}
Suppose that the convergence conditions in Assumptions~\ref{ass:span}-\ref{ass:rates} hold and assume the existence of a consistent preliminary estimator $\tilde{\theta}_\ell \xrightarrow{P} \theta_0$. Then Assumption~\ref{ass:CASF_conditions} guarantees that Assumptions~\ref{ass:sq_convergence}, \ref{ass:general_erro}, \ref{ass:linear_approx}, \ref{ass:gen_reg}, and \ref{ass:jacobian} are satisfied for the moment condition identifying the CASF.
\end{prop}
\begin{comment}
\subsection{CASF}
\subsubsection{Data Generating Process and Estimators}
This subsection describes the Monte Carlo simulation to evaluate the finite sample properties of the CASF estimator proposed in this paper. The DGP is
\begin{equation*}
(Z,U,V)\sim N\left(0,
\begin{bmatrix}
\operatorname{Id}_{6} & 0 & 0 \\
0 & 1 & 1/2 \\
0 & 1/2 & 1
\end{bmatrix}
\right),
\end{equation*}
where $\operatorname{Id}_6$ denotes the $6\times6$ Identity Matrix. Therefore, $Z$
is a 6-dimensional random vector. The correlation between $U$ and $V$ is $1/2
$. Note that the fact that $Z \perp U$ and $Z \perp V$ guarantees that the
Control Function Assumption is satisfied. Here $X=(D,Z_{1},...,Z_{5})$.
Both $D$ and $Y$ are generated by the following linear models:
\begin{align*}
Y &= \sum_{k=1}^5 Z_k + 2D + U \text{ and } \\
D &= \sum_{k=1}^6 Z_k + V.
\end{align*}
So $Z_6$ is excluded from the structural equation (i.e., it does not
directly affect $Y$) and may be used as an instrument.
We estimate the CASF for the following counterfactual distribution $F_X^*$:
(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 first step is $
g_{0}(z)=\sum_{k=1}^6 z_k$, and the second step is $h_{0}(x,v)=\sum_{k=1}^5
z_k+2d+v/2$.
We note here that, even if the model considered is linear, the second-step
correction nuisance parameter is highly non-linear. Letting $s\equiv
\sum_{k=1}^5z_k$, the Riesz representer is
\begin{equation*}
r_{2}(z_1,\dots,z_5,d,v)= C \cdot \exp\left(-\frac14 -\frac{s}{2} +\frac{d}{2
}+\frac{s^2}{4}+\frac{v^2}{2} +\frac{d^2}{4} -\frac{sd}{2}+sv-dv\right),
\end{equation*}
for a constant $C$. The function $\alpha_{02}$ is the orthogonal projection
of $r_{2}$ onto $\Delta_{2}(g_{0})$.
We display results for three different estimators of the CASF:
\begin{itemize}
\item The naive plug-in estimator: $\hat{\theta}_{PI}\equiv
n^{-1}\sum_{i=1}^n m(W_i, \hat{g}, \hat{h})$, where $m$ is given in equation~
\eqref{eq:CASF_moment}.
\item A cross-fitted Doubly Robust debiased estimator that only corrects for
the effect of pluging-in $\hat{h}$: $\hat{\theta}_{DR}$ is as in equation~
\eqref{eq:CASF_estimator} but with $\phi(W_{i},\hat{g}_{\ell}, \hat{h}_\ell,
\hat{\alpha}_{\ell},\tilde{\theta}_{\ell})$ replaced by $\hat{\alpha}
_{2\ell}(X_i, \hat{V}_{i\ell})\cdot(Y_i-\hat{h}_\ell(X_i, \hat{V}_{i\ell}))$
. That is, the correction term for the first step is omitted.
\item The cross-fitted fully Locally Robust debiased estimator: $\hat{\theta}
_{LR}$ as in equation~\eqref{eq:CASF_estimator}. That is, the estimator is
based on the fully debiased moment condition in equation~\eqref{eq:dmon_i}.
\end{itemize}
Numerical integration, with a sample of size $S=10^7$, is used to compute
the integrals w.r.t. $F_X^*$. The estimators for the nuisance parameters $g_0
$, $h_0$, $\alpha_{01}$, and $\alpha_{02}$ are Lasso with three
dictionaries: one that includes linear terms, another including linear and
quadratic terms, and a last one including linear, quadratic, and interaction
terms. The number of splits for cross-fitting is $L=5$ for every sample size.
To perform inference with each estimator, we present results that parallel
common practice. The fully debiased estimator uses the correct asymptotic
variance, the one accounting for first and second step estimation. This is
given by equation~\eqref{eq:asymp_variance}. The estimator $\hat{\theta}_{DR}
$ only accounts for the second step when computing the asymptotic variance
(as it does for estimation). Its asymptotic variance can be constructed by
replacing $\phi(W_{i},\hat{g}_{\ell}, \hat{h}_\ell, \hat{\alpha }_{\ell},
\tilde{\theta}_{\ell})$ by $\hat{\alpha}_{2\ell}(X_i, \hat{V}
_{i\ell})\cdot(Y_i-\hat{h}_\ell(X_i, \hat{V}_{i\ell}))$ in the second step
IF. To emphasize that plug-in estimation leads to an asymptotic bias
problem, confidence intervals for the plug-in estimator are built with
correct asymptotic variance (the one in equation~\eqref{eq:asymp_variance}).
\subsubsection{Results}
The next tables report results for a Monte Carlo simulation with $B=1098$
replications. Each table gives results for a different dictionary: linear,
quadratic or the one which also includes interaction terms.
\begin{table}[!htb]
\centering
\begin{tabular}{r|ccc|ccc|ccc}
\multicolumn{1}{l|}{} & \multicolumn{3}{c|}{Mean Absolute Bias} &
\multicolumn{3}{c|}{Standard Error} & \multicolumn{3}{c}{Coverage (95\%)} \\
\cline{2-10}
n & PI & DR & LR & PI & DR & LR & PI & DR & LR \\ \hline
100 & 0.2285 & 0.1463 & 0.1482 & 0.1649 & 0.1812 & 0.1733 & 0.6388 & 0.8681
& 0.8626 \\
500 & 0.1425 & 0.0594 & 0.0516 & 0.0685 & 0.0692 & 0.0645 & 0.3876 & 0.8954
& 0.9208 \\
1000 & 0.1236 & 0.0435 & 0.0376 & 0.049 & 0.0488 & 0.0462 & 0.2266 & 0.8744
& 0.9272 \\
5000 & 0.0852 & 0.0234 & 0.0169 & 0.0227 & 0.0212 & 0.0202 & 0.0227 & 0.7925
& 0.9290 \\
10000 & 0.0746 & 0.0181 & 0.0123 & 0.017 & 0.0148 & 0.0142 & 0.0018 & 0.7489
& 0.9163 \\ \hline
\end{tabular}
\caption{CASF results for the dictionary including linear terms.}
\label{tab:linear}
\end{table}
Tables \ref{tab:linear} and \ref{tab:quadratic} present results for the
linear and quadratic dictionaries, respectively. Correcting for the second
step already reduces a large amount of the bias of the plug-in estimator.
Adding the first-step correction further decreases bias. As shown in the
tables, however, the estimator accounting only for the second step fails to
keep coverage at the nominal 95\% level as the sample size increases.
\begin{figure}[!htb]
\minipage{0.32\textwidth} \includegraphics[width=
\linewidth]{Plots/L_MODEL_gquad_hquad_pi.png} \endminipage\hfill
\minipage{0.32\textwidth} \includegraphics[width=
\linewidth]{Plots/L_MODEL_gquad_hquad_db2.png} \endminipage\hfill
\minipage{0.32\textwidth} \includegraphics[width=
\linewidth]{Plots/L_MODEL_gquad_hquad_dbf.png} \endminipage
\caption{Distribution of the CASF estimators using a quadratic dictionary
for a sample size of $n=10000$.}
\label{fig:quad_n10000}
\end{figure}
The tables highlight that the plug-in estimator suffers from severe
asymptotic bias issues: coverage decreases rapidly, even if the confidence
interval is constructed with correct standard errors. Indeed, Figure~\ref
{fig:quad_n10000} shows that, for a sample size of $n=10000$, the
distribution of plug-in estimators has almost zero mass near the true
parameter $\theta_0=2$.
\begin{table}[!htb]
\centering
\begin{tabular}{r|ccc|ccc|ccc}
\multicolumn{1}{l|}{} & \multicolumn{3}{c|}{Mean Absolute Bias} &
\multicolumn{3}{c|}{Standard Error} & \multicolumn{3}{c}{Coverage (95\%)} \\
\cline{2-10}
n & PI & DR & LR & PI & DR & LR & PI & DR & LR \\ \hline
100 & 0.2764 & 0.1898 & 0.1957 & 0.1766 & 0.2003 & 0.2 & 0.5532 & 0.7534 &
0.7561 \\
500 & 0.1462 & 0.0603 & 0.0527 & 0.0692 & 0.0709 & 0.0663 & 0.3794 & 0.8963
& 0.9327 \\
1000 & 0.1214 & 0.0446 & 0.0374 & 0.0491 & 0.0488 & 0.0465 & 0.2329 & 0.8717
& 0.929 \\
5000 & 0.079 & 0.0259 & 0.0161 & 0.0236 & 0.0212 & 0.0202 & 0.0437 & 0.737 &
0.9354 \\
10000 & 0.0694 & 0.0209 & 0.0115 & 0.0172 & 0.0148 & 0.0143 & 0.0036 & 0.6533
& 0.9327 \\ \hline
\end{tabular}
\caption{CASF results for the dictionary including linear and quadratic
terms.}
\label{tab:quadratic}
\end{table}
Table~\ref{tab:interaction} displays results for the dictionary that also
includes interaction terms. The results are striking, as the Doubly Robust
estimator that only accounts for the second step performs well. It is able
to keep coverage at nominal levels, outperforming the fully debiased
estimator. Nevertheless, the decrease in coverage is small: 1-2\% for
intermediate sample sizes ($n=500$ and $1000$) and 4-5\% for large samples ($
n=5000$ and $10000$). We believe that this fact rests on the dictionaries
performing well to estimate the second step correction, but notably worst to
estimate the more complex fist step correction. This result suggest that the
complexity of the dictionary must be increased faster when accounting for
the first step.
\begin{table}[!htb]
\centering
\begin{tabular}{r|ccc|ccc|ccc}
\multicolumn{1}{l|}{} & \multicolumn{3}{c|}{Mean Absolute Bias} &
\multicolumn{3}{c|}{Standard Error} & \multicolumn{3}{c}{Coverage (95\%)} \\
\cline{2-10}
n & PI & DR & LR & PI & DR & LR & PI & DR & LR \\ \hline
100 & 0.3447 & 0.3583 & 0.3857 & 0.179 & 0.3606 & 0.4054 & 0.9016 & 0.7969 &
0.806 \\
500 & 0.1707 & 0.0997 & 0.1085 & 0.0693 & 0.1213 & 0.1359 & 0.7925 & 0.9481
& 0.9227 \\
1000 & 0.1367 & 0.0635 & 0.0674 & 0.0488 & 0.0776 & 0.0849 & 0.5883 & 0.9372
& 0.9262 \\
5000 & 0.0846 & 0.0259 & 0.0283 & 0.0229 & 0.0312 & 0.0349 & 0.1383 & 0.9372
& 0.8926 \\
10000 & 0.0729 & 0.0173 & 0.0205 & 0.017 & 0.0207 & 0.0228 & 0.0337 & 0.9399
& 0.8926 \\ \hline
\end{tabular}
\caption{CASF results for the dictionary including linear, quadratic, and
interaction terms.}
\label{tab:interaction}
\end{table}
\end{comment}
\section{Conclusion}
\label{sec: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.
\newpage \addcontentsline{toc}{section}{References} \makeatletter
\makeatother
\bibliographystyle{apalike}
\bibliography{references}
\cleardoublepage