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.
115,824 characters
Locally Robust Inference with Outcomes Predicted by Machine Learning
\maketitle
\begin{abstract} Many parameters of interest in economics depend on outcomes generated by machine learning (ML).These generated outcomes are subsequently used in standard econometric estimation and typically treated as observed, a practice that threatens the validity of downstream inference. Empirically relevant cases include dynamic treatment effects under sequential exogeneity, long-term effects with surrogates, the partially linear model with an ML-predicted outcome, productivity persistence with an unrestricted production function and constructed instruments, and outcomes residualized of confounding variation. This paper develops a three-step locally robust GMM framework for valid inference when the dependent variable is generated by ML. We construct a debiased moment condition that augments the identifying moment with first- and second-step influence functions, making inference on the parameter of interest locally insensitive to regularization and model selection biases at both steps. We establish consistency and asymptotic normality of the proposed
estimator and illustrate its finite-sample performance through Monte
Carlo simulations. We also apply the method to U.S.\ mortgage data to
study conditional Black--White disparities in a race-blind ML
default-risk score. The score predicts subsequent default well overall,
but the conditional disparity embedded in the score is much smaller
than the corresponding disparity in realized serious delinquency;
accounting for first-step estimation also raises the naive
out-of-sample estimate of the score disparity.
\vspace{1em}
\noindent\textbf{Keywords:} Machine learning, generated outcomes, debiased GMM. \\
\textbf{JEL Classification:} C13, C14. \end{abstract}
\newpage
\section{Introduction}
The rapid growth in rich, high-dimensional data availability has created new opportunities for empirical economic research. Examples include unstructured text, images, and satellite imagery, from which machine learning (ML) methods have become increasingly common. In particular, ML is often used to construct key variables, whether regressors or outcomes, that are subsequently employed using standard econometric techniques. However, these constructed variables, known as generated variables, depend fundamentally on first-stage estimation procedures, threatening valid inference.\par
Generated dependent variables arise structurally in many applications in economics. In dynamic treatment effects under sequential exogeneity, the counterfactual outcome under a fixed treatment path is unobserved and must be generated at each stage. Empirical examples include the impact of sequential training programmes for the unemployed in West Germany \citep{lechner2010identification} and the evaluation of vocational training sequences for disadvantaged youth under the U.S. Job Corps program \citep{bodory2022evaluating}. In long-term treatment effects with surrogates, a long-run outcome is unavailable in the experimental sample and is instead predicted from short-term proxies observed alongside the treatment, as in the effect of early childhood classroom assignment on predicted adult earnings \citep{chetty2011does} and the long-run employment effects of welfare-to-work programs \citep{athey2025surrogate}.\par
Generated outcomes also arise in empirical research, where the parameter of interest depends not on a directly observed outcome but on one constructed in a preliminary estimation step. A growing part of the empirical literature already relies on ML to construct outcomes of interest. For instance, researchers use remote sensing and other high-dimensional data to construct variables such as consumption, GDP, and agricultural yields \citep{jean2016combining, henderson2012measuring, burke2017satellite, graetz2018mapping}. In such contexts, the parameter of interest is typically the effect of some key variable on the constructed outcome, controlling for confounders. In other settings, practitioners rely on parametric specifications to construct the outcome. For inequality of opportunity measurement, e.g., \cite{brunori2025inequality}, income is purged of life-cycle variation or cohort fixed effects before estimating how much of its remaining variation is attributable to circumstances beyond individual control, such as parental background or place of birth. In partial insurance models, e.g., \cite{blundell2008consumption}, consumption is predicted from auxiliary data before being used to study how households smooth income shocks. Whether the outcome is constructed by machine learning or by parametric methods, its preliminary estimation threatens the validity of inference. \par
Machine learning offers an appealing choice for constructing generated outcomes. Specifically, it handles high-dimensional covariates flexibly and avoids the misspecification bias that parametric specifications introduce. However, ML estimation leads to regularization and model selection biases that standard inference procedures do not account for, threatening the validity of inference when ML-generated outcomes are used as dependent variables. The empirical literature exploiting ML's flexibility typically plugs ML-generated outcomes directly into standard estimators without correcting for regularization bias, yielding biased point estimates and confidence intervals that fail to achieve nominal coverage. Relying instead on parametric specifications avoids regularization bias at the cost of misspecification bias that propagates through all subsequent stages. This trade-off motivates the development of a framework that preserves the flexibility of machine learning while delivering valid inference despite the regularization bias it introduces.\par
This paper develops a framework for locally robust inference with outcomes predicted by machine learning. We consider a three-step setting in which the dependent variable is generated in a first step via machine learning regression of an observed outcome on a vector of covariates. In a second step, a nuisance function is constructed using possibly another machine learning projection of the generated outcome. Finally, the parameter of interest is identified in a third step through a Generalized Method of Moments (GMM) criterion that integrates the outputs of the first two steps with the observed data. By accommodating ML-generated outcomes, our framework enables researchers to exploit the flexibility of machine learning while delivering valid inference.\par
A closely related literature studies inference when machine-learning
predictions are used as outcomes in subsequent statistical analysis.
\cite{wang2020postprediction} show that treating predicted outcomes as
observed can distort downstream inference and propose post-prediction
corrections based on separate training, testing, and validation samples.
\cite{angelopoulos2023prediction} develop prediction-powered inference,
which combines labeled observations with predictions from an arbitrary
machine-learning system to conduct inference on features of the
underlying observed outcome distribution. These approaches establish
that prediction error must be accounted for in downstream inference.
Our problem is complementary. The population prediction function
$g_0$ is itself an unknown nuisance entering the parameter of interest,
and it may also generate a subsequent nuisance $h_0(g_0)$. We derive
orthogonal GMM moments that remove the first-order effects of estimating
these objects by machine learning. In particular, ordinary sample
splitting or cross-fitting makes predictions out of sample but does not
by itself eliminate the first-order effect of estimating $g_0$ on a
downstream parameter defined through that function.\par
Our paper also builds on the literature on Neyman-orthogonal inference
with high-dimensional nuisance parameters. \cite{chernozhukov2018double}
and \cite{chernozhukov2022locally} provide general methods for valid
inference with estimated nuisance functions, while
\cite{escanciano2023automatic} studies a three-step setting with
generated regressors. \cite{chernozhukov2023automatic} develop recursive
orthogonalization for nested functionals, including dynamic treatment
effects. Their recursive framework covers our first dynamic-treatment
example, which we use as a benchmark. Our contribution is a
generated-outcome GMM formulation in which the first-step prediction may
affect the identifying moment both through a downstream nuisance and
directly. The resulting total-derivative representation separates these
indirect and direct channels and yields their Riesz corrections in a
common framework. It also accommodates restricted population projection
spaces, as illustrated by the constructed-instrument application to
productivity persistence.
\par
Our proposal consists of constructing a locally robust/debiased moment condition that augments the original identifying moment with the first- and second-step influence functions. In doing so, we account for both the effect of the ML-generated outcome and the nonparametric nuisance estimation on inference. Specifically, making the moment locally robust to the second-step nuisance annihilates its indirect dependence on the generated outcome, so that the first-step correction needs only to address the remaining direct effect of the generated outcome on the moment condition. The resulting moment condition satisfies Neyman orthogonality with respect to both nuisance parameters, ensuring inference on the parameter of interest is locally insensitive to the regularization and model selection biases introduced by their ML estimation.\par
We illustrate our framework through three examples. We begin with dynamic treatment effects, where the parameter of interest depends on the generated outcome only through the second-step nuisance. In this setting, no direct first-step correction is required, and the first-step effect is only through the indirect effect on the second step. This example falls under the nested functional setting of \cite{chernozhukov2023automatic}, which arises here as a special case of our framework. We then extend the partially linear model to allow for an ML-predicted
dependent variable, where the generated outcome enters the identifying
moment directly and a nonzero first-step correction is generally
required. As a concrete example, we consider conditional group
disparities in algorithmic credit scoring: a protected characteristic
is excluded from an ML default-risk model, but the remaining borrower
information may proxy for group membership. In this case the first-step
Riesz representer measures the incremental ability of the credit-model
information set to predict group membership beyond the controls used in
the downstream comparison. Finally, we study the estimation of the persistence of unobserved
productivity with an unrestricted production function. Rather than
assuming the availability of a valid instrument, we construct
instruments by orthogonalizing observable period-1 productivity signals
with respect to additive functions of period-1 and period-2 production
inputs. Under primitive timing restrictions, this construction removes
the unrestricted production function and predictable components of
productivity from the identifying moment. We characterize when relevant
instruments of this form exist and then apply our locally robust
framework to account for ML estimation of both residualized productivity
and the constructed instrument.\par
We establish consistency and asymptotic normality of the proposed estimator under mild conditions. We study its finite-sample behavior through Monte Carlo simulations for the partially linear model. Our simulations suggest our Debiased ML (DML) estimator exhibits satisfactory finite-sample performance, with negligible bias that vanishes as the sample size increases and coverage rates close to nominal levels. In contrast, the naive plug-in estimator, which treats the ML-generated outcome as observed and applies no correction for regularization bias, exhibits sizable bias and coverage rates far below nominal, confirming that ignoring the generated nature of the outcome leads to invalid inference.\par
Motivated by the credit-scoring example, we provide an empirical
application using U.S.\ mortgage data. We estimate a race-blind ML
model of subsequent serious delinquency and study conditional
Black--White disparities after controlling for conventional
underwriting characteristics. The second-step outcome regression
is constructed from out-of-sample first-step predictions. In the
preferred analysis-weighted XGBoost specification, the locally
robust estimate of the predicted-risk disparity is $0.082$
percentage points, compared with $1.017$ percentage points for
realized serious delinquency. Treating the out-of-sample
predictions as observed yields $0.064$ percentage points, so the
locally robust adjustment raises the estimate by approximately
$28\%$. These estimates describe conditional disparities rather
than causal effects.
\par
The remainder of this paper proceeds as follows.
Section~\ref{sec:set} presents the three-step setting and the general
construction of locally robust moments in the presence of a generated
dependent variable. Section~\ref{sec:illustrate} illustrates the
framework through three examples. Section~\ref{sec:estim} develops the
debiased GMM estimator, establishes its consistency and asymptotic
normality, and examines its finite-sample performance through Monte Carlo
simulations. Section~\ref{sec:empirical} applies the method to conditional
group disparities in mortgage default risk, and
Section~\ref{sec:conc} concludes. Proofs are provided in the Appendix.
\section{A General Framework for Valid Inference with ML Generated Outcomes}\label{sec:set}
\subsection{Three-Step Setting}
We observe data $W = (Y, X)$ with cumulative distribution function (cdf) $F_0$, where $Y$ is an observed outcome and $X$ collects all available covariates. The parameter of interest $\theta_0$ is identified by a moment condition that depends on an outcome $\tilde{Y}$ not directly observed, but instead generated from $Y$ through a regression function $g_0$ of some components of $X$, estimated via machine learning. We develop locally robust inference for $\theta_0$ that remains valid despite the regularization and model selection biases introduced by ML estimation of $g_0$. We consider the case where the generated outcome is one-dimensional, although our framework naturally extends to the multivariate case by applying the construction to each component.\par
We describe our general three-step setting as follows:\par \bigskip
\textbf{First step.}
Let $X_1$ denote some components of the available covariates $X$, used to construct the generated outcome. There is a first step $g_0(X_1)$ solving the moment restrictions
\begin{align}
\label{eq:gen}
\mathbb{E}[\delta_1(X_1)\left(Y - g_0(X_1)\right)] = 0 \quad \text{for all } \delta_1 \in \Delta_1,
\end{align}
where $Y - g_0(X_1)$ is the prediction error for $g_0 \in \Delta_1$, and $\Delta_1$ is a linear and closed subspace of $L_2(X_1)$. 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$. If $\Delta_1 = L_2(X_1)$, then $g_0(X_1) = \mathbb{E}[Y|X_1]$, as when the researcher predicts an economic outcome from satellite imagery or constructs a counterfactual continuation value in a dynamic treatment regime. Since $\text{dim}(X_1)$ is allowed to be high, this nonparametric first step may not be feasible to implement directly. We can then consider a high-dimensional additive regression model with the same error with $\Delta_1 = \sum_{j=1}^{\text{dim}(X_1)}\Delta_{1,j}$, where $\Delta_{1,j}$ is a subset of $L_2(X_{1j})$ for $X_{1j}$ the $j$-th component of $X_1$. When $\Delta_1$ is the mean square limit of linear combinations $\sum_{k=1}^{K}\beta_{0k}c_k(X_1)$ for $K \in \mathbb{N}$, a sequence of real numbers $(\beta_{0k})_{k=1}^{\infty}$, and a dictionary $(c_k)_{k=1}^{\infty}$ of functions in $L_2(X_1)$, this setting covers high-dimensional linear regression, commonly estimated via Lasso or other ML methods.\par
The first step $g_0$ in \eqref{eq:gen} is used to construct the generated outcome
\begin{align}
\label{eq:varphi}
\tilde{Y} \equiv \tilde{Y}(g_0) \equiv \varphi(W, g_0),
\end{align}
where $\varphi$ is a known mapping from observed variables $W$ and the unknown function $g_0$ to the generated outcome $\tilde{Y}$. In the leading case $\varphi(W, g_0) = g_0(X_1)$, the generated outcome is directly the ML prediction of $Y$ given $X_1$. Alternatively, $\varphi(W, g_0) = Y - g_0(X_1)$ covers settings where the outcome of interest is residualized of confounding variation, such as income or consumption purged of life-cycle variation.\par\bigskip
\textbf{Second step.}
Let $S \equiv S(g_0)$ denote some components of $(\tilde{Y}(g_0), W)$, and $X_2$ denote some components of $X$ onto which $S$ is projected. The second step is an orthogonal projection of $S$ on $X_2$, identified by
\begin{align}
\label{eq:orth2}
\mathbb{E}[\delta_2(X_2)(S(g_0) - h_0(X_2))] = 0 \quad \text{for all } \delta_2 \in \Delta_2,
\end{align}
where, if $S(g_0)$ is $d_S$-dimensional,
$\Delta_2$ denotes the corresponding product Hilbert space
$\Delta_2^{(1)}\times\cdots\times\Delta_2^{(d_S)}$,
with each $\Delta_2^{(j)}$ a closed linear subspace of
$L_2(X_2)$. In the scalar case, this reduces to a single
closed linear subspace of $L_2(X_2)$. Here, $S(g_0)$ may be the generated outcome itself, a vector containing it alongside observed regressors, or a transformed version of it (as when partialling out confounders from an ML-predicted dependent variable or projecting a counterfactual outcome onto period-1 state variables). This dependence of $h_0$ on $g_0$ through $S(g_0)$ gives the problem a compositional structure that we exploit explicitly in the construction of locally robust moments and in the corresponding asymptotic analysis. \par \bigskip
\textbf{Third step.}
Let $\theta_0 \in \Theta \subseteq \mathbb{R}^{p}$ denote the $p$-dimensional parameter of interest, identified by $q \geq p$ moment conditions. The moment function $m\colon \mathbb{R}^{\text{dim}(W)}\times \Delta_1 \times \Delta_2 \times \Theta \rightarrow \mathbb{R}^{q}$ identifies $\theta_0$ through the GMM moment condition
\begin{align}
\label{eq:main}
\mathbb{E}[m(W,g_{0},h_{0},\theta_{0})]=0,
\end{align}
where $\theta_0$ is assumed to be the unique solution to $\mathbb{E}[m(W,g_{0},h_{0},\theta)]=0$ over $\theta \in \Theta$. The moment $m$ may depend on $g_0$ directly, as in the partially linear model where the dependent variable is an ML prediction, or only through $h_0$, as in dynamic treatment effects where the parameter of interest is identified by a marginal structural moment condition.\par
\subsection{Locally Robust Moment Construction}
The moment in equation (\ref{eq:main}) identifies the parameter $\theta_0$ as a function of the observed data $W$ and the nuisance parameters $g_0$ and $h_0$. Accordingly, we can consider the nuisance parameter $\gamma(F)\equiv\left(g(F), h(F,g(F))\right)$, where $h(F,g)$ denotes the second-step nuisance when $W$ follows the distribution $F$ and the first-step is fixed at $g$, so that the moment identifying the parameter of interest $\theta_0$ in equation (\ref{eq:main}) can be expressed as
\begin{align}
\label{eq:main2}\bar{m}(\gamma_0,\theta_0)\equiv\mathbb{E}[m(W,\gamma_0,\theta_0)]=0,
\end{align}
where $\bar{m}(\gamma,\theta)$ denotes the expected moment function, $\mathbb{E}\left[\cdot\right]$ is the expectation under $F_0$, $g_0\equiv g(F_0)$, $h_0\equiv h(F_0,g_0)$, and $\gamma_0\equiv(g_0, h_0)$ is the probability limit (plim) of $\hat{\gamma}$ under $F_0$.\par
The challenge of estimating $\theta_0$ using equation (\ref{eq:main2}) lies in the sensitivity of the identifying moment to estimation of $\gamma_0$. Formally, for an alternative distribution function $H$ and a regular path $F_\tau=(1-\tau)F_0+\tau H$ for $\tau \in [0,1]$, such that $\gamma(F_\tau)\equiv(g(F_\tau), h(F_\tau, g(F_\tau)))$ exists for $\tau$ small enough, the identifying moment is not generally locally robust to estimation of $\gamma_0$, i.e.,
\begin{align*}
\frac{d}{d\tau}\bar{m}(\gamma(F_\tau),\theta)\neq 0,
\end{align*}
where, henceforth, $d/d\tau$ is the derivative with respect to $\tau$ from the right and evaluated at $\tau=0$.\par
This motivates constructing a moment condition whose expectation is insensitive to local perturbations of $\gamma$ around $\gamma_0$, that is, one whose derivative with respect to $\tau$ vanishes at $\tau = 0$, a property known as Neyman orthogonality \citep{neyman1959optimal,chernozhukov2018double, chernozhukov2022locally}, or equivalently, local robustness with respect to $(g_0, h_0)$.\par
We derive the following standard form for the debiased/locally robust moment function:
\begin{align}
\label{eq:debiased}
\psi(w, g_0, h_0, \alpha_0, \theta)= m(w, g_0, h_0, \theta)+ \underbrace{\alpha_{01}(x_1;\theta) \cdot (y - g_0(x_1))}_{\equiv \phi_1(w,g_0,\alpha_{01},\theta)}
+ \underbrace{\alpha_{02}(x_2;\theta)' \left(s - h_0(x_2)\right)}_{\equiv \phi_2(w,h_0,\alpha_{02},\theta)},
\end{align}
where $\alpha_0\equiv(\alpha_{01},\alpha_{02})$ denotes the Riesz representers and \[ \phi\equiv\phi_1+\phi_2 \] denotes the total influence-function correction, with $\phi_1$ and $\phi_2$ the first- and second-step corrections, respectively. The first-step influence function corresponds to $\phi_1$ in equation (\ref{eq:debiased}) and captures the effect of the generated variable on the identifying moment. In general, it is nonzero, implying that inference procedures that ignore generated outcomes are invalid. Importantly, $\phi_1$ depends on the first-step Riesz representer $\alpha_{01}$ and the corresponding identifying error $\epsilon(w, g_0)=y - g_0(x_1)$. The second-step influence function, $\phi_2$, accounts for the effect of estimating $h_0$ on the moment condition and similarly depends on $\alpha_{02}$. When there are finitely many first-step regressions, as in the
production-function example, we interpret $\Delta_1$ as the corresponding
product Hilbert space. The first-step correction is then the sum of the
componentwise corrections, which follows directly from linearity of the
Hadamard derivative and the Riesz representation on the product space. \par
To construct locally robust/debiased estimators, we show that Neyman orthogonal moments are constructed by sequentially adding to the original identifying moment the second- and first-step IFs. This sequential construction exploits the idea of downstream local robustness of \cite{escanciano2023automatic}: making the moment locally robust to $h$ first annihilates the indirect effect of $g$ through $h$, so that the first-step IF needs only correct for the remaining direct effect of $g$ on the moment condition. To fix ideas, notice that the impact of $\gamma(F_\tau)$ on $\bar{m}$ is given by:
\begin{align*}
\frac{d}{d\tau}\bar{m}(\gamma(F_\tau),\theta) &= \frac{d}{d\tau}\bar{m}(g(F_\tau), h(F_0, g(F_\tau)),\theta) + \frac{d}{d\tau}\bar{m}(g_0, h(F_\tau, g_0),\theta) \nonumber
\\&= \underbrace{\frac{d}{d\tau}\bar{m}(g(F_\tau), h_0,\theta)}_{\text{direct effect of }g} + \underbrace{\frac{d}{d\tau}\bar{m}(g_0, h(F_0, g(F_\tau)),\theta)}_{\text{indirect effect of }g\text{ through }h} + \underbrace{\frac{d}{d\tau}\bar{m}(g_0, h(F_\tau, g_0),\theta)}_{\text{direct effect of }h} \nonumber
\\&= \frac{d}{d\tau}\bar{m}(g(F_\tau), h_0,\theta) + \frac{d}{dh}\bar{m}(g_0, h_0,\theta)\cdot\left[\frac{d}{d\tau}h(F_0, g(F_\tau)) + \frac{d}{d\tau}h(F_\tau, g_0)\right],
\end{align*}
where the last equality follows by the chain rule, under Hadamard differentiability of $\bar{m}(g,h,\theta)$ with respect to $h$ and of $h(F_0,g)$ with respect to $g$; see Proposition~3.1 of \citet{escanciano2023automatic}.\par
The last display highlights that if the moment condition is locally robust to $h$, meaning $$\frac{d}{dh}\bar{m}(g_0, h_0,\theta)=0,$$ then there is no indirect effect of $g$ through $h$. This is the idea of the downstream local robustness in \cite{escanciano2023automatic} applied to generated outcomes. Exploiting this insight, we proceed sequentially: first, we augment $m$ with the influence function
of $h$,
\begin{align*}
m^*(w, g_0, h_0, \alpha_{02}, \theta)\equiv m(w, g_0, h_0,\theta)+\phi_2(w,h_0,\alpha_{02},\theta),
\end{align*}
and then we augment $m^*$ with the influence function of $g$, yielding
\begin{align*}
\psi(w, g_0, h_0, \alpha_0, \theta)&\equiv m^*(w, g_0, h_0, \alpha_{02}, \theta)+\phi_1(w,g_0,\alpha_{01},\theta)
\\&=m(w, g_0, h_0,\theta)+\phi_1(w,g_0,\alpha_{01},\theta)+\phi_2(w,h_0,\alpha_{02},\theta).
\end{align*}
Notice that the locally robust moment depends on the Riesz representers
$\alpha_0 \equiv (\alpha_{01}, \alpha_{02})$, which are identified by the
orthogonality conditions (with the notation $h_0(g)\equiv h(F_0,g)$)
\begin{align*}
\mathbb{E}[\alpha_{01}(X_1;\theta)\delta_1(X_1)]
&=
\frac{d}{d\tau}
\mathbb{E}\!\left[
m\left(W, g_0+\tau\delta_1,
h_0(g_0+\tau\delta_1), \theta\right)
\right]\bigg|_{\tau=0}
\quad \text{for all } \delta_1 \in \Delta_1,
\\
\mathbb{E}[\alpha_{02}(X_2;\theta)'\delta_2(X_2)]
&=
\frac{d}{d\tau}
\mathbb{E}\!\left[
m\left(W,g_0,h_0(g_0)+\tau\delta_2,\theta\right)
\right]\bigg|_{\tau=0}
\quad \text{for all } \delta_2 \in \Delta_2.
\end{align*}
Here $\alpha_{01}$ represents the Riesz representer of the total
Gateaux derivative of the expected moment with respect to the first-step
nuisance $g$, including both the direct effect of $g$ on the moment and
the indirect effect through the induced change in $h_0(g)$. By contrast,
$\alpha_{02}$ represents the Riesz representer of the Gateaux derivative
with respect to the second-step nuisance $h$, holding $g_0$ fixed.
The following assumption imposes pathwise differentiability and Riesz representability conditions ensuring that the orthogonal corrections are well-defined. We begin by recalling the notion of Hadamard tangential differentiability.
\begin{definition}
Let $\mathbb{D}$ and $\mathbb{E}$ be normed vector spaces, let
$\mathbb{D}_0\subseteq\mathbb{D}$, and let
$\phi:\mathbb{D}_\phi\subseteq\mathbb{D}\to\mathbb{E}$.
The map $\phi$ is Hadamard differentiable at $\eta\in\mathbb{D}_\phi$
tangentially to $\mathbb{D}_0$ if there exists a continuous linear map
$D\phi(\eta):\mathbb{D}_0\to\mathbb{E}$ such that
\[
\frac{
\phi(\eta+t_n h_n)-\phi(\eta)
}{t_n}
\to
D\phi(\eta)[h]
\]
for every sequence $t_n\downarrow0$ and every sequence
$h_n\to h\in\mathbb{D}_0$ satisfying
$\eta+t_n h_n\in\mathbb{D}_\phi$ for all sufficiently large $n$.
The map $D\phi(\eta)$ is called the Hadamard derivative of $\phi$ at
$\eta$.
\end{definition}
\begin{assumption}
\label{ass:lr}
For each $\theta\in\Theta$, let
\[
\bar m(g,h,\theta)\equiv \mathbb{E}[m(W,g,h,\theta)]
\]
and define $h_0(g)\equiv h(F_0,g)$, with
$h_0\equiv h_0(g_0)$. The map
\[
(g,h)\mapsto \bar m(g,h,\theta)
\]
is Hadamard differentiable at $(g_0,h_0)$ tangentially to
$\Delta_1\times\Delta_2$. Moreover, the generated outcome map
$g\mapsto S(g)$, viewed as a map from $\Delta_1$ into the appropriate
$L_2(F_0)$ product space, is Hadamard differentiable at $g_0$
tangentially to $\Delta_1$.
\end{assumption}
Under Assumption \ref{ass:lr}, the second-step map $g\mapsto h_0(g)$
is also Hadamard differentiable. Indeed, since $h_0(g)$ is the $L_2$
projection of $S(g)$ onto $\Delta_2$, its derivative is the
corresponding projection of $D_gS(g_0)[\delta_1]$. Thus, denoting
\[
\dot h_g(\delta_1)\equiv D_gh_0(g_0)[\delta_1],
\]
the derivative $\dot h_g(\delta_1)\in\Delta_2$ is characterized by
\[
\mathbb{E}\left[
\delta_2(X_2)
\left\{
D_gS(g_0)[\delta_1]
-
\dot h_g(\delta_1)(X_2)
\right\}
\right]=0
\quad
\text{for all } \delta_2\in\Delta_2.
\]
Consequently, by the chain rule, the composite map
$g\mapsto \bar m(g,h_0(g),\theta)$ is Hadamard differentiable at $g_0$.
Define the partial derivative functionals
\[
\dot m_g(\delta_1;\theta)
\equiv
D_g\bar m(g_0,h_0,\theta)[\delta_1],
\qquad
\dot m_h(\delta_2;\theta)
\equiv
D_h\bar m(g_0,h_0,\theta)[\delta_2],
\]
and the total first-step derivative
\[
\dot m_{g,\mathrm{tot}}(\delta_1;\theta)
\equiv
\dot m_g(\delta_1;\theta)
+
\dot m_h(\dot h_g(\delta_1);\theta).
\]
Equivalently,
\[
\dot m_{g,\mathrm{tot}}(\delta_1;\theta)
=
\frac{d}{d\tau}
\mathbb{E}\!\left[
m\{W,g_0+\tau\delta_1,h_0(g_0+\tau\delta_1),\theta\}
\right]\bigg|_{\tau=0}.
\]
When $q>1$, the following Riesz representations are understood
componentwise, so that each component of the $q$-dimensional
moment has its own first- and second-step representers.
The Riesz representers are defined by
\begin{align}
\mathbb{E}[\alpha_{01}(X_1;\theta)\delta_1(X_1)]
&=
\dot m_{g,\mathrm{tot}}(\delta_1;\theta),
\quad
\text{for all } \delta_1\in\Delta_1,
\label{eq:riesz1}
\\
\mathbb{E}[\alpha_{02}(X_2;\theta)'\delta_2(X_2)]
&=
\dot m_h(\delta_2;\theta),
\quad
\text{for all } \delta_2\in\Delta_2.
\label{eq:riesz2}
\end{align}
\begin{theorem}
\label{thm:lr}
Suppose Assumption \ref{ass:lr} holds. Then there exist unique Riesz
representers $\alpha_0=(\alpha_{01},\alpha_{02})$, with
$\alpha_{01}\in\Delta_1^q$ and $\alpha_{02}\in\Delta_2^q$,
satisfying \eqref{eq:riesz1}--\eqref{eq:riesz2} componentwise. Define
\[
\psi(w,g,h,\alpha_0,\theta)
=
m(w,g,h,\theta)
+
\alpha_{01}(x_1;\theta)(y-g(x_1))
+
\alpha_{02}(x_2;\theta)'(s(g)-h(x_2)),
\]
where $s(g)$ denotes the realization of $S(g)$. Then
\[
\mathbb{E}[\psi(W,g_0,h_0,\alpha_0,\theta_0)]=0
\]
whenever
\[
\mathbb{E}[m(W,g_0,h_0,\theta_0)]=0.
\]
Moreover, $\psi$ is Neyman orthogonal: for every regular path
\[
g_\tau=g_0+\tau\delta_1+o(\tau),
\qquad
h_\tau=h_0+\tau\delta_2+o(\tau),
\]
with $\delta_1\in\Delta_1$ and $\delta_2\in\Delta_2$,
\[
\frac{d}{d\tau}
\mathbb{E}\!\left[
\psi(W,g_\tau,h_\tau,\alpha_0,\theta)
\right]\bigg|_{\tau=0}
=
0.
\]
\end{theorem}
\begin{remark}
\label{rem:linear}
In the leading case where $\Delta_2=L_2(X_2)$, the derivative of the second-step projection satisfies
\[
D_gh_0(g_0)[\delta_1](X_2)
=
\mathbb{E}[D_gS(g_0)[\delta_1]\mid X_2].
\]
Let $a_g(X_1;\theta)$ denote the Riesz representer of the direct
partial derivative $\dot m_g(\delta_1;\theta)$. Since
$\alpha_{01}$ represents the total derivative
\[
\dot m_{g,\mathrm{tot}}(\delta_1;\theta)
=
\dot m_g(\delta_1;\theta)
+
\dot m_h(\dot h_g(\delta_1);\theta),
\]
it follows that
\[
\alpha_{01}(X_1;\theta)
=
a_g(X_1;\theta)
+
\mathbb{E}[\alpha_{02}(X_2;\theta)\mid X_1]
\]
when $S(g)=g(X_1)$, while
\[
\alpha_{01}(X_1;\theta)
=
a_g(X_1;\theta)
-
\mathbb{E}[\alpha_{02}(X_2;\theta)\mid X_1]
\]
when $S(g)=Y-g(X_1)$. In particular, if $\alpha_{02}=0$, then
$\alpha_{01}=a_g(X_1;\theta)$, corresponding to the absence of an
indirect effect through the second-step nuisance, a special case of
downstream local robustness.
\end{remark}
\begin{remark}
For simplicity of notation, we suppress the dependence of the Riesz
representers $\alpha_{01}(\cdot;\theta)$ and
$\alpha_{02}(\cdot;\theta)$, and other nuisance objects on $\theta$ whenever no confusion may arise.
\end{remark}
\begin{comment}
A key simplification relative to \cite{escanciano2023automatic} is that $\Delta_2$, the second-step projection space defined in equation
\eqref{eq:orth2}, does not depend on $g_0$, so that $\alpha_{01}$ and
$\alpha_{02}$ can be identified and estimated separately. By contrast, in the case of ML generated regressors, the dependence of $\Delta_2(g)$
on $g_0$ requires a joint characterization of both representers. In the PLM, the structure of the moment condition is sufficiently simple that the Riesz representers can be derived analytically, as we show in the following section. In general, however, analytic expressions for $\alpha_{01}$ and $\alpha_{02}$ may not be available or may be cumbersome to derive. Section \ref{sec:estim} provides automatic estimators of $\alpha_{01}$ and $\alpha_{02}$.
\end{comment}
\section{Three Leading Examples}\label{sec:illustrate}
This section illustrates our general framework through three leading
cases. We begin with dynamic treatment effects, where the parameter of
interest depends on the generated outcome only through the second-step
nuisance, so that no direct first-step correction is required ($a_g\equiv0$ in Remark \ref{rem:linear}). This shows that the nested functional setting of \cite{chernozhukov2023automatic} arises as a special case of our
framework. We then extend the partially linear model to allow for an
ML-predicted dependent variable, where the generated outcome enters the
moment condition directly, and a nonzero first-step correction is
required. Finally, we study the estimation of the persistence of unobserved productivity with an unrestricted production function. This application features two first-step regressions: one generates residualized productivity and enters the identifying moment through a downstream conditional expectation, while the other generates a constructed instrument and enters the identifying moment directly.
\subsection{Dynamic Treatment Effects under Sequential Exogeneity}
\label{subsec:dte}
We start by considering dynamic treatment effects with two periods under
sequential exogeneity. This example is closely related to the nested regression formulation in \cite{chernozhukov2023automatic} and \cite{rutgers2026generated}, and earlier work on dynamic treatment regimes (e.g., \citealp{robins1986new}; \citealp{murphy2003optimal}). We do not claim
any methodological innovation relative to
\cite{chernozhukov2023automatic} in this example. Rather, we revisit
this setting to highlight how \emph{downstream local robustness}
provides a unified and transparent interpretation of the recursive
construction of the locally robust moment.\par
Let $W=(Y,S_1,D_1,S_2,D_2)$, where $D_t$ is the treatment at time
$t=1,2$, and $S_t$ are state variables summarizing the relevant
histories. For a fixed treatment path $d=(d_1,d_2)$, the parameter of
interest for a potential outcome $Y^{(d_1,d_2)}$ is
\[
\theta_0(d)=\mathbb{E}\left[Y^{(d_1,d_2)}\right].
\]
Under sequential exogeneity and overlap, with $\perp$ denoting statistical independence,
\begin{align*}
Y^{(d_1,d_2)} &\perp D_1 \mid S_1,
\\
Y^{(D_1,d_2)} &\perp D_2 \mid S_2,
\\
P(D_t=d_t\mid S_t) &\geq c>0
\quad a.s., \quad t=1,2,
\end{align*}
the parameter admits the nested representation
\begin{align*}
g_0(s_2,d_2)
&=
\mathbb{E}[Y\mid S_2=s_2,D_2=d_2],
\\
h_0(s_1,d_1;d_2)
&=
\mathbb{E}[g_0(S_2,d_2)\mid S_1=s_1,D_1=d_1],
\\
\theta_0(d)
&=
\mathbb{E}[h_0(S_1,d_1;d_2)].
\end{align*}
\paragraph{Three-step representation.}
This problem fits our framework as follows.\par
\textbf{First step.}
Let $X_1=(S_2,D_2)$ and define $g_0\in L_2(X_1)$ by
\[
\mathbb{E}[\delta_1(X_1)\left(Y-g_0(X_1)\right)]
=
0
\quad
\text{for all } \delta_1\in L_2(X_1).
\]
The generated outcome is
\[
\tilde Y
=
\varphi(W,g_0)
=
g_0(S_2,d_2),
\]
that is, the continuation value obtained by evaluating the period-2
regression function at the counterfactual treatment level $d_2$.\par
\textbf{Second step.}
Let $X_2=(S_1,D_1)$ and define
\[
S(g_0;d_2)=\tilde Y(g_0;d_2)
=
g_0(S_2,d_2),
\]
and
\[
h_0(X_2;d_2)
=
\mathbb{E}[S(g_0;d_2)\mid X_2]
=
\mathbb{E}[g_0(S_2,d_2)\mid S_1,D_1].
\]
\textbf{Third step.}
The parameter $\theta_0(d)$ is identified by
\[
\mathbb{E}[m(W,g_0,h_0,\theta_0)]
=
0,
\qquad
m(W,g,h,\theta)
=
h(S_1,d_1;d_2)-\theta.
\]
\paragraph{Local robustness with respect to $h_0$.}
Define $\alpha_{02}$ as the Riesz representer of the linear functional
\[
\delta_2
\mapsto
\mathbb{E}[\delta_2(S_1,d_1)]
\]
on $L_2(S_1,D_1)$. Then
\[
\mathbb{E}[\alpha_{02}(S_1,D_1)\delta_2(S_1,D_1)]
=
\mathbb{E}[\delta_2(S_1,d_1)].
\]
Under discrete treatments and overlap,
\[
\alpha_{02}(S_1,D_1)
=
\frac{\mathbf{1}\{D_1=d_1\}}
{P(D_1=d_1\mid S_1)}.
\]
The locally robust moment for the second step is
\begin{align*}
m^*(W,g,h,\alpha_{02},\theta)
&=
h(S_1,d_1;d_2)-\theta
\\
&\quad
+
\alpha_{02}(S_1,D_1)
\big(g(S_2,d_2)-h(S_1,D_1;d_2)\big).
\end{align*}
\paragraph{Downstream local robustness and the first-step correction.}
The parameter depends on $g_0$ only through $h_0$, that is,
\[
g_0 \longrightarrow h_0 \longrightarrow \theta_0.
\]
Hence, once the moment is locally robust to $h_0$, the first-order
indirect effect of $g_0$ through $h_0$ vanishes. This is the key idea
of \emph{downstream local robustness}. As a result, to construct the
first-step correction it suffices to differentiate $m^*$ with respect to
$g$, holding $h_0$ fixed.\par
The only remaining dependence of $m^*$ on $g$ is through
\[
\alpha_{02}(S_1,D_1)
\big(g(S_2,d_2)-h_0(S_1,D_1;d_2)\big).
\]
Thus $\alpha_{01}$ is the Riesz representer of
\[
\delta_1
\mapsto
\mathbb{E}\big[
\alpha_{02}(S_1,D_1)\delta_1(S_2,d_2)
\big]
\]
on $L_2(S_2,D_2)$, that is,
\[
\mathbb{E}[
\alpha_{01}(S_2,D_2)\delta_1(S_2,D_2)
]
=
\mathbb{E}\big[
\alpha_{02}(S_1,D_1)\delta_1(S_2,d_2)
\big].
\]
Under overlap,
\[
\alpha_{01}(S_2,D_2)
=
\frac{\mathbf{1}\{D_2=d_2\}}
{P(D_2=d_2\mid S_2)}
\,
\mathbb{E}[
\alpha_{02}(S_1,D_1)\mid S_2].
\]
If $S_2$ contains $(S_1,D_1)$, then
\[
\alpha_{01}(S_2,D_2)
=
\frac{\mathbf{1}\{D_1=d_1\}}
{P(D_1=d_1\mid S_1)}
\frac{\mathbf{1}\{D_2=d_2\}}
{P(D_2=d_2\mid S_2)}.
\]
\paragraph{Locally robust moment.}
The resulting locally robust moment is
\begin{align*}
\psi(W,g_0,h_0,\alpha_0,\theta)
&=
h_0(S_1,d_1;d_2)-\theta
\\
&\quad
+
\alpha_{02}(S_1,D_1)
\big(
g_0(S_2,d_2)-h_0(S_1,D_1;d_2)
\big)
\\
&\quad
+
\alpha_{01}(S_2,D_2)
\big(
Y-g_0(S_2,D_2)
\big),
\end{align*}
where $\alpha_0=(\alpha_{01},\alpha_{02})$.
\paragraph{Discussion.}
For $T=2$, this moment coincides with the recursive debiased moment in
\cite{chernozhukov2023automatic}. The equivalence is not accidental: it
follows from the fact that the parameter depends on $g_0$ only through
$h_0$. The contribution of our framework in this example is therefore
not a new estimator, but rather the interpretation of the recursive
construction through downstream local robustness. This perspective
extends naturally beyond nested functionals to settings where the moment
depends directly on the first-step nuisance, as in the partially linear
model with an ML-generated outcome.
\subsection{The Partially Linear Model}\label{subsec:plm}
In this context, we observe $W = (Y, X)$, where $X$ includes a vector
$X_1$ used to predict $Y$ in a first step, a vector of controls $X_2$
onto which the generated outcome is projected in a second step, and a
low-dimensional vector of regressors of interest $X_3$ whose effect on
the ML prediction is the parameter of interest in a third step. Although
$X_2$ and $X_3$ are disjoint, we allow $X_1$ to overlap with either or
contain additional covariates in $X$. For instance, $X_3$ may be a
policy variable or a socioeconomic characteristic of interest, $X_2$
demographic and geographic controls, and $X_1$ may collect
satellite-derived predictors of an economic outcome unavailable at
scale, potentially augmented by $X_2$.\par
Let $\Delta_1\subseteq L_2(X_1)$ be a specified closed linear
prediction space, and let $\Pi_{\Delta_1}$ denote the
$L_2(F_0)$ orthogonal projection onto that space. The generated
outcome is
\[
\tilde Y=g_0(X_1),
\qquad
g_0=\Pi_{\Delta_1}Y.
\]
Thus, $g_0$ is the best mean-square predictor within $\Delta_1$,
which need not equal $\mathbb E[Y\mid X_1]$. We estimate $g_0$
using machine learning and use its prediction as the dependent
variable in the model
\begin{align}
\label{eq:plm}
\tilde{Y} = X_3'\theta_0 + r_0(X_2) + U,
\qquad
\mathbb{E}[U|X_2,X_3] = 0,
\end{align}
which identifies $\theta_0$ as the partially linear projection
coefficient of the ML prediction on $X_3$ controlling for $X_2$. \par
Following \cite{robinson1988root}, the identifying moment condition for
$\theta_0$ is
\begin{align}
\label{eq:robinsmoment}
\mathbb{E}\left[
\left(
\tilde{Y}
-
h_{0Y}(X_2)
-
\theta_0'(X_3-h_{0X}(X_2))
\right)
(X_3-h_{0X}(X_2))
\right]
=
0,
\end{align}
where
\[
h_{0Y}(X_2)=\mathbb{E}[\tilde{Y}|X_2],
\qquad
h_{0X}(X_2)=\mathbb{E}[X_3|X_2].
\]
However, since $\tilde{Y}$ is itself estimated by machine learning,
this moment condition does not account for the regularization bias
introduced in the first step. To construct the locally robust moment,
we start by mapping the problem into our three-step framework.
\paragraph{Three-step representation.}
This problem fits our framework as follows.\par
\bigskip \noindent
\textbf{First step.}
The first-step function $g_0\in\Delta_1$ is identified by
\[
\mathbb E\left[
\delta_1(X_1)\{Y-g_0(X_1)\}
\right]=0
\quad
\text{for all }\delta_1\in\Delta_1.
\]
The generated outcome is
\[
\tilde Y=\varphi(W,g_0)=g_0(X_1).
\]
\bigskip \noindent
\textbf{Second step.}
Let $X_2$ denote the controls onto which the generated outcome and the
regressors of interest are projected, and define
\[
S(g_0)\equiv(\tilde{Y}(g_0),X_3)',
\qquad
h_0(X_2)
=
\mathbb{E}[S(g_0)|X_2]
=
(h_{0Y}(X_2),h_{0X}(X_2))'.
\]
\bigskip \noindent
\textbf{Third step.}
The parameter $\theta_0$ is identified by the moment condition
\[
\mathbb{E}[m(W,g_0,h_0,\theta_0)]
=
0,
\]
where
\[
m(W,g,h,\theta)
=
\big(
g(X_1)-h_Y(X_2)-\theta'(X_3-h_X(X_2))
\big)
(X_3-h_X(X_2)),
\]
with $h=(h_Y,h_X)'$.
The following proposition derives the locally robust moment for
$\theta_0$.
\begin{proposition}
\label{prop:plm}
Let $g_0=\Pi_{\Delta_1}Y$, and define
\[
h_{0Y}(X_2)=\mathbb E[g_0(X_1)\mid X_2],
\qquad
h_{0X}(X_2)=\mathbb E[X_3\mid X_2],
\qquad
\nu=X_3-h_{0X}(X_2).
\]
Suppose $B=\mathbb E[\nu\nu']$ is nonsingular and let
\[
\theta_0
=
B^{-1}\mathbb E[
\{g_0(X_1)-h_{0Y}(X_2)\}\nu
].
\]
The first- and second-step Riesz representers are
\[
\alpha_{01}=\Pi_{\Delta_1}\nu,
\qquad
\alpha_{02}=0,
\]
where the projection is applied componentwise. A locally robust
moment for $\theta_0$ is
\begin{align}
\psi(W,g_0,h_0,\alpha_0,\theta)
&=
\{g_0(X_1)-h_{0Y}(X_2)-\theta'\nu\}\nu
\nonumber\\
&\quad
+\alpha_{01}(X_1)\{Y-g_0(X_1)\}.
\label{eq:lr_plm_general_main}
\end{align}
If $\Delta_1=L_2(X_1)$, then
\[
g_0(X_1)=\mathbb E[Y\mid X_1],
\qquad
\alpha_{01}(X_1)=\mathbb E[\nu\mid X_1].
\]
If, additionally, $X_2\subseteq X_1$ and $X_3\subseteq X_1$,
then $\alpha_{01}=\nu$ and
$h_{0Y}(X_2)=\mathbb E[Y\mid X_2]$, so
\begin{align}
\psi(W,h_0,\theta)
=
\{Y-h_{0Y}(X_2)-\theta'\nu\}\nu.
\label{eq:lr_plm_main}
\end{align}
\end{proposition}
\begin{remark}[Partially linear projection]
\label{rem:plm_projection}
The conditional mean restriction in \eqref{eq:plm} is stronger than
needed for the moment construction. Let
\[
\nu\equiv X_3-h_{0X}(X_2),\qquad
B\equiv\mathbb{E}[\nu\nu'],
\]
and define the partially linear projection coefficient
\[
\theta_0
=
B^{-1}\mathbb{E}\big[
(\tilde Y-h_{0Y}(X_2))\nu
\big].
\]
Setting
\[
U
=
\tilde Y-h_{0Y}(X_2)-\theta_0'\nu
\]
gives
\[
\mathbb{E}[U\mid X_2]=0,
\qquad
\mathbb{E}[U\nu]=0.
\]
These conditions are sufficient for the orthogonality calculations
in Proposition~\ref{prop:plm}. Hence the locally robust moments
\eqref{eq:lr_plm_general_main}--\eqref{eq:lr_plm_main}, and the
asymptotic results below, remain valid when $\theta_0$ is interpreted
as this partially linear projection coefficient.
\end{remark}
\begin{remark}[Algorithmic credit scoring and conditional group disparities]
\label{rem:credit_scoring}
A concrete application of the general PLM arises in algorithmic credit
scoring. Let $Y\in\{0,1\}$ indicate future default and let
\[
g_0(X_1)
=
\mathbb{P}(Y=1\mid X_1)
\]
denote the population default-risk score associated with the information
$X_1$ used by a machine-learning credit model. Let $A\in\{0,1\}$ denote
a protected group characteristic, such as race or gender, which is
excluded from $X_1$, and let $X_2$ collect characteristics on which the
researcher wants borrowers to be compared. The parameter of interest is
the partially linear projection coefficient of the ML risk score on
group membership controlling for $X_2$, as defined in
Remark~\ref{rem:plm_projection}. It measures conditional group
disparities in predicted default risk and need not have a causal
interpretation.
This setting is empirically relevant because credit-scoring algorithms
are commonly constructed without directly using protected
characteristics, while the remaining borrower information may still
contain substantial information about group membership
\citep{fuster2022predictably,hurlin2026fairness}.
Consider the empirically natural case
\[
X_2\subseteq X_1,
\qquad
A\notin X_1.
\]
Setting $X_3=A$, define
\[
h_{0A}(X_2)
=
\mathbb{E}[A\mid X_2],
\qquad
\nu
=
A-h_{0A}(X_2).
\]
Then Proposition~\ref{prop:plm} gives
\[
\alpha_{01}(X_1)
=
\mathbb{E}[\nu\mid X_1]
=
\mathbb{E}[A\mid X_1]
-
\mathbb{E}[A\mid X_2].
\]
Hence the locally robust moment is
\begin{align*}
\psi
&=
\left(
g_0(X_1)
-h_{0Y}(X_2)
-\theta\nu
\right)\nu
\\
&\quad
+
\left\{
\mathbb{E}[A\mid X_1]
-
\mathbb{E}[A\mid X_2]
\right\}
\left\{
Y-g_0(X_1)
\right\},
\end{align*}
where
\[
h_{0Y}(X_2)
=
\mathbb{E}[g_0(X_1)\mid X_2].
\]
Here we consider the conditional-mean case
$\Delta_1=L_2(X_1)$ of Proposition~\ref{prop:plm}.
At the population truth, $X_2\subseteq X_1$ implies
\[
h_{0Y}(X_2)=\mathbb E[Y\mid X_2].
\]
Nevertheless, we retain the sequential implementation: estimate $g_0$
first and then estimate $h_{0Y}$ by regressing out-of-sample first-step
predictions on $X_2$. The motivation is to keep the second-step
regression aligned with the prediction used in the identifying moment,
rather than to rely on an identity that holds at the true first-step
function.
To make this distinction precise, for a candidate prediction function
$g$ define
\[
h_Y(g)(X_2)=\mathbb E[g(X_1)\mid X_2].
\]
Then $h_{0Y}=h_Y(g_0)$, but for a possibly misspecified first-step
limit $g^\dagger$,
\[
h_Y(g^\dagger)(X_2)
\neq
\mathbb E[Y\mid X_2]
\]
in general. Regressing generated predictions on $X_2$ targets
$h_Y(g^\dagger)$, provided the second-step conditional regression is
consistently estimated, and preserves the conditional centering
\[
\mathbb E[
g^\dagger(X_1)-h_Y(g^\dagger)(X_2)
\mid X_2
]=0.
\]
This centering maintains orthogonality of the Robinson moment with
respect to the group regression $h_{0A}$, even at the misspecified
first-step limit. Regressing $Y$ directly on $X_2$ need not preserve
this property.
This alignment does not, by itself, establish valid inference under
persistent first-step misspecification. Such inference also requires
controlling the effect of estimating the Riesz representer, including
the drift term
\[
\mathbb E\left[
\{\hat\alpha_{01}(X_1)-\alpha_{01}(X_1)\}
\{g_0(X_1)-\hat g(X_1)\}
\right],
\]
where fitted functions are held fixed when taking expectations.
Thus, the sequential construction has a misspecification-robustness
motivation, while root-$n$ inference remains subject to the relevant
nuisance convergence and remainder conditions. The first-step Riesz representer therefore has a direct interpretation:
it measures the incremental ability of the rich information used by the
credit-scoring algorithm to predict group membership beyond the
comparison variables $X_2$. In particular, if
\[
\mathbb{E}[A\mid X_1]
=
\mathbb{E}[A\mid X_2],
\]
then $\alpha_{01}=0$ and estimation of the credit-risk score has no
first-order effect on inference for the conditional group disparity.
Otherwise, treating the estimated score as an observed dependent
variable ignores a first-order estimation effect.
Sample splitting alone does not generally remove this problem. Let
$\hat g$ be estimated on an independent training sample and let $n$
denote the size of the evaluation sample. Conditional on the training
sample, the naive downstream moment contains the first-order term
\[
\mathbb{E}\left[
\alpha_{01}(X_1)
\{\hat g(X_1)-g_0(X_1)\}
\,\middle|\,
\text{training sample}
\right].
\]
For the contribution from estimating $g_0$ to be negligible at the
root-$n$ scale, this term must be $o_p(n^{-1/2})$. For example, if the
training sample has size $N$ and
\[
\|\hat g-g_0\|_2=O_p(N^{-a}),
\]
Cauchy--Schwarz gives the sufficient condition
\[
\sqrt n\,N^{-a}\to0.
\]
At the parametric rate $a=1/2$, this sufficient rate-separation
condition becomes $n/N\to0$. It need not be necessary when the
prediction error has additional orthogonality to $\alpha_{01}$.
Sample splitting makes predictions out of sample but does not
automatically make this first-step contribution negligible. The
locally robust correction removes its linear contribution when the
Riesz representer is known; with an estimated representer, the
remaining drift must satisfy the appropriate remainder conditions.
\end{remark}
The cancellation in \eqref{eq:lr_plm_main} is a special case.
For a restricted first-step prediction space,
\[
\mathbb E[Y\mid X_2]-h_{0Y}(X_2)
=
\mathbb E[Y-g_0(X_1)\mid X_2],
\]
which need not vanish even when $X_2\subseteq X_1$. Thus,
$h_{0Y}$ generally must be estimated from the generated outcome
rather than by regressing $Y$ directly on $X_2$. The identity
$h_{0Y}(X_2)=\mathbb E[Y\mid X_2]$ holds under the additional
condition $\mathbb E[Y-g_0(X_1)\mid X_2]=0$, for example when
$L_2(X_2)\subseteq\Delta_1$. We do not impose this condition in
the general PLM formulation.
Although the second-step correction is zero, the second-step
nuisance remains generated: its estimation depends on the
first-step prediction. Orthogonality removes its first-order
effect on the expected moment, not its role in constructing
the estimator.
\begin{remark}[Outcomes residualized of confounding variation]
\label{rem:resid}
This PLM setting extends directly to a common practice in empirical work:
residualizing an observed outcome of confounding variation before using
it in subsequent analysis. In the measurement of inequality of
opportunity, income is purged of life-cycle or cohort variation before
its remaining dispersion is attributed to circumstances beyond
individual control \citep{brunori2025inequality}; in partial insurance
models, consumption is purged of predictable components before
household smoothing is studied \citep{blundell2008consumption}. The
object entering the subsequent analysis is then a generated outcome,
corresponding to the choice
$\tilde{Y}\equiv\varphi(W,g_0)=Y-g_0(X_1)$.
Robinson's moment remains orthogonal with respect to $h_0$, so
$\alpha_{02}=0$. Perturbing $g$ enters with the opposite sign
to the predicted-outcome case, giving
\[
\alpha_{01}
=
-\Pi_{\Delta_1}\nu,
\qquad
\nu=X_3-h_{0X}(X_2),
\]
where the projection is applied componentwise. The locally robust
moment is therefore
\[
\psi
=
\left(
\tilde Y-h_{0Y}(X_2)-\theta'\nu
\right)\nu
+
\alpha_{01}(X_1)
\left(
Y-g_0(X_1)
\right),
\]
where
\[
h_{0Y}(X_2)=\mathbb E[\tilde Y\mid X_2],
\qquad
\tilde Y=Y-g_0(X_1).
\]
Thus, the first-step identifying error coincides with the generated
outcome itself. In the conditional-mean case
$\Delta_1=L_2(X_1)$, the representer simplifies to
\[
\alpha_{01}(X_1)
=
-\mathbb E[\nu\mid X_1].
\]
The same construction applies when $\theta_0$ is interpreted as
the partially linear projection coefficient in
Remark~\ref{rem:plm_projection}.
\end{remark}
\begin{comment}
\begin{remark}[The surrogate index as an illustration]\label{rem:surrogate}
\cite{athey2025surrogate} study the estimation of long-term average treatment effects by combining two samples: an experimental sample, indexed by $P_i = E$, in which treatment $W_i \in \{0,1\}$ and short-term surrogate outcomes $S_i$ are observed but the long-term outcome $Y_i$ is missing; and an observational sample, indexed by $P_i = O$, in which $Y_i$ and $S_i$ are observed but $W_i$ is not. Pre-treatment covariates $X_i$ are observed in both samples. Under a surrogacy assumption, requiring that the causal path from treatment to the long-term outcome runs entirely through the surrogates, and a comparability assumption across samples, \citet{athey2025surrogate} establish that the average treatment effect $\tau \equiv \mathbb{E}[Y_i(1) - Y_i(0) \mid P_i = E]$ is nonparametrically identified, and that four distinct moment representations of $\tau$ in terms of
observables are available. We use these four representations to illustrate the scope of our framework.\par
Define the surrogate index $\mu(s,x,O) \equiv \mathbb{E}[Y_i \mid S_i{=}s,X_i{=}x, P_i{=}O]$, the propensity score $\rho(x) \equiv \Pr(W_i{=}1 \mid X_i{=}x, P_i{=}E)$, the surrogate score $\rho(s,x) \equiv \Pr(W_i{=}1 \mid S_i{=}s, X_i{=}x, P_i{=}E)$, and the sampling score $\pi(s,x) \equiv \Pr(P_i{=}E \mid S_i{=}s, X_i{=}x)$, with $\pi \equiv \Pr(P_i{=}E)$. Also define the conditional mean of the surrogate index in the experimental sample $\mu(w,x) \equiv \mathbb{E}[\mu(S_i,X_i,O) \mid W_i{=}w, X_i{=}x, P_i{=}E]$.
\medskip
\noindent\textit{First representation: nested functional.}
The statistical estimand
\begin{align*} \tau^{*} \;\equiv\; &\mathbb{E}\Bigl[\Bigl\{ \mathbb{E}\bigl[\mathbb{E}[Y_i \mid S_i, X_i, P_i{=}O] \mid W_i{=}1, X_i, P_i{=}E\bigr]\\-&\mathbb{E}\bigl[\mathbb{E}[Y_i \mid S_i, X_i, P_i{=}O] \mid W_i{=}0, X_i, P_i{=}E\bigr] \Bigr\} \Bigm| P_i{=}E\Bigr] \end{align*} involves a nested functional: the moment depends on the generated outcome $\tilde{Y} \equiv \mu(S_i,X_i,O)$ only through a subsequent conditional expectation, not directly. Downstream local robustness to $\mu(w,x)$ automatically annihilates the indirect first-step effect of $g_0 = \mu(\cdot,\cdot,O)$, so no direct first-step correction is required. This case falls under the nested-functional setting of \cite{chernozhukov2023automatic}, which arises here as a special case of our
framework.\par
\medskip
\noindent\textit{Second representation: generated outcome enters directly.}
The surrogate-index inverse probability weighting estimand,
\begin{equation*} \tau^{E} \;\equiv\; \mathbb{E}\!\left[ \mu(S_i,X_i,O)\cdot\frac{W_i}{\rho(X_i)}- \mu(S_i,X_i,O)\cdot\frac{1-W_i}{1-\rho(X_i)}\;\bigg|\; P_i{=}E\right], \end{equation*}
has a different structure: the surrogate index $\mu(S_i,X_i,O)$ appears in the identifying moment \emph{directly} as the generated dependent variable, with no further conditioning. This is precisely the setting of the present section, with $\mu(S_i,X_i,O)$ in the role of $\tilde{Y} \equiv g_0(X_1)$. In this representation, the surrogate index enters the identifying moment directly rather than only through a downstream conditional expectation. Theorem~\ref{thm:lr} therefore provides a direct locally robust construction for this moment representation, correcting for the regularization bias introduced by ML estimation of the surrogate index.
\medskip
\noindent\textit{Third representation: no generated outcome.}
The surrogate-score estimand,
\begin{equation*} \tau^{O} \;\equiv\; \mathbb{E}\!\left[ Y_i \cdot \frac{\rho(S_i,X_i)\,\pi(S_i,X_i)\,(1-\pi)}{\rho(X_i)\,(1-\pi(S_i,X_i))\,\pi}- Y_i \cdot\frac{(1-\rho(S_i,X_i))\,\pi(S_i,X_i)\,(1-\pi)}{(1-\rho(X_i))\,(1-\pi(S_i,X_i))\,\pi}\;\bigg|\; P_i{=}O\right],
\end{equation*} involves only the observed outcome $Y_i$ in the observational sample, reweighted by the surrogate and sampling scores. No generated outcome appears, so our framework is not needed: the standard two-step locally robust construction of \cite{chernozhukov2022locally} applies directly.
\medskip
\noindent\textit{Fourth representation: the efficient doubly-robust moment.}
The estimand $$\tau^{O,E} = \mathbb{E}[\psi(Y_i,S_i,W_i,X_i,P_i) + \tau],$$ with
influence function
\begin{align*} \psi(y,s,w,x,p) \;=\;& \frac{\mathbf{1}_{p=E}}{\pi}\left( \frac{w\cdot(\mu(s,x,O)-\mu(1,x))}{\rho(x)}- \frac{(1-w)\cdot(\mu(s,x,O)-\mu(0,x))}{1-\rho(x)} \right) \\ +\;& \frac{\mathbf{1}_{p=E}}{\pi} \bigl(\mu(1,x) - \mu(0,x) - \tau\bigr) \\+\;& \frac{\mathbf{1}_{p=O}}{\pi}\cdot\frac{\pi(s,x)}{1-\pi(s,x)} \cdot\frac{(y-\mu(s,x,O))(\rho(s,x)-\rho(x))}{\rho(x)(1-\rho(x))}, \end{align*}
combines both samples and achieves the semiparametric efficiency bound. This is precisely the kind of moment our framework is designed to construct: a locally robust moment obtained by correcting an identifying moment that originally depended on the generated outcome $\tilde{Y} \equiv \mu(S_i,X_i,O)$.
\end{remark}
\end{comment}
\subsection{Productivity Persistence with Nonparametric Production Functions}
Our framework accommodates structural environments in which the
production function is left unrestricted. We consider identification
and inference on the persistence of latent productivity. Rather than
requiring an externally specified instrument, we construct valid
instruments from observable period-1 productivity signals using
orthogonal projections. This construction separates instrument validity
from relevance and fits directly into our general locally robust
framework.
Productivity is a key driver of firm performance, yet its latent nature
requires inferring its dynamics from data on inputs and outputs. We
observe
\[
W=\big((Y_t,K_t,L_t,I_t)_{t=1,2}\big),
\]
where $Y_t$ is output, $(K_t,L_t)$ are capital and labor inputs, and
$I_t$ is investment, all in logs. Our parameter of interest is the
persistence of latent productivity, denoted by $\theta_0$.
\paragraph{Structural model and timing.}
Assume output satisfies
\begin{align}
Y_t
=
\tilde g_{0t}(K_t,L_t)+\omega_t+U_t,
\qquad t=1,2,
\label{eq:prod_struct}
\end{align}
where $\tilde g_{0t}$ is an unrestricted, possibly period-specific
production function and $\omega_t$ is latent productivity. Following
the proxy-variable approach
\citep{olley1996dynamics,levinsohn2003estimating,ackerberg2015identification},
assume that investment is strictly monotone in productivity conditional
on capital and labor, so that
\begin{align}
\omega_1
=
\tilde h_0(K_1,L_1,I_1)
\label{eq:prod_proxy}
\end{align}
for some unknown function $\tilde h_0$.
Productivity evolves according to
\begin{align}
\omega_2
=
\theta_0\omega_1+\xi_2.
\label{eq:prod_ar1_struct}
\end{align}
Let $\mathcal F_1$ denote the firm's information set after the
period-2 production inputs $(K_2,L_2)$ have been chosen but before the
period-2 productivity innovation $\xi_2$ and transitory output shock
$U_2$ are realized. Assume that
\[
(K_1,L_1,I_1,K_2,L_2)\subseteq\mathcal F_1
\]
and
\begin{align}
\mathbb E[U_1\mid K_1,L_1,I_1]&=0,
\label{eq:prod_timing1}
\\
\mathbb E[\xi_2+U_2\mid\mathcal F_1]&=0.
\label{eq:prod_timing2}
\end{align}
Thus, $I_1$ may respond to period-1 productivity, whereas the
period-2 inputs used in the instrument construction are predetermined
with respect to the innovation $\xi_2$ and the transitory shock $U_2$.
This timing structure is appropriate when these inputs are committed
before the period-2 innovation is observed; it is stronger than
production-function models in which flexible inputs respond
contemporaneously to productivity.
\paragraph{Reduced-form productivity.}
To map the problem into our general framework, define
\[
X_{11}\equiv(K_1,L_1),
\qquad
X_{12}\equiv(K_2,L_2),
\qquad
X_{21}\equiv(K_1,L_1,I_1).
\]
The first reduced-form regression projects output on the period-1 inputs
\[
g_{01}(X_{11})
\equiv
\mathbb E[Y_1\mid X_{11}],
\]
which by \eqref{eq:prod_struct} and \eqref{eq:prod_timing1} decomposes as
\[
g_{01}(X_{11})
=
\tilde g_{01}(K_1,L_1)
+
\mathbb E[\omega_1\mid K_1,L_1].
\]
Define the residualized generated outcome
\[
\tilde Y_1
\equiv
Y_1-g_{01}(X_{11})
\]
and its conditional mean
\begin{align}
h_0(X_{21})
&\equiv
\mathbb E[\tilde Y_1\mid X_{21}]
\nonumber\\
&=
\omega_1-\mathbb E[\omega_1\mid K_1,L_1].
\label{eq:prod_h0}
\end{align}
Hence $h_0$ recovers the component of period-1 productivity that is
not predictable from period-1 production inputs.
\paragraph{Constructed instruments.}
The remaining challenge is to construct an instrument that is
correlated with $h_0$ but orthogonal to the unrestricted production
function and to the predictable component of period-1 productivity.
Let
\[
\mathcal H
\equiv
\overline{
\left\{
a(X_{11})+b(X_{12}):
a\in L_2(X_{11}),\
b\in L_2(X_{12})
\right\}
}^{\,L_2}
\]
denote the closed additive space generated by period-1 and period-2
production inputs.
Let
\[
R\equiv r(X_{21})\in L_2
\]
be an observable period-1 productivity signal, where the transformation
$r$ is chosen by the researcher. Define
\begin{align}
q_0(X_{11},X_{12})
&\equiv
\Pi_{\mathcal H}R,
\label{eq:prod_q}
\\
Z_0
&\equiv
R-q_0(X_{11},X_{12}).
\label{eq:prod_Z}
\end{align}
Thus, $q_0$ is the $L_2$ projection of $R$ onto $\mathcal H$ and
$Z_0$ is its orthogonal residual. By construction,
\begin{align}
\mathbb E[Z_0 a(X_{11})]&=0
\quad
\text{for all }a\in L_2(X_{11}),
\label{eq:prod_orth11}
\\
\mathbb E[Z_0 b(X_{12})]&=0
\quad
\text{for all }b\in L_2(X_{12}).
\label{eq:prod_orth12}
\end{align}
Using \eqref{eq:prod_h0} and the productivity law of motion,
\begin{align}
Y_2-\theta_0 h_0(X_{21})
&=
\tilde g_{02}(K_2,L_2)
+
\theta_0
\mathbb E[\omega_1\mid K_1,L_1]
+
\xi_2+U_2.
\label{eq:prod_reduced}
\end{align}
The first two terms on the right-hand side belong to $\mathcal H$.
The constructed instrument $Z_0$ therefore eliminates both without
requiring either $\tilde g_{02}$ or
$\mathbb E[\omega_1\mid K_1,L_1]$ to be separately estimated.
\begin{proposition}[Constructed-instrument identification]
\label{prop:app}
Suppose equations
\eqref{eq:prod_struct}--\eqref{eq:prod_timing2} hold and let
$Z_0$ be defined by
\eqref{eq:prod_q}--\eqref{eq:prod_Z}. Then
\begin{align}
\mathbb E\left[
Z_0
\left\{
Y_2-\theta_0 h_0(X_{21})
\right\}
\right]
=0.
\label{eq:prod_moment}
\end{align}
If
\[
\mathbb E[Z_0h_0(X_{21})]\neq0,
\]
then $\theta_0$ is identified as
\[
\theta_0
=
\frac{\mathbb E[Z_0Y_2]}
{\mathbb E[Z_0h_0(X_{21})]}.
\]
Moreover, consider the class of constructed instruments
\[
\mathcal Z
\equiv
\left\{
r(X_{21})-\Pi_{\mathcal H}r(X_{21}):
r\in L_2(X_{21})
\right\}.
\]
There exists a relevant instrument in $\mathcal Z$ if and only if
\begin{align}
h_0(X_{21})\notin\mathcal H.
\label{eq:prod_relevance}
\end{align}
In particular, when \eqref{eq:prod_relevance} holds,
\[
Z_0^*
\equiv
h_0(X_{21})-\Pi_{\mathcal H}h_0(X_{21})
\]
satisfies
\[
\mathbb E[Z_0^*h_0(X_{21})]
=
\mathbb E[(Z_0^*)^2]
>0.
\]
\end{proposition}
Condition~\eqref{eq:prod_relevance} separates instrument validity from
relevance. Validity follows from the timing restrictions and the
orthogonalization in \eqref{eq:prod_Z}. Relevance requires the
investment proxy to contain variation in residualized period-1
productivity that cannot be reconstructed by adding an arbitrary
function of period-1 production inputs to an arbitrary function of
period-2 production inputs.
Importantly, proxy invertibility and instrument relevance are distinct
requirements. Strict monotonicity of investment implies that period-1
productivity can be recovered from $X_{21}$, but it does not imply
condition~\eqref{eq:prod_relevance}. If the residualized productivity
component $h_0(X_{21})$ can be reconstructed from an additive function
of $X_{11}$ and $X_{12}$, then every instrument obtained by
orthogonalizing a function of $X_{21}$ against $\mathcal H$ is
irrelevant. Hence identification requires the period-1 proxy to contain
productivity information that is not exhausted by current and future
production inputs.
The choice of the signal $R$ is important for relevance. Raw investment
need not itself generate a useful instrument. For example, under a
deterministic capital-accumulation equation, investment may belong to
the additive space $\mathcal H$, in which case its orthogonal residual
is identically zero. We therefore consider a dictionary
$\{r_j(X_{21})\}_{j=1}^J$ of nonlinear transformations and interactions
of $(K_1,L_1,I_1)$ and construct
\[
Z_{0j}
=
r_j(X_{21})
-
\Pi_{\mathcal H}r_j(X_{21}),
\qquad j=1,\ldots,J.
\]
For a fixed finite dictionary, the resulting vector
$Z_0=(Z_{01},\ldots,Z_{0J})'$ can be used directly in the
overidentified GMM framework developed above. For scalar $\theta_0$,
relevance requires
\[
\mathbb E[Z_0h_0(X_{21})]\neq0.
\]
At the population level, if the linear span of a sequence of candidate
signals is dense in $L_2(X_{21})$, condition~\eqref{eq:prod_relevance}
implies that some sufficiently rich finite subset contains relevant
instruments. In implementation, we use a fixed finite set of instruments chosen
ex ante from the candidate dictionary. Our asymptotic theory maintains strong identification. For a fixed
finite vector of constructed instruments, this corresponds here to
the nonsingularity condition in Assumption~\ref{ass:normality}(ii);
weak-instrument asymptotics are beyond the scope of the paper.
The projection on $\mathcal H$ must preserve its additive structure.
In particular, $q_0$ is the $L_2$ projection of $R$ onto
$L_2(X_{11})+L_2(X_{12})$, rather than the unrestricted conditional
expectation $\mathbb E[R\mid X_{11},X_{12}]$. This distinction is
important for relevance. Orthogonality to the additive space is
sufficient to remove the unrestricted period-2 production function and
the predictable component of period-1 productivity, while projecting
onto the larger space $L_2(X_{11},X_{12})$ may unnecessarily remove
additional variation in the period-1 proxy that is informative about
$h_0$.
\paragraph{Three-step representation.}
The construction fits our general framework with two first-step
regressions.
\bigskip \noindent
\textbf{First step.}
The reduced-form regression $g_{01}\in L_2(X_{11})$ is identified by
\[
\mathbb E\left[
\delta_{11}(X_{11})
\{Y_1-g_{01}(X_{11})\}
\right]
=0
\quad
\text{for all }\delta_{11}\in L_2(X_{11}).
\]
The instrument projection $q_0\in\mathcal H$ is identified by
\[
\mathbb E\left[
\delta_q(X_{11},X_{12})
\{R-q_0(X_{11},X_{12})\}
\right]
=0
\quad
\text{for all }\delta_q\in\mathcal H.
\]
Thus, $q_0$ is a first-step regression of the form considered in
Section~\ref{sec:set}, with the closed linear regression space
$\Delta_{1q}=\mathcal H$ rather than the unrestricted space
$L_2(X_{11},X_{12})$.
Consequently, using the product-space convention introduced in
Section~\ref{sec:set}, the first-step nuisance is
\[
g_0=(g_{01},q_0).
\]
The first component generates the residualized outcome
\[
\tilde Y_1(g_{01})
=
Y_1-g_{01}(X_{11}),
\]
while the second generates the instrument
\[
Z(q)=R-q(X_{11},X_{12}).
\]
\bigskip \noindent
\textbf{Second step.}
The second-step nuisance is
\[
h_0(X_{21})
=
\mathbb E[
\tilde Y_1(g_{01})\mid X_{21}
],
\]
or equivalently
\[
\mathbb E\left[
\delta_2(X_{21})
\left\{
\tilde Y_1(g_{01})-h_0(X_{21})
\right\}
\right]
=0
\quad
\text{for all }\delta_2\in L_2(X_{21}).
\]
\bigskip \noindent
\textbf{Third step.}
The identifying moment is
\[
m(W,g,h,\theta)
=
Z(q)
\left\{
Y_2-\theta h(X_{21})
\right\},
\]
so that at the truth Proposition~\ref{prop:app} gives
\[
\mathbb E[
m(W,g_0,h_0,\theta_0)
]
=0.
\]
This representation illustrates two distinct forms of
orthogonalization. The first is economic: $Z_0$ is constructed to be
orthogonal to arbitrary functions of period-1 and period-2 production
inputs. The second is statistical: the identifying moment must be made
locally robust to ML estimation of $g_{01}$, $q_0$, and $h_0$.
To derive the locally robust moment, define
\[
a_0(X_{21})
\equiv
\mathbb E[Z_0\mid X_{21}]
\]
and let
\[
\pi_0(X_{11},X_{12};\theta)
\equiv
\Pi_{\mathcal H}
\left[
Y_2-\theta h_0(X_{21})
\right].
\]
Both objects are well-defined functions of observables and nuisance
parameters. The function $a_0$ is the conditional mean of the
constructed instrument, while $\pi_0(\cdot;\theta)$ is the additive
projection of the third-step residual on $\mathcal H$.
\begin{proposition}\label{prop:prod_lr}
Consider the identifying moment in \eqref{eq:prod_moment}. A locally
robust moment for $\theta_0$ is
\begin{align}
\psi(W,g_0,h_0,\alpha_0,\theta)
&=
Z_0
\left\{
Y_2-\theta h_0(X_{21})
\right\}
\notag\\
&\quad
-\theta a_0(X_{21})
\left\{
Y_1-g_{01}(X_{11})-h_0(X_{21})
\right\}
\notag\\
&\quad
-\pi_0(X_{11},X_{12};\theta)
\left\{
R-q_0(X_{11},X_{12})
\right\}.
\label{eq:prod_lr}
\end{align}
The first-step Riesz representer is
$\alpha_{01}=(\alpha_{011},\alpha_{01q})$, where
\begin{align*}
\alpha_{011}(X_{11};\theta)&=0,
\\
\alpha_{01q}(X_{11},X_{12};\theta)
&=
-\pi_0(X_{11},X_{12};\theta),
\end{align*}
and the second-step Riesz representer is
\[
\alpha_{02}(X_{21};\theta)
=
-\theta a_0(X_{21}).
\]
\end{proposition}
The zero representer $\alpha_{011}$ has an economic interpretation.
Estimation of $g_{01}$ affects the downstream nuisance $h_0$, but the
resulting first-order effect on the identifying moment is proportional
to
\[
\mathbb E[Z_0\delta_{11}(X_{11})],
\]
which is zero by construction of the instrument. Thus, orthogonalizing
the instrument against $L_2(X_{11})$ also eliminates the first-order
effect of estimating the regression used to construct residualized
productivity.
The two nonzero corrections have distinct roles.
The term involving $a_0$ corrects for estimation of the downstream
productivity nuisance $h_0$, whereas the term involving $\pi_0$
corrects for estimation of the projection $q_0$ used to construct the
instrument. Hence the application combines a constructed instrument
with locally robust estimation of that instrument and of latent
productivity.
The procedure is directly implementable with cross-fitting. On each
training sample, estimate $g_{01}$ by regressing $Y_1$ on $X_{11}$,
estimate $q_0$ by an additive regression of $R$ on $X_{11}$ and
$X_{12}$, thereby enforcing the projection space $\mathcal H$, and form $\hat Z=R-\hat q$. Next, estimate $h_0$
by regressing $Y_1-\hat g_{01}(X_{11})$ on $X_{21}$ and estimate
$a_0(X_{21})=\mathbb E[Z_0\mid X_{21}]$. Finally, estimate the additive
projection $\pi_0(\cdot;\theta)$. Since
\[
\pi_0(X_{11},X_{12};\theta)
=
\Pi_{\mathcal H}Y_2
-
\theta
\Pi_{\mathcal H}h_0(X_{21}),
\]
the two projections on the right-hand side can be estimated once and
combined for each value of $\theta$. The resulting cross-fitted moment
in \eqref{eq:prod_lr} can then be used in the GMM estimator of
Section~\ref{sec:estim}.
The residualized outcome
\[
\tilde Y_1
=
Y_1-g_{01}(X_{11})
\]
is an instance of the generated-outcome map
$\varphi(W,g_0)=Y-g_0(X_1)$ of Remark~\ref{rem:resid}, here applied to
latent productivity rather than to income or consumption. Two features
distinguish the present application. First, the second-step correction
no longer vanishes: unlike in the partially linear case, where
Robinson's moment is already orthogonal to $h_0$, the moment in
\eqref{eq:prod_moment} depends on $h_0$ through the identifying
equation, so that
$\alpha_{02}(X_{21};\theta)=-\theta a_0(X_{21})\neq0$. Second, the
identifying instrument is itself generated by a first-step projection,
and locally robust inference must therefore correct for two distinct
first steps.
\section{Estimation, Inference, and Finite-Sample Performance}\label{sec:estim}
\subsection{Estimation and Inference}
The orthogonal moment $\psi$ in Theorem~\ref{thm:lr} requires consistent estimators of the nuisance parameters $g_0$, $h_0$, and the Riesz representers $\alpha_0\equiv(\alpha_{01},\alpha_{02})$. We suppress the dependence of $\alpha_0(\theta)$ and
$\hat{\alpha}_\ell(\theta)$ on $\theta$ whenever no confusion arises. In our examples, the orthogonality conditions \eqref{eq:riesz1}--\eqref{eq:riesz2} admit analytical solutions, so consistent plug-in estimators are available directly from consistent estimators of the nuisance parameters they depend on. In settings where analytical expressions for the Riesz representers are unavailable or cumbersome to derive, automatic estimation \citep{chernozhukov2022locally} can be used in their place, without any modification to the locally robust moment construction of Theorem~\ref{thm:lr}.
To avoid own-observation bias and circumvent the need to verify Donsker-type conditions, we use cross-fitting: the sample
$(W_i)_{i=1}^n$ is partitioned into $L$ folds
$\mathcal I_\ell$, $\ell=1,\ldots,L$, and for each fold we obtain
estimators $\hat g_\ell$, $\hat h_\ell$, and
$\hat\alpha_\ell\equiv(\hat\alpha_{1\ell},\hat\alpha_{2\ell})$
using only observations outside $\mathcal I_\ell$.
Let
\[
\hat S_{i\ell}\equiv S(\hat g_\ell)(W_i)
\]
denote the cross-fitted generated outcome evaluated at observation
$i$. The debiased sample moment function is then
\begin{align}
\label{eq:sample_moment}
\hat\psi(\theta)
\equiv
\frac{1}{n}
\sum_{\ell=1}^{L}
\sum_{i\in\mathcal I_\ell}
\hat\psi_{i\ell}(\theta),
\end{align}
where
\begin{align*}
\hat\psi_{i\ell}(\theta)
\equiv
m(W_i,\hat g_\ell,\hat h_\ell,\theta)
+
\hat\alpha_{1\ell}(X_{1i})'\bigl(Y_i-\hat g_\ell(X_{1i})\bigr)
+
\hat\alpha_{2\ell}(X_{2i})'\bigl(\hat S_{i\ell}-\hat h_\ell(X_{2i})\bigr).
\end{align*}
When there is more than one moment condition, each component of $m$ is debiased by its own Riesz representers. The debiased GMM estimator is
\begin{align}
\label{eq:gmm}
\hat{\theta} = \operatorname*{arg\,min}_{\theta \in \Theta}\, \hat{\psi}(\theta)'\hat{\Upsilon}\,\hat{\psi}(\theta),
\end{align}
where $\hat{\Upsilon}$ is a positive semi-definite weighting matrix of dimension $q \times q$. \par
\bigskip
\textbf{Example: PLM with generated dependent variable, cont.}
In the partially linear model, the moment condition is exactly
identified and linear in $\theta_0$, so the locally robust moment in
equation~\eqref{eq:lr_plm_general_main} gives $\theta_0$ in closed
form,
\begin{align*}
\theta_0
&=
B^{-1}C,
\nonumber
\\
B
&\coloneqq
\mathbb{E}\left[
(X_3-h_{0X}(X_2))
(X_3-h_{0X}(X_2))'
\right],
\nonumber
\\
C
&\coloneqq
\mathbb E\left[
(\tilde Y-h_{0Y}(X_2))
(X_3-h_{0X}(X_2))
+
\alpha_{01}(X_1)(Y-\tilde Y)
\right],
\end{align*}
and the debiased GMM estimator in equation~\eqref{eq:gmm} reduces to
the sample analogue
\begin{align}
\label{eq:plm_estimator}
\hat\theta
=
\hat B^{-1}\hat C,
\end{align}
where
\begin{align*}
\hat B
&\coloneqq
\frac{1}{n}
\sum_{\ell=1}^L
\sum_{i\in\mathcal I_\ell}
\bigl(
X_{3i}-\hat h_{X\ell}(X_{2i})
\bigr)
\bigl(
X_{3i}-\hat h_{X\ell}(X_{2i})
\bigr)',
\\
\hat C
&\coloneqq
\frac{1}{n}
\sum_{\ell=1}^L
\sum_{i\in\mathcal I_\ell}
\Big[
\bigl(
\hat g_\ell(X_{1i})
-
\hat h_{Y\ell}(X_{2i})
\bigr)
\bigl(
X_{3i}-\hat h_{X\ell}(X_{2i})
\bigr)
\\
&\hspace{4cm}
+
\hat\alpha_{01\ell}(X_{1i})
\bigl(
Y_i-\hat g_\ell(X_{1i})
\bigr)
\Big].
\end{align*}
Here
\[
\hat h_{Y\ell}(X_2)
=
\hat{\mathbb E}_\ell
[
\hat g_\ell(X_1)\mid X_2
],
\qquad
\hat h_{X\ell}(X_2)
=
\hat{\mathbb E}_\ell
[
X_3\mid X_2
],
\]
where $\hat{\mathbb E}_\ell$ denotes ML estimators of the indicated
conditional means. The estimator $\hat\alpha_{01\ell}$ targets
\[
\alpha_{01}
=
\Pi_{\Delta_1}\{X_3-h_{0X}(X_2)\},
\]
with each fitted component belonging to $\Delta_1$.
All nuisance functions are estimated using observations outside
fold $\mathcal I_\ell$ and satisfy the convergence conditions
specified below.
We now establish consistency and asymptotic normality of $\hat{\theta}$ in the general case of (\ref{eq:gmm}). Both results follow the proof strategy of \cite{chernozhukov2022locally}, Theorem A3 for consistency and Theorem 9 for asymptotic normality, adapted to the compositional setting in which $\hat h_\ell$ is constructed from $S(\hat g_\ell)$. This dependence is handled through joint nuisance and remainder conditions rather than by treating the two nuisance estimators as independent first steps.\par
Consistency follows from standard M-estimation theory adapted to the locally robust framework. The following regularity conditions are required.
\begin{assumption}
\label{ass:consistency}\leavevmode
\begin{enumerate}[(i)]
\item The parameter space $\Theta \subset \mathbb{R}^p$ is compact and $\theta_0$ is the unique solution to the moment condition $\mathbb{E}[\psi(W, g_0, h_0, \alpha_0, \theta)] = 0$ over $\theta \in \Theta$.
\item $\mathbb{E}[\|\psi(W, g_0, h_0, \alpha_0, \theta_0)\|] < \infty$ and $\int \|m(w, \hat{g}_\ell, \hat{h}_\ell, \theta) - m(w, g_0, h_0, \theta)\| dF_0(w) \overset{p}{\to} 0$ for all $\theta \in \Theta$.
\item There exist $C>0$ and $d(W)$ with
$\mathbb{E}[d(W)]<\infty$ such that, almost surely,
for all $\theta,\tilde{\theta}\in\Theta$,
\[
\|\psi(W,g_0,h_0,\alpha_0,\theta)
-\psi(W,g_0,h_0,\alpha_0,\tilde{\theta})\|
\leq d(W)\|\theta-\tilde{\theta}\|^{1/C}.
\]
Moreover, for each fold $\ell$,
\[
\int\sup_{\theta\in\Theta}
\|\psi(w,\hat{g}_\ell,\hat{h}_\ell,\hat{\alpha}_\ell,\theta)
-\psi(w,g_0,h_0,\alpha_0,\theta)\|\,dF_0(w)
\overset{p}{\to}0.
\]
\item For each fold $\ell$, the estimated Riesz representers satisfy \[ \hat{\alpha}_{1\ell}\in\Delta_1^q, \qquad \hat{\alpha}_{2\ell}\in\Delta_2^q \] with probability approaching one. Moreover, with $\phi=\phi_1+\phi_2$,
\begin{align*}
\int \|\phi(w, \hat{g}_\ell, h_0, \alpha_0, \theta_0) - \phi(w, g_0, h_0, \alpha_0, \theta_0)\|^2 dF_0(w) &\overset{p}{\to} 0,\\
\int \|\phi(w, g_0, \hat{h}_\ell, \alpha_0, \theta_0) - \phi(w, g_0, h_0, \alpha_0, \theta_0)\|^2 dF_0(w) &\overset{p}{\to} 0,\\
\int \|\phi(w, g_0, h_0, \hat{\alpha}_\ell, \theta_0) - \phi(w, g_0, h_0, \alpha_0, \theta_0)\|^2 dF_0(w) &\overset{p}{\to} 0.
\end{align*}
\item For each fold $\ell$,
$\frac{1}{\sqrt{n}}\sum_{i\in\mathcal{I}_\ell}
\hat{\Delta}_\ell(W_i)\overset{p}{\to}0$, where
\begin{align*}
\hat{\Delta}_\ell(w)\coloneqq{}&
\psi(w,\hat{g}_\ell,\hat{h}_\ell,\hat{\alpha}_\ell,\theta_0)
-\psi(w,\hat{g}_\ell,\hat{h}_\ell,\alpha_0,\theta_0)\\
&-\psi(w,g_0,h_0,\hat{\alpha}_\ell,\theta_0)
+\psi(w,g_0,h_0,\alpha_0,\theta_0).
\end{align*}
\item $\hat{\Upsilon} \overset{p}{\to} \Upsilon$ positive definite.
\end{enumerate}
\end{assumption}
The identification condition in $(i)$ is standard. The $L_1$ continuity condition in $(ii)$ ensures the difference between the moment evaluated at estimated and true nuisance parameters becomes negligible as $n \to \infty$. Condition~$(iii)$ ensures continuity of the population moment
and uniform approximation of the estimated score by its
true-nuisance counterpart, yielding the uniform convergence
needed for consistency. Conditions~$(iv)$ and $(v)$ are standard mean-square
consistency and remainder conditions, and condition~$(vi)$
imposes convergence of the weighting matrix. The dependence of
$\hat h_\ell$ on $\hat g_\ell$ through the generated outcome is
handled by these joint nuisance and remainder conditions rather
than by imposing a separate product-rate restriction at the
consistency stage.
\begin{lemma}[Consistency]
\label{lem:consistency}
Under Assumptions~\ref{ass:lr} and \ref{ass:consistency}, the debiased GMM estimator $\hat{\theta}$ in equation \eqref{eq:gmm} satisfies
\[
\hat{\theta} \overset{p}{\to} \theta_0.
\]
\end{lemma}
Having established consistency, asymptotic normality requires additional structure on the moment integrability and Jacobian stability beyond Assumption~\ref{ass:consistency}. Let $\bar{\psi}$ denote the expected debiased moment evaluated at estimated nuisances, $\bar{\psi}(\hat{g}_\ell, \hat{h}_\ell, \alpha_0, \theta_0) \equiv \mathbb{E}[\psi(W, \hat{g}_\ell, \hat{h}_\ell, \alpha_0, \theta_0)]$.
\begin{assumption}
\label{ass:normality}\leavevmode
\begin{enumerate}[(i)]
\item $\mathbb{E}[\|\psi(W, g_0, h_0, \alpha_0, \theta_0)\|^2] < \infty$.
\item $\theta_0\in\operatorname{int}(\Theta)$. The Jacobian
\[
G\coloneqq
\mathbb{E}\left[
\partial_\theta
\psi(W,g_0,h_0,\alpha_0,\theta_0)
\right]
\]
exists and $G'\Upsilon G$ is nonsingular.
\item Let $\eta=(g,h,\alpha)$ and
$\eta_0=(g_0,h_0,\alpha_0)$. For all
$\|\eta-\eta_0\|$ sufficiently small,
$\psi(W,\eta,\theta)$ is differentiable with respect to $\theta$
with probability approaching one, and there exist $C>0$ and
$d(W)$ such that, for all $\theta\in\Theta$,
\[
\left\|
\frac{\partial\psi(W,\eta,\theta)}{\partial\theta}
-
\frac{\partial\psi(W,\eta,\theta_0)}{\partial\theta}
\right\|
\leq
d(W)\|\theta-\theta_0\|^{1/C},
\qquad
\mathbb{E}[d(W)]<C.
\]
\item For each fold $\ell$,
\[
\int
\left\|
\frac{\partial
\psi(w,\hat g_\ell,\hat h_\ell,\hat\alpha_\ell,\theta_0)}
{\partial\theta}
-
\frac{\partial
\psi(w,g_0,h_0,\alpha_0,\theta_0)}
{\partial\theta}
\right\|
dF_0(w)
\overset{p}{\to}0.
\]
\item For each fold $\ell$,
\begin{align*}
\int \|m(w,\hat{g}_\ell,\hat{h}_\ell,\theta_0)
-m(w,g_0,h_0,\theta_0)\|^2\,dF_0(w)
&\overset{p}{\to}0,\\
\int \|\hat{\Delta}_\ell(w)\|^2\,dF_0(w)
&\overset{p}{\to}0.
\end{align*}
\item For each fold $\ell$, one of the following holds: either $(a)$ there exists $C > 0$ such that, with probability tending to one,
\[
\left|\bar{\psi}(\hat{g}_\ell, \hat{h}_\ell, \alpha_0, \theta_0)\right| \leq C\left(\|\hat{g}_\ell - g_0\|_2^2 + \|\hat{h}_\ell - h_0\|_2^2\right),
\]
and $\sqrt{n}\left(\|\hat{g}_\ell - g_0\|_2^2 + \|\hat{h}_\ell - h_0\|_2^2\right) = o_p(1)$; or $(b)$
\[
\sqrt{n}\,\bar{\psi}(\hat{g}_\ell, \hat{h}_\ell, \alpha_0, \theta_0) \overset{p}{\to} 0.
\]
\end{enumerate}
\end{assumption}
Condition~$(i)$ strengthens the integrability requirement in Assumption~\ref{ass:consistency}$(ii)$ from $L_1$ to $L_2$, needed to establish a well-defined asymptotic variance and apply a central limit theorem. Condition~$(ii)$ guarantees the GMM problem is well-posed and the asymptotic variance is nonsingular. Since the influence adjustments have mean zero for every admissible $\alpha_0(\theta)$,
\[ G = \mathbb{E}\left[ \partial_\theta \psi(W,g_0,h_0,\alpha_0,\theta_0) \right] = \mathbb{E}\left[ \partial_\theta m(W,g_0,h_0,\theta_0) \right]. \] Conditions~$(iii)$ and~$(iv)$ control local stability of the Jacobian with respect to $\theta$ and the nuisance parameters respectively, ensuring that first-stage estimation error does not distort the asymptotic behavior of $\hat{\theta}$. Condition~$(v)$ supplies mean-square control of the identifying
moment and the interaction remainder, complementing
Assumption~\ref{ass:consistency}$(iv)$ on the influence corrections. Condition~$(vi)$ corresponds to Assumption~6.8 of \cite{escanciano2023automatic} and controls the first-order bias of the debiased moment when nuisance functions are estimated jointly: it requires either a quadratic remainder bound on the expected moment at estimated nuisances, or directly that $\sqrt{n}\bar{\psi}(\hat{g}_\ell, \hat{h}_\ell) \overset{p}{\to} 0$. Either condition ensures the $\sqrt{n}$-normalized expected remainder vanishes, which is not needed for Lemma~\ref{lem:consistency} but is required for asymptotic normality. Hadamard differentiability alone does not imply a quadratic
remainder bound. For the PLM, an exact expansion yields such
a bound, and condition~$(a)$ follows from
Assumption~\ref{ass:plm}$(iv)$; see the proof of
Corollary~\ref{cor:plm_normality}.
\begin{theorem}[Asymptotic Normality]
\label{thm:normality}
Under Assumptions~\ref{ass:lr}--\ref{ass:normality}, the debiased GMM estimator $\hat{\theta}$ in equation \eqref{eq:gmm} satisfies
\[
\sqrt{n}(\hat{\theta} - \theta_0) \xrightarrow{d} \mathcal{N}(0, V),
\] where $V = (G'\Upsilon G)^{-1}(G'\Upsilon\Psi\Upsilon G)(G'\Upsilon G)^{-1}$, $\Psi = \mathbb{E}[\psi(W, g_0, h_0, \alpha_0, \theta_0)\psi(W, g_0, h_0, \alpha_0, \theta_0)']$, and $\Upsilon$ is the probability limit of $\hat{\Upsilon}$. Moreover, the plug-in estimator is consistent: $\hat{V} \overset{p}{\to} V$.
\end{theorem}
For inference, let
\[
\hat G
\equiv
\frac{\partial\hat\psi(\hat\theta)}{\partial\theta},
\qquad
\hat\Psi
\equiv
\frac{1}{n}\sum_{\ell=1}^L\sum_{i\in\mathcal I_\ell}
\hat\psi_{i\ell}(\hat\theta)\hat\psi_{i\ell}(\hat\theta)'.
\]
The plug-in variance estimator is
\[
\hat V
=
(\hat G'\hat\Upsilon\hat G)^{-1}
\hat G'\hat\Upsilon\hat\Psi\hat\Upsilon\hat G
(\hat G'\hat\Upsilon\hat G)^{-1}.
\]
We now verify that Assumptions~\ref{ass:consistency} and \ref{ass:normality} are satisfied for the PLM estimator in equation~\eqref{eq:plm_estimator} under the following primitive conditions, and establish its asymptotic normality as a corollary of Theorem~\ref{thm:normality}.\par
\begin{assumption}
\label{ass:plm}
The following conditions hold for the partially linear model:
\begin{enumerate}[(i)]
\item $B = \mathbb{E}[(X_3 - h_{0X}(X_2))(X_3 - h_{0X}(X_2))']$ exists and is nonsingular.
\item
$\mathbb E[|Y|^4]<\infty$,
$\mathbb E[\|X_3\|^4]<\infty$,
$\mathbb E[|g_0(X_1)|^4]<\infty$, and
$\mathbb E[\|\alpha_{01}(X_1)\|^4]<\infty$.\item For each fold $\ell$,
$\hat g_\ell\in\Delta_1$ and
$\hat\alpha_{01\ell}\in\Delta_1^p$
with probability approaching one. For each fold $\ell$,
\[
\|\hat{g}_\ell-g_0\|_4
+\|\hat{h}_\ell-h_0\|_4
+\|\hat{\alpha}_{01\ell}-\alpha_{01}\|_4
\overset{p}{\to}0,
\]
where $\|\cdot\|_4$ denotes the $L_4(F_0)$ norm,
with fitted functions held fixed when integrating.
\item For each fold $\ell$,
\begin{align*}
\sqrt{n}\left(\|\hat{g}_\ell-g_0\|_2^2
+\|\hat{h}_\ell-h_0\|_2^2\right)&=o_p(1),\\
\sqrt{n}\,
\|\hat{\alpha}_{01\ell}-\alpha_{01}\|_2
\|\hat{g}_\ell-g_0\|_2&=o_p(1).
\end{align*}
\end{enumerate}
\end{assumption}
\begin{corollary}[Asymptotic Normality for the PLM]\label{cor:plm_normality}
Under Assumption~\ref{ass:plm}, the estimator $\hat{\theta}$ in equation~\eqref{eq:plm_estimator} satisfies
\[\sqrt{n}(\hat{\theta} - \theta_0) \xrightarrow{d} \mathcal{N}(0, V),\]
where $V = B^{-1}\Psi B^{-1}$ and $\Psi = \mathbb{E}[\psi(W, g_0, h_0, \alpha_0, \theta_0)\psi(W, g_0, h_0, \alpha_0, \theta_0)']$. Moreover, the plug-in estimator $\hat{V} \overset{p}{\to} V$.
\end{corollary}\subsection{Monte Carlo Simulations}
This section studies the finite-sample behavior of the proposed estimator
through Monte Carlo simulations for the partially linear model with an
ML-generated dependent variable. We begin with a first design in which
the covariates used to construct the generated outcome contain both the
controls and the regressor of interest. This design illustrates the
cancellation result in equation~\eqref{eq:lr_plm_main}: although the
parameter of interest is defined through an ML-generated outcome, the
locally robust moment eliminates the first-step regularization bias and
can be implemented using the observed outcome directly.
\paragraph{Design 1: generated outcome using all partially linear covariates.}
This simulation design mimics the following empirical scenario. A researcher seeks to estimate the effect of a binary variable of interest $X_3$ on an economic outcome, controlling for confounders $X_2$. The observed outcome $Y$ is noisy, and the researcher exploits additional high-dimensional features $Z$, such as unstructured text, images, or satellite-derived measurements, to construct a denoised prediction of the outcome using machine learning. The parameter of interest is the effect of $X_3$ on this ML-predicted outcome, controlling for $X_2$.
Specifically, we consider the data-generating process
\[Y = \theta_0 X_3+r_0(X_2)+s_0(Z)+\varepsilon, \qquad \varepsilon\sim\mathcal N(0,1),\]
where
\[r_0(X_2)=X_{21}^2+0.5X_{22}, \qquad s_0(Z)=0.1\sum_{j=1}^5 Z_j.\]
The confounders satisfy
\[X_{2j}\overset{\mathrm{i.i.d.}}{\sim}\mathcal N(0,1), \qquad j=1,\ldots,5,\]
the additional features satisfy
\[Z_j\overset{\mathrm{i.i.d.}}{\sim}\mathcal N(0,1), \qquad j=1,\ldots,10,\]
and the binary variable of interest is generated according to
\[X_3 = \mathbf 1\left(0.25X_{21}+0.25X_{22}+\eta>0\right), \qquad \eta\sim\mathcal N(0,1).\]
Thus $X_3$ is correlated with $X_2$, making the partialling-out step nontrivial. The true parameter is $\theta_0=1$.
In Design 1, the generated outcome is constructed using
\[X_1=(X_2,X_3,Z)\in\mathbb R^{16}.\]
Hence
\[\tilde Y \equiv g_0(X_1) \equiv \mathbb E[Y\mid X_1] = \theta_0X_3+r_0(X_2)+s_0(Z).\]
Since $X_2\subseteq X_1$ and $X_3\subseteq X_1$, the first-step Riesz representer simplifies to
\[\alpha_{01}(X_1) = X_3-h_{0X}(X_2), \qquad \alpha_{02}=0.\]
Therefore, the locally robust score reduces to Robinson's partialling-out moment evaluated at the observed outcome $Y$ rather than at the generated outcome $g(X_1)$, as established in equation \eqref{eq:lr_plm_main}. Thus, although the parameter of interest is defined through the ML-generated outcome $\tilde Y$, the locally robust estimator can be implemented without directly using $g$ in the final score.
The nuisance functions are estimated using $L=5$-fold cross-fitting. Let $\mathcal I_\ell$, $\ell=1,\ldots,L$, denote the folds. For each fold $\ell$, the nuisance estimators $\hat g_\ell$, $\hat h_{0Y,\ell}$, and $\hat h_{0X,\ell}$ are trained using only observations outside $\mathcal I_\ell$. Specifically, $\hat g_\ell$ is obtained by regressing $Y$ on $X_1=(X_2,X_3,Z)$ using LightGBM, while
\[\hat h_{0Y,\ell}(X_2) = \hat{\mathbb E}_\ell[Y\mid X_2], \qquad \hat h_{0X,\ell}(X_2) = \hat{\mathbb E}_\ell[X_3\mid X_2],\]
are obtained by regressing $Y$ and $X_3$ on $X_2$ using LightGBM. Estimating $\hat h_{0Y,\ell}$ directly from $Y$ rather than from $\hat g_\ell$ is valid since
\[\mathbb E[Y\mid X_2] = \mathbb E[g_0(X_1)\mid X_2]\]
by the law of iterated expectations.
The simulations compare three estimators.
\paragraph{Oracle estimator.}
The infeasible oracle estimator uses the population nuisance functions
\[h_{0Y}(X_2)=\mathbb E[Y\mid X_2], \qquad h_{0X}(X_2)=\mathbb E[X_3\mid X_2],\]
and solves
\[\hat\theta_{\mathrm{oracle}} = \left(\frac1n \sum_{i=1}^n (X_{3i}-h_{0X}(X_{2i}))^2\right)^{-1} \frac1n \sum_{i=1}^n (Y_i-h_{0Y}(X_{2i}))(X_{3i}-h_{0X}(X_{2i})).\]
\paragraph{Locally robust estimator.}
The feasible locally robust estimator uses the cross-fitted nuisance estimates evaluated only on held-out observations and solves
\[\hat\theta_{\mathrm{LR}} = \hat B^{-1}\hat C,\]
where
\[\hat B = \frac1n \sum_{\ell=1}^L \sum_{i\in\mathcal I_\ell} \left(X_{3i}-\hat h_{0X,\ell}(X_{2i})\right)^2,\]
and
\[\hat C = \frac1n \sum_{\ell=1}^L \sum_{i\in\mathcal I_\ell} \left(Y_i-\hat h_{0Y,\ell}(X_{2i})\right) \left(X_{3i}-\hat h_{0X,\ell}(X_{2i})\right).\]
\paragraph{Naive plug-in estimator.}
The naive plug-in estimator ignores the generated nature of the dependent variable and instead applies Robinson's partialling-out score directly to the cross-fitted ML-generated outcome:
\[\hat\theta_{\mathrm{plugin}} = \hat B^{-1}\hat C_{\mathrm{plugin}},\]
where
\[\hat C_{\mathrm{plugin}} = \frac1n \sum_{\ell=1}^L \sum_{i\in\mathcal I_\ell} \left(\hat g_\ell(X_{1i})-\hat h_{0Y,\ell}(X_{2i})\right) \left(X_{3i}-\hat h_{0X,\ell}(X_{2i})\right).\]
This estimator ignores the regularization bias introduced by the first-step ML estimation of $g_0$.
Table~\ref{tab:sims} reports bias, root mean square error (RMSE), and empirical coverage of nominal 95\% confidence intervals across $R=1,000$ Monte Carlo replications for sample sizes $n\in\{500,1000,2000,5000\}$. The locally robust estimator exhibits sound finite-sample performance. Its bias is negligible across all sample sizes, fluctuating around zero, and its RMSE is very similar to the oracle's. Coverage ranges from $0.934$ to $0.951$ and tracks that of the infeasible oracle. In contrast, the naive plug-in estimator exhibits a sizable bias that declines more slowly than $n^{-1/2}$, with coverage collapsing from $0.078$ to $0.003$.
\begin{table}[H]
\centering \resizebox{\textwidth}{!}{
\begin{threeparttable} \begin{tabular}{l ccc ccc ccc} \hline & \multicolumn{3}{c}{Oracle} & \multicolumn{3}{c}{Locally Robust} & \multicolumn{3}{c}{Naive Plug-in} \\ \cmidrule(lr){2-4} \cmidrule(lr){5-7} \cmidrule(lr){8-10}
$n$ & Bias & RMSE & Coverage & Bias & RMSE & Coverage & Bias & RMSE & Coverage \\
\hline\hline
500 & 0.003 & 0.095 & 0.936 & 0.010 & 0.111 & 0.951 & -0.277 & 0.304 & 0.078 \\
1,000 & 0.000 & 0.069 & 0.947 & -0.004 & 0.077 & 0.943 & -0.209 & 0.224 & 0.035 \\
2,000 & 0.001 & 0.048 & 0.939 & -0.002 & 0.053 & 0.934 & -0.147 & 0.157 & 0.028 \\
5,000 & -0.002 & 0.029 & 0.956 & -0.006 & 0.031 & 0.952 & -0.100 & 0.105 & 0.003 \\
\hline
\end{tabular}
\begin{tablenotes}
\item Results based on $R=1,000$ Monte Carlo replications. The true parameter is $\theta_0=1$. Nominal coverage is $0.95$. All nuisance functions are estimated via LightGBM with $L=5$-fold cross-fitting.
\end{tablenotes}
\end{threeparttable}}
\caption{Finite-sample performance of the oracle, locally robust, and naive plug-in estimators in Design 1.}
\label{tab:sims}
\end{table}
\paragraph{Design 2: residualized generated outcome with a nested second step.}
This design mimics a second empirical scenario. A researcher is interested in the association of a binary variable of interest $X_3$, say sex, with lifetime income, but the outcome is first purged of confounding variation, such as year and cohort dummies, before projecting the residualized outcome on $X_3$ controlling for $X_2$. The parameter of interest is therefore the partially linear projection coefficient of the residualized outcome, as in Remark~\ref{rem:plm_projection}. This is also the general case of our framework, in which the first-step correction does not collapse and the second-step nuisance depends on the first step.\par
We mantain the outcome equation of Design 1 and the distributions of the confounders $X_2$, the shock $\eta$ and the treatment $X_3$, with $\beta_0=1$, and modify four features:
\[Z_1 = \rho\,\eta+\sqrt{1-\rho^2}\,\omega, \quad \rho=0.6, \qquad s_0(Z)=0.1\sum_{j=2}^{6}Z_j, \qquad \dim(Z)=50,\]
\[X_1=(X_{21},X_{22},Z)\in\mathbb R^{52}, \qquad r_0(X_2)=X_{21}^2+0.5X_{22}+X_{23},\]
with $\omega\sim\mathcal N(0,1)$ and $Z_j\overset{\mathrm{i.i.d.}}{\sim}\mathcal N(0,1)$ for $j=2,\ldots,50$. The proxy $Z_1$ makes $X_1$ partially informative about the treatment's identifying variation $\nu\equiv X_3-h_{0X}(X_2)$, so that $\alpha_{01}$ neither vanishes nor equals $\nu$; excluding $Z_1$ from $s_0$ keeps it free of any direct effect on the outcome; the higher dimension prevents the regularization bias in $\hat g$ from vanishing too quickly; and excluding $X_{23}$ from $X_1$ makes the second step nontrivial.
The generated outcome is the observed outcome residualized of the variation explained by $X_1$,
\[\tilde Y \equiv Y-g_0(X_1), \qquad g_0(X_1)\equiv\mathbb E[Y\mid X_1],\]
as in Remark \ref{rem:resid}, so that $\alpha_{02}=0$ and $\alpha_{01}(X_1)=-\mathbb E[\nu\mid X_1]$, and
\[\psi(W,g_0,h_0,\alpha_0,\theta) = \left(\tilde Y-h_{0Y}(X_2)-\theta\nu\right)\nu + \alpha_{01}(X_1)\left(Y-g_0(X_1)\right).\]
Two features distinguish this design from Design 1. First, $\alpha_{01}$ must be estimated, since $X_3\not\subseteq X_1$. Second, $h_{0Y}(X_2)=X_{23}\neq0$ cannot be recovered from $Y$ alone, so the second step must be estimated from the first. This makes Design~2 a genuinely nested case of the framework in Section~\ref{sec:set}, in which estimation error from the first step propagates into the second-step nuisance and is handled through the joint remainder conditions of our asymptotic theory. The estimand is
\[\theta_0 = \beta_0\left(1-R^2_{\nu\mid X_1}\right), \qquad R^2_{\nu\mid X_1} \equiv \frac{\mathbb E[\alpha_{01}^2]}{\mathbb E[\nu^2]},\]
which admits a closed form in this design and equals $0.23$ for $\rho=0.6$, so that $\theta_0=0.77$. Unlike Design 1, the partially linear projection coefficient differs from $\beta_0$: residualizing removes the component of the identifying variation $\nu$ that is explained by $X_1$.
Nuisances are estimated as in Design 1, with $L=5$-fold cross-fitting and two differences. The Riesz representer is estimated by $\hat\alpha_{01,\ell}=-\big(\hat m_\ell(X_1)-\hat{\mathbb E}_\ell[\hat h_{0X,\ell}(X_2)\mid X_1]\big)$, with both conditional expectations of the binary $X_3$ obtained by a logistic objective and an inner cross-fit for the values of $\hat h_{0X}$ entering the projection; projecting onto $X_1$ makes the representer $X_1$-measurable by construction, without imposing any condition on the data-generating process. The second-step nuisance $\hat h_{0Y,\ell}$ is obtained by regressing $Y-\hat g$ on $X_2$, again using an inner cross-fit so that the first-step fits are out of sample.
The simulations compare three estimators.
\paragraph{Oracle estimator.}
The infeasible oracle estimator uses the population nuisance functions and solves
\[\hat\theta_{\mathrm{oracle}} = \hat B_{\mathrm{oracle}}^{-1}\hat C_{\mathrm{oracle}},\]
where\fontsize{10.5}{12.5}\selectfont
\begin{align*}
\hat B_{\mathrm{oracle}} = \frac1n \sum_{i=1}^n \nu_i^2, \quad \hat C_{\mathrm{oracle}} = \frac1n \sum_{i=1}^n \Big[\left(Y_i-g_0(X_{1i})-h_{0Y}(X_{2i})\right)\nu_i + \alpha_{01}(X_{1i})\left(Y_i-g_0(X_{1i})\right)\Big].
\end{align*} \normalsize
\paragraph{Locally robust estimator.}
The feasible locally robust estimator uses the cross-fitted nuisance estimators evaluated only on held-out observations and solves
\[\hat\theta_{\mathrm{LR}} = \hat B^{-1}\hat C,\]
where
\[\hat B = \frac1n \sum_{\ell=1}^L \sum_{i\in\mathcal I_\ell} \hat\nu_{i\ell}^2, \qquad \hat\nu_{i\ell} \equiv X_{3i}-\hat h_{0X,\ell}(X_{2i}),\]
and
\[\hat C = \frac1n \sum_{\ell=1}^L \sum_{i\in\mathcal I_\ell} \Big[\left(Y_i-\hat g_\ell(X_{1i})-\hat h_{0Y,\ell}(X_{2i})\right)\hat\nu_{i\ell} + \hat\alpha_{01,\ell}(X_{1i})\left(Y_i-\hat g_\ell(X_{1i})\right)\Big].\]
\paragraph{Naive plug-in estimator.}
The naive plug-in estimator ignores the first-step correction term and applies Robinson's partialling-out score directly to the cross-fitted generated outcome:
\[\hat\theta_{\mathrm{plugin}} = \hat B^{-1}\hat C_{\mathrm{plugin}},\]
where
\[\hat C_{\mathrm{plugin}} = \frac1n \sum_{\ell=1}^L \sum_{i\in\mathcal I_\ell} \left(Y_i-\hat g_\ell(X_{1i})-\hat h_{0Y,\ell}(X_{2i})\right)\hat\nu_{i\ell}.\]
This estimator ignores the regularization bias introduced by the first-step ML estimation of $g_0$.\par
Table~\ref{tab:sims2} reports the same statistics for Design 2. The locally robust estimator again exhibits sound finite-sample performance, with negligible bias across all sample sizes and RMSE approaching that of the oracle. However, the coverage ranges from $0.918$ to $0.934$, some two percentage points below the oracle's. Two sources could account for this shortfall: the interval may be off-centre, or it may be the wrong width. To tell them apart, note that coverage depends on both through a single expression. Let $sd$ denote the Monte Carlo standard deviation of the estimator and write $\lambda \equiv \text{mean}(\hat{se})/sd$ and $b\equiv\text{bias}/sd$. Standardising $|\hat\theta-\theta_0|\le 1.96\,\hat{se}$ by $sd$ and using the asymptotic normality established in Theorem~\ref{thm:normality}, so that $(\hat\theta-\theta_0)/sd$ is approximately $\mathcal N(b,1)$, gives
\[\text{Coverage} \approx \Phi(1.96\lambda-b)-\Phi(-1.96\lambda-b),\]
which equals the nominal $0.95$ when $\lambda=1$ and $b=0$. The two departures can therefore be switched off separately.\par
\begin{table}[H]
\centering \resizebox{\textwidth}{!}{
\begin{threeparttable} \begin{tabular}{l ccc ccc ccc} \hline & \multicolumn{3}{c}{Oracle} & \multicolumn{3}{c}{Locally Robust} & \multicolumn{3}{c}{Naive Plug-in} \\ \cmidrule(lr){2-4} \cmidrule(lr){5-7} \cmidrule(lr){8-10}
$n$ & Bias & RMSE & Coverage & Bias & RMSE & Coverage & Bias & RMSE & Coverage \\
\hline\hline
500 & 0.002 & 0.099 & 0.946 & 0.023 & 0.138 & 0.924 & 0.118 & 0.171 & 0.829 \\
1,000 & 0.001 & 0.071 & 0.947 & 0.000 & 0.089 & 0.924 & 0.094 & 0.123 & 0.792 \\
2,000 & -0.001 & 0.050 & 0.950 & -0.006 & 0.061 & 0.918 & 0.073 & 0.093 & 0.695 \\
5,000 & -0.002 & 0.031 & 0.952 & -0.003 & 0.035 & 0.934 & 0.055 & 0.064 & 0.582 \\
\hline
\end{tabular}
\begin{tablenotes}
\item Results based on $R=1,000$ Monte Carlo replications. The true parameter is $\theta_0=\beta_0(1-R^2_{\nu\mid X_1})=0.77$. Nominal coverage is $0.95$. All nuisance functions are estimated via LightGBM with $L=5$-fold cross-fitting.
\end{tablenotes}
\end{threeparttable}}
\caption{Finite-sample performance of the oracle, locally robust, and naive plug-in estimators in Design 2.}
\label{tab:sims2}
\end{table}
Table~\ref{tab:coverage_decomp} reports this decomposition for both designs. The standardised bias is negligible at every sample size, so that setting $b=0$ leaves coverage essentially unchanged; setting $\lambda=1$ instead, and leaving the bias as it is, restores coverage of $0.946$--$0.950$ in every cell. The shortfall in Design 2 is therefore attributable primarily to the reported standard error rather than to the estimator. The pattern is consistent with higher-order variability from nuisance estimation that is not captured by the first-order asymptotic variance approximation. The contrast between the two designs is instructive in this respect. Design 1 estimates three nuisance functions, and $\lambda$ is close to unity throughout, so that the cost of width never exceeds one percentage point and coverage is essentially nominal. Design 2 estimates five, including a Riesz representer and a second-step nuisance that is itself constructed from the first-step fit, and $\lambda$ falls to roughly $0.90$, with a corresponding cost of two to three percentage points. The gap closes as the sample grows ($\lambda$ rises to $0.947$ at $n=5,000$), consistent with the omitted contribution being of smaller order.
\begin{table}[H]
\centering
\begin{threeparttable} \begin{tabular}{l cc ccc c} \hline
$n$ & $\lambda$ & $b$ & Observed & Predicted & At $\lambda=1$ & Cost of width \\
\hline\hline
\multicolumn{7}{l}{\textit{Design 1}} \\\hline
500 & 0.994 & 0.094 & 0.951 & 0.947 & 0.949 & 0.002 \\
1,000 & 0.956 & -0.054 & 0.943 & 0.939 & 0.950 & 0.011 \\
2,000 & 0.953 & -0.032 & 0.934 & 0.938 & 0.950 & 0.012 \\
5,000 & 1.028 & -0.197 & 0.952 & 0.952 & 0.946 & -0.006 \\
\hline
\multicolumn{7}{l}{\textit{Design 2}} \\\hline
500 & 0.906 & 0.171 & 0.924 & 0.920 & 0.947 & 0.027 \\
1,000 & 0.915 & 0.002 & 0.924 & 0.927 & 0.950 & 0.023 \\
2,000 & 0.893 & -0.092 & 0.918 & 0.919 & 0.949 & 0.031 \\
5,000 & 0.947 & -0.091 & 0.934 & 0.936 & 0.949 & 0.014 \\
\hline
\end{tabular}
\begin{tablenotes}
\item Decomposition of the coverage attained by the locally robust estimator, based on $R=1,000$ Monte Carlo replications. Predicted coverage evaluates $\Phi(1.96\lambda-b)-\Phi(-1.96\lambda-b)$ at the measured $\lambda$ and $b$. The column at $\lambda=1$ sets the interval to the estimator's sampling standard deviation, leaving the bias unchanged; it is a counterfactual, not a feasible procedure. The final column is the difference between the two.
\end{tablenotes}
\end{threeparttable}
\caption{Decomposition of coverage into interval width and bias.}
\label{tab:coverage_decomp}
\end{table}
Monte Carlo simulations illustrate that our framework delivers valid inference with a machine-learning generated dependent variable. The locally robust estimator is essentially unbiased in both designs, its root mean square error converges on that of the infeasible oracle, and its confidence intervals attain close to nominal coverage, with the residual shortfall traceable to the variance formula rather than to the estimator itself. The naive plug-in estimator, by contrast, carries a bias that vanishes too slowly for valid inference, and its coverage deteriorates as the sample grows. The framework thus offers a practical toolkit for valid inference in the presence of outcome predicted by ML.
\section{Empirical Application: Conditional Group Disparities in Mortgage Default Risk}
\label{sec:empirical}
\subsection{Empirical question and data}
We apply our method to study conditional group disparities in mortgage
default risk. The empirical question is whether a race-blind
machine-learning model of future mortgage performance assigns different
predicted default risk to Black and White borrowers who are comparable
in conventional underwriting characteristics, and how this disparity
compares with the corresponding disparity in realized default. The
exercise is descriptive: the resulting coefficients measure conditional
group disparities and do not have a causal interpretation.
We use the Public Use File of the National Survey of Mortgage
Originations (NSMO), combining borrower and loan characteristics measured
at origination with subsequent mortgage performance. The sample covers
mortgages originated between 2013 and 2023Q3. We define the outcome as
an indicator that the mortgage is reported as 90 or more days delinquent
or in a more severe adverse status at any of the eight quarterly
post-origination observations. We restrict the sample to non-Hispanic
borrowers who identify as Black only or White only and for whom this
post-origination performance history is sufficiently observed to
classify the outcome. Our final sample contains $43{,}070$ mortgages, of which
$3{,}183$ are to Black borrowers and $39{,}887$ to White borrowers.
There are $280$ observed serious-delinquency/default events over the
eight post-origination quarterly observations. Using the NSMO analysis
weights, the corresponding event frequencies are $2.57\%$ for Black
borrowers and $0.60\%$ for White borrowers.
Let $Y$ denote the indicator for subsequent serious delinquency/default
and let $A$ equal one for a Black borrower and zero for a White borrower.
We estimate the race-blind default-risk function
\[
g_0(X_1)=\mathbb{P}(Y=1\mid X_1),
\]
where $X_1$ contains rich information available at origination but
excludes race. The comparison variables $X_2\subset X_1$ contain a
smaller set of conventional underwriting characteristics: origination
credit score, loan-to-value ratio, debt-to-income ratio, origination
timing, number of borrowers, loan type, adjustable-rate status, and loan
purpose.
Our main parameter is the partially linear projection coefficient
$\theta_0$ of the predicted default risk $g_0(X_1)$ on $A$, controlling
flexibly for $X_2$, as in Remark~\ref{rem:credit_scoring}. We also
estimate the analogous coefficient $\theta_Y$ using the realized
outcome $Y$. Comparing $\theta_0$ with $\theta_Y$ therefore asks how
much of the conditional Black--White disparity in realized mortgage
performance is reflected in the race-blind ML risk score.
Moreover, because $X_2\subseteq X_1$ and
$g_0(X_1)=\mathbb{E}[Y\mid X_1]$,
\[
\theta_Y-\theta_0
=
\frac{
\mathbb{E}\left[
\{Y-g_0(X_1)\}
\{A-\mathbb{E}[A\mid X_1]\}
\right]
}{
\mathbb{E}\left[
\{A-\mathbb{E}[A\mid X_2]\}^2
\right]
}.
\]
Thus, the difference between the realized-default disparity and the
predicted-risk disparity reflects the component of the group disparity
in realized default that remains after conditioning on the rich
race-blind information set $X_1$.
\subsection{Estimation}
We retain the conditional-mean estimand
$g_0(X_1)=\mathbb E[Y\mid X_1]$ and implement the second step
sequentially, using generated outcomes. Specifically, the nuisance
functions are
\[
h_{0Y}(X_2)=\mathbb E[g_0(X_1)\mid X_2],
\qquad
h_{0A}(X_2)=\mathbb E[A\mid X_2],
\]
and the first-step Riesz representer is
\[
\alpha_{01}(X_1)
=
\mathbb E[\nu\mid X_1],
\qquad
\nu=A-h_{0A}(X_2).
\]
Although $h_{0Y}(X_2)=\mathbb E[Y\mid X_2]$ at the population
truth, we estimate $h_{0Y}$ from first-step predictions rather
than directly from $Y$. As explained in
Remark~\ref{rem:credit_scoring}, this construction keeps the
second-step regression aligned with the prediction entering the
identifying moment. It does not, by itself, establish valid
inference under persistent first-step misspecification.
We use five outer folds. For each evaluation fold, the first-step
risk model is fitted using only the complementary training sample.
Within that training sample, five-fold inner cross-fitting produces
out-of-sample predictions of the generated outcome. We regress
these predictions on $X_2$ to estimate $h_{0Y}$ and evaluate the
resulting regression on the outer evaluation fold. To estimate
$\alpha_{01}$, we similarly construct inner-cross-fitted residuals
$A-\widehat h_{0A}(X_2)$ and regress them on $X_1$.
Our main specification uses XGBoost. The models for $g_0$ and
$h_{0A}$ use a binary logistic objective, whereas the regressions
for $h_{0Y}$ and $\alpha_{01}$ use squared error. All tuning and
preprocessing exclude the corresponding evaluation observations.
We also report the sequential implementation using elastic net.
Let
\[
\widehat\nu_i
=
A_i-\widehat h_{0A}(X_{2i}),
\]
and denote by $\mathbb{E}_n^w[\cdot]$ the empirical average using the
NSMO analysis weights, normalized to sum to one. Define
\[
\widehat B
=
\mathbb{E}_n^w[\widehat\nu_i^2].
\]
The naive estimator treats the cross-fitted predicted default risk as
an observed outcome:
\[
\widehat\theta_{\mathrm{naive}}
=
\widehat B^{-1}
\mathbb{E}_n^w
\left[
\left\{
\widehat g(X_{1i})
-
\widehat h_{0Y}(X_{2i})
\right\}
\widehat\nu_i
\right].
\]
Our locally robust estimator instead uses
\[
\widehat\theta_{\mathrm{LR}}
=
\widehat B^{-1}
\mathbb{E}_n^w
\left[
\left\{
\widehat g(X_{1i})
-
\widehat h_{0Y}(X_{2i})
\right\}
\widehat\nu_i
+
\widehat\alpha_{01}(X_{1i})
\left\{
Y_i-\widehat g(X_{1i})
\right\}
\right].
\]
For the realized-default benchmark, we separately estimate
\[
\mu_{0Y}(X_2)=\mathbb E[Y\mid X_2]
\]
directly from observed outcomes, using the same outer folds, and
compute
\[
\widehat\theta_Y
=
\widehat B^{-1}\mathbb E_n^w
\left[
\{Y_i-\widehat\mu_{0Y}(X_{2i})\}\widehat\nu_i
\right].
\]
Although $\mu_{0Y}=h_{0Y}$ at the population truth, their estimators
use different training responses: observed outcomes for
$\widehat\mu_{0Y}$ and generated predictions for $\widehat h_{0Y}$.
Standard errors for $\widehat\theta_{\mathrm{LR}}$ and
$\widehat\theta_Y$ use the corresponding cross-fitted scores under
the maintained nuisance convergence and remainder conditions.
The standard error for their difference uses the joint scores.
The naive standard error treats predicted outcomes as observed;
the reported uncertainty for the LR correction is a plug-in
diagnostic. These calculations do not establish inference robust
to arbitrary persistent first-step misspecification. Our preferred
specification uses the NSMO analysis weights, and we also report
unweighted estimates.
\subsection{Results}
The race-blind XGBoost model has an out-of-sample AUC of $0.889$
in the analysis-weighted sample. Table~\ref{tab:mortgage_main}
reports the estimates from the sequential implementation. The
conditional Black--White disparity in realized serious delinquency
is $1.017$ percentage points. The locally robust estimate of the
disparity in predicted default risk is $0.082$ percentage points,
with a standard error of $0.034$ percentage points. The latter is
approximately $8\%$ of the estimated disparity in realized default.
Treating the cross-fitted predictions as observed yields an estimate
of $0.064$ percentage points. The locally robust correction is
$0.018$ percentage points, increasing the naive estimate by
approximately $28\%$. Its plug-in diagnostic interval includes
zero. We report this comparison as an illustration of the
estimation adjustment, rather than as a formal test of robustness
to first-step misspecification.
The difference between the estimated disparities in realized
default and predicted risk is $0.935$ percentage points, with a
standard error of $0.321$ percentage points. Thus, under the
maintained conditions for inference, the estimated conditional
disparity represented by the race-blind risk function is
substantially smaller than the corresponding disparity in
realized mortgage performance.
\begin{table}[t]
\centering
\begin{threeparttable}
\caption{Conditional Black--White Disparities in Mortgage Default Risk}
\label{tab:mortgage_main}
\begin{tabular}{lccc}
\hline
& Estimate & Standard error & $p$-value \\
\hline
Realized default: $\widehat\theta_Y$
& 1.017 & 0.324 & 0.002 \\
ML score: naive out-of-sample
& 0.064 & 0.031 & 0.037 \\
ML score: locally robust
& 0.082 & 0.034 & 0.016 \\
LR correction
& 0.018 & 0.015 & 0.244 \\
$\widehat\theta_Y-\widehat\theta_{\mathrm{LR}}$
& 0.935 & 0.321 & 0.004 \\
\hline
\end{tabular}
\begin{tablenotes}
\footnotesize
\item Notes: Estimates and standard errors are in percentage points.
Results use XGBoost, five outer folds, five inner folds for generated
training responses, and the NSMO analysis weights. The second-step
outcome regression is estimated from first-step predictions. The
ML default-risk model excludes race. The naive estimator treats
cross-fitted predictions as observed, whereas the locally robust
estimator adds the first-step correction. Standard errors for the
LR and realized-default estimators use their cross-fitted scores;
the difference uses joint scores. The naive and correction rows
report plug-in uncertainty diagnostics, not misspecification-robust
inference. The estimates are descriptive and have no causal
interpretation.
\end{tablenotes}
\end{threeparttable}
\end{table}
In the unweighted XGBoost specification, the locally robust estimate
is $0.092$ percentage points, compared with $0.895$ percentage
points for realized default. The sequential elastic-net
implementation gives a smaller analysis-weighted score-disparity
estimate of $0.010$ percentage points, with a standard error of
$0.040$ percentage points, compared with $1.076$ percentage
points for realized default. Thus, the magnitude and precision
of the estimated score disparity are learner-sensitive, although
the contrast with the substantially larger realized-default
disparity persists.
The fitted Riesz channel is small. In the preferred XGBoost
specification,
\[
\frac{\mathbb E_n^w[\widehat\alpha_{01}^2]}
{\mathbb E_n^w[\widehat\nu^2]}
=0.0015.
\]
The out-of-sample AUC for predicting race increases from $0.744$
using $X_2$ to $0.747$ using $X_1$. These diagnostics suggest
that the fitted models extract little additional information
about group membership from $X_1$ beyond the conventional
underwriting characteristics in $X_2$, consistent with a
first-step correction that is small in absolute terms.
\section{Conclusions}\label{sec:conc}
In this paper, we provide a framework for locally robust inference with outcomes predicted by machine learning. We consider a three-step setting in which, in the first step, an outcome is generated by ML; in the second step, a nuisance function is constructed using nonparametric methods that involve the generated outcome; and in the third step, the parameter of interest is identified through a GMM criterion. To attain valid inference, we construct a moment that is Neyman-orthogonal to the first- and second-step nuisances, exploiting the idea of downstream local robustness: insensitivity to the second-step nuisance annihilates the indirect effect of the first step through the second, so that only the direct effect of the first step on the moment condition needs to be corrected. By replacing restrictive parametric first steps with flexible machine learning methods, researchers can reduce sensitivity to functional-form misspecification, while our debiased moment condition corrects for the first-order regularization bias that ML estimation introduces, enabling valid inference under the stated conditions.\par
We establish the consistency and asymptotic normality of the resulting estimator and illustrate its finite-sample behavior through Monte Carlo simulations. The locally robust estimator exhibits small bias and coverage close to nominal levels, whereas the naive estimator that treats the ML-generated outcome as observed can display substantial bias and severe undercoverage. We also apply the method to conditional Black--White disparities
in U.S.\ mortgage default risk, estimating the second-step outcome
regression from generated first-step predictions. In the preferred
specification, the estimated disparity in predicted risk is
approximately $8\%$ of the corresponding disparity in realized
serious delinquency. The locally robust adjustment raises the
naive estimate of the predicted-risk disparity by approximately
$28\%$. The application illustrates the substantive interest of
parameters defined through ML-generated outcomes and the role of
sequential nuisance estimation. The reported inference remains
subject to the stated nuisance convergence and remainder
conditions.\par
\bibliography{cites}\clearpage