Automatic Debiased Machine Learning of Causal and Structural Effects

Automatic Debiased Machine Learning of Causal and Structural Effects

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.

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

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

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:

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

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.

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

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

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,

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.$

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

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).$

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

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.$

We focus on $m(w,\gamma)$ where there exists a function $\alpha_{0}(X)$ with $\mathrm{E}[\alpha_{0}(X)^{2}]<\infty$ and

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) summarizes the functional $m(w,\gamma)$ and the Rr in each of the examples:

Equation ((ref)) 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)) implies that the effect of interest can be represented in three different ways, as

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

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

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

which satisfies equation ((ref)) 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)) 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)) 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

where $\bar{\gamma}(X)$ is the solution to equation ((ref)) 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),

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

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)) with $F=F_{0}$ are that $\mathrm{E}[\alpha(X)\{Y-\bar{\gamma} (X)\}]=0$ for all $\alpha\in\Gamma.$ Therefore

where the second equality follows by equation ((ref)) and the third equality by the necessary and sufficient condition for equation ((ref)) 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

does not depend on $\gamma$. Also equation ((ref)) 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)). 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)) 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

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:

\

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. $

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. $

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)) 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

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].$

For brevity we omit further discussion of Example 4 from the paper and refer the interested reader to Chernozhukov, Hausman, and Newey (2021).

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 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

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)), 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

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

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)) 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)). It is also notationally convenient to drop the $0$ mean normalization of $b(x)$ and consider $\hat{\alpha}$ having the form

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.$

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. $

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)). 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

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), 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)) and $ \bar{\alpha}$ satisfying equation ((ref)), which imply that for $ m(w,b)=(m(w,b_{1}),...,m(w,b_{p}))^{\prime}$,

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

be unbiased estimators of $M$ and $G.$ The coefficient estimator is given by

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)) 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)) 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))

Replacing the expectation by a sample average and $\bar{\alpha}(X)$ by $ b(X)^{\prime}\rho$ gives

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)) without the penalty term. The first-order conditions for equation ((ref)) 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.

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)). 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$,

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}).$

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$,

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}).$

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

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

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

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,$

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)) 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.

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}]}$

Assumption 2: There exists $C>1,$ $\xi>0$\ such that for each positive integer $s\leq C\varepsilon_{n}^{-2/(2\xi+1)}$ \ there is $\bar{\rho}$\ with $s$\ nonzero elements such that

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 .$

Assumption 3: $G=\mathrm{E}[b(X)b(X)^{\prime}]$ has largest eigenvalue bounded uniformly in $n$\ and there is $C,c>0$\ such that for all $s\approx C\varepsilon_{n}^{-2}$\ with probability approaching one

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.

Assumption 4: There is $C>0$\ such that with probability one $\sup_{j}|b_{j}(X)|\leq C.$

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.$

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}.$

We also hypothesize a convergence rate for $\hat{M}.$

Assumption 6: $\left\Vert \hat{M}-M\right\Vert _{\infty}=O_{p}(\varepsilon_{n})$ for $\varepsilon_{n}\longrightarrow0.$

We use this condition to accommodate $\hat{M}$ that can depend on the regression learner $\hat{\gamma}$ as needed for Section 5.

Theorem 1: If Assumptions 1 - 6 are satisfied then for all $c>0,$

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:

Assumption 7: $\bar{\alpha}(X)=\sum_{j=1}^{\infty}\rho_{j0}b_{j}(X)$ , $\sum_{j=1}^{\infty}\left\vert \rho_{j0}\right\vert <\infty $ , and for $C>0$\ and $\bar{s}=C\sqrt{n}$\ the $ b_{j}(x)$\ corresponding to the largest $\bar{s}$\textit{\ values of }$\left\vert \rho_{j0}\right\vert $\textit{\ are included in }$b(x).$

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} .$

Theorem 2: If Assumptions 1 and 5-7 are satisfied and $ \left\Vert \hat{G}-G\right\Vert _{\infty}=O_{p}(\varepsilon_{n})$\ then for all $c>0,$

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.

Assumption 8: There is $C>0$\ such that with probability one $\max_{j\leq p}|m(W,b_{j}))|\leq C.$

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.

Assumption 9:\ $\mathrm{E}[\{Y-\bar{\gamma}(X)\}^{2}|X]$ \ and $\bar{\alpha}(X)$\ are bounded.

We impose this condition for simplicity; it could be weakened. We also impose the following condition.

Assumption 10: $\mathrm{E}[m(W,\gamma_{0})^{2}]<\infty$ and $\int[m(w,\hat{\gamma})-m(w,\bar{\gamma})]^{2}F_{W}(dw)\overset{p}{ \longrightarrow}0.$

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.

Assumption 11: With probability approaching one $\hat { \gamma}_{\ell}\in\Gamma$\ and there is $d_{\gamma}>0$\ such that $\Vert\hat{\gamma}-\bar{\gamma}\Vert=O_{p}(n^{-d_{\gamma}})$\ and either Assumptions 2 and 3 are satisfied with

or Assumption 7 is satisfied and $d_{\gamma}>1/4.$

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

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.

Theorem 3: If Assumptions 1-5, and 8-11 are satisfied then $\sqrt{n}(\hat{\theta}-\bar{\theta})\overset{d}{\longrightarrow}N(0,V).$ If in addition Assumption 7 is satisfied then $\hat{V}\overset{p}{ \longrightarrow}V$.

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.

Corollary 4: If Assumptions 1-5, and 8-11 are satisfied, $ \mathrm{E}[m(W,\gamma )^{2}]\leq C\left\Vert \gamma \right\Vert^{2} $ for all $\gamma \in \Gamma ,$\ and $\bar{\alpha}(X)\neq 0$\ then $ \sqrt{n}(\tilde{\theta}-\bar{\theta})\overset{d}{\longrightarrow }N(0,V).$

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.

Corollary 5 (Example 1): If Assumptions 1-5, 9, and 11 are satisfied and there is $C>0$\ such that $\left\vert [f_{1}(x)-f_{0}(x)]/f(x)\right\vert \leq C$ then $\sqrt{n}(\hat{ \theta}-\bar{\theta})\overset{d}{\longrightarrow }N(0,V).$\ If in addition Assumption 7 is satisfied then $\hat{V}\overset{p}{\longrightarrow }V$.

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.

Corollary 6 (Example 2): If Assumptions 1-5, 9, and 11 are satisfied and there is $C>0$\ such that $\left\vert S(u)\right\vert \leq C$, $f(D|Z)^{-1}\omega (D)\leq C$\ then $\sqrt{ n}(\hat{\theta}-\bar{\theta})\overset{d}{\longrightarrow }N(0,V).$\ If in addition Assumption 7 is satisfied then $\hat{V}\overset{p}{ \longrightarrow }V$.

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.

Corollary 7 (Example 3): If Assumptions 1, 4-5, 9, and 11 are satisfied and there is $C>0$\ with $\pi _{0}(Z)\in \lbrack C,1-C]$\ then $\sqrt{n}(\hat{\theta}-\bar{\theta})\overset{d}{ \longrightarrow }N(0,V).$\ If in addition Assumption 7 is satisfied then $\hat{V}\overset{p}{\longrightarrow }V$.

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.

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)). Specifically it is

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

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)) 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)) rather than the one in equation ( (ref)).

The $\hat{M}_{k\ell j}$ given here generalizes equation ((ref)) 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)). 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)) is automatic, in the same way as the estimator of equation ((ref)), 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})$.

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

The causal mediation effect of Imai, Keele, and Tingley (2010, Theorem 1) is

This effect, or parameter, has the form $\theta_{0}(d,d^{\prime})=\mathrm{E} [m(W,\gamma)]$ for $W=(Y,D,Q,Z)$ and

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

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

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,$

Then we have

We can then compute $\hat{\alpha}_{k\ell}(x)$ as in equation ((ref) ) and $\hat{\theta}$ for $Y_{ki}=1(Q_{i}=k),$ $(k=1,...,K-1)$ and $ Y_{Ki}=Y_{i}$ as in equation ((ref)).

The orthogonal moment function corresponding to this estimator is

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},$

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},$

for any $\gamma _{K}\in \Gamma _{K}$ and $\gamma _{k}\in \Gamma _{1},$ $ (k\leq K-1).$

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

be the Gateaux derivative of $m(W,\gamma)$ with respect to $\gamma_{k}$ when it exists. Comparing this definition with equation ((ref)) 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.

Assumption 12: There are $C,$\ $\varepsilon >0,$ \ $a_{kj}(w),$\ and $A_{k}(w,\gamma)$\ 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$

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).$

Lemma 8: 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)$, and Assumption 12 is satisfied then

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.

Assumption 13: i) For $(k=1,...,K)$\ there is $ \bar{\alpha}_{k}\in\Gamma_{k}$\ such that for all $ \gamma_{k}\in\Gamma_{k}$, $\mathrm{E}[D_{k}(W,\gamma_{k},\bar{\gamma })]=\mathrm{E}[\bar{\alpha }_{k}(X_{k})\gamma_{k}(X_{k})];$\ 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,$

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.

Assumption 14: There is $1/4<d_{\gamma}<1/2$\ such that $\left\Vert \hat{\gamma}_{k}-\bar{\gamma}_{k}\right\Vert =O_{p}(n^{-d_{\gamma}}),$\ $(k=1,...,K)$\ and for $\bar{ \alpha}=\bar{\alpha}_{k}$\ 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{\ }

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

Here $\bar{\theta}$ will be the object estimated by $\hat{\theta}$ for $\bar{ \gamma}=$plim$(\hat{\gamma}).$

Theorem 9: If for $\Gamma =\Gamma _{k}$, $ b(x)=b_{k}(x_{k})$, $r=r_{k}$\ for $(k=1,...,K)$\ 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$.

Example 6: It is straightforward to specify regularity conditions for causal mediation that are sufficient for the conditions of Theorem 9 to hold.

Assumption 15: $\bar{\gamma}_{k}(X_{k})$ is bounded $ (k=1,...,K),$\ there is $C>0$\ such that $\Pr(D=d,Q=q|Z)>C$ \ for all $d\in\{0,1\},$\ $q\in\{1,...,K-1\},$ \textit{and }$\mathrm{E}[\{Y-\bar{\gamma}_{K}(D,Q,Z)\}^{2}|D,Q,Z]\leq C.$

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

Corollary 10: If for $\Gamma =\Gamma _{k}$, $ b(x)=b_{k}(x_{k})$, $r=r_{k}$\ for $(k=1,...,K)$\ 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$.

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.

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.

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

where $\gamma_{0}(D,Z)=\mathrm{E}[Y|D,Z]$. We will focus here on the response effect

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

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.$

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

The orthogonal moment function for this object is

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

if either $\bar{\gamma}(X)=\mathrm{E}[Y|X]$ or $\alpha_{0}(X)\in\Gamma$.

An Auto-DML is given by

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

Then by block diagonality of $\hat{G}_{\ell}$ and the first block of $\hat {M }_{\ell}$ being zero

The first order conditions for the Lasso coefficients $\hat{\rho}_{\ell2}$ are

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).$

Corollary 11: If i) there is $C>0$ with $\pi _{0}(Z)<1-C$\ 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)$ \ and $\psi (W)=\Pr (D=1)^{-1}\{D[Y-\bar{\gamma}(0,Z)-\bar{\theta}]- \bar{\alpha}(X)[Y-\bar{\gamma}(X)]\},$

If Assumption 7 is also satisfied then $\hat{V}\overset{p}{ \longrightarrow }V.$

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.

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). 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), (ref), and (ref) 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).

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

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

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

Combining ((ref)), ((ref)), and ((ref)), we summarize the correlated random effects model as follows.

In summary, the choice of $K$-dimensional dictionary $b_{1}(X_{it})$ in the demand model ((ref)) 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)). 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})$.

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

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

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) for details.

For completeness, we present $\hat{M}_{\ell }$ for average derivative using the panel data dictionary $b_{it}$.

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) 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). We use $L=5$ folds in cross-fitting. We calculate clustered standard errors by delta method, as described in Appendix (ref).

Table (ref) 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.

For further comparison, we report results from the plug-in approach in Table (ref). 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.

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.