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.
123,209 characters
Automatic Debiased Machine Learning of Causal and Structural Effects
\title[Auto DML for Causal Effects]{Automatic Debiased Machine Learning of
Causal and Structural Effects}
\thanks{The present paper formed the basis of the Fisher-Schultz Lecture given by
Victor Chernozhukov at the 2019 European Meeting of the Econometric Society
in Manchester. This research was supported by NSF grants 1559172 and 1757140. Helpful
comments were provided by the editor G. Imbens, three referees, J. Robins,
Y. Zhu, and participants at a 2016 demand workshop at Boston College and a
2018 machine learning and statistical inference workshop at the Banff
International Research Station.}
\author{Victor Chernozhukov}
\author{Whitney K. Newey}
\author{Rahul Singh}
\date{October 7, 2021}
\begin{abstract}
Many causal and structural effects depend on regressions. Examples include
policy effects, average derivatives, regression decompositions, average
treatment effects, causal mediation, and parameters of economic structural
models. The regressions may be high dimensional, making machine learning
useful. Plugging machine learners into identifying equations can lead to
poor inference due to bias from regularization and/or model selection. This
paper gives automatic debiasing for linear and nonlinear functions of
regressions. The debiasing is automatic in using Lasso and the function of
interest without the full form of the bias correction. The debiasing can be
applied to any regression learner, including neural nets, random forests,
Lasso, boosting, and other high dimensional methods. In addition to
providing the bias correction we give standard errors that are robust to
misspecification, convergence rates for the bias correction, and primitive
conditions for asymptotic inference for estimators of a variety of
estimators of structural and causal effects. The automatic debiased machine
learning is used to estimate the average treatment effect on the treated for
the NSW job training data and to estimate demand elasticities from Nielsen
scanner data while allowing preferences to be correlated with prices and
income.
Keywords: Debiased machine learning, causal parameters, structural
parameters, regression effects, Lasso, Riesz representation.
\end{abstract}
\maketitle
\section{Introduction}
Many causal and structural parameters of economic interest depend on
regressions, i.e. on conditional expectations or least squares projections.
Examples include policy effects, average derivatives, regression
decompositions, average treatment effects, causal mediation, and parameters
of economic structural models. Often, regressions may be high dimensional,
depending on many variables. There may be many covariates for policy
effects, average derivatives, and treatment effects, or many prices and
covariates in the economic demand for some commodity. This paper is about
estimating economic and causal parameters that depend on high dimensional
regressions.
Machine learning is a collection of modern, adaptive statistical learning
methods for estimating regression functions and other statistical objects.
These methods exploit structured parsimony restrictions (such as approximate
sparsity) on regressions, together with various forms of regularization and
model selection, to enable high quality prediction in high dimensional
settings. Key methods include neural nets (deep learning), random forests,
and Lasso. The goal of this paper is to deploy these methods to infer causal
and structural parameters that depend on regression functions, including
policy, derivative, decomposition, and treatment effects as well as economic
structural parameters.
Machine learning is different than other methods in ways that are useful in
high dimensional settings. For example, Lasso has good properties with very
many potential regressors (possibly many more than sample size) when
relatively few important regressors give a good approximation but the
identity of those few is not known (i.e. the regression is approximately
sparse). In contrast, series regression is based on relatively few
regressors, often many fewer than the sample size. Lasso and series
regression are similar in that they both depend on a few regressors giving a
good approximation. They differ in that series regression requires that the
identity of the important regressors is known, while with Lasso their
identity need not be known. For Lasso, the important regressors just need to
be included somewhere among the many potential regressors. This difference
is useful in high dimensional settings, where there are potentially very
many regressors needed to approximate a function of many variables.
Typically, economics and statistics provide little guidance about which
regressors are important. With Lasso, such information is not needed, since
very many terms can be included among the potential regressors. Other
machine learning methods, such as random forests and neural nets, are also
well suited to high dimensional regression.
Machine learners provide remarkably good predictions in a variety of
settings but are inherently biased. The bias arises from using
regularization and/or model selection to control the variance of the
prediction. To obtain small mean squared prediction errors, machine learners
regularize and/or select among models so that variance and squared bias are
approximately equal. Although such equality is good for prediction, it is
not good for inference. Confidence intervals based on estimators with
approximately equal variance and squared bias will tend to have poor
coverage. This inference problem can be even worse when machine learners are
plugged into a formula for a causal or structural effect. These formulae
often involve averaging over regressor values which reduces variance without
affecting as much the bias. Variance could also potentially also be a
problem but machine learners control that for prediction purposes.
For causal and structural estimators that plug-in regularized machine
learners, the squared bias can shrink slower than the variance, leading to
extremely poor confidence interval coverage and estimators that are not
root-n consistent. Chernozhukov et al. (2017, 2018) give Lasso and random
forest examples respectively and Chernozhukov et al. (2020) shows that Lasso
plug-in estimators are not root-n consistent. Model selection inherent in
machine learners also creates inference problems. Model selection creates
bias from incorrect model choice under local alternatives, making the usual
asymptotic confidence intervals invalid over local alternatives, as shown by
Leeb and Potscher (2008a,b). Estimators of parameters of interest obtained
by plugging in machine learners can inherit this problem, as pointed out by
Belloni, Chernozhukov, and Kato (2015) and Chernozhukov, Hansen, and
Spindler (2015) and shown in Chernozhukov et al. (2020).
To reduce regularization and model selection bias we use a Neyman orthogonal
moment function where there is no first-order effect of the regression on
the expected moment function. The orthogonal moment function is constructed
by adding to an identifying moment the nonparametric influence function of
the regression on the identifying moment function. This construction is
model free, nonparametric, and based on the probability limit of the
regression learner for any distribution, as in Chernozhukov et al. (2016,
2020). As a result the orthogonality property is model free, meaning that
regression learners have no first order effect on the moments for
unrestricted, possibly misspecified, nonparametric distributions.
Consequently the standard errors are robust to misspecification because they
are constructed from the orthogonal moments while ignoring the presence of
the regression learners.
The orthogonal moment function depends on another unknown function $\bar{
\alpha}$ in addition to the regression. We develop a Lasso minimum distance
learner of $\bar{\alpha}$ that is automatic and nonparametric, in the sense
that it depends only on the identifying moment function and not on the form
of $\bar{\alpha}$. The structure of the identifying moment function is used
to approximate $\bar{\alpha}$ as a linear combination of a dictionary (i.e.
basis) of known functions. We use the Lasso learner of $\bar{\alpha}$ and a
regression learner in the orthogonal moment functions to construct an
automatic debiased machine learner (Auto-DML) of parameters of interest. We
introduce debiased machine learning estimators for a wide variety of
effects, including policy effects, average derivatives, bounds on average
equivalent variation, and any other linear function of a regression where
debiased machine learners were not previously available. We also allow for
the identifying moment functions to be nonlinear in regressions. In addition
we give novel estimators of average treatment effects, causal mediation, and
regression decomposition.
We allow any regression learner, including neural nets, random forests,
Lasso, and other high dimensional learners to be used in the orthogonal
moment function. The primary requirement of the regression learner is that
the product of mean-square convergence rates for the learner of $\bar{\alpha}
$ and the regression learner is faster than $n^{-1/2}.$ Under this condition
and a few other regularity conditions we show root-n consistency and
asymptotic normality of the estimator of the parameter of interest. We give
convergence rates for the Lasso learner of $\bar{\alpha}$ and combine them
with existing convergence rates for regressions to verify conditions for
particular estimators. A learner of $\bar{\alpha}$ and large sample theory
is given for parameters that depend nonlinearly on regressions as well as
parameters that are linear in a regression.
The large sample theory in this paper takes the probability limit of the
regression learner and $\bar{\alpha}$ to be fixed. It would be
straightforward to extend the results to allow the regression limit and $
\bar{\alpha}$ to change with sample size. Such a change would allow us to
accommodate sparse specifications where number of nonzero coefficients in
the true regression grows with the sample size but would complicate notation
and detail. We choose to work with a fixed regression for simplicity while
accommodating high dimensional regressions via approximate sparsity.
We give an application to estimating the treatment effect on the treated of
job training from the National Supported Work Demonstration (NSW). For many
large sets of covariates, we find similar estimates based on neural net,
random forest, and Lasso regressions with the automatic bias correction for
each. We also give an application to estimating price elasticities from
scanner panel data while allowing endogeneity of prices. We estimate the
elasticities from Auto-DML of an average derivative that includes many
covariates that account for correlated random effects. We find price
elasticities that are much smaller than cross-section elasticities,
consistent with though larger t than fixed effects elasticities found in
Chernozhukov, Hausman, and Newey (2021). We also find that plug in estimates
are similar to the cross-section elasticity estimates, so that debiasing is
important in this application.
The estimators of parameters of interest use cross-fitting, as in
Chernozhukov et al. (2018), where orthogonal moment functions are averaged
over groups of observations, the regression and $\bar{\alpha}$ learners use
all observations not in the group, and each observation is included in the
average over one group. Cross-fitting removes a source of bias and
eliminates any need for Donsker conditions for the regression learner. Early
work by Bickel (1982), Schick (1986), and Klaassen (1987) used similar
sample splitting ideas.
Auto-DML for a general linear functional of a regression, convergence rates,
and asymptotic normality results for a Dantzig selector of $\bar{\alpha}$
and the regression were given in Chernozhukov, Newey, and Robins (2018).
Chernozhukov, Newey, and Singh (2018) gave Auto-DML for any regression
learner, for nonlinear functions of a regression, and convergence rates for
a Lasso learner of $\bar{\alpha}.$ The current paper is a revised version of
Chernozhukov, Newey, and Singh (2018) with a different title. Chernozhukov,
Newey, and Singh (2019) is a revised version of Chernozhukov, Newey, and
Robins (2018) and is distinguished from the current paper and previous work
in giving and analyzing Auto-DML for local (nonparametric) effects as well
as focusing on the Dantzig selector for $\bar{\alpha}$ and the regression
for global effects. All of these papers make use of model free orthogonal
moment functions for regression learners given in Chernozuhkov et al. (2016)
and the automatic debiasing in Chernozhukov et al. (2020) builds on this
paper. The combined use of cross-fitting and orthogonal moment functions for
debiased machine learning is like Chernozhukov et al. (2018). The Auto-DML
in Chernozhukov, Newey, and Robins (2018), Chernozhukov, Newey, and Singh
(2018), and here innovates by not requiring an explicit formula for the bias
correction that is required in Chernozhukov et al. (2018) and earlier papers.
This work builds upon ideas in classical semi- and nonparametric learning
theory with low-dimensional regressions using traditional smoothing methods
(Van Der Vaart, 1991; Bickel et al., 1993; Newey 1994; Robins and Rotnitzky,
1995; Van der Vaart, 1998), that do not apply to the current
high-dimensional setting. The orthogonal moment functions developed in
Chernozhukov et al. (2016) and used here build on previous work on model
free orthogonal moment functions. Hasminskii and Ibragimov (1979) and Bickel
and Ritov (1988) suggest such estimators for functionals of a density. Newey
(1994) develops such scores for densities and regressions from computation
of the semiparametric efficiency bound for regular functionals. Doubly
robust estimating equations for treatment effects as in Robins, Rotnitzky,
and Zhao (1995) and Robins and Rotnitzky (1995) constitute model based
orthogonal moment functions and have motivated much subsequent work. Newey,
Hsieh, and Robins (1998, 2004) extend model free orthogonal moment functions
to any functional of a density or distribution in a low dimensional setting.
Model free, orthogonal moments for any learner are given and their general
properties derived in Chernozhukov et al. (2016, 2020). We use those model
free, orthogonal moment functions for regressions.
This paper also builds upon and contributes to the literature on modern
orthogonal/debiased estimation and inference, including Zhang and Zhang
(2014), Belloni et al. (2012, 2014a,b), Robins et al. (2013), van der Laan
and Rose (2011), Javanmard and Montanari (2014a,b, 2015), Van de Geer et al.
(2014), Farrell (2015), Ning and Liu (2017), Chernozhukov et al. (2015),
Neykov et al. (2018), Ren et al. (2015), Jankova and Van De Geer (2015,
2016a, 2016b), Bradic and Kolar (2017), Zhu and Bradic (2017a,b). This prior
work is about regression coefficients, treatment effects, and semiparametric
likelihood models. The objects of interest we consider are different than
those analyzed in Cai and Guo (2017). The continuity properties of
functionals we consider provide additional structure that we exploit, namely
the $\bar{\alpha}\,$, an object that is not considered in Cai and Guo
(2017).
Targeted maximum likelihood was developed by Scharfstein, Rotnitzky, Robins
(1999) and Van Der Laan and Rubin (2006). The use of machine learning for
these estimators was proposed by Van der Laan and Rose (2011) and large
sample theory given by Luedtke and Van Der Laan (2016), Toth and van der
Laan (2016), and Zheng et al. (2016). In this paper we give a targeted
version of Auto-DML with automatic debiasing that we refer to as Auto-TML.
This estimator differs from previous ones in the objects we consider and the
use of automatic debiasing in Auto-TML.
Various papers have considered direct estimation of $\bar{\alpha}$ for
treatment effects, where $\bar{\alpha}$ is a Riesz representer that depends
on inverse propensity scores. Our work is the first to present a framework
for direct estimation of the Riesz representer of a broad class of linear
and nonlinear functionals, in a high-dimensional setting, without requiring
strong Donsker class assumptions. The earliest reference of which we know is
Robins et al. (2007), which gives a linear estimator for $\bar{\alpha}$ for
only the average treatment effect. Vermeulen and Vansteelandt (2015) base
parametric propensity score and regression estimators on double robustness
conditions for the average treatment effect. We differ in using a linear
approximation to $\bar{\alpha}$, which is restrictive in a parametric
setting but is general in high dimensional and/or nonparametric settings.
Newey and Robins (2018) present and analyze estimators based on regression
splines, while we present and analyze sparse methods for the
high-dimensional setting. The Lasso minimum distance learner of $\bar{\alpha}
$ given in Chernozhukov, Newey, and Singh (2018) and here is a direct
estimator of the Riesz representer for a broad class of linear and nonlinear
functionals that can be interpreted as being based on orthogonality of the
moment functions. Chernozhukov et al. (2020) extends this learner of $\bar{
\alpha}$ to functions of high dimensional regression quantiles and other
objects.
In independent work on treatment effects Avagyan and Vansteelandt (2017)
give a model assisted estimator based on regularized first order conditions
and Tan (2020) developed a model assisted, multistep method of doubly robust
estimation with Lasso type regression learners having standard errors that
are robust to misspecification of the regression or propensity score.
Smucler, Rotnitzky, and Robins (2019) extended that approach to the linear
functionals of a regression considered in Chernozhukov, Newey, and Singh
(2018). For treatment effects the estimator we give is single step, allows
for any regression learner (e.g. neural nets), is model free, and has
correct standard errors if either or both the regression and the propensity
score are misspecified. Farrell, Liang, and Misra. (2021) gave a neural nets
and model based estimator of the average treatment effect and Wooldridge and
Zhu (2020) give a Lasso based debiased machine learner for panel data with
correlated random effects that depend on high dimensional regressions. Our
results also allow for a neural net regression learner but are model free
with specification robust standard error.
Chernozhukov, Newey, and Robins (2018) gave Auto-DML for linear functionals
using the Dantzig selector. More recently Hirshberg and Wager (2018) gave estimators for linear functionals based on minimax estimation of sample weights that are consistent for realizations of $\bar{\alpha}$ in sample mean square error, rather than a linear approximation to the $\bar{\alpha}$ function, in the low dimensional case, using the same orthogonal moment functions considered here.
The objects considered by Chernozhukov, Newey, and Robins (2018) include
average derivatives. More recently Hirshberg and Wager (2020) gave an
average derivative estimator based on debiasing a Lasso regression learner
of a single index high dimensional regression and Rothenhausler and Yu
(2019) gave an average derivative estimator using debiased Lasso regression.
Singh and Sun (2019) extend the present work to the instrumental variable
setting and present estimators of the local average treatment effect,
average complier characteristics, and complier counter factual
distributions. Previous to the current version of this paper Farbmacher et
al. (2020) gave DML (debiased machine learning) for causal mediation. We
propose an Auto-DML for causal mediation analysis as an example in Section 5.
In summary, contributions of the paper include the construction of DML for a
wide range of interesting policy effects and structural parameters where DML
was not previously available. This construction is based on a Lasso minimum
distance learner of $\bar{\alpha}$ we propose. The debiasing and inference
is model free and robust to misspecification and carried out in a single
step, unlike previous estimators of average treatment effects. For average
treatment and other effects we construct DML for a variety of regression
learners, such as neural nets, random forests, or high dimensional methods.
In Section 2 we describe the objects of interest we consider and associated
orthogonal moment functions. In Section 3 we give the Lasso learner of $\bar{
\alpha},$ the Auto-DML and Auto-TML estimators, and a consistent estimator
of their asymptotic variance. Section 4 derives mean square convergence
rates for the Lasso learner of $\bar{\alpha}$ and conditions for root-n
consistency and asymptotic normality of Auto-DML and Auto-TML including
primitive conditions in examples. Section 5 gives Auto-DML for nonlinear
functionals of multiple regressions and as an example develops Auto-DML for
causal mediation analysis. Section 6 gives Auto-DML for regression
decomposition and estimates the average treatment on the treated for the NSW
experiment. Section 7 gives Auto-DML estimates of price elasticities that
allow for correlated random effects in scanner panel data. Section 8 offers
some conclusions and possible extensions.
\section{Average Linear Effects and Orthogonal Moment Functions}
For expositional purposes, in this Section we first consider parameters that
depend linearly on a single conditional expectation. To describe such an
object, let $W$ denote a data observation, and consider a subvector $
(Y,X^{\prime})^{\prime}$ where $Y$ is a scalar outcome with finite second
moment and $X$ is a covariate vector. Denote the conditional expectation of $
Y$ given $X\in\mathcal{X}$ as
\begin{equation*}
\gamma_{0}(x)=\mathrm{E}[Y|X=x].
\end{equation*}
Let $m(w,\gamma)$ denote a function of the function $\gamma$ (i.e. a
functional of $\gamma),$ where $\gamma$ denotes a possible conditional
expectation function $\gamma:\mathcal{X}\longrightarrow\mathbb{R}$, that
depends on a data observation $w$ and is linear in $\gamma.$ We will
consider effects of the form
\begin{equation*}
\theta_{0}=\mathrm{E}[m(W,\gamma_{0})].
\end{equation*}
The parameter of interest $\theta_{0}$ is an expectation of some known
formula $m(W,\gamma)$ of a data observation $W$ and a regression $\gamma.$
We also give results in later Sections for important parameters having more
general forms. In Section 5 we allow $m(W,\gamma)$ to be nonlinear in
multiple regressions and propose an estimator of causal effects with
mediation. In Section 6 we give estimators of regression decompositions and
their properties. These important examples extend the framework of this
Section to parameters that are nonlinear in multiple regressions
Several important examples of linear effects are:
\bigskip
\textsc{Example 1:} (Average Policy Effect). An average effect of a counter
factual shift in the distribution of regressors from a known $F_{0}$ to
another known $F_{1}$, when $\gamma _{0}$ does not vary with the
distribution of $X$, is
\begin{equation*}
\theta _{0}=\int \gamma _{0}(x)d\mu (x);\text{ }\mu (x)=F_{1}(x)-F_{0}(x).
\end{equation*}
Here $m(w,\gamma )=\int \gamma (x)d\mu (x)$ which does not depend on $w.$
This policy effect builds on but is different than Stock (1989) in comparing
averages over two known distributions rather than the empirical distribution.
\bigskip
\textsc{Example 2:} (Weighted Average Derivative). Here $X=(D,Z)$ for a
continuously distributed random variable $D,$ $\gamma_{0}(x)=
\gamma_{0}(d,z), $ $\omega(d)$ is a pdf, and
\begin{equation*}
\theta_{0}=\mathrm{E} \left[ \int\omega(u)\frac{\partial\gamma_{0}(u,Z)}{
\partial d}du\right] =\mathrm{E} \left[ \int S(u)\gamma_{0}(u,Z)\omega(u)du
\right] =\mathrm{E}[S(U)\gamma_{0}(U,Z)],
\end{equation*}
where $S(u)=-\omega(u)^{-1}\partial\omega(u)/\partial u$ is the negative
score for the pdf $\omega(u),$ the second equality follows by integration by
parts, and $U$ is a random variable that is independent of $Z$ with pdf $
\omega(u).$ This $U$ could be thought of as one simulation draw from the pdf
$\omega(u).$ Here $m(w,\gamma)=S(u)\gamma(u,x)$ where $W$ includes $U.$
This $\theta _{0}$ can be interpreted as an average treatment effect on $Y$
of a continuous treatment $D$ in a model where $Y=Y(D)$ for a potential
outcome stochastic process $Y(d)$ that is independent of $D$ conditional on
covariates $Z.$ By conditional independence
\begin{equation*}
\mathrm{E}[\gamma _{0}(u,Z)]=\int \mathrm{E}[Y(D)|D=u,Z=z]F_{Z}(dz)=\int
\mathrm{E}[Y(u)|Z=z]F_{Z}(dz)=\mathrm{E}[Y(u)],
\end{equation*}
for $\omega (u)>0$ assuming that the joint pdf of $(D,Z)$ is positive where $
\omega (D)>0$, as in Chamberlain (1984), Wooldridge (2002), and Blundell and
Powell (2004). The $\mathrm{E}[Y(u)]$ is the average outcome at $D=u$ and is
sometimes referred to as the average structural function. Assuming that we
can interchange the order of differentiation and integration,
\begin{equation*}
\theta _{0}=\int \omega (u)\frac{\partial \mathrm{E}[\gamma _{0}(u,Z)]}{
\partial u}du=\int \frac{\partial \mathrm{E}[Y(u)]}{\partial u}\omega
(u)du=\int \mathrm{E}\left[ \frac{\partial Y(u)}{\partial u}\right] \omega
(u)du,
\end{equation*}
similarly to Imbens and Newey (2009) and Rothenh{\"{a}}usler and Yu (2019),
which build on but are different than Powell, Stock, and Stoker (1989).
Regarding $\mathrm{E}[\partial Y(u)/\partial u]$ as the average treatment
effect at $u$ we see that $\theta _{0}$ is a weighted average treatment
effect. Alternatively, $\theta _{0}$ can be regarded as an average
derivative of the average structural function. The averaging over a known
pdf $\omega (u)$ helps fulfill regularity conditions for the Auto-DML
developed here that can be used to estimate $\theta _{0}$ for high
dimensional covariates $Z.$
\bigskip
\textsc{Example 3: }(Average Treatment Effect). In this example $X=(D,Z)$
and $\gamma_{0}(x)=\gamma_{0}(d,z)$, where $D\in\{0,1\}$ is the treatment
indicator and $Z$ are covariates. The object of interest is
\begin{equation*}
\theta_{0}=\mathrm{E}[\gamma_{0}(1,Z)-\gamma_{0}(0,Z)].
\end{equation*}
If potential outcomes are mean independent of treatment $D$ conditional on
covariates $Z$, then $\theta_{0}$ is the average treatment effect (Rosenbaum
and Rubin, 1983). Here $m(w,\gamma)=\gamma(1,z)-\gamma(0,z).$
\bigskip
\textsc{Example 4:} (Average Equivalent Variation Bound). An economic
example is a bound on average equivalent variation for heterogenous demand.
Here $Y$ is the share of income spent on a commodity and $X=(P_{1},Z),$
where $P_{1}$ is the price of the commodity and $Z$ includes income $Z_{1}$,
prices of other goods, and other observable variables affecting utility. Let
$\check{p}_{1}<\bar{p}_{1}$ be lower and upper prices over which the price
of the commodity can change, $\kappa$ a bound on the income effect, $
\omega(z)$ some weight function, and $U$ a random variable that is uniformly
distributed over $(\check{p}_{1},\bar{p}_{1})$ and independent of $(Y,X).$ $
U $ can be thought of as one simulation draw from a uniform distribution on $
(\check{p}_{1},\bar{p}_{1}).$ The object of interest is
\begin{equation*}
\theta_{0}=\mathrm{E} \left[ \Lambda(U,Z)\gamma_{0}(U,Z)\right] ,\text{ }
\Lambda(u,z)=\omega(z)1(\check{p}_{1}<u<\bar{p}_{1})(\bar{p}_{1}-\check{p}
_{1})\frac{z_{1}}{u}\exp(-\kappa\lbrack u-\check{p}_{1}]).
\end{equation*}
If individual heterogeneity in consumer preferences is independent of $X$
and $\kappa$ is a lower (upper) bound on the derivative of consumption with
respect to income for all individuals, then $\theta_{0}$ is an upper (lower)
bound on the weighted average over consumers of equivalent variation for a
change in the price of the first good from $\check{p}_{1}$ to $\bar{p}_{1}$;
see Hausman and Newey (2016). Here $m(w,\gamma)=\Lambda(u,z)\gamma(u,z),$
where $W$ includes $U.$
\bigskip
We focus on $m(w,\gamma)$ where there exists a function $\alpha_{0}(X)$ with
$\mathrm{E}[\alpha_{0}(X)^{2}]<\infty$ and
\begin{equation}
{\mathrm{E}}[m(W,\gamma)]={\mathrm{E}}[\alpha_{0}(X)\gamma(X)]\text{ \ for
all }\gamma\text{ such that }{\mathrm{E}}[\gamma(X)^{2}]<\infty.
\label{Riesz rep}
\end{equation}
By the Riesz representation theorem, existence of such a $\alpha_{0}(X)$ is
equivalent to $\mathrm{E}[m(W,\gamma)]$ being a mean-square continuous
functional of $\gamma,$ i.e. $\mathrm{E}[m(W,\gamma)]\leq C\left\Vert
\gamma\right\Vert $ for all $\gamma$, where $\left\Vert \gamma\right\Vert =
\sqrt{\mathrm{E}[\gamma(X)^{2}]}$ and $C>0.$ We will refer to this $
\alpha_{0}(X)$ as the Riesz representer (Rr). Existence of the Rr is
equivalent to the semiparametric variance bound for $\theta_{0}$ being
finite, as stated in Newey (1994) and shown in Hirshberg and Wager (2018)
for conditional expectations and in Chernozhukov, Newey, and Singh (2019)
more generally for least squares projections. Thus, in assuming existence of
$\alpha_{0}(X)$ we are just assuming that $\theta_{0}$ has a finite
semiparametric variance bound.
Each of Examples 1-4 has such a Rr. Let $f(x)$ denote the pdf of $X$ in
Example 1, $f(d|z)$ the pdf of $D$ conditional on $Z$ in Example 2, $\pi
_{0}(z)=\Pr(D=1|Z=z)$ the propensity score in Example 3, and $f(p_{1}|z)$
the pdf of $P_{1}$ conditional on $Z$ in Example 4. Table~\ref{tab:RR} summarizes the
functional $m(w,\gamma)$ and the Rr in each of the examples:
\begin{table}[ptb]
\centering
\begin{tabular}{c|c|c}
\hline\hline
Effect & $m(W,\gamma)$ & Riesz Representer \\ \hline
Policy Effect & $\int\gamma(x)[f_{1}(x)-f_{0}(x)]dx$ & $
f(X)^{-1}[f_{1}(X)-f_{0}(X)]$ \\
Weighted Average Derivative & $S(U)\gamma(U,Z)$ & $f(D|Z)^{-1}\omega(D)S(D)$
\\
Average Treatment Effect & $\gamma(1,Z)-\gamma(0,Z)$ & $\pi_{0}(Z)^{-1}D-(1-
\pi_{0}(Z))^{-1}(1-D)$ \\
Equivalent Variation Bound & $\Lambda(U,Z)\gamma(U,Z)$ & $(\bar{p}_{1}-
\check{p}_{1})^{-1}f(P_{1}|Z)^{-1}\Lambda(P_{1},Z)$ \\
\hline\hline
\end{tabular}
\caption{$m$ and Rr for Examples 1-4}
\label{tab:RR}
\end{table}
Equation (\ref{Riesz rep}) follows in Example 1 by multiplying and dividing
by $f(x)$ inside the integral, in Example 2 by integration and multiplying
and dividing by $f(d|z)$, in Example 3 in a standard way for average
treatment effects, and in Example 4 by multiplying and dividing by $
f(p_{1}|z)$. For $\mathrm{E}[\alpha_{0}(X)^{2}]<\infty$ to hold the
denominator must not be too small relative to the numerator in each $
\alpha_{0}(X)$, on average. For instance Example 3 must have $\mathrm{E}
[\{\pi_{0}(Z)(1-\pi_{0}^{{}}(Z))\}^{-1}]<\infty.$
Equation (\ref{Riesz rep}) implies that the effect of interest can be
represented in three different ways, as
\begin{equation*}
\theta_{0}=\mathrm{E}[m(W,\gamma_{0})]=\mathrm{E}[\alpha_{0}(X)
\gamma_{0}(X)]=\mathrm{E}[\alpha _{0}(X)Y],
\end{equation*}
where the last equality follows by iterated expectations. Any of these three
expressions could be used to estimate $\theta_{0}$. We could estimate $
\theta_{0}$ from the first expression using a learner (estimator) of $
\gamma_{0}$. We could also estimate $\theta_{0}$ from the last expression
using a learner of $\alpha_{0}(X).$ In addition we could use learners of
both $\gamma_{0}$ and $\alpha_{0}$ to estimate $\theta_{0}$ from the middle
expression. We focus here on using a learner of $\gamma_{0}$, though $
\alpha_{0}$ will be important for the bias correction to follow.
We rely on a regression learner (estimator) $\hat{\gamma}$ of $\gamma_{0}$
to estimate $\theta_{0}.$ The $\hat{\gamma}$ can be any of a variety of
machine learners including neural nets, random forests, Lasso, and other
high dimensional methods. All we require is that $\hat{\gamma}$ converge in
mean square at a sufficiently fast rate, as specified in Section 4.
Whatever the choice of $\hat{\gamma},$ estimating $\theta_{0}$ by plugging $
\hat{\gamma}$ into $m(W,\gamma)$ and averaging over observations on $W$ can
lead to large biases when $\hat{\gamma}$ involves regularization and/or
model selection, as discussed in the Introduction. For that reason we use an
orthogonal moment function for $\theta_{0}$, where the regression learner $
\hat{\gamma}$ has no first-order effect on the moments. We follow
Chernozhukov et al. (2016, 2020) in basing the orthogonal moment function on
the probability limit (plim) $\gamma(F)$ of $\hat{\gamma}$ when one
observation $W$ has CDF $F,$ where $F$ is unrestricted except for regularity
conditions. Here $\gamma(F)$ can be thought of as the plim of $\hat{\gamma}$
under general misspecification, where $\gamma(F)$ need not be the
conditional expectation $\mathrm{E}_{F}[Y|X]$.
The plim $\gamma(F)$ of $\hat{\gamma}$ depends on the learner. For example
Lasso, the Dantzig selector, boosting, and other high dimensional methods
are based on a sequence of potential regressors $X=(X_{1},X_{2},...)$. These
learners have the form
\begin{equation*}
\hat{\gamma}(x)=\sum_{j=1}^{\infty}\hat{\beta}_{j}x_{j}\text{, }\hat{\beta }
_{j^{\prime}}\neq0\text{ for a finite number of }j^{\prime}\text{,}
\end{equation*}
where $x=(x_{1},x_{2},...)$ denotes a possible realization of $X$. Because
each $\hat{\gamma}(X)$ is a linear combination of $X=(X_{1},X_{2},...)$ the
plim $\gamma(F)$ of $\hat{\gamma}$ will also be a linear combination of $X$,
or at least will be approximated by such a linear combination. Define $
\Gamma $ to be the mean square closure of the set of finite linear
combinations of $X$, i.e. $\Gamma$ is the set of $\gamma(X)$ such that $
\mathrm{E}[\gamma(X)^{2}]<\infty$ and for every $\varepsilon>0$ there exists
$(\beta_{j}^{\varepsilon})_{j=1}^{\infty}$ such that $\beta_{j^{\prime}}^{
\varepsilon}\neq0$ for a finite number of $j^{\prime}$ and $\mathrm{E}
[\{\gamma(X)-\sum_{j=1}^{\infty}\beta
_{j}^{\varepsilon}X_{j}\}^{2}]<\varepsilon.$ It will be the case that $
\gamma(F)\in\Gamma.$ Because Lasso and other high dimensional methods are
being used for least squares prediction of $Y$ it will also be the case that
\begin{equation}
\gamma(F)=\arg\min_{\gamma\in\Gamma}\mathrm{E}_{F}[\{Y-\gamma(X)\}^{2}],
\label{BLP}
\end{equation}
This $\gamma(F)$ minimizes population least squares criteria over the (mean
square closure of) linear combinations of $X,$ i.e. it is the best linear
predictor of $Y$ by linear combinations of $X.$ Here $\gamma(F)$ is the
infinite dimensional linear regression that is nonparametrically estimated
by Lasso and other high dimensional methods.
Neural nets and random forests may have a different $\gamma(F)$. A neural
net or random forest is often a nonparametric regression estimator for a
finite (but high) dimensional $X$. In that case
\begin{equation*}
\gamma(F)=\mathrm{E}_{F}[Y|X],
\end{equation*}
which satisfies equation (\ref{BLP}) when $\Gamma$ is the set of all
(measurable) functions of $X$ with finite second moment. The plim of Lasso
and other high dimensional methods will also be this $\gamma(F)$ if $
X=(X_{1},X_{2},...)$ can approximate any function of a fixed set of
regressors, but otherwise will not. A third type of learner $\hat{\gamma}$
is one that imposes additivity restrictions on $\hat{\gamma}$, such as $\hat{
\gamma}(X)=\hat{\gamma}_{1}(X_{1})+\hat{\gamma}_{2}(X_{2})$, allowing for
nonparametric learners $\hat{\gamma}_{1}(X_{1})$ and $\hat{\gamma}
_{2}(X_{2}).$ In that case $\gamma(F)$ will be satisfy equation (\ref{BLP})
where $\Gamma$ is the mean square closure of functions that are additive in $
X_{1}$ and $X_{2}$.
We use the orthogonal moment function from Chernozhukov et al. (2016, 2020)
for a regression learner $\hat{\gamma}$ having plim $\gamma(F)$ satisfying
equation (\ref{BLP}) for any linear, closed $\Gamma.$ The orthogonal moment
function is constructed by adding to the identifying moment function $
m(w,\gamma)-\theta$ the nonparametric influence function of of $\mathrm{E}
[m(W,\gamma (F))].$ As shown in Newey (1994) the nonparametric influence
function of $\mathrm{E}[m(W,\gamma(F))]$ is
\begin{equation*}
\bar{\alpha}(X)[Y-\bar{\gamma}(X)],
\end{equation*}
where $\bar{\gamma}(X)$ is the solution to equation (\ref{BLP}) for $F=F_{0}$
and $\bar{\alpha}\in\Gamma$ satisfies $\mathrm{E}[m(W,\gamma)]=\mathrm{E}[
\bar{\alpha}(X)\gamma(X)]$ for all $\gamma\in\Gamma.$ As in Chernozhukov,
Newey, and Singh (2019),
\begin{equation}
\bar{\alpha}=\arg\min_{\alpha\in\Gamma}\mathrm{E}[\{\alpha_{0}(X)-\alpha(X)
\}^{2}]. \label{BLP alpha}
\end{equation}
This $\bar{\alpha}$ can be thought of as the Riesz representer for the
linear functional $\mathrm{E}[m(W,\gamma)]$ with domain $\Gamma.$ Evaluating
the nonparametric influence function at possible values $\gamma$ and $\alpha$
of $\bar{\gamma}$ and $\bar{\alpha}$ and adding it to the the identifying
moment function gives the orthogonal moment function
\begin{equation}
\psi(w,\theta,\gamma,\alpha)=m(w,\gamma)-\theta+\alpha(x)[y-\gamma(x)].
\label{Debiased mom}
\end{equation}
The moment function $\psi(w,\theta,\gamma,\alpha)$ depends on a possible
value $\alpha$ of the unknown function $\bar{\alpha}$ as well as a possible
value $\gamma$ of the plim $\bar{\gamma}$ of the regression learner. A
learner $\hat{\alpha}$ of $\bar{\alpha}$ is needed to use this orthogonal
moment function to estimate $\theta_{0}.$ In Section 3 we will describe how
to construct $\hat{\alpha}.$ In Chernozhukov et al. (2016, 2020) $\psi
(w,\theta,\gamma,\alpha)$ is shown to be orthogonal without being specific
about the form of $\hat{\alpha}.$ For exposition we repeat that
demonstration here. Consider any $\gamma,\alpha\in\Gamma$, representing
possible realizations of learners $\hat{\gamma}$ and $\hat{\alpha}$ that are
in $\Gamma.$ The well known necessary and sufficient conditions for equation
(\ref{BLP}) with $F=F_{0}$ are that $\mathrm{E}[\alpha(X)\{Y-\bar{\gamma}
(X)\}]=0$ for all $\alpha\in\Gamma.$ Therefore
\begin{align}
\mathrm{E}[\psi(W,\theta,\gamma,\alpha)-\psi(W,\theta,\bar{\gamma},\bar{
\alpha})] & =\mathrm{E}[m(W,\gamma)]-\mathrm{E}[m(W,\bar{\gamma})]+\mathrm{E}
[\alpha(X)\{Y-\gamma (X)\}] \label{2nd order} \\
& =\mathrm{E}[\alpha_{0}(X)\{\gamma(X)-\bar{\gamma}(X)\}]+\mathrm{E}
[\alpha(X)\{Y-\gamma (X)\}] \notag \\
& =\mathrm{E}[\bar{\alpha}(X)\{\gamma(X)-\bar{\gamma}(X)\}]+\mathrm{E}
[\alpha(X)\{\bar{\gamma }(X)-\gamma(X)\}] \notag \\
& =-\mathrm{E}[\{\alpha(X)-\bar{\alpha}(X)\}\{\gamma(X)-\bar{\gamma}(X)\}],
\notag
\end{align}
where the second equality follows by equation (\ref{Riesz rep}) and the
third equality by the necessary and sufficient condition for equation (\ref
{BLP alpha}) that $\mathrm{E}[\{\alpha_{0}(X)-\bar{\alpha}(X)\}\gamma(X)]=0$
for all $\gamma\in\Gamma.$ Here we see that $\psi(w,\theta,\gamma,\bar{\alpha
})$ "partials out" $\gamma$ in the sense that
\begin{equation*}
\mathrm{E}[m(W,\gamma)+\bar{\alpha}(X)\{Y-\gamma(X)\}]=\mathrm{E}[m(W,\bar{
\gamma})]
\end{equation*}
does not depend on $\gamma$. Also equation (\ref{2nd order}) gives an
explicit formula showing that the effect of $\gamma$ and $\alpha$ on $
\mathrm{E}[\psi (W,\theta,\gamma,\alpha)]$ is second order and hence $
\psi(W,\theta ,\gamma,\alpha)$ is orthogonal.
The orthogonality property of $\psi (W,\theta ,\gamma ,\alpha )$ only
depends on $\gamma ,$ $\alpha \in \Gamma $ and $\bar{\gamma}$ satisfying
equation (\ref{BLP}). In particular orthogonality does not depend on either $
\bar{\gamma}$ being $\mathrm{E}[Y|X]$ or on $\bar{\alpha}=\alpha _{0}.$ In
this sense orthogonality of $\psi (W,\theta ,\gamma ,\alpha )$ is model
free, i.e. nonparametric. Consequently the estimator of $\theta $ will be
asymptotically normal and standard errors consistent even if either $\bar{
\gamma}\neq \gamma _{0}$ or $\bar{\alpha}\neq \alpha _{0}$ or both, which is
possible when neither $\gamma _{0}(X)=\mathrm{E}[Y|X]$ nor $\alpha _{0}(X)$
satisfying equation (\ref{Riesz rep}) is an element of $\Gamma .$ This
robustness of the standard errors results from the orthogonality of the
moments only depending on the $\bar{\gamma}$ limit of the regression
estimator, so that the sample average of the estimated orthogonal moment
function will be asymptotically equivalent to the sample average at the
truth, without any model assumptions.
The orthogonal moment function could also be viewed as the efficient
influence function of $\mathrm{E}[m(W,\bar{\gamma})]$ which clarifies that
the Auto-DML is an efficient semiparametric estimator of $\mathrm{E}[m(W,
\bar{\gamma})]$. Viewing $\psi(w,\theta,\gamma,\alpha)$ in this way is not
useful for debiasing because the results of Chernozhukov et. al. (2016,
2020) already imply model free orthogonality.
The moment function $\psi(w,\theta,\gamma,\alpha)$ is doubly robust for
estimation of the true parameter $\theta_{0}.$ Evaluating at $\theta_{0},
\bar{\gamma},\bar{\alpha}$ and taking the expectation gives
\begin{align}
\mathrm{E}[\psi(W,\theta_{0},\bar{\gamma},\bar{\alpha})] & =\mathrm{E}[m(W,
\bar{\gamma })]-\theta_{0}+\mathrm{E}[\bar{\alpha}(X)\{Y-\bar{\gamma}(X)\}]
\label{Double rob} \\
& =\mathrm{E}[\alpha_{0}(X)\{\bar{\gamma}(X)-\gamma_{0}(X)\}]+\mathrm{E}[
\bar{\alpha }(X)\{\gamma_{0}(X)-\bar{\gamma}(X)\}] \notag \\
& =-\mathrm{E}[\{\bar{\alpha}(X)-\alpha_{0}(X)\}\{\bar{\gamma}(X)-\gamma
_{0}(X)\}], \notag
\end{align}
which is zero for $\bar{\gamma}=\gamma_{0}$ or $\bar{\alpha}=\alpha_{0}.$
Thus $\mathrm{E}[\psi(W,\theta_{0},\bar{\gamma},\bar{\alpha})]=0$, so that
the orthogonal moment condition identifies $\theta_{0},$ when either $\bar{
\gamma}(X)=\mathrm{E}[Y|X]$ or $\alpha_{0}(X)\in\Gamma.$ These conditions
both hold when the regression learner is nonparametric so that $\Gamma$ is
the set of all functions of $X$ with finite second moment. For high
dimensional regressions where $\Gamma$ is the closed linear span of $
X=(X_{1},X_{2},...)$ the plim of the learner $\hat{\gamma}$ may not be $
\mathrm{E}[Y|X]$ but the orthogonal moment function still identifies $
\theta_{0}$ when $\alpha_{0}(X)\in\Gamma.$ That is, $\theta_{0}$ is
identified when $\alpha_{0}(X)$ can be approximated arbitrarily well in mean
square by a linear combination of $X.$ This robustness condition can be
interpreted in each of Examples 1-4:
\bigskip\
\textsc{Example 1:} For high dimensional $\hat{\gamma},$ where $\Gamma$ is
the mean square closure of linear combinations of $X,$ $\mathrm{E}
[\psi(W,\theta_{0},\bar{\gamma},\bar{\alpha})]=0$ even when $\bar{\gamma}
(X)\neq \mathrm{E}[Y|X]$ if $\alpha_{0}(X)=[f_{1}(X)-f_{0}(X)]/f(X)\in
\Gamma. $
\bigskip
\textsc{Example 2:} For high dimensional $\hat{\gamma},$ where $\Gamma$ is
the mean square closure of linear combinations of $X,$ $\mathrm{E}
[\psi(W,\theta_{0},\bar{\gamma},\bar{\alpha})]=0$ even when $\bar{\gamma}
(X)\neq \mathrm{E}[Y|X]$ if $\alpha_{0}(X)=f(D|Z)^{-1}\omega(D)S(D)\in
\Gamma. $
\bigskip
\textsc{Example 3: }For the average treatment effect where $\Gamma$ is
nonparametric, so that $\bar{\gamma}(X)=\mathrm{E}[Y|X]$ and $\bar{\alpha}
(X)=\alpha_{0}(X),$ the orthogonal moment function in equation (\ref
{Debiased mom}) corresponds to the seminal doubly robust moment function of
Robins, Rotnitzky, and Zhao (1995). When $\hat{\gamma}$ is high dimensional,
with say $X=(DZ,(1-D)\tilde{Z})$ for sequences $Z=(Z_{1},Z_{2},...)$ and $
\tilde{Z}=(\tilde{Z}_{1},\tilde{Z}_{2},...)$, with each $\tilde{Z}_{j}$ a
function of $Z,$ the orthogonal moment function is
\begin{equation*}
\psi(W,\theta,\bar{\gamma},\bar{\alpha})=\bar{\gamma}(1,Z)-\bar{\gamma }
(1,0)-\theta+\bar{\alpha}(X)[Y-\bar{\gamma}(X)].
\end{equation*}
This orthogonal moment function is different than those previously
considered in $\bar{\alpha}(X)$ being the projection of $\alpha_{0}(X)$ on $
\Gamma$ rather than $\alpha_{0}(X)$. Here $\mathrm{E}[\psi(W,\theta_{0},\bar{
\gamma},\bar{\alpha})]=0$ if linear combinations of $Z$ and $\tilde{Z}$ can
approximate abitrarily well $\pi_{0}(Z)^{-1}$ and $[1-\pi_{0}(Z)]^{-1}$
respectively, even when $\bar{\gamma}(X)\neq \mathrm{E}[Y|X].$
\bigskip
For brevity we omit further discussion of Example 4 from the paper and refer
the interested reader to Chernozhukov, Hausman, and Newey (2021).
\bigskip
\section{Estimation}
To estimate (learn) $\theta_{0}$ we use cross-fitting where the orthogonal
moment function $\psi(w,\gamma,\alpha,\theta)$ is averaged over observations
different than used to estimate $\bar{\gamma}$ and $\bar{\alpha}.$ We assume
that the data $W_{i},$ $(i=1,...,n)$ are i.i.d.. Let $I_{\ell},$ $
(\ell=1,...,L)$, be a partition of the observation index set $\{1,...,n\}$
into $L$ distinct subsets of about equal size. In practice $L=5$ (5-fold) or
$L=10$ (10-fold) cross-fitting is often used. Let $\hat{\gamma}_{\ell}$ and $
\hat{\alpha}_{\ell}$ be estimators constructed from the observations that
are \textit{not} in $I_{\ell}.$ We construct the estimator $\hat{\theta}$ by
setting the sample average of $\psi(W_{i},\theta,\hat{\gamma}_{\ell},\hat{
\alpha}_{\ell})$ to zero and solving for $\theta.$ This $\hat{\theta}$ and
an associated asymptotic variance estimator $\hat{V}$ have explicit forms
\begin{align}
\hat{\theta} & =\frac{1}{n}\sum_{\ell=1}^{L}\sum_{i\in I_{\ell}}\{m(W_{i},
\hat{\gamma}_{\ell})+\hat{\alpha}_{\ell}(X_{i})[Y_{i}-\hat{\gamma }
_{\ell}(X_{i})]\}, \label{Estimator} \\
\hat{V} & =\frac{1}{n}\sum_{\ell=1}^{L}\sum_{i\in I_{\ell}}\hat{\psi}
_{i\ell}^{2},\text{ }\hat{\psi}_{i\ell}=m(W_{i},\hat{\gamma}_{\ell})-\hat{
\theta}+\hat{\alpha}_{\ell}(X_{i})[Y_{i}-\hat{\gamma}_{\ell}(X_{i})], \notag
\end{align}
Any regression learner $\hat{\gamma}_{\ell }$ can be used here as long as
its mean-square convergence rate is a power of $1/n,$ as assumed in Section
4. Such a convergence rate is available for neural nets (Chen and White,
1999, Schmidt-Heiber, 2020, Farrell, Liang, and Misra, 2021), random forests
(Syrgkanis and Zampetakis, 2020), Lasso (Bickel, Ritov, and Tsybakov, 2009),
boosting (Luo and Spindler, 2016), and other high dimensional methods. As a
result any of these regression learners can be used to construct an Auto-DML
$\hat{\theta}$ from equation (\ref{Estimator}), in conjunction with a
learner $\hat{\alpha}_{\ell }$ of $\bar{\alpha}.$
The correctness of $\hat{V}$ relies on consistency of the regression learner
$\hat{\gamma}_{\ell }.$ It would be interesting to investigate whether the
finite sample approximation could be improved by using a variance estimator
that allowed $\hat{\gamma}_{\ell }$ to not be consistent because the
dimension of the regression grows as fast as the sample size, e.g. as in
Cattaneo, Jansson, and Newey (2018).
An alternative estimator of $\theta _{0}$ can be constructed that extends
the targeted maximum likelihood approach of Scharfstein, Rotnitzky, and J.M.
Robins (1999) and van der Laan and Rubin (2006) to the objects we consider.
This Auto-TML estimator is a plug-in estimator based on a regression learner
that has been debiased in a direction specific to the object of interest.
This estimator is given by
\begin{equation}
\tilde{\theta}=\frac{1}{n}\sum_{\ell =1}^{L}\sum_{i\in I_{\ell }}m(W_{i},
\tilde{\gamma}_{\ell }),\text{ }\tilde{\gamma}_{\ell }(x)=\hat{\gamma}_{\ell
}(x)+\frac{\sum_{i\in I_{\ell }}\hat{\alpha}_{\ell }(X_{i})[Y_{i}-\hat{\gamma
}_{\ell }(X_{i})]}{\sum_{i\in I_{\ell }}\hat{\alpha}_{\ell }(X_{i})^{2}}\hat{
\alpha}_{\ell }(x). \label{TMLE}
\end{equation}
As with other targeted estimators the plug-in form of Auto-TML allows
imposition of constraints through $m(W,\gamma )$. In Section 4 we show that
this estimator is asymptotically equivalent to $\hat{\theta}$.
To describe $\hat{\alpha}_{\ell}$ let $b(x)=(b_{1}(x),...,b_{p}(x))$ be a $
p\times1$ dictionary of functions of $x,$ where $p$ can be large, with each $
b_{j}(x)$ standardized to have mean $0$ and standard deviation $1,$ to be
further discussed in this Section. For convenience we ignore dependence of $
b(x)$ on the data in the notation. The learner $\hat{\alpha}_{\ell}$ given
here is
\begin{align}
\hat{\alpha}_{\ell}(x) & =\frac{1}{n-n_{\ell}}\sum_{i\notin
I_{\ell}}m(W_{i},1)+b(x)^{\prime}\hat{\rho}_{\ell},\text{ }\hat{\rho}
_{\ell}=\arg \min_{\rho}\{-2\hat{M}_{\ell}^{\prime}\rho+\rho^{\prime}\hat{G}
_{\ell}\rho+2r\sum_{j=1}^{J}\left\vert \rho_{j}\right\vert \},
\label{CFLasso} \\
\hat{M}_{\ell} & =\frac{1}{n-n_{\ell}}\sum_{i\notin I_{\ell}}m(W_{i},b),
\text{ }\hat{G}_{\ell}=\frac{1}{n-n_{\ell}}\sum_{i\notin
I_{\ell}}b(X_{i})b(X_{i})^{\prime}, \notag
\end{align}
where $n_{\ell}$ is the number of observations in $I_{\ell}$ and $r>0$ is a
positive scalar. This $\hat{\alpha}_{\ell}$ is used in equation (\ref
{Estimator}) to construct $\hat{\theta}$ and $\hat{V}$.
To explain and motivate $\hat{\alpha}_{\ell}$ it is notationally convenient
to drop the $\ell$ subscript, with the understanding that $\hat{\alpha}
_{\ell}$ is computed using only observations not in $I_{\ell}$ for each $
\ell,$ as in equation (\ref{CFLasso}). It is also notationally convenient to
drop the $0$ mean normalization of $b(x)$ and consider $\hat{\alpha}$ having
the form
\begin{equation}
\hat{\alpha}(x)=b(x)^{\prime}\hat{\rho}\text{,} \label{Riesz est}
\end{equation}
where $\hat{\rho}$ is a vector of estimated coefficients.
The $\hat{\alpha}$ depends on the choice of dictionary $b(x)$ and penalty
degree $r.$ For the dictionary we require that each $b_{j}(x)$ belongs to
the set $\Gamma$ of possible plims of $\hat{\gamma}(x)$ discussed in Section
2 and that linear combinations of the dictionary "span" $\Gamma.$
\bigskip
\textsc{Assumption 1:} $b(x)=(b_{1}(x),...,b_{p}(x))^{\prime}$ where i) $
b_{j}\in\Gamma$ for all $j$ and ii) for any $\alpha\in\Gamma$ and $
\varepsilon>0$ there is $p$ and $\rho\in
\mathbb{R}
^{p}$ such that $\mathrm{E}[\{\alpha(X)-b(X)^{\prime}\rho\}^{2}]<
\varepsilon. $
\bigskip
One key feature of this condition is that each $b_{j}\in\Gamma.$ This
feature allows us to use $m(w,\gamma)$ to construct $\hat{\alpha}$ and will
guarantee that $\hat{\alpha}\in\Gamma$, as required for the orthogonality
shown in equation (\ref{2nd order}). Another key feature is that linear
combinations of $b(x)$ can approximate anything that belongs to $\Gamma.$
This feature will lead to $\hat{\alpha}$ estimating $\bar{\alpha}.$ The link
imposed by Assumption 1, between the regression learner $\hat{\gamma}$ and
the dictionary $b(x)$ used to construct $\hat{\alpha},$ is important for the
orthogonality property of $\psi(w,\gamma,\alpha,\theta)$ and hence for $\hat{
\theta}$ to be asymptotically normal and $\hat{V}$ to be a consistent
estimator of the asymptotic variance under general misspecification.
Assumption 1 requires that linear combinations of $b(x)$ must be able to
approximate any $\gamma$ in the set of possible plims of $\hat{\gamma}$ and
that each $b_{j}$ must be a possible plim of $\hat{\gamma}$. For Lasso and
other high dimensional regression learners where $X=(X_{1},X_{2},...)$
Assumption 1 will be satisfied for
\begin{equation}
b(x)=(x_{1},...,x_{p})^{\prime}. \label{high dimensional dict}
\end{equation}
Evidently each element $b_{j}(X)=X_{j}$ is an element of $\Gamma$ and the
spanning condition is satisfied because any linear combination of $X$ with a
finite number of nonzero coefficients will also be a linear combination of $
b(x)$ for $p$ large.
We emphasize that $b(X)$ is required to approximate only the projection $
\bar{\alpha}(X)$ and not $\alpha _{0}(X).$ For instance, in the average
treatment effect example $\bar{\alpha}(X)$ is the projection of the
difference of inverse propensity scores on the space spanned by $
X=(X_{1},X_{2},...)$ which is naturally approximated by linear combinations
of $X=(X_{1},...,X_{p}).$ Assumption 1 does not require that this $b(X)$
approximate the inverse propensity score.
For neural nets, random forests, and other learners that nonparametrically
estimate $\mathrm{E}[Y|X],$ Assumption 1 will require that a linear
combination of $b(X)$ can approximate any function of $X$ for large enough $
p.$ Such a $b(x)$ can be formed from low order multivariate powers of
components of $x$, with a full set of approximating functions included as $p$
grows. In applications one may use a variety of nonlinear functions
including powers of transformations of $X.$
The learner $\hat{\alpha}$ also depends on the choice of penalty degree $r.$
An important, useful feature of Lasso is that $r=A\sqrt{\ln(p)/n}$ for a
constant $A$ gives the fastest possible mean square convergence rate for
Lasso, that optimally trades off bias and variance. In Appendix~\ref
{sec:computing}, we describe cross-validation and theoretical methods for
choosing the choosing $r$ based on data that have proven stable across
several different applications. We also provide R code, available upon
request, for the construction of $\hat{\alpha}(x)$ and $\hat{\theta}$.
We can motivate $\hat{\rho}$ in $\hat{\alpha}(x)=b(x)^{\prime}\hat{\rho}$ as
being based on the Riesz representation in equation (\ref{Riesz rep}) and $
\bar{\alpha}$ satisfying equation (\ref{BLP alpha}), which imply that for $
m(w,b)=(m(w,b_{1}),...,m(w,b_{p}))^{\prime}$,
\begin{equation}
M:=\mathrm{E}[m(W,b)]=\mathrm{E}[\alpha_{0}(X)b(X)]=\mathrm{E}[\bar{\alpha}
(X)b(X)], \label{Rr basis}
\end{equation}
where the last equality is satisfied by $b_{j}\in\Gamma,$ which implies $
\mathrm{E}[b_{j}(X)\{\alpha_{0}(X)-\bar{\alpha}(X)\}]=0$ for each $j$. We
see that the cross moments $M$ between the true, unknown $\bar{\alpha}(x)$
and the dictionary $b(x)$ are equal to the expectation of the known vector
of functions $m(w,b).$ Also, the second moment matrix $G=\mathrm{E}
[b(X)b(X)^{\prime}]$ of the dictionary is an expectation of a known function
of the data. Estimating $M$ and $G$ enables learning coefficients $\rho$ of
the least squares regression of $\bar{\alpha}(X)$ on $b(X),$ satisfying $
M=G\rho.$ We learn $\rho$ using a Lasso minimum distance objective function
to allow for large $p$. Let
\begin{equation*}
\hat{M}=\frac{1}{n}\sum_{i=1}^{n}m(W_{i},b),\text{ }\hat{G}=\frac{1}{n}
\sum_{i=1}^{n}b(X_{i})b(X_{i})^{\prime},
\end{equation*}
be unbiased estimators of $M$ and $G.$ The coefficient estimator is given by
\begin{equation}
\hat{\rho}=\arg\min_{\rho}\{-2\hat{M}^{\prime}\rho+\rho^{\prime}\hat{G}
\rho+2r\left\Vert \rho\right\Vert _{1}\},\text{ }\left\Vert \rho\right\Vert
_{1}=\sum_{j=1}^{p}|\rho_{j}|. \label{RRLasso}
\end{equation}
The estimator $\hat{\rho}$ can be interpreted as a minimum distance version
of Lasso. Here $\hat{M}$ is analogous to $\sum_{i=1}^{n}Y_{i}b(X_{i})/n$ in
Lasso. The objective function in equation (\ref{RRLasso}) can be thought of
as the Lasso objective with $\sum_{i=1}^{n}Y_{i}b(X_{i})/n$ replaced by $
\hat{M}$ and $\sum_{i=1}^{n}Y_{i}^{2}/n$ dropped. In this way the objective
function is a penalized approximation to the least squares regression of $
\alpha_{0}(x)$ on $b(x),$ where $2r\left\Vert \rho\right\Vert _{1}$ is the
penalty. We refer to this as minimum distance Lasso because $\hat{M}$ does
not have the product form of Lasso regression.
The learner $\hat{\alpha}(x)$ of $\bar{\alpha}(x)$ is automatic in being
based on $\hat{M}$ and $\hat{G}$, neither of which requires knowledge of the
form of $\bar{\alpha}.$ In particular, $\hat{\alpha}(x)=b(x)^{\prime }\hat{
\rho}$ does not depend on plugging in nonparametric estimates of components
of $\bar{\alpha}(x).$ Instead, $b(x)^{\prime }\hat{\rho}$ is linear in the
dictionary $b(x)$ and uses the known functional $m(w,\gamma )$ in the
construction of $\hat{M}$ to obtain the learner $\hat{\rho}.$ This automatic
nature of $\hat{\alpha}(x)$ is especially useful for Lasso and other high
dimensional regression learners where $b(x)$ can be taken to be the first $p$
elements of $x=(x_{1},x_{2},...),$ and where $\bar{\alpha}(x)$ is a least
squares projection of $\alpha _{0}(X)$ on $\Gamma ,$ as in Section 2. The
projection $\bar{\alpha}(x)$ will generally not have a simple form that can
be learned by plugging in nonparametric learners to an explicit formula. For
instance, in the average treatment effect example the projection of the
inverse propensity score on the high dimensional regressors $
(X_{1},X_{2},...)$ does not have a closed form but is naturally approximated
by a linear combination of the first $p$ regressors where $
b(X)=(X_{1},...,X_{p})^{\prime }.$
The learner $\hat{\alpha}(x)=b(x)^{\prime }\hat{\rho}$ also avoids inverting
a learner of a conditional probability or pdf. The finite sample properties
of methods that rely on inverses of learners can be poor; see Singh and Sun
(2019) for recent examples. Instead, $\hat{\alpha}$ approximates and learns $
\bar{\alpha}$ by a linear combination of functions. In this way the $\hat{
\alpha}$ that we propose here avoids potential instability from inverting a
high dimensional estimator. The inverse of a conditional probability or
density is present in $\alpha _{0}(x)$ in all of the examples in this paper.
We anticipate that this feature is present quite generally for causal and
structural models involving shifts in regressors, because the Rr equation (
\ref{Riesz rep}) involves an expectation with respect to the data
distribution rather than the shifted distribution. Thus absence of an
inverse of a machine learner in $\hat{\alpha}$ may prove to be widely
useful. In some economic structural models the linearity of $\hat{\alpha}$
in $b(x)$ may not be quite as appealing, because inverse densities can have
a parametric form and so mitigate the problem of inverting a high
dimensional learner. An example is the dynamic discrete choice learner of
Chernozhukov et al. (2016, 2020). Also there is more work to be done to see
whether this approach has better properties than previously proposed ones in
practical settings.
This learner $\hat{\alpha}(x)$ can be thought of as being based on
orthogonality of the moment function with respect to $\gamma.$ Let $\tau$
denote a scalar and $b_{j}(x)$ an element of $b(x)$. Then by equation (\ref
{Rr basis})
\begin{equation*}
\frac{\partial}{\partial\tau}\mathrm{E}[\psi(W,\theta,\gamma+\tau b_{j},\bar{
\alpha })]=\mathrm{E}[m(W,b_{j})-\bar{\alpha}(X)b_{j}(X)]=0,\text{ }
(j=1,...,p).
\end{equation*}
Replacing the expectation by a sample average and $\bar{\alpha}(X)$ by $
b(X)^{\prime}\rho$ gives
\begin{equation*}
\frac{1}{n}\sum_{i=1}^{n}\{m(W_{i},b_{j})-[b(X_{i})^{\prime}
\rho]b_{j}(X_{i})\}=e_{j}^{\prime}(\hat{M}-\hat{G}\rho),
\end{equation*}
where $e_{j}$ is the jth columin of a $p$ dimensional identity matrix. This
sample average is a scaled version of the derivative of objective function
in equation (\ref{RRLasso}) without the penalty term. The first-order
conditions for equation (\ref{RRLasso}) will set $\hat{\rho}$ so that this
object is close to zero, subject to the penalty, i.e. will solve penalized
versions of a moment equation. Thus, the Lasso minimum distance learner can
be thought of as a method that uses orthogonality of $\psi(W,\theta,\gamma,
\alpha)$ with respect to $\gamma$ to learn $\bar{\alpha}$ while penalizing
to facilitate high dimensional estimation. In Section 6 we use an extension
of this approach to construct an Auto-DML when $m(W,\gamma)$ is nonlinear in
$\gamma.$
To illustrate $\hat{\alpha}$ we consider the choice of dictionary and the
form of $\hat{\alpha}$ for Examples 1-3.
\bigskip
\textsc{Example 1:} If the regression learner $\hat{\gamma}$ is
nonparametric the dictionary $b(X)$ should also be nonparametric while if $
\hat{\gamma}$ is a high dimensional regression the dictionary should be
chosen as in equation (\ref{high dimensional dict}). Here $m(w,b)=\int
b(x)[f_{1}(x)-f_{0}(x)]dx$ does not depend on the data observation $w$ and
the first order conditions for $\hat{\rho}$ imply that for each $j$,
\begin{equation*}
\left\vert \int b_{j}(x)[f_{1}(x)-f_{0}(x)]dx-\frac{1}{n-n_{\ell}}
\sum_{i\notin I_{\ell}}b_{j}(X_{i})\hat{\alpha}_{\ell}(X_{i})\right\vert
\leq r.
\end{equation*}
Here $\hat{\alpha}_{\ell}(X_{i})$ acts to approximately re-weight so that
the integral of the basis function $b_{j}(x)$ over the policy shift is
approximately equal to the sample average of the re-weighted basis function $
b_{j}(X_{i})\hat{\alpha}_{\ell}(X_{i}).$
\bigskip
\textsc{Example 2:} The dictionary $b(X)$ should be chosen as in Example 1.
Also by $m(w,b)=S(u)\gamma(u,z)$ the first order conditions for $\hat{\rho}$
imply that for each $j$,
\begin{equation*}
\left\vert \frac{1}{n-n_{\ell}}\sum_{i\notin
I_{\ell}}\{S(U_{i})b_{j}(U_{i},Z_{i})-b_{j}(X_{i})\hat{\alpha}
_{\ell}(X_{i})\}\right\vert \leq r.
\end{equation*}
Here $\hat{\alpha}_{\ell}(X_{i})$ acts approximately as a re-weighting
scheme, making the sample average of the score $S(U_{i})$ times the basis
function $b_{j}(U_{i},Z_{i})$ be approximately equal to the sample average
of the re-weighted basis function $b_{j}(X_{i})\hat{\alpha}_{\ell}(X_{i}).$
\bigskip
\textsc{Example 3: }The dictionary should be chosen similarly to Example 1.
For instance suppose that $X=(DZ,(1-D)Z)$, where $Z=(Z_{1},Z_{2},...)$ is a
sequence or possible covariates. Then the dictionary
\begin{equation}
b(x)=(dq(z)^{\prime},(1-d)q(z)^{\prime})^{\prime},\text{ }
q(z)=(z_{1},...,z_{p/2})^{\prime}, \label{ATE dic}
\end{equation}
would satisfy Assumption 1. The estimator $\hat{\alpha}_{\ell}$ has an
interesting form for this dictionary. Note that $
m(w,b)=b(1,z)-b(0,z)=(q(z)^{\prime},0^{\prime})^{\prime}-(0^{\prime},q\left(
z\right) ^{\prime})^{\prime}=(q(z)^{\prime},-q(z)^{\prime})$. Then
\begin{equation*}
\hat{M}_{\ell}=\left(
\begin{array}{c}
\bar{q}_{\ell} \\
-\bar{q}_{\ell}
\end{array}
\right) ,\text{ }\bar{q}_{\ell}=\frac{1}{n-n_{\ell}}\sum_{i\notin
I_{\ell}}q(Z_{i}).
\end{equation*}
Let $\hat{\rho}_{\ell}^{1}$ be the estimated coefficients of $dq(z)$ and $
\hat{\rho}_{\ell}^{0}$ be the estimated coefficients of $(1-d)q(z)$. Then
the learner of $\bar{\alpha}(X_{i})$ is
\begin{equation*}
\hat{\alpha}_{\ell}(X_{i})=D_{i}\hat{\omega}_{\ell i}^{1}-(1-D_{i})\hat {
\omega}_{\ell i}^{0},\text{ }\hat{\omega}_{\ell i}^{1}=q(Z_{i})^{\prime}\hat{
\rho}_{\ell}^{1},\text{ }\hat{\omega}_{\ell i}^{0}=-q(Z_{i})^{\prime}\hat{
\rho}_{\ell}^{0},
\end{equation*}
where $\hat{\omega}_{\ell i}^{1}$ and $\hat{\omega}_{\ell i}^{0}$ might be
thought of as \textquotedblleft weights.\textquotedblright\ These weights
sum to one if $q(z)$ includes a constant but may be negative. The first
order conditions for $\hat{\alpha}$ are that for each $j,$
\begin{equation}
\left\vert \frac{1}{n-n_{\ell}}\sum_{i\notin I_{\ell}}q_{j}(Z_{i})[1-D_{i}
\hat{\omega}_{\ell i}^{1}]\right\vert \leq r,\text{ }\left\vert \frac {1}{
n-n_{\ell}}\sum_{i\notin I_{\ell}}q_{j}(Z_{i})[1+(1-D_{i})\omega_{\ell
i}^{0}]\right\vert \leq r. \label{ATE balance}
\end{equation}
Here $\hat{\rho}_{\ell}$ sets the weights $\hat{\omega}_{\ell i}^{1}$ and $
\hat{\omega}_{\ell i}^{0}$ to approximately \textquotedblleft
balance\textquotedblright\ the overall sample average with the treated and
untreated averages for each element of the dictionary $q(z).$ The
constraints of equation (\ref{ATE balance}) are like the balancing
conditions of Zubizarreta (2015) and Athey, Imbens, and Wager (2018). The
source of these constraints is regularized least squares approximation of $
\bar{\alpha }(x)=proj(\pi_{0}(z)^{-1}d-[1-\pi_{0}(z)]^{-1}(1-d)|Z)$ by a
linear combination of the dictionary $b(x)$. The approach of this paper
shows that this type of balancing is sufficient to debias any regression
learner under regularity conditions in Section 4.
\section{Large Sample Inference}
In this Section, we give mean square convergence rates for the Lasso minimum
distance learner of $\hat{\alpha}$ and root-n consistency and asymptotic
normality results for the learner $\hat{\theta}$ of the object of interest
and its asymptotic variance estimator $\hat{V}$. Let $\varepsilon_{n}$
denote a sequence that converges to zero no faster than $\sqrt{\ln(p)/n}$
and for a random variable $a(W)$ let $\left\Vert a\right\Vert =\sqrt{\mathrm{
E}[a(W)^{2}]}$
\bigskip
\textsc{Assumption 2:} \textit{There exists }$C>1,$ $\xi>0$\textit{\ such
that for each positive integer }$s\leq C\varepsilon_{n}^{-2/(2\xi+1)}$
\textit{\ there is }$\bar{\rho}$\textit{\ with }$s$\textit{\ nonzero
elements such that}
\begin{equation*}
\left\Vert \bar{\alpha}-b^{\prime}\bar{\rho}\right\Vert \leq C(s)^{-\xi}.
\end{equation*}
\bigskip
Here $\left\Vert \bar{\alpha}-b^{\prime}\bar{\rho}\right\Vert $ is the mean
square approximation error from using the linear combination $b^{\prime}\bar{
\rho}$ to approximate $\bar{\alpha}.$ This approximate sparsity condition
specifies that there is a sparse $\bar{\rho}$, having only $s$ nonzero
elements, so that the approximation error is bounded by $C(s)^{-\xi}.$ Note
that it is not required that $\bar{\alpha}$ be equal to linear combination
of $s$ terms, i.e. it is not required that $\bar{\alpha}$ be strictly
sparse. Assumption 2 does allow unknown identity of the elements of $b(x)$
that give the approximation rate $s^{-\xi}$. In this way this condition
allows for high dimensional $x$ where statistics and economics do not
provide much guidance on which elements of $b(x)$ are important.
The $\varepsilon_{n}$ in this condition represents a convergence rate for $
\hat{M}$ and $\hat{G}$ that will be no faster than $\sqrt{\ln(p)/n}$ under
the conditions given in the rest of this Section. When $s$ is chosen to be
approximately $C\varepsilon_{n}^{-2/(2\xi+1)},$ which is the largest $s$
allowed by Assumption 2, $s$ will grow no faster than $(\sqrt{n/\ln (p)}
)^{2/(2\xi+1)}\leq n^{1/(2\xi+1)},$ which grows slower than $n.$ Because $
p\geq s$ is implicitly required by this condition, Assumption 2 puts a quite
a weak restriction on $p.$ An important feature of Assumption 2 is that the
sparse approximation is based on functions included in the $p\times1$
dictionary $b(x)$. Thus larger values of $p$ give more flexibility and will
help Assumption 2 to be satisfied.
Our results will require a convergence rate for $\hat{\alpha}$ that is
faster than some power of $n.$ Assumption 2 is a natural condition that
leads to such a rate. Sufficient conditions for Assumption 2 are well known
from the approximation literature when $\bar{\alpha}(x)$ belongs to a Besov
or Holder class of function and linear combinations of $b(x)$ can
approximate any function of $x$.
We will also make use of a sparse eigenvalue condition as considered in much
of the Lasso literature. Let $\rho$ denote a $p\times1$ vector, $\rho_{J}$ a
$J\times1$ subvector of $\rho,$ and $\rho_{J^{c}}$ the vector consisting of
components of $\rho$ that are not in $\rho_{J}$. Also for a matrix $A$ let $
\left\Vert A\right\Vert _{1}=\sum_{i,j}\left\vert a_{ij}\right\vert .$
\bigskip
\textsc{Assumption 3: }$G=\mathrm{E}[b(X)b(X)^{\prime}]$ \textit{has largest
eigenvalue bounded uniformly in }$n$\textit{\ and there is }$C,c>0$\textit{\
such that for all }$s\approx C\varepsilon_{n}^{-2}$\textit{\ with
probability approaching one}
\begin{equation*}
\min_{J\leq s}\min_{\left\Vert \rho_{J^{c}}\right\Vert _{1}\leq3\left\Vert
\rho_{J}\right\Vert _{1}}\frac{\rho^{\prime}\hat{G}\rho}{\rho_{J}^{\prime}
\rho_{J}}\geq c
\end{equation*}
\bigskip
This is a sparse eigenvalue condition that is familiar from the Lasso
literature, including Bickel, Ritov, Tsybakov (2009), Belloni and
Chernozhukov (2013), and Rudelson and Zhou (2013).
We will work with a dictionary $b(X)$ with elements that are uniformly
bounded.
\bigskip
\textsc{Assumption 4}: \textit{There is }$C>0$\textit{\ such that with
probability one }$\sup_{j}|b_{j}(X)|\leq C.$
\bigskip
This condition implies a convergence rate of $\sqrt{\ln(p)/n}$ for $
\left\Vert \hat{G}-G\right\Vert _{\infty},$ where $\left\Vert A\right\Vert
_{\infty}=\max_{i,j}\left\vert a_{ij}\right\vert $ for a matrix $A=[a_{ij}]$.
Lasso mean square convergence rates are often stated in terms of finite
sample bounds. Because the focus of this paper is root-n consistency for $
\hat {\theta}$ and for that we only need convergence at certain powers of $n$
we can simplify the statement of convergence rates without affecting the
conditions for $\hat{\theta}$ by allowing the Lasso regularization value $r$
to shrink slightly slower than $\varepsilon_{n}.$ This does lead to
approximate sparseness conditions that are strict inequalities on the size
of $\xi$ but Bradic et al. (2019) have shown that strict inequalities are
necessary for root-n consistent estimation, meaning that there is no loss of
generality in these conditions. We also limit the growth of $p$ to be slower
than some power of $n.$
\bigskip
\textsc{Assumption 5:} $\varepsilon_{n}=o(r),$ $r=o(n^{c}\varepsilon_{n})$
for all $c>0$, and there exists $C>0$ such that $p\leq Cn^{C}.$
\bigskip
We also hypothesize a convergence rate for $\hat{M}.$
\bigskip
\textsc{Assumption 6}: $\left\Vert \hat{M}-M\right\Vert
_{\infty}=O_{p}(\varepsilon_{n})$ for $\varepsilon_{n}\longrightarrow0.$
\bigskip
We use this condition to accommodate $\hat{M}$ that can depend on the
regression learner $\hat{\gamma}$ as needed for Section 5.
\bigskip
\textsc{Theorem 1:} \textit{If Assumptions 1 - 6 are satisfied then for all }
$c>0,$
\begin{equation*}
\Vert\hat{\alpha}-\bar{\alpha}\Vert=o_{p}(n^{c}\varepsilon_{n}^{2\xi/(2\xi
+1)}).
\end{equation*}
\bigskip
This theorem is based on extending Lemmas of Bradic et al. (2019) to allow $
\varepsilon_{n}$ to shrink slower than $\sqrt{\ln(p)/n}.$ The extension will
be used in Section 5 to obtain convergence rates when $\hat{M}$ depends on a
nonparametric estimator.
The sparse eigenvalue condition of Assumption 3 seems strong in some
settings. It is possible to drop Assumption 3 and Assumption 2 if the
following condition is satisfied:
\bigskip
\textsc{Assumption 7: }$\bar{\alpha}(X)=\sum_{j=1}^{\infty}\rho_{j0}b_{j}(X)$
\textit{, }$\sum_{j=1}^{\infty}\left\vert \rho_{j0}\right\vert <\infty $
\textit{, and for }$C>0$\textit{\ and }$\bar{s}=C\sqrt{n}$\textit{\ the }$
b_{j}(x)$\textit{\ corresponding to the largest }$\bar{s}$\textit{\ values
of }$\left\vert \rho_{j0}\right\vert $\textit{\ are included in }$b(x).$
\bigskip
This condition allows us to drop Assumption 2 because absolute summability
of the coefficients $\rho_{0j}$ implies a sparse approximation rate of $
\xi=1/2.$ It also allows $\hat{G}$ to converge at a rate slower $
\varepsilon_{n}$ in order to accommodate nonparametric estimation in $\hat{G}
.$
\bigskip
\textsc{Theorem 2:} \textit{If Assumptions 1 and 5-7 are satisfied and }$
\left\Vert \hat{G}-G\right\Vert _{\infty}=O_{p}(\varepsilon_{n})$\textit{\
then for all }$c>0,$
\begin{equation*}
\Vert\hat{\alpha}-\bar{\alpha}\Vert=o_{p}(n^{c}\sqrt{\varepsilon_{n}}).
\end{equation*}
\bigskip
This result extends Chatterjee and Javarov (2015) to allow $\varepsilon_{n}$
to shrink slower than $\sqrt{\ln(p)/n}.$ When $\varepsilon_{n}=\sqrt{\ln
(p)/n}$ in Assumption 6 this result gives a mean square convergence rate for
$\hat{\alpha}$ that is faster than $n^{-1/4+c}$ for all $c>0,$ without a
sparse eigenvalue condition.
We now use these results to obtain root-n consistency and asymptotic
normality for the Auto-DML $\hat{\theta}$ and consistency of its asymptotic
variance estimator $\hat{V}.$ We impose some additional regularity
conditions.
\bigskip
\textsc{Assumption 8}: \textit{There is }$C>0$\textit{\ such that with
probability one }$\max_{j\leq p}|m(W,b_{j}))|\leq C.$
\bigskip
Under this condition Assumption 6 will be satisfied with $\varepsilon _{n}=
\sqrt{\ln(p)/n}.$ This condition will be satisfied under by Assumption 4 in
each of Examples 1-3 under conditions of Corollaries 4-6 to follow.
\bigskip
\textsc{Assumption 9:}\textit{\ }$\mathrm{E}[\{Y-\bar{\gamma}(X)\}^{2}|X]$
\textit{\ and }$\bar{\alpha}(X)$\textit{\ are bounded.}
\bigskip
We impose this condition for simplicity; it could be weakened. We also
impose the following condition.
\bigskip
\textsc{Assumption 10:} $\mathrm{E}[m(W,\gamma_{0})^{2}]<\infty$ \textit{and
}$\int[m(w,\hat{\gamma})-m(w,\bar{\gamma})]^{2}F_{W}(dw)\overset{p}{
\longrightarrow}0.$
\bigskip
This condition will be implied by existence of $C>0$ with $\left\vert
\mathrm{E}[m(W,\gamma)^{2}]\right\vert \leq C\left\Vert \gamma\right\Vert
^{2}$ for all $\gamma$, which will be satisfied in the examples we consider
under regularity conditions to be specified.
\bigskip
\textsc{Assumption 11: }\textit{With probability approaching one }$\hat {
\gamma}_{\ell}\in\Gamma$\textit{\ and there is }$d_{\gamma}>0$\textit{\ such
that }$\Vert\hat{\gamma}-\bar{\gamma}\Vert=O_{p}(n^{-d_{\gamma}})$\textit{\
and either Assumptions 2 and 3 are satisfied with}
\begin{equation}
\frac{\xi}{2\xi+1}+d_{\gamma}>\frac{1}{2}, \label{rate dr}
\end{equation}
\textit{or Assumption 7 is satisfied and }$d_{\gamma}>1/4.$
\bigskip
This assumption allows $\hat{\gamma}$ to be any learner that converges in
mean square at a rate that is some power of $n.$ By Theorem 1, the mean
square convergence rate for $\hat{\alpha}$ is as close as desired to $
n^{-\xi /(2\xi+1)}.$ Thus Assumption 11 requires that the product of
convergence rates for $\hat{\alpha}$ and $\hat{\gamma}$ must go to zero
faster than $1/\sqrt {n}.$ This is a rate double robustness condition that
appears in earlier low dimensional and high dimensional literatures cited in
the introduction. Under Assumptions 2 and 3 a full trade-off in rates
between $\hat{\alpha}$ and $\hat{\gamma}$ is permitted, since Assumption 11
is satisfied for any $\xi$ if $d_{\gamma}$ is large enough and for any $
d_{\gamma}$ if $\xi$ is large enough. Under Assumption 7 this trade-off is
not present, since $d_{\gamma }>1/4$ is required by Assumption 11. Assumption 11 can be dropped if $\alpha_0(X)$ is known and is used in place of $\hat\alpha(X)$ in the construction of $\hat\theta$ in equation (3.1). In that case only mean square consistency of $\hat\gamma$ will be required for root-n consistency and asymptotic normality of $\hat\theta$.
The following gives the large sample inference results for $\hat{\theta}$
and $\hat{V}.$ Define
\begin{equation*}
\bar{\theta}=\mathrm{E}[m(W,\bar{\gamma})],\text{ }\psi(w)=m(w,\bar{\gamma})-
\bar{\theta}+\bar{\alpha}(x)[y-\bar{\gamma}(x)],\text{ }V=\mathrm{E}
[\psi(W)^{2}].
\end{equation*}
Here $\bar{\theta}$ will be the object estimated by $\hat{\theta}$ when
neither of the double robustness conditions $\bar{\gamma}(X)=\mathrm{E}[Y|X]$
nor $\bar{\alpha}(X)\in\Gamma$ is satisfied.
\bigskip
\textsc{Theorem 3}: \textit{If Assumptions 1-5, and 8-11 are satisfied then }
$\sqrt{n}(\hat{\theta}-\bar{\theta})\overset{d}{\longrightarrow}N(0,V).$
\textit{If in addition Assumption 7 is satisfied then }$\hat{V}\overset{p}{
\longrightarrow}V$.
\bigskip
It is possible to construct a consistent estimator of $V$ without Assumption
7 by using a trimmed version of $\hat{\alpha}_{\ell }(x)$ but we omit that
demonstration to avoid further complicating $\hat{V}$. The conclusion of
Theorem 3 implies that asymptotic test statistics and confidence intervals
can be formed in the usual manner from $\hat{\theta}$ and $\hat{V}.$ Theorem
3 is proven by using the convergence rate results of Theorem 1 and Theorem 2
to show that the hypotheses of Lemma 15 of Chernozhukuv et al. (2020) are
satisfied.
The asymptotic variance $V$ is fixed rather than varying with $n$ because we
have chosen to work with i.i.d. data and an approximately sparse regression
for simplicity. It would be straightforward to extend the results to allow
the regression to change with sample size in order to accomodate sparse
regressions and corresponding variances that change with $n$.
Under similar conditions as Theorem 3 Auto-TML is also consistent and
asymptotically normal.
\bigskip
\textsc{Corollary 4:} \textit{If Assumptions 1-5, and 8-11 are satisfied, }$
\mathrm{E}[m(W,\gamma )^{2}]\leq C\left\Vert \gamma \right\Vert^{2} $ \textit{for all
}$\gamma \in \Gamma ,$\textit{\ and }$\bar{\alpha}(X)\neq 0$\textit{\ then }$
\sqrt{n}(\tilde{\theta}-\bar{\theta})\overset{d}{\longrightarrow }N(0,V).$
\bigskip
Most of the conditions of Theorem 3 are quite general, with only Assumptions
8 and 10 pertaining to a particular $m(w,\gamma )$. It is straightforward to
specify conditions under which Assumptions 8 and 10 are satisfied for
Examples 1-3.
\bigskip
\textsc{Corollary 5 (Example 1):} \textit{If Assumptions 1-5, 9, and 11 are
satisfied and there is }$C>0$\textit{\ such that }$\left\vert
[f_{1}(x)-f_{0}(x)]/f(x)\right\vert \leq C$ \textit{then }$\sqrt{n}(\hat{
\theta}-\bar{\theta})\overset{d}{\longrightarrow }N(0,V).$\textit{\ If in
addition Assumption 7 is satisfied then }$\hat{V}\overset{p}{\longrightarrow
}V$\textit{.}
\bigskip
The specific regularity condition for the policy effect in Corollary 5 is
that the Rr $\alpha _{0}(X)=[f_{1}(X)-f_{0}(X)]/f(x)$ be bounded.
\bigskip
\textsc{Corollary 6 (Example 2):} \textit{If Assumptions 1-5, 9, and 11 are
satisfied and there is }$C>0$\textit{\ such that }$\left\vert
S(u)\right\vert \leq C$, $f(D|Z)^{-1}\omega (D)\leq C$\textit{\ then }$\sqrt{
n}(\hat{\theta}-\bar{\theta})\overset{d}{\longrightarrow }N(0,V).$\textit{\
If in addition Assumption 7 is satisfied then }$\hat{V}\overset{p}{
\longrightarrow }V$\textit{.}
\bigskip
The regularity conditions for the weighted average derivative in Corollary 6
are that the score $S(u)$ is bounded and the Rr $\alpha
_{0}(X)=f(D|Z)^{-1}\omega (D)S(D)$ is also bounded.
\bigskip
\textsc{Corollary 7 (Example 3): }\textit{If Assumptions 1, 4-5, 9, and 11
are satisfied and there is }$C>0$\textit{\ with }$\pi _{0}(Z)\in \lbrack
C,1-C]$\textit{\ then }$\sqrt{n}(\hat{\theta}-\bar{\theta})\overset{d}{
\longrightarrow }N(0,V).$\textit{\ If in addition Assumption 7 is satisfied
then }$\hat{V}\overset{p}{\longrightarrow }V$\textit{.}
\bigskip
The additional condition in Corollary 7 is that the propensity score is
bounded away from $0$ and $1$, an overlap condition that is common in
asymptotic theory for estimators of the average treatment effect. Together
Corollaries 5--7 demonstrate how simple primitive conditions involving $
m(w,\gamma )$ can be specified so that the Auto-DML $\hat{\theta}$ of an
object of interest will be asymptotically normal and the asymptotic variance
estimator $\hat{V}$ consistent.
\section{Nonlinear Effects of Multiple Regressions}
Some important effects of interest are expectations of nonlinear functions
of multiple regressions. Causal mediation analysis is an important example
that we consider in this Section. The regression decomposition in Section 6
is another important example. In this Section we give Auto-DML for such
effects. Such effects have the form $\theta_{0}=\mathrm{E}[m(W,\gamma_{0})]$
where $m(w,\gamma)$ is nonlinear in a possible value $\gamma$ of multiple
regressions $(\gamma _{1}(X_{1}),...,\gamma_{K}(X_{K}))^{\prime}$ with
regressors $X_{k}$ specific to each regression $\gamma_{k}(X_{k})$. The
corresponding orthogonal moment functions are like those discussed in
Section 3 except that the bias correction is a sum of $K$ terms with the $
k^{th}$ term being the bias correction for the learner of $\gamma_{k}$, as
in Newey (1994, p. 1357). The estimated bias corrections are like those of
Section 4 with the $k^{th}$ term being the product of a Lasso learner $\hat{
\alpha}_{k\ell}(X_{k})$ and the residual $Y_{k}-\hat{\gamma}_{k\ell}(X_{k}).$
Each $\hat{\alpha}_{k\ell}(X_{k})$ differs from Section 3 in the
corresponding $\hat{M}_{k\ell}$ being a derivative evaluated at a
preliminary estimator of $\bar{\gamma}$. Because the construction of $\hat{
\theta}$ is so closely related to that in Section 3 we proceed immediately
with its description here and fill in details concerning the orthogonal
moment function below.
The Auto-DML of a nonlinear effect is similar to equation (\ref{Estimator}).
Specifically it is
\begin{align}
\hat{\theta}& =\frac{1}{n}\sum_{\ell =1}^{L}\sum_{i\in I_{\ell }}\{m(W_{i},
\hat{\gamma}_{\ell })+\sum_{k=1}^{K}\hat{\alpha}_{k\ell }(X_{ki})[Y_{ki}-
\hat{\gamma}_{k\ell }(X_{ki})]\}, \label{nonlin est} \\
\hat{V}& =\frac{1}{n}\sum_{\ell =1}^{L}\sum_{i\in I_{\ell }}\hat{\psi}
_{i\ell }^{2},\text{ }\hat{\psi}_{i\ell }=m(W_{i},\hat{\gamma}_{\ell })-\hat{
\theta}+\sum_{k=1}^{K}\hat{\alpha}_{k\ell }(X_{ki})[Y_{ki}-\hat{\gamma}
_{k\ell }(X_{ki})], \notag
\end{align}
where each $\hat{\alpha}_{k\ell }(X_{ki})$ is obtained as follows: For each $
k$ let $b_{k}(x_{k})=(b_{k1}(x_{k}),....,b_{kp}(x_{k}))^{\prime }$ be a $
p\times 1$ dictionary vector specific to the $k^{th}$ regression $\gamma
_{k}(x_{k})$ and let $\hat{\gamma}_{\ell ,\ell ^{\prime }}$ be the vector of
regressions computed from all observations not in either $I_{\ell }$ or $
I_{\ell ^{\prime }}$. Also let $\tau $ denote a scalar, and $e_{k}$ the $
k^{th}$ column of the $K$ dimensional identity matrix. Then
\begin{align}
\hat{\alpha}_{k\ell }(X_{ki})& =b_{k}(X_{ki})^{\prime }\hat{\rho}_{k\ell },
\text{ }\hat{\rho}_{k\ell }=\arg \min_{\rho }\{-2\hat{M}_{k\ell }^{\prime
}\rho +\rho ^{\prime }\hat{G}_{k\ell }\rho +2r_{k}\left\Vert \rho
\right\Vert _{1}\},\text{ }\left\Vert \rho \right\Vert
_{1}=\sum_{j=1}^{p}|\rho _{j}|, \label{nonlin Rr} \\
\hat{M}_{k\ell }& =(\hat{M}_{k\ell 1},...,\hat{M}_{k\ell p})^{\prime },\text{
}\hat{G}_{k\ell }=\left( \frac{1}{n-n_{\ell }}\right) \sum_{i\notin I_{\ell
}}b_{k}(X_{ki})b_{k}(X_{ki})^{\prime }, \notag \\
\hat{M}_{k\ell j}& =\left. \frac{d}{d\tau }\left( \frac{1}{n-n_{\ell }}
\right) \sum_{\ell ^{\prime }\neq \ell }\sum_{i\in I_{\ell ^{\prime
}}}m(W_{i},\hat{\gamma}_{\ell ,\ell ^{\prime }}+\tau e_{k}b_{kj})\right\vert
_{\tau =0},\text{ }(j=1,...,p). \notag
\end{align}
where $b_{kj}$ denotes the $j^{th}$ element of the dictionary $b_{k}(x_{k})$
as a function of $x_{k}.$ Thus the $\hat{\alpha}_{k\ell }(X_{i})$ in
equation (\ref{nonlin est}) is a Lasso minimum distance estimator like that
of Section 3 that is specific to $\hat{\gamma}_{k}$ and uses the $\hat{M}
_{k\ell }$ from equation (\ref{nonlin Rr}) rather than the one in equation (
\ref{CFLasso}).
The $\hat{M}_{k\ell j}$ given here generalizes equation (\ref{CFLasso}) to
allow for nonlinearity of $m(w,\gamma)$ in $\gamma.$ The derivative with
respect to the scalar $\tau$ in $\hat{M}_{k\ell j}$ is generally simple to
compute analytically using the chain rule of calculus, as we will illustrate
for causal mediation analysis. When $m(w,\gamma)$ is linear in a single $
\gamma$ this derivative just evaluates $m(W_{i},\gamma)$ at $\gamma=b_{j}$,
giving the $\hat{M}_{\ell j}$ of equation (\ref{CFLasso}). As with linear $
m(w,\gamma)$ the $\hat{M}_{k\ell j}$ and the rest of the $\hat{\theta}$
depends just on $m(w,\gamma)$ and the first step. Thus the $\hat{\theta}$ in
equation (\ref{nonlin est}) is automatic, in the same way as the estimator
of equation (\ref{Estimator}), in only requiring $m(w,\gamma)$ and the
regression residuals $Y_{{}}$for its construction.
The $\hat{M}_{k\ell j}$ given here does depend on a cross-fit regression
learner $\hat{\gamma}_{\ell ,\ell ^{\prime }}$ in order to allow for the
nonlinearity of $m(w,\gamma )$ in $\gamma .$ The cross-fitting will make the
sample average used in the construction of $\hat{M}_{k\ell j}$ independent
of the regression learner $\hat{\gamma}_{\ell ,\ell ^{\prime }}$ used in its
construction. This independence helps $\hat{M}_{k\ell j}$ to be uniformly
consistent over $j=1,...,p$ for large $p$ with only mean square convergence
convergence rates for $\hat{\gamma}_{\ell ,\ell ^{\prime }}.$ This feature
of the theory helps $\hat{\theta}$ to be root-n consistent and
asymptotically normal for a wide variety of regression learners $\hat{\gamma}
_{\ell ,\ell ^{\prime }}.$ This $\hat{M}_{k\ell j}$ was given in
Chernozhukov, Newey, and Singh (2018, p. 17). Multiple cross-fitting has
also been used in Newey and Robins (2018) and Kennedy (2020).
The dictionary $b_{k}(x_{k})$ used in the construction of $\hat{\alpha}
_{k\ell}(x_{k})$ should be chosen analogously to the $b(x)$ in Section 3.
Each $b_{kj}$ should be an element of the set $\Gamma_{k}$ of possible
plim's of $\hat{\gamma}_{k}$. Also linear combinations of $b_{k}(x_{k})$
should be able to approximate any element of $\Gamma_{k}$ arbitrarily well
in mean square. That is, Assumption 1 should be satisfied with $\Gamma_{k}$
and $b_{k}(x)$ replacing $\Gamma$ and $b(x)$ respectively. In particular if $
\hat{\gamma}_{k}$ is a high dimensional regression then $
b(x)=(x_{k1},...,x_{kp})^{\prime }$ will do. If $\hat{\gamma}_{k}$ is a
nonparametric estimator then $b_{k}(x_{k})$ should be chosen so that linear
combinations can approximate any function of $x_{k}$.
An important difference between the Lasso minimum distance learner in
Section 3 and each $\hat{\alpha}_{k\ell}(x_{k})$ here is that the penalty
size $r_{k}$ must be chosen to be larger than $\sqrt{\ln(p)/n}$ when $
m(w,\gamma)$ depends nonlinearly on $\gamma.$ The reason for larger $r_{k}$
is that $\hat{M}_{k\ell}$ depends on the machine learner $\hat{\gamma}
_{\ell,\ell^{\prime}}$ and so will converge at a slower rate, leading to a
requirement that $r_{k}$ converge to zero slightly slower than the mean
square convergence rate of $\hat{\gamma}_{\ell,\ell^{\prime}}.$ A choice of $
r_{k}$ proportional to $n^{-1/4}$ will generally suffice for this purpose,
since $\hat{\gamma}_{\ell,\ell^{\prime}}$ will be required to converge
faster than $n^{-1/4}$.
This estimator will not be doubly robust due to the nonlinearity of $
m(w,\gamma)$ in $\gamma$; see Chernozhukov et al. (2016). Nevertheless it
will have zero first order bias and so be root-n consistent and
asymptotically normal under sufficient regularity conditions. It has zero
first order bias because $\hat{\alpha}_{k\ell}(x_{k})$ will consistently
estimate $\bar{\alpha }_{k}(x_{k})$ such that $\sum_{k=1}^{K}\bar{\alpha}
_{k}(x)[y_{k}-\bar{\gamma }_{k}(x_{k})]$ is the influence function for $
\mathrm{E}[m(W,\gamma(F))]$ at $\gamma(F)=\bar{\gamma}$ where $\gamma(F)=$
plim$(\hat{\gamma})$.
\bigskip
\textsc{Example 5:} (Causal Mediation Analysis) Causal mediation analysis
provides an interesting example of a nonlinear function of multiple
regressions. This effect allows for intermediate variables, called
mediators, that lie between treatment and outcome. In this example there is
an outcome variable $Y$, a treatment indicator $D\in\{0,1\},$ and covariates
$Z$ similar to the average treatment effect in Example 3. In addition there
is a mediation variable that we will denote by $Q,$ where we assume that $
Q\in\{1,...,K-1)$ for an integer $K\geq3.$ Let
\begin{equation*}
\gamma_{K0}(D,Q,Z)=\mathrm{E}[Y|D,Q,Z]\text{, }\gamma_{k0}(D,Z)=\Pr
(Q=k|D,Z)=\mathrm{E}[1(Q=k)|D,Z],\text{ }(k=1,...,K-1).
\end{equation*}
The causal mediation effect of Imai, Keele, and Tingley (2010, Theorem 1) is
\begin{equation*}
\theta_{0}(d,d^{\prime})=\mathrm{E}[\sum_{k=1}^{K-1}\gamma_{K0}(d,k,Z)\gamma
_{k0}(d^{\prime},Z)].
\end{equation*}
This effect, or parameter, has the form $\theta_{0}(d,d^{\prime})=\mathrm{E}
[m(W,\gamma)]$ for $W=(Y,D,Q,Z)$ and
\begin{equation*}
m(W,\gamma)=\sum_{k=1}^{K-1}\gamma_{K}(d,k,Z)\gamma_{k}(d^{\prime},Z).
\end{equation*}
In this example we have $X_{k}=(D,Z)$, $(k=1,...,K-1)$ and $X_{K}=(D,Q,Z).$
To construct the Auto-DML $\hat{\theta}$ we need to choose the dictionaries $
b_{k}(X_{k})$ for each $k$. We choose
\begin{equation*}
b_{K}(X_{K})=(b_{K1}(D,Q,Z),...,b_{Kp}(D,Q,Z))^{\prime}
\end{equation*}
to be a nonparametric dictionary if $\hat{\gamma}_{K}$ is a nonparametric
estimator such as a neural net or random forest or choose $b_{K}(D,Q,Z)$ to
be the leading $p$ regressors used in a high dimension regression learner $
\hat{\gamma}_{K}.$ For $k\leq K-1$ we choose the same dictionary $
b_{k}(X_{k})=b_{1}(D,Z)$ with
\begin{equation*}
b_{1}(D,Z)=(b_{11}(D,Z),...,b_{1p}(D,Z))^{\prime},
\end{equation*}
for each $k\leq K-1.$ We specify $b_{1}(D,Z)$ to be a nonparametric
dictionary if each $\hat{\gamma}_{k}$ is a nonparametric estimator such as a
neural net or random forest or choose $b_{1}(D,Z)$ to be the leading $p$
regressors used in a high dimension regression learner for each $\hat{\gamma}
_{k}.$
It is straightforward to compute each $\hat{M}_{k\ell j}.$ Note that for $
k\leq K-1,$
\begin{align*}
\left. \frac{d}{d\tau }m(W,\gamma +\tau e_{k}b_{kj})\right\vert _{\tau =0}&
=\left. \frac{d}{d\tau }\gamma _{K}(d,k,Z)\{\gamma _{k}(d^{\prime },Z)+\tau
b_{1j}(d^{\prime },Z)\}\right\vert _{\tau =0}=\gamma
_{K}(d,k,Z)b_{1j}(d^{\prime },Z), \\
\left. \frac{d}{d\tau }m(W,\gamma +\tau e_{K}b_{Kj})\right\vert _{\tau =0}&
=\left. \frac{d}{d\tau }\{\sum_{k=1}^{K-1}\{\gamma _{K}(d,k,Z)+\tau
b_{Kj}(d,k,Z)\}\gamma _{k}(d^{\prime },Z)]\right\vert _{\tau =0} \\
& =\sum_{k=1}^{K-1}b_{Kj}(d,k,Z)\gamma _{k}(d^{\prime },Z).
\end{align*}
Then we have
\begin{align*}
\hat{M}_{k\ell j} & =\frac{1}{n-n_{\ell}}\sum_{\ell^{\prime}\neq\ell}\sum_{i
\in I_{\ell^{\prime}}}\hat{\gamma}_{K\ell,\ell^{
\prime}}(d,k,Z_{i})b_{1j}(d^{\prime},Z_{i}),\text{ }(k=1,...,K-1), \\
\hat{M}_{K\ell j} & =\frac{1}{n-n_{\ell}}\sum_{\ell^{\prime}\neq\ell}\sum_{i
\in I_{\ell^{\prime}}}\sum_{k=1}^{K-1}b_{Kj}(d,k,Z_{i})\hat{\gamma }
_{k\ell,\ell^{\prime}}(d^{\prime},Z_{i}),\text{ }(j=1,...,p).
\end{align*}
We can then compute $\hat{\alpha}_{k\ell}(x)$ as in equation (\ref{nonlin Rr}
) and $\hat{\theta}$ for $Y_{ki}=1(Q_{i}=k),$ $(k=1,...,K-1)$ and $
Y_{Ki}=Y_{i}$ as in equation (\ref{nonlin est}).
The orthogonal moment function corresponding to this estimator is
\begin{align*}
\psi (W,\gamma ,\alpha ,\theta )& =\sum_{k=1}^{K-1}\gamma _{K}(d,k,Z)\gamma
_{k}(d^{\prime },Z)-\theta +\alpha _{K}(D,Q,Z)[Y-\gamma _{K}(D,Q,Z)] \\
+\sum_{k=1}^{K-1}\alpha _{k}(D,Z)[1(Q& =k)-\gamma _{k}(D,Z)],\text{ }\gamma
_{K},\alpha _{K}\in \Gamma _{K},\text{ }\gamma _{k},\alpha _{k}\in \Gamma
_{1},\text{ }(k\leq K-1).
\end{align*}
where $\Gamma _{K}$ is the set of possible plims of $\hat{\gamma}_{K}$ and $
\Gamma _{1}$ is the set of plims of $\hat{\gamma}_{k}$ for $k\leq K-1$. This
moment function differs from the multiply robust moment function of Tchetgen
Tchetgen and Shipster (2012) in imposing the constraint that each $\gamma
_{k}$ and $\alpha _{k}$ are contained in the set $\Gamma _{k}$ of possible
plim's of $\hat{\gamma}_{k}.$ For example, when $\hat{\gamma}_{K}$ is a high
dimensional regression estimator $\gamma _{K}$ and $\alpha _{K}$ must be
elements of the mean square span of $(X_{1},X_{2},...)$ similarly to Section
2. It has the multiple robustness feature that for $\bar{\theta}=\mathrm{E}
[m(W,\bar{\gamma})]$ and any $\alpha =(\alpha _{1},...,\alpha _{K})\in \Pi
_{k=1}^{K}\Gamma _{k},$
\begin{equation*}
\mathrm{E}[\psi (W,\bar{\gamma},\alpha ,\bar{\theta})]=0,
\end{equation*}
shown in Chernozhukov et al. (2020) to be a general feature of orthogonal
moment functions constructed from the influence function of $\mathrm{E}
[m(W,\gamma (F))]$. It also has other multiple robustness features. For $
\alpha _{k0},$ $(k=1,...,K)$ given in the proof of Corollary 10 in the
Appendix, when $\alpha _{k0}\in \Gamma _{1},$ $(k\leq K-1)$ and $\alpha
_{K0}\in \Gamma _{K},$
\begin{equation*}
\mathrm{E}[\psi (W,\gamma _{10},...,\gamma _{K-1,0},\gamma _{K},\alpha
_{0},\theta _{0})]=0,\text{ }\mathrm{E}[\psi (W,\gamma _{1},...,\gamma
_{K-1},\gamma _{K0},\alpha _{0},\theta _{0})]=0,\text{ }
\end{equation*}
for any $\gamma _{K}\in \Gamma _{K}$ and $\gamma _{k}\in \Gamma _{1},$ $
(k\leq K-1).$
\bigskip
We now return to the general learner $\hat{\theta}$ and give regularity
conditions for asymptotic normality and consistent estimation of the
asymptotic variance of $\hat{\theta}$. For $\tilde{\gamma}=(\tilde{\gamma}
_{1},...,\tilde{\gamma}_{K})^{\prime}\in\Pi_{k=1}^{K}\Gamma_{k}$ and $
\gamma_{k}\in\Gamma_{k}$ let
\begin{equation*}
D_{k}(W,\gamma_{k},\tilde{\gamma}):=\left. \frac{\partial m(W,\tilde{\gamma }
+e_{k}\tau\gamma_{k})}{\partial\tau}\right\vert _{\tau=0}
\end{equation*}
be the Gateaux derivative of $m(W,\gamma)$ with respect to $\gamma_{k}$ when
it exists. Comparing this definition with equation (\ref{nonlin Rr}) we see
that each $\hat{M}_{kj\ell}$ is an average of values of this Gateaux
derivative. We impose the following condition on these derivatives.
\bigskip
\textsc{Assumption 12: }\textit{There are }$C,$\textit{\ }$\varepsilon >0,$
\textit{\ }$a_{kj}(w),$\textit{\ and }$A_{k}(w,\gamma)$\textit{\ such that
for all }$\gamma$\textit{\ with }$\Vert\gamma-\bar{\gamma}\Vert\leq
\varepsilon,$\textit{\ }$D_{k}(W,b_{kj},\gamma)$\textit{\ exists and for }$
k=1,...,K$
\begin{align*}
D_{k}(W,b_{kj},\gamma) & =a_{kj}(W)A_{k}(W,\gamma),\text{ }\max_{j\leq
p}\left\vert \mathrm{E}[a_{kj}(W)\{A_{k}(W,\gamma)-A_{k}(W,\bar{\gamma}
)\}]\right\vert \leq C\left\Vert \gamma-\bar{\gamma}\right\Vert , \\
\max_{j\leq p}\left\vert a_{kj}(W)\right\vert & \leq C,\text{ }\mathrm{E}
[A_{k}(W,\gamma)^{2}]\leq C.
\end{align*}
\bigskip
This condition and the use of the cross-fit $\hat{\gamma}_{\ell,\ell^{
\prime}}$ in $\hat{M}_{k\ell}$ lead to a convergence rate for $\hat{M}
_{k\ell}.$ Let $M_{kj}=\mathrm{E}[D_{k}(W,b_{kj},\bar{\gamma})]\,$\ and $
M_{k}=(M_{k1},...,M_{kp}),$ $(j=1,...,p;k=1,...,K).$
\bigskip
\textsc{Lemma 8:} \textit{If there is }$0<d_{\gamma }<1/2$ such that $
\left\Vert \hat{\gamma}_{k\ell ,\ell ^{\prime }}-\bar{\gamma}_{k\ell ,\ell
^{\prime }}\right\Vert =O_{p}(n^{-d_{\gamma }}),$ $(k=1,...,K;\ell ,\ell
^{\prime }=1,...L)$, \textit{and Assumption 12 is satisfied then}
\begin{equation*}
\left\Vert \hat{M}_{k\ell }-M_{k}\right\Vert _{\infty }=O_{p}(n^{-d_{\gamma
}}).
\end{equation*}
\bigskip
This result can be utilized to obtain mean square convergence rates for $
\hat{\alpha}_{k}$ from Theorems 1 and 2. As for linear functionals the limit
$\bar{\alpha}_{k}$ of the estimators $\hat{\alpha}_{k}$ are important for
the properties of $\hat{\theta}$. Here the $\bar{\alpha}_{k}$ are associated
with the Gateaux derivatives $D_{k}(W,\gamma_{k},\bar{\gamma}),$ $
(k=1,...,K).$ The following condition specifies each $\bar{\alpha}_{k}$ and
specifies the size of the remainder in a linearization using the Gateaux
derivatives.
\bigskip
\textsc{Assumption 13: }\textit{i) For }$(k=1,...,K)$\textit{\ there is }$
\bar{\alpha}_{k}\in\Gamma_{k}$\textit{\ such that for all }$
\gamma_{k}\in\Gamma_{k}$\textit{, }$\mathrm{E}[D_{k}(W,\gamma_{k},\bar{\gamma
})]=\mathrm{E}[\bar{\alpha }_{k}(X_{k})\gamma_{k}(X_{k})];$\textit{\ ii) }$
\bar{\alpha}_{k}(X_{k})$ \textit{and }$\mathrm{E}[\{Y_{k}-\bar{\gamma}
_{k}(X_{k})\}^{2}|X_{k}]$\textit{\ are bounded;} \textit{iii) there are }$
\varepsilon,$\textit{\ }$C>0$\textit{\ such that for all }$
\gamma\in\Pi_{k=1}^{K}\Gamma_{k}$\textit{\ with }$\Vert \gamma-\bar{\gamma}
\Vert<\varepsilon,$
\begin{equation*}
|\mathrm{E}[m(W,\gamma)-m(W,\bar{\gamma})-\sum_{k=1}^{K}D_{k}(W,\gamma_{k}-
\bar{\gamma },\bar{\gamma})]|\leq C\Vert\gamma-\bar{\gamma}\Vert^{2}.
\end{equation*}
\bigskip
Here each $\bar{\alpha}_{k}$ is specified as the Riesz representer for the
linear functional $\mathrm{E}[D_{k}(W,\gamma_{k},\bar{\gamma})]$ on $
\gamma_{k}\in\Gamma_{k}$ as in Newey (1994, equation 4.4). Here the
linearization $\mathrm{E}[D_{k}(W,\gamma_{k},\bar{\gamma})]$ has the role
that was fulfilled by the linear functional $\mathrm{E}[m(W,\gamma)]$
earlier. Indeed when $m(W,\gamma)$ is linear then $m(W,\gamma)$ will be its
Gateaux derivative.
From Lemma 8 we see that the convergence rate for each $\hat{M}_{k\ell }$ is
the convergence rate $n^{-d_{\gamma }}$ of $\hat{\gamma}$ rather than $\sqrt{
\ln (p)/n}.$ Consquently conditions for root-n consistency are different in
the nonlinear $m(W,\gamma )$ case than in the linear one. The following
condition imposes the rate conditions for a nonlinear functional.
\bigskip
\textsc{Assumption 14:} \textit{There is }$1/4<d_{\gamma}<1/2$\textit{\ such
that }$\left\Vert \hat{\gamma}_{k}-\bar{\gamma}_{k}\right\Vert
=O_{p}(n^{-d_{\gamma}}),$\textit{\ }$(k=1,...,K)$\textit{\ and for }$\bar{
\alpha}=\bar{\alpha}_{k}$\textit{\ and }$b(x)=b_{k}(x_{k})$\textit{, either
i) Assumptions 2 and 3 are satisfied and }$d_{\gamma}(1+4\xi)/(1+2\xi )>1/2$
\textit{\ or ii) Assumption 7 is satisfied and }$d_{\gamma}>1/3.$\textit{\ }
\bigskip
The requirement $d_{\gamma}>1/4$ given here is familiar for estimators that
depend nonlinearly on unknown functions, e.g. Newey (1994).. Condition i)
allows $d_{\gamma}$ to be any rate greater than $1/4$ if $\xi$ is large
enough. Condition ii), which drops the sparse eigenvalue assumption but
requires absolute summability of the coefficients of each $\bar{\alpha}_{k},$
requires $d_{\gamma}>1/3.$
The following gives the large sample inference results for $\hat{\theta}$
and $\hat{V}.$ Define
\begin{equation*}
\bar{\theta}=\mathrm{E}[m(W,\bar{\gamma})],\text{ }\psi(w)=m(w,\bar{\gamma})-
\bar{\theta}+\sum_{k=1}^{K}\bar{\alpha}_{k}(x_{k})[y_{k}-\bar{\gamma}
_{k}(x)],\text{ }V=\mathrm{E}[\psi(W)^{2}].
\end{equation*}
Here $\bar{\theta}$ will be the object estimated by $\hat{\theta}$ for $\bar{
\gamma}=$plim$(\hat{\gamma}).$
\bigskip
\textsc{Theorem 9}: \textit{If for }$\Gamma =\Gamma _{k}$\textit{, }$
b(x)=b_{k}(x_{k})$\textit{, }$r=r_{k}$\textit{\ for }$(k=1,...,K)$\textit{\
and }$\varepsilon _{n}=n^{-d_{\gamma }}$ Assumptions 1, 4, 5, 10, and 12-14
are satisfied \textit{then }$\sqrt{n}(\hat{\theta}-\bar{\theta})\overset{d}{
\longrightarrow }N(0,V).$ \textit{If in addition Assumption 7 is satisfied
for }$\bar{\alpha}=\bar{\alpha}_{k}$\textit{\ and }$b=b_{k}$ \textit{for
each }$(k=1,...,K)$\textit{\ then }$\hat{V}\overset{p}{\longrightarrow }V$.
\bigskip
\textsc{Example 6: }It is straightforward to specify regularity conditions
for causal mediation that are sufficient for the conditions of Theorem 9 to
hold.
\bigskip
\textsc{Assumption 15:} $\bar{\gamma}_{k}(X_{k})$ \textit{is bounded }$
(k=1,...,K),$\textit{\ there is }$C>0$\textit{\ such that }$\Pr(D=d,Q=q|Z)>C$
\textit{\ for all }$d\in\{0,1\},$\textit{\ }$q\in\{1,...,K-1\},$ \textit{and
}$\mathrm{E}[\{Y-\bar{\gamma}_{K}(D,Q,Z)\}^{2}|D,Q,Z]\leq C.$
\bigskip
This condition is used to guarantee that $\bar{\alpha}_{k}(X_{k})$ is
bounded for each $k.$ For brevity the form of $\bar{\alpha}_{k}(X_{k})$ and $
\psi(w)$ is given in the Appendix
\bigskip
\textsc{Corollary 10:} \textit{If for }$\Gamma =\Gamma _{k}$\textit{, }$
b(x)=b_{k}(x_{k})$\textit{, }$r=r_{k}$\textit{\ for }$(k=1,...,K)$\textit{\
and }$\varepsilon _{n}=n^{-d_{\gamma }}$\textit{\ Assumptions 1, 4, 5, 14,
and 15 are satisfied and there is }$C>0$\textit{\ such that }$\left\vert
\hat{\gamma}_{k}(x_{k})\right\vert \leq C$\textit{\ for all }$x_{k}$ \textit{
then }$\sqrt{n}(\hat{\theta}-\bar{\theta})\overset{d}{\longrightarrow }
N(0,V).$\textit{\ If in addition Assumption 7 is satisfied for }$\bar{\alpha}
=\bar{\alpha}_{k}$\textit{\ and }$b=b_{k}$ \textit{for each }$(k=1,...,K)$
\textit{\ then }$\hat{V}\overset{p}{\longrightarrow }V$.
\bigskip
The conditions of this result are simple relative to the general regularity
conditions in Assumptions 12 and 13. This simplicity is facilitated by $
m(W,\gamma )$ being quadratic in $\gamma .$ The condition that $\left\vert
\hat{\gamma}_{k}(x_{k})\right\vert \leq C$ is not strong for $k=1,...,K-1$
because $Y_{ki}\in \{0,1\}$. For $k=K$ this restriction could be imposed by
truncating $\hat{\gamma}_{k}(x)$ for some $C$ larger than a known bound on $
\gamma _{K}(X_{k})$ without affecting Assumption 14. In this way Corollary
10 provides a quite simple set of conditions for Auto-DML\ of causal
mediation effects.
\section{Regression Decomposition and the Average Treatment Effect on the
Treated}
In this Section we consider regression decompositions and the average
treatment effect on the treated (ATET). We also give an empirical
application of the ATET using Auto-DML.
\bigskip
\textsc{Example 6:} (Regression Decomposition and ATET): The effect of some
dummy variable $D\in\{0,1\}$ on an outcome variable $Y$ is often of
interest. Regression analysis can be used to decompose the unconditional
effect into an effect conditional on covariates and an effect from a shift
in the covariate distribution when $D$ shifts. One such decomposition takes
the form
\begin{align*}
\mathrm{E}[Y|D & =1]-\mathrm{E}[Y|D=0]=\Delta_{response}+
\Delta_{composition}, \\
\Delta_{response} & =\mathrm{E}[Y|D=1]-\frac{\mathrm{E}[D\gamma_{0}(0,Z)]}{
\Pr(D=1)},\text{ }\Delta_{composition}=\frac{\mathrm{E}[D\gamma_{0}(0,Z)]}{
\Pr(D=1)}-\mathrm{E}[Y|D=0],
\end{align*}
where $\gamma_{0}(D,Z)=\mathrm{E}[Y|D,Z]$. We will focus here on the
response effect
\begin{equation*}
\theta_{0}=\Delta_{response}=\frac{\mathrm{E}[D\gamma_{0}(1,Z)]-\mathrm{E}
[D\gamma_{0}(0,Z)]}{\Pr(D=1)}=\frac{\mathrm{E}[D\{\gamma_{0}(1,Z)-
\gamma_{0}(0,Z)\}]}{\Pr(D=1)}.
\end{equation*}
This $\theta_{0}$ is the average effect of changing $D$ on the outcome $Y$
conditional on $Z,$ averaged over the subpopulation with $D=1.$ One could
also consider a corresponding effect on the subpopulation with $D=0.$ That
could also be estimated using Auto-DML similarly to $\theta_{0}$ but for
brevity we omit this discussion.
This $\theta_{0}$ is also the ATET when $D$ is a treatment indicator and
potential outcomes are mean independent of treatment conditional on
covariates $Z$. Thus the estimator $\hat{\theta}$ and the asymptotic
variance estimator $\hat{V}$ we give could be applied for inference for the
ATET. We do so in the application given later in this Section.
The key regression functional of interest for $\theta_{0}$ is
\begin{align}
\mathrm{E}[D\gamma_{0}(0,Z)] & =\mathrm{E}[\pi_{0}(Z)\gamma_{0}(0,Z)]=
\mathrm{E}[\pi_{0}(Z)\frac {1-D}{1-\pi_{0}(Z)}\gamma_{0}(0,Z)]
\label{Rr ATET} \\
& =\mathrm{E}[\alpha_{0}(X)\gamma_{0}(X)],\text{ }\alpha_{0}(X)=\frac{
(1-D)\pi_{0}(Z)}{1-\pi_{0}(Z)}. \notag
\end{align}
Here $\alpha_{0}(X)$ is the Rr of a linear effect as in Section 2 with $
m(w,\gamma)=d\gamma(0,z).$ The condition $\mathrm{E}[\alpha_{0}(X)^{2}]<
\infty$ for a finite semiparametric variance bound is $\mathrm{E}
[1/\{1-\pi_{0}(Z)\}]<\infty.$
\bigskip
The effect $\theta_{0}=\Delta_{response}=ATET$ is a special case of the
nonlinear effect in Section 5 where $\gamma=(\gamma_{1},\gamma_{2})$, $
Y_{1}=Y,$ $X_{1}=(D,Z),$ $Y_{2}=D,$ $X_{2}=1$, and
\begin{equation*}
m(w,\gamma)=\frac{y_{2}}{\gamma_{2}}[y_{1}-\gamma_{1}(0,z)].
\end{equation*}
The orthogonal moment function for this object is
\begin{equation*}
\psi(w,\gamma,\gamma_{2},\alpha,\theta)=\frac{1}{\gamma_{2}}\{d[y-\gamma
(0,z)-\theta]-\alpha(x)[y-\gamma(x)]\},
\end{equation*}
where for notational convenience we let $y_{1}=y,$ $y_{2}=d,$ and $\gamma
_{1}=\gamma$. Similarly to Section 2 this moment function is doubly robust
in that
\begin{equation*}
\mathrm{E}[\psi(W,\bar{\gamma},\gamma_{20},\bar{\alpha},\theta_{0})]=0
\end{equation*}
if either $\bar{\gamma}(X)=\mathrm{E}[Y|X]$ or $\alpha_{0}(X)\in\Gamma$.
An Auto-DML is given by
\begin{align}
\hat{\theta} & =\frac{1}{n_{D}}\{\sum_{\ell=1}^{L}\sum_{i\in
I_{\ell}}\{D_{i}[Y_{i}-\hat{\gamma}_{\ell}(0,Z_{i})]-\hat{\alpha}
_{\ell}(X_{i})[Y_{i}-\hat{\gamma}_{\ell}(X_{i})]\}\}, \\
\hat{V} & =\frac{1}{n}\sum_{\ell=1}^{L}\sum_{i\in I_{\ell}}\hat{\psi}
_{i\ell}^{2},\text{ }\hat{\psi}_{i\ell}=(\frac{n}{n_{D}})\{D_{i}[Y_{i}-\hat{
\gamma}_{\ell}(0,Z_{i})-\hat{\theta}]-\hat{\alpha}_{\ell}(X_{i})[Y_{i}-\hat{
\gamma}_{\ell}(X_{i})]\}, \notag
\end{align}
where $n_{D}$ is the number of treated observations and $\hat{\alpha}_{\ell
}(x)$ is the Lasso learner of the Rr for $m(w,\gamma)=d\gamma(0,z).$
Similarly to the ATE\ in Example 3 we specify the dictionary to be $
b(x)=[dq(z)^{\prime },(1-d)q(z)^{\prime}]^{\prime},$ where $
q(z)=(z_{1},...,z_{p/2})^{\prime}$ when $\hat{\gamma}_{\ell}$ is high
dimensional and $q(z)$ is a vector of approximating functions when $\hat{
\gamma}_{\ell}$ is nonparametric. Then $m(w,b_{j})=d\cdot
b_{j}(0,z)=d\cdot1(j>p/2)q_{j-p/2}(z),$ so that
\begin{equation*}
\hat{M}_{\ell}=\frac{1}{n-n_{\ell}}\sum_{i\notin I_{\ell}}m(W_{i},b)=\left(
\begin{array}{c}
0 \\
\bar{q}_{\ell}
\end{array}
\right) ,\text{ }\bar{q}_{\ell}=\frac{1}{n-n_{\ell}}\sum_{i\notin
I_{\ell}}D_{i}q(Z_{i}).
\end{equation*}
Then by block diagonality of $\hat{G}_{\ell}$ and the first block of $\hat {M
}_{\ell}$ being zero
\begin{align*}
\hat{\alpha}_{\ell}(x) & =(1-d)q(z)^{\prime}\hat{\rho}_{\ell2},\text{ }\hat{
\rho}_{\ell2}=\arg\min_{\rho}\{-2\bar{q}_{\ell}^{\prime}\rho_{2}+\rho
_{2}^{\prime}\hat{G}_{\ell2}\rho_{2}+2r\left\Vert \rho_{2}\right\Vert _{1}\},
\\
\hat{G}_{\ell2} & =\frac{1}{n-n_{\ell}}\sum_{i\notin
I_{\ell}}(1-D_{i})q(Z_{i})q(Z_{i})^{\prime}.
\end{align*}
The first order conditions for the Lasso coefficients $\hat{\rho}_{\ell2}$
are
\begin{equation}
\left\vert \frac{1}{n-n_{\ell}}\sum_{i\notin
I_{\ell}}q_{j}(Z_{i})[D_{i}-(1-D_{i})\hat{\omega}_{\ell i}]\right\vert \leq
r,\text{ }\hat{\omega }_{\ell i}=q(Z_{i})^{\prime}\hat{\rho}_{\ell2},\text{ }
(j=1,...,p/2).
\end{equation}
The $\hat{\alpha}_{\ell}$ learner sets the \textquotedblleft
weights\textquotedblright\ $\hat{\omega}_{\ell i}$ to approximately
\textquotedblleft balance\textquotedblright\ the treated and untreated
averages for each element of $q(z).$
\bigskip
\textsc{Corollary 11}: \textit{If i) there is }$C>0$ \textit{with }$\pi
_{0}(Z)<1-C$\textit{\ and ii) Assumptions 1, 4, 5, and 9, 11 are satisfied
then for }$\bar{\theta}=\mathrm{E}[D\{Y-\bar{\gamma}(0,Z)\}]/\Pr (D=1)$
\textit{\ and }$\psi (W)=\Pr (D=1)^{-1}\{D[Y-\bar{\gamma}(0,Z)-\bar{\theta}]-
\bar{\alpha}(X)[Y-\bar{\gamma}(X)]\},$
\begin{equation*}
\sqrt{n}(\hat{\theta}-\theta _{0})\overset{d}{\longrightarrow }N(0,V),\text{
}V=\mathrm{E}[\psi (W)^{2}].
\end{equation*}
\textit{If Assumption 7 is also satisfied then} $\hat{V}\overset{p}{
\longrightarrow }V.$
\bigskip
As an empirical application, we use the Auto-DML of the ATET to estimate the
effect of job training in the National Supported Work Demonstration (NSW), a
job training program for disadvantaged workers that operated in the
mid-1970s. We follow the empirical strategy of LaLonde (1986) and Dehijia
and Wahba (1999), who compare the difference-in-means estimator applied to
an experimental data set with various econometric estimators applied to
\textquotedblleft quasi-experimental\textquotedblright\ data sets. The
experimental data set consists of the treatment and control groups from a
field experiment. A quasi-experimental data set consists of the treatment
group from a field experiment and a comparison group from an unrelated
national survey.
We use sample selection and variable construction as in Dehijia and Wahba
(1999) and Farrell (2015). The outcome $Y$ is earnings in 1978. The
treatment $D$ is an indicator of participation in job training. We consider
three specifications of covariates $Z$. We impose common support of the
propensity score for the treated and untreated groups based on covariates $Z$
as in Farrell (2015). Specifically, we calculate the range of propensity
scores for the treated group, and drop observations in the untreated group
whose propensity scores lie outside this range. We implement this procedure
for each of the three specifications (inducing three different propensity
scores), and ultimately keep the untreated observations that pass all three
tests. In estimation, we consider the fully-interacted dictionary $
b(D,Z)=(1,D,Z,DZ)$ for all three specifications of $Z$.
The covariate specifications are as follows.
\begin{enumerate}
\item Demographics and earnings, with quadratic terms of continuous
variables. In particular, the covariates are: age, education, black
indicator, Hispanic indicator, married indicator, 1974 earnings, 1975
earnings, age squared, education squared, 1974 earnings squared, and 1975
earnings squared. This specification is moderately flexible. It is one that
an analyst may reasonably implement without knowing the experimental
benchmark ex ante. Here $dim(Z)=11$ and $p=dim(b(D,Z))=24$.
\item Demographics and earnings, with quadratic terms of continuous
variables and constructed indicators. In particular, the covariates are:
those in specification 1; unemployed in 1974 indicator, unemployed in 1975
indicator, and no degree indicator. This specification includes some domain
knowledge about which signals employers may respond to while making hiring
decisions. Note that it does not include conveniently hand-crafted basis
functions to get closer to the experimental benchmark. Here $dim(Z)=14$ and $
p=dim(b(D,Z))=30$.
\item A high dimensional specification where the covariates are: those in
specification 2; all possible first order interactions, and all polynomials
up to order five of the continuous variables (age, education, 1974 earnings,
1975 earnings). This specification was introduced by Farrell (2015). Here $
dim(Z)=171$ and $p=dim(b(D,X))=344$.
\end{enumerate}
We estimate the Rr with Lasso minimum distance, and the regression with
Lasso minimum distance, random forests (RF), or neural networks (NN). For
Lasso minimum distance, we use the tuning procedure described in Appendix~
\ref{sec:computing}. We use the same settings of random forest as
Chernozhukov et al. (2018). We implement a neural network with two hidden
layers of eight units each and linear activation. We use $L=5$ folds in
cross-fitting.
Tables~\ref{tab:ATT_nsw2},~\ref{tab:ATT_psid2}, and~\ref{tab:ATT_cps2}
summarize results for the NSW, PSID, and CPS data sets, respectively. For
comparison, LaLonde (1986) reports $1794$ $(633)$ by difference-in-means
applied to the NSW data, which is the experimental benchmark. Farrell (2015)
reports $1737$ $(869)$ by group Lasso applied to the PSID data using
specification 3. Our corresponding estimate is $1763$ $(1026)$, and our
other results are broadly consistent. To validate the robustness of our
results with respect to the choice of tuning procedure, we report analogous
tables using cross validated regularization in Appendix~\ref
{sec:additional_empirics}.
\begin{table}[ptb]
\centering
\begin{tabular}{c|c|c|c|c|c|c|c|c}
\hline\hline
spec. & treated & untreated & Lasso ATET & Lasso SE & RF ATET & RF SE & NN
ATET & NN SE \\ \hline
1 & 185 & 172 & 3022.84 & 1278.54 & 3106.55 & 1327.02 & 2585.19 & 1183.85 \\
2 & 185 & 172 & 2959.72 & 1253.13 & 3077.26 & 1318.67 & 2606.43 & 1020.56 \\
3 & 185 & 172 & 2289.65 & 836.19 & 2785.13 & 819.17 & 2504.83 & 770.24 \\
\hline\hline
\end{tabular}
\caption{ATET using NSW treatment and NSW control, by Auto-DML}
\label{tab:ATT_nsw2}
\end{table}
\begin{table}[ptb]
\centering
\begin{tabular}{c|c|c|c|c|c|c|c|c}
\hline\hline
spec. & treated & untreated & Lasso ATET & Lasso SE & RF ATET & RF SE & NN
ATET & NN SE \\ \hline
1 & 185 & 727 & 900.58 & 873.62 & 1521.92 & 977.08 & 197.53 & 946.30 \\
2 & 185 & 727 & 1466.35 & 882.67 & 1336.66 & 956.22 & 1447.65 & 980.73 \\
3 & 185 & 727 & 1763.20 & 1026.09 & 2010.53 & 987.73 & 2698.55 & 1036.24 \\
\hline\hline
\end{tabular}
\caption{ATET using NSW treatment and PSID comparison, by Auto-DML}
\label{tab:ATT_psid2}
\end{table}
\begin{table}[ptb]
\centering
\begin{tabular}{c|c|c|c|c|c|c|c|c}
\hline\hline
spec. & treated & untreated & Lasso ATET & Lasso SE & RF ATET & RF SE & NN
ATET & NN SE \\ \hline
1 & 185 & 5904 & 703.21 & 583.23 & 1639.95 & 616.08 & 1686.77 & 611.81 \\
2 & 185 & 5904 & 971.46 & 583.48 & 1584.12 & 616.33 & 1094.86 & 590.31 \\
3 & 185 & 5904 & 1358.46 & 614.56 & 1906.62 & 651.77 & 2235.09 & 742.86 \\
\hline\hline
\end{tabular}
\caption{ATET using NSW treatment and CPS comparison, by Auto-DML}
\label{tab:ATT_cps2}
\end{table}
\section{Panel Average Derivative and Demand Elasticities}
In this Section, we apply Auto-DML to estimating demand elasticities while
allowing for individual preferences that are correlated with prices and
total expenditure. Specifically, we estimate own-price elasticity in a panel
data model with correlated random slopes. We apply this approach to Nielsen
scanner data.
A panel data model requires double indexing. Let $Y_{it}$, $
(t=1,...,T_{i},i=1,...,n)$, denote the share of total expenditure on some
good for household $i$ in time period $t$. Let $X_{it}$ be a vector of log
prices, log expenditure, and covariates. Let $\tilde{X}_{i}=(X_{i1}^{
\prime},...,X_{i,T_{i}}^{\prime})^{\prime}$ collect observations over all
time periods for individual $i$ into one vector. We allow for an unbalanced
panel where different households may have different numbers of observations $
T_{i}$ as in Wooldridge (2019).
Consider the demand model of Chernozhukov, Hausman, and Newey (2021) given
by
\begin{equation}
\mathrm{E}[Y_{it}|\tilde{X}_{i},B_{it}]=b_{1}(X_{it})^{\prime }B_{it}.
\label{eq:demand}
\end{equation}
The $K$-dimensional dictionary $b_{1}(X_{it})$ is a vector of functions of $
X_{it}$ that includes a constant and, for example, powers of log price and
log expenditure. $B_{it}$ represents household specific preferences that may
vary over time and that may be correlated with regressors from each time
period. We assume the conditional mean of $B_{it}$ is time stationary with
\begin{equation}
\mathrm{E}[B_{it}|\tilde{X}_{i}]=[I_{K}\otimes \tilde{H}_{i}]^{\prime }\pi
_{0},\quad \pi _{0}=(\pi _{10}^{\prime },...,\pi _{K,0}^{\prime })^{\prime },
\label{eq:stationary}
\end{equation}
where $I_{K}$ is a $K$-dimensional identity matrix. $\tilde{H}_{i}$ is a
vector of functions of $\tilde{X}_{i}$ with length that does not depend on $
T_{i}$. This panel model is like that of Chamberlain (1982, 1992),
Chernozhukov et al. (2013b), Graham and Powell (2012), and Wooldridge
(2019), as further discussed in Chernozhukov, Hausman, and Newey (2021).
We will consider identifying and estimating transformations of $\beta _{0}=
\mathrm{E}[B_{it}]$. $\beta_{0}$ is interpretable as the average marginal
effect of changing $b_{1}(X_{it})$. The transformations we consider will be
interpretable as average income, own-price, and cross-price elasticities. By
law of iterated expectations, our model implies
\begin{equation}
\beta_{0}=\mathrm{E}[B_{it}]=[I_{K}\otimes \mathrm{E}[\tilde{H}
_{i}]]^{\prime}\pi_{0}. \label{eq:LIE}
\end{equation}
Combining~(\ref{eq:demand}),~(\ref{eq:stationary}), and~(\ref{eq:LIE}), we
summarize the correlated random effects model as follows.
\begin{align}
\gamma_{0}(\tilde{X}_{i}) & =\mathrm{E}[Y_{it}|\tilde{X}_{i}] \notag \\
& =b_{1}(X_{it})^{\prime}\{\beta_{0}+\mathrm{E}[B_{it}|\tilde{X}_{i}]-\beta
_{0}\} \notag \\
& =b_{1}(X_{it})^{\prime}\{\beta_{0}+[I_{K}\otimes\tilde{H}
_{i}]^{\prime}\pi_{0}-[I_{K}\otimes \mathrm{E}[\tilde{H}_{i}]]^{\prime}
\pi_{0}\} \notag \\
& =b_{1}(X_{it})^{\prime}\beta_{0}+[b_{1}(X_{it})\otimes(\tilde{H}_{i}-
\mathrm{E}[\tilde{H}_{i}])]^{\prime}\pi_{0}. \label{eq:corrRE}
\end{align}
In summary, the choice of $K$-dimensional dictionary $b_{1}(X_{it})$ in the
demand model~(\ref{eq:demand}) induces a $p$-dimensional dictionary $
b_{it}=b(\tilde{X}_{i})=(b_{1}(X_{it})^{\prime},[b_{1}(X_{it})\otimes (
\tilde{H}_{i}-\mathrm{E}[\tilde{H}_{i}])]^{\prime})^{\prime}$ in the
correlated random effects model~(\ref{eq:corrRE}). In practice, we replace $
\mathrm{E}[\tilde{H}_{i}]$ with $\frac{1}{n}\sum_{i=1}^{n} \tilde{H}_{i}$
and set $\tilde{H}_{i}=\frac {1}{T_{i}}\sum_{t=1}^{T_{i}} b_{1}(X_{it})$.
\bigskip
\textsc{Example 10:} Demand elasticities. Denote $X_{it}=(D_{it},Z_{it})$
where $D_{it}$ is log own price. By the derivation in Chernozhukov et al.
(2019) for budget share regressions, an average own-price elasticity is
\begin{equation*}
\theta_{0}^{\ast}=\frac{\theta_{0}}{\mathrm{E}[Y]}-1,\quad\theta_{0}=\mathrm{
E} \left[ \frac{\partial\gamma_{0}(\tilde{X}_{i})}{\partial d}\right] .
\end{equation*}
Own-price elasticity $\theta_{0}^{\ast}$ is a smooth transformation of a
linear effect $\theta_{0}$, which in this case is average derivative.
Auto-DML of own-price elasticity is then given by
\begin{equation*}
\hat{\theta}^{\ast}=\frac{\hat{\theta}}{\frac{1}{n\sum_{i=1}^{n}T_{i}}
\sum_{i=1}^{n}\sum_{t=1}^{T_{i}}Y_{it}}-1
\end{equation*}
where $\hat{\theta}$ is the Auto-DML of average derivative from Example 4.
Income elasticity and cross-price elasticity have a similar structure; see
Appendix~\ref{sec:delta} for details.
For completeness, we present $\hat{M}_{\ell }$ for average derivative using
the panel data dictionary $b_{it}$.
\begin{equation*}
\hat{M}_{\ell }=\frac{1}{\sum_{i=1}^{n}T_{i}-\sum_{i\in I_{\ell }}T_{i}}
\sum_{i\not\in I_{\ell }}\sum_{t=1}^{T_{i}}\frac{\partial b_{it}}{\partial d}
=\frac{1}{\sum_{i=1}^{n}T_{i}-\sum_{i\in I_{\ell }}T_{i}}\sum_{i\not\in
I_{\ell }}\sum_{t=1}^{T_{i}}\left(
\begin{array}{c}
\frac{\partial b_{1}(X_{it})}{\partial d} \\
0
\end{array}
\right) .
\end{equation*}
Recall Theorem 4 provides consistency and asymptotic normality guarantees
for Auto-DML $\hat{\theta}$. A more sophisticated estimator $\hat{V}$ of the
asymptotic variance of $\hat{\theta}$ is required that accounts for
clustering of observations by household. See the Appendix~\ref{sec:delta}
for details. Importantly, the cluster structure is also preserved in
cross-fitting. Clustering methods for DML were previously used by Chiang et
al. (2019) and Chernozhukov, Hausman, and Newey (2021). The consistency of
the own-price elasticity $\hat{\theta}^{\ast }$ follows from the continuous
mapping theorem, and the asymptotic normality of $\hat{\theta}^{\ast }$
follows from delta method.
As an empirical application, we apply Auto-DML to estimate own-price
elasticity of milk and soda with Nielsen scanner data. The empirical work
here is the researchers' own analyses calculated (or derived) based in part
on data from Nielsen Consumer LLC and marketing databases provided through
the NielsenIQ Datasets at the Kilts Center for Marketing Data Center at The
University of Chicago Booth School of Business. The conclusions drawn from
the NielsenIQ data are those of the researchers and do not reflect the views
of NielsenIQ. NielsenIQ is not responsible for, had no role in, and was not
involved in analyzing and preparing the results reported herein.
The data we use are a subset of the Nielsen Homescan Panel as in Burda,
Harding, and Hausman (2008, 2012). The data include 1483 households from the
Houston-area zip codes for the years 2004-2006. The number of monthly
observations for each household ranges from 12 to 36, with some households
being added and taken away throughout the three years covered. 609
households are included the entire time. Expenditures include all purchases
of the household in each month. The original data had time stamps for
purchases. If a household purchased a good more than once in a month, the
\textquotedblleft monthly price\textquotedblright\ is the average price that
the household paid (i.e. total amount spent on good/total quantity
purchased). We include observations with zero expenditure share as justified
in Chernozhukov, Hausman, and Newey (2021). For those observations, $
Y_{it}=0 $ and own price is imputed in the ways described in Chernozhukov,
Hausman, and Newey (2021).
We consider 15 groups of goods: bread, butter, cereal, chips, coffee,
cookies, eggs, ice cream, milk, orange juice, salad, soda, soup, water, and
yogurt. As in Burda, Harding, and Hausman (2008, 2012), we choose these
groups because they make up a relatively large proportion of total food
expenditure. We consider budget share regressions for two of these goods:
milk and soda. $Y_{it}$ is share of expenditure spent on milk (soda) by
household $i$ in month $t$. We take as $b_{1}(X_{it})$ the concatenation of
the following variables: fourth order polynomial of log expenditure; fourth
order polynomial of log price for milk (soda); up to fourth order
interactions thereof; and log price of other goods. For $\tilde{H}_{i}$, we
use the time averages of $b_{1}(X_{it})$. Note that $K=dim(b_{1}(X_{it}))=42$
and $p=1521$.
We estimate own-price elasticity according to the procedure outlined
previously in this Section. We estimate both the Rr and the regression with
Lasso minimum distance. For Lasso minimum distance, we use the tuning
procedure described in Appendix~\ref{sec:computing}. We use $L=5$ folds in
cross-fitting. We calculate clustered standard errors by delta method, as
described in Appendix~\ref{sec:delta}.
Table~\ref{tab:dml_elasticity} summarizes results for the milk and soda
own-price elasticities using Auto-DML. For comparison, the cross sectional
estimates for milk and soda elasticities are $-1.27$ ($0.0163$) and $-0.859$
($0.00485$), respectively (Table 1 of Chernozhukov, Hausman, and Newey 2021)
and the corresponding fixed effects estimates are $-.739$ $(.0197)$ and $
-.853$ ($.00517$). Our results show that allowing for correlated random
coefficients lowers these elasticity estimates, especially the milk
elasticity. These results confirm the finding in Table 5 of Chernozhukov,
Hausman, and Newey (2019), that panel elasticity estimates allowing for
correlation of preferences with prices and total expenditure are much
smaller than cross-section estimates for milk. Our own-price elasticity
estimates are not as small as their slope fixed effect estimates, which for
milk are between $-0.626$ $(0.00849)$ and $-0.496$ $(0.0479)$ and for soda
are between $-0.805$ $(0.00830)$ and $-0.780$ $(0.0235)$ depending on choice
of regularization parameter.
\begin{table}[ptb]
\centering
\begin{tabular}{c|c|c}
\hline\hline
good & elasticity & SE \\ \hline
milk & -.645 & .00649 \\
soda & -.826 & .00379 \\
\hline\hline
\end{tabular}
\vspace{8pt}
\caption{Average own-price elasticity, by Auto-DML}
\label{tab:dml_elasticity}
\end{table}
For further comparison, we report results from the plug-in approach in Table~
\ref{tab:dml_elasticity_bias}. The plug-in elasticity estimates are much
closer to the cross-section estimates than the Auto-DML estimates. The
results of this table confirm the importance of debiasing in this
application, with debiased estimates differing from plug-in estimates by
much more than the associated standard errors.
\begin{table}[ptb]
\centering
\begin{tabular}{c|c|c}
\hline\hline
good & elasticity & SE \\ \hline
milk & -.863 & .00255 \\
soda & -.863 & .00305 \\
\hline\hline
\end{tabular}
\vspace{8pt}
\caption{Average own-price elasticity, by plug-in}
\label{tab:dml_elasticity_bias}
\end{table}
\section{Conclusions}
In this paper we have given an automatic method of debiasing a machine
learner of a parameter of interest that depends on a high dimensional and/or
nonparametric regression. The method only requires the form of the object of
interest. The regression learners are allowed to be anything that converges
in mean square at a fast enough rate. We have shown root-n consistency and
asymptotic normality and given a consistent asymptotic variance estimator
for a wide variety of causal and structural estimators, including nonlinear
functionals of regression. We have applied these methods to estimate the
average treatment effect on the treated in a job training experiment and
have found similar results for Lasso, neural nets, and random forests
regressions. We also have also estimated a correlated random slopes
specification for consumer demand from scanner data and found estimates that
are similar to fixed slope effect elasticities.\newpage