Double/debiased machine learning for treatment and structural parameters

Summary

We revisit the classic semi‐parametric problem of inference on a low‐dimensional parameter θ₀ in the presence of high‐dimensional nuisance parameters η₀. We depart from the classical setting by allowing for η₀ to be so high‐dimensional that the traditional assumptions (e.g. Donsker properties) that limit complexity of the parameter space for this object break down. To estimate η₀, we consider the use of statistical or machine learning (ML) methods, which are particularly well suited to estimation in modern, very high‐dimensional cases. ML methods perform well by employing regularization to reduce variance and trading off regularization bias with overfitting in practice. However, both regularization bias and overfitting in estimating η₀ cause a heavy bias in estimators of θ₀ that are obtained by naively plugging ML estimators of η₀ into estimating equations for θ₀. This bias results in the naive estimator failing to be $N^{-1/2}$ consistent, where N is the sample size. We show that the impact of regularization bias and overfitting on estimation of the parameter of interest θ₀ can be removed by using two simple, yet critical, ingredients: (1) using Neyman‐orthogonal moments/scores that have reduced sensitivity with respect to nuisance parameters to estimate θ₀; (2) making use of cross‐fitting, which provides an efficient form of data‐splitting. We call the resulting set of methods double or debiased ML (DML). We verify that DML delivers point estimators that concentrate in an $N^{-1/2}$‐neighbourhood of the true parameter values and are approximately unbiased and normally distributed, which allows construction of valid confidence statements. The generic statistical theory of DML is elementary and simultaneously relies on only weak theoretical requirements, which will admit the use of a broad array of modern ML methods for estimating the nuisance parameters, such as random forests, lasso, ridge, deep neural nets, boosted trees, and various hybrids and ensembles of these methods. We illustrate the general theory by applying it to provide theoretical properties of the following: DML applied to learn the main regression parameter in a partially linear regression model; DML applied to learn the coefficient on an endogenous variable in a partially linear instrumental variables model; DML applied to learn the average treatment effect and the average treatment effect on the treated under unconfoundedness; DML applied to learn the local average treatment effect in an instrumental variables setting. In addition to these theoretical applications, we also illustrate the use of DML in three empirical examples.

1. Introduction and Motivation

Motivation

We develop a series of simple results for obtaining root‐N consistent estimation, where N is the sample size, and valid inferential statements about a low‐dimensional parameter of interest, θ₀, in the presence of a high‐dimensional or ‘highly complex’ nuisance parameter, η₀. The parameter of interest will typically be a causal parameter or treatment effect parameter, and we consider settings in which the nuisance parameter will be estimated using machine learning (ML) methods, such as random forests, lasso or post‐lasso, neural nets, boosted regression trees, and various hybrids and ensembles of these methods. These ML methods are able to handle many covariates and they provide natural estimators of nuisance parameters when these parameters are highly complex. Here, highly complex formally means that the entropy of the parameter space for the nuisance parameter is increasing with the sample size in a way that moves us outside the traditional framework considered in the classical semi‐parametric literature where the complexity of the nuisance parameter space is taken to be sufficiently small. The main contribution of this paper is to offer a general and simple procedure for estimating and to perform inference on θ₀ that is formally valid in these highly complex settings.  

Example 1.1. Partially Linear Regression

As a lead example, consider the following partially linear regression (PLR) model as in Robinson (1988):  

$$\begin{array}{ccc}& Y=D{\theta }_0+g_0(X)+U,& E[U\mid X,D]=0,\hfill \end{array}$$

(1.1)

 

$$\begin{array}{ccc}& D=m_0(X)+V,& E[V\mid X]=0.\hfill \end{array}$$

(1.2)

Here, Y is the outcome variable, D is the policy/treatment variable of interest, vector  

$$X=(X_1,\dots ,X_p)$$

consists of other controls, and U and V are disturbances.¹ The first equation is the main equation, and θ₀ is the main regression coefficient that we would like to infer. If D is exogenous conditional on controls X, θ₀ has the interpretation of the treatment effect parameter or ‘lift’ parameter in business applications. The second equation keeps track of confounding, namely the dependence of the treatment variable on controls. This equation is not of interest per se but it is important for characterizing and removing regularization bias. The confounding factors X affect the policy variable D via the function $m_0(X)$ and the outcome variable via the function $g_0(X)$⁠. In many applications, the dimension p of vector X is large relative to N. To capture the feature that p is not vanishingly small relative to the sample size, modern analyses then model p as increasing with the sample size, which causes traditional assumptions that limit the complexity of the parameter space for the nuisance parameters ${\eta }_0=(m_0,g_0)$ to fail.

Regularization Bias

A naive approach to estimation of θ₀ using ML methods would be, for example, to construct a sophisticated ML estimator $D{\widehat{\theta }}_0+{\widehat{g}}_0(X)$ for learning the regression function $D{\theta }_0+g_0(X)$⁠.² Suppose, for the sake of clarity, that we randomly split the sample into two parts: a main part of size n, with observation numbers indexed by $i\in I$⁠, and an auxiliary part of size $N-n$⁠, with observations indexed by $i\in I^c$⁠. For simplicity, we take $n=N/2$ for the moment and we turn to more general cases that cover unequal split‐sizes, using more than one split, and achieving the same efficiency as if the full sample were used for estimating θ₀ in the formal development in Section 3. Suppose ${\widehat{g}}_0$ is obtained using the auxiliary sample and that, given this ${\widehat{g}}_0$⁠, the final estimate of θ₀ is obtained using the main sample:  

$${\displaystyle {\widehat{\theta }}_0={\left(\frac{1}{n}\sum _{i\in I}{D}_i^2\right)}^{-1}\frac{1}{n}\sum _{i\in I}{D}_i(Y_i-{\widehat{g}}_0(X_i)).}$$

(1.3)

The estimator ${\widehat{\theta }}_0$ will generally have a slower than $1/\sqrt{n}$ rate of convergence, namely,  

$$|\sqrt{n}({\widehat{\theta }}_0-{\theta }_0)|\stackrel{p}{\to }\infty .$$

(1.4)

As detailed below, the driving force behind this ‘inferior’ behaviour is the bias in learning g₀.

To heuristically illustrate the impact of the bias in learning g₀, we can decompose the scaled estimation error in ${\widehat{\theta }}_0$ as  

$$\sqrt{n}({\widehat{\theta }}_0-{\theta }_0)=\underset{:=a}{\underbrace{{\displaystyle {\left(\frac{1}{n}\sum _{i\in I}{D}_i^2\right)}^{-1}\frac{1}{\sqrt{n}}\sum _{i\in I}{D}_i{U}_i}}}+\underset{:=b}{\underbrace{{\displaystyle {\left(\frac{1}{n}\sum _{i\in I}{D}_i^2\right)}^{-1}\frac{1}{\sqrt{n}}\sum _{i\in I}{D}_i(g_0(X_i)-{\widehat{g}}_0(X_i))}}}.$$

The first term is well behaved under mild conditions, obeying $a⇝N(0,\overline{\mathrm{\Sigma }})$ for some $\overline{\mathrm{\Sigma }}$⁠. Term b is the regularization bias term, which is not centred and diverges in general. Indeed, we have  

$${\displaystyle b={(E[{D_i}^2])}^{-1}\frac{1}{\sqrt{n}}\sum _{i\in I}{m}_0(X_i)(g_0(X_i)-{\widehat{g}}_0(X_i))+o_P(1)}$$

to the first order. Heuristically, b is the sum of n terms that do not have mean zero, $m_0(X_i)(g_0(X_i)-{\widehat{g}}_0(X_i))$⁠, divided by $\sqrt{n}$⁠. These terms have non‐zero mean because, in high‐dimensional or otherwise highly complex settings, we must employ regularized estimators – such as lasso, ridge, boosting or penalized neural nets – for informative learning to be feasible. The regularization in these estimators keeps the variance of the estimator from exploding but also necessarily induces substantive biases in the estimator ${\widehat{g}}_0$ of g₀. Specifically, the rate of convergence of (the bias of) ${\widehat{g}}_0$ to g₀ in the root‐mean‐squared error sense will typically be $n^{-{\varphi }_g}$ with ${\varphi }_g<1/2$⁠. Hence, we expect b to be of stochastic order $\sqrt{n}{n}^{-{\varphi }_g}\to \infty$ as $D_i$ is centred at $m_0(X_i)\ne 0$⁠, which then implies 1.4.

Overcoming Regularization Biases Using Orthogonalization

Now consider a second construction that employs an orthogonalized formulation obtained by directly partialling out the effect of X from D to obtain the orthogonalized regressor $V=D-m_0(X)$⁠. Specifically, we obtain $\widehat{V}=D-{\widehat{m}}_0(X)$⁠, where ${\widehat{m}}_0$ is an ML estimator of m₀ obtained using the auxiliary sample of observations. We are now solving an auxiliary prediction problem to estimate the conditional mean of D given X, so we are doing ‘double prediction’ or ‘double machine learning’.

After partialling the effect of X out from D and obtaining a preliminary estimate of g₀ from the auxiliary sample as before, we can formulate the following debiased ML (DML) estimator for θ₀ using the main sample of observations:³ 

$${\displaystyle {\check{\theta }}_0={\left(\frac{1}{n}\sum _{i\in I}{\widehat{V}}_i{D}_i\right)}^{-1}\frac{1}{n}\sum _{i\in I}{\widehat{V}}_i(Y_i-{\widehat{g}}_0(X_i)).}$$

(1.5)

By approximately orthogonalizing D with respect to X and approximately removing the direct effect of confounding by subtracting an estimate of g₀, ${\check{\theta }}_0$ removes the effect of regularization bias that contaminates 1.3. The formulation of ${\check{\theta }}_0$ also provides direct links to both the classical econometric literature, as the estimator can clearly be interpreted as a linear instrumental variable (IV) estimator, and to the more recent literature on debiased lasso in the context where g₀ is taken to be well approximated by a sparse linear combination of pre‐specified functions of X; see, e.g., Belloni et al. (2013, 2014a,b), Javanmard and Montanari (2014b), van de Geer et al. (2014) and Zhang and Zhang (2014).⁴

To illustrate the benefits of the auxiliary prediction step and the estimation of θ₀ with ${\check{\theta }}_0$⁠, we sketch the properties of ${\check{\theta }}_0$ here. We can decompose the scaled estimation error of ${\check{\theta }}_0$ into three components:  

$$\sqrt{n}({\check{\theta }}_0-{\theta }_0)=a^{\ast }+b^{\ast }+c^{\ast }.$$

The leading term, $a^{\ast }$⁠, will satisfy  

$${\displaystyle a^{\ast }={(E[V^2])}^{-1}\frac{1}{\sqrt{n}}\sum _{i\in I}{V}_i{U}_i⇝N(0,\mathrm{\Sigma })}$$

under mild conditions. The second term, $b^{\ast }$⁠, captures the impact of regularization bias in estimating g₀ and m₀. Specifically, we will have  

$${\displaystyle b^{\ast }={(E[V^2])}^{-1}\frac{1}{\sqrt{n}}\sum _{i\in I}({\widehat{m}}_0(X_i)-m_0(X_i))({\widehat{g}}_0(X_i)-g_0(X_i)),}$$

which now depends on the product of the estimation errors in ${\widehat{m}}_0$ and ${\widehat{g}}_0$⁠. Because this term depends only on the product of the estimation errors, it can vanish under a broad range of data‐generating processes. Indeed, this term is upper‐bounded by $\sqrt{n}{n}^{-({\varphi }_m+{\varphi }_g)}$⁠, where $n^{-{\varphi }_m}$ and $n^{-{\varphi }_g}$ are respectively the rates of convergence of ${\widehat{m}}_0$ to m₀ and ${\widehat{g}}_0$ to g₀; this upper bound can clearly vanish even though both m₀ and g₀ are estimated at relatively slow rates. Verifying that ${\check{\theta }}_0$ has good properties then requires that the remainder term, $c^{\ast }$⁠, is sufficiently well behaved. Sample splitting will play a key role in allowing us to guarantee that $c^{\ast }=o_P(1)$ under weak conditions as outlined below and discussed in detail in Section 3.

Figure 1 provides a numerical illustration of the negative impact of regularization bias and the benefit of orthogonalization. The left panel shows the behaviour of a conventional (non‐orthogonal) ML estimator, ${\widehat{\theta }}_0$⁠, in the partially linear model in a simple simulation experiment where we learn g₀ using a random forest. The g₀ in this experiment is a very smooth function of a small number of variables, so the experiment is seemingly favourable to the use of random forests a priori. The histogram shows the simulated distribution of the centred estimator, ${\widehat{\theta }}_0-{\theta }_0$⁠. The estimator is badly biased, shifted much to the right relative to the true value θ₀. The distribution of the estimator (approximated by the blue histogram) is substantively different from a normal approximation (shown by the red curve) derived under the assumption that the bias is negligible. The right panel shows the behaviour of the orthogonal, DML estimator, ${\check{\theta }}_0$⁠, in the partially linear model in a simple experiment where we learn nuisance functions using random forests. Note that the simulated data are exactly the same as those underlying in the left panel. The simulated distribution of the centred estimator, ${\check{\theta }}_0-{\theta }_0$ (given by the blue histogram) illustrates that the estimator is approximately unbiased, concentrates around θ₀, and is well approximated by the normal approximation obtained in Section 3 (shown by the red curve).

Figure 1.

Comparison of the conventional and double ML estimators. [Color figure can be viewed at wileyonlinelibrary.com]

The Role of Sample Splitting in Removing Bias Induced By Overfitting

Our analysis makes use of sample splitting, which plays a key role in establishing that remainder terms, such as $c^{\ast }$⁠, vanish in probability. In the partially linear model, we find that the remainder $c^{\ast }$ contains terms such as  

$${\displaystyle \frac{1}{\sqrt{n}}\sum _{i\in I}{V}_i({\widehat{g}}_0(X_i)-g_0(X_i)),}$$

(1.6)

which involve $1/\sqrt{n}$‐normalized sums of products of structural unobservables from model 1.1–1.2 with estimation errors in learning the nuisance functions g₀ and m₀. The use of sample splitting allows simple and tight control of such terms. To see this, assume that observations are independent and recall that ${\widehat{g}}_0$ is estimated using only observations in the auxiliary sample. Then, conditioning on the auxiliary sample and recalling that $E[V_i|X_i]=0$⁠, it is easy to verify that term 1.6 has mean zero and variance of order  

$${\displaystyle \frac{1}{n}\sum _{i\in I}{({\widehat{g}}_0(X_i)-g_0(X_i))}^2\stackrel{p}{\to }0.}$$

Thus, the term 1.6 vanishes in probability by Chebyshev's inequality.

While sample splitting allows us to deal with remainder terms such as $c^{\ast }$⁠, its direct application does have the drawback that the estimator of the parameter of interest only makes use of the main sample, which can result in a substantial loss of efficiency as we are only making use of a subset of the available data. However, we can flip the role of the main and auxiliary samples to obtain a second version of the estimator of the parameter of interest. By averaging the two resulting estimators, we can regain full efficiency. Indeed, the two estimators will be approximately independent, so simply averaging them offers an efficient procedure. We call this sample‐splitting procedure – where we swap the roles of main and auxiliary samples to obtain multiple estimates and then average the results – ‘cross‐fitting’. We formally define this procedure and discuss a K‐fold version of cross‐fitting in Section 3.

Without sample splitting, terms such as 1.6 might not vanish and can lead to poor performance of estimators of θ₀. The difficulty arises because model errors, such as $V_i$⁠, and estimation errors, such as ${\widehat{g}}_0(X_i)-g_0(X_i)$⁠, are generally related because the data for observation i are used in forming the estimator ${\widehat{g}}_0$⁠. The association can then lead to poor performance of an estimator of θ₀ that makes use of ${\widehat{g}}_0$ as a plug‐in estimator for g₀ even when this estimator converges at a very favourable rate, say $N^{-1/2+\epsilon }$⁠.

As an artificial but illustrative example of the problems that can result from overfitting, let ${\widehat{g}}_0(X_i)=g_0(X_i)+(Y_i-g_0(X_i))/N^{1/2-\epsilon }$ for any i in the sample used to form estimator ${\widehat{g}}_0$⁠, and note that the second term provides a simple model that captures overfitting of the outcome variable within the estimation sample. This estimator is excellent in terms of rates. If $U_i$ and $D_i$ are bounded, ${\widehat{g}}_0$ converges uniformly to g₀ at the nearly parametric rate $N^{-1/2+\epsilon }$⁠. Despite this fast rate of convergence, term $c^{\ast }$ now explodes if we do not use sample splitting. For example, suppose that the full sample is used to estimate both ${\widehat{g}}_0$ and ${\check{\theta }}_0$⁠. A simple calculation then reveals that term $c^{\ast }$ becomes  

$${\displaystyle \frac{1}{\sqrt{N}}\sum _{i=1}^N{V}_i({\widehat{g}}_0(X_i)-g_0(X_i))\propto N^{\epsilon }\to \infty .}$$

This bias due to overfitting is illustrated in the left panel of Figure 2. The histogram in the figure gives a simulated distribution for the studentized $\check{\theta }$ resulting from using the full sample and the contrived estimator $\widehat{g}(X_i)$ given above. We can see that the histogram is shifted markedly to the left demonstrating substantial bias resulting from overfitting. The right panel of Figure 2 also illustrates that this bias is completely removed by sample splitting. The results in the right panel of Figure 2 make use of the twofold cross‐fitting procedure discussed above using the estimator $\check{\theta }$ and the contrived estimator $\widehat{g}(X_i)$ exactly as in the left panel. The difference is that $\widehat{g}(X_i)$ is formed in one half of the sample and then $\check{\theta }$ is estimated using the other half of the sample. This procedure is then repeated, swapping the roles of the two samples, and the results are averaged. We can see that the substantial bias from the full sample estimator has been removed and that the spread of the histogram corresponding to the cross‐fitting estimator is roughly the same as that of the full sample estimator, clearly illustrating the bias‐reduction property and efficiency of the cross‐fitting procedure.

Figure 2.

Comparison of full‐sample and cross‐fitting procedures. [Color figure can be viewed at wileyonlinelibrary.com]

A less contrived example that highlights the improvements brought by sample splitting is the sparse high‐dimensional IV model analysed in Belloni et al. (2012). Specifically, they consider the IV model  

$$Y=D{\theta }_0+\epsilon ,$$

where $E[\epsilon |D]\ne 0$ but instruments Z exist such that $E[D|Z]$ is not a constant and $E[\epsilon |Z]=0$⁠. Within this model, Belloni et al. (2012) focus on the problem of estimating the optimal instrument, ${\eta }_0(Z)=E[D|Z]$⁠, using lasso‐type methods. If ${\eta }_0(Z)$ is approximately sparse in the sense that only s terms of the dictionary of series transformations $B(Z)=(B_1(Z),\dots ,B_p(Z))$ are needed to approximate the function accurately, Belloni et al. (2012) require that $s^2\ll n$ to establish their asymptotic results when sample splitting is not used. However, they show that these results continue to hold under the much weaker requirement that $s\ll n$ if one employs sample splitting. We note that this example provides a prototypical example where Neyman orthogonality holds and ML methods can usefully be adopted to aid in learning structural parameters of interest. We also note that the weaker conditions required when using sample splitting would also carry over to sparsity‐based estimators in the partially linear model cited above. We discuss this in more detail in Section 4.

While we find substantial appeal in using sample splitting, one can also use empirical process methods to verify that biases introduced due to overfitting are negligible. For example, consider the problematic term in the partially linear model described previously, ${\displaystyle (1/\sqrt{n})\sum _{i\in I}{V}_i({\widehat{g}}_0(X_i)-g_0(X_i))}$⁠. This term is clearly bounded by  

$${\displaystyle \underset{g\in {\mathcal{G}}_N}{sup}\left|\frac{1}{\sqrt{n}}\sum _{i\in I}{V}_i(g(X_i)-g_0(X_i))\right|,}$$

(1.7)

where ${\mathcal{G}}_N$ is the smallest class of functions that contains estimators of g₀, $\widehat{g}$⁠, with high probability. In conventional semi‐parametric statistical and econometric analysis, the complexity of ${\mathcal{G}}_N$ is controlled by invoking Donsker conditions, which allow verification that terms such as 1.7 vanish asymptotically. Importantly, Donsker conditions require that ${\mathcal{G}}_N$ has bounded complexity, specifically a bounded entropy integral. Because of the latter property, Donsker conditions are inappropriate in settings using ML methods where the dimension of X is modelled as increasing with the sample size and estimators necessarily live in highly complex spaces. For example, Donsker conditions rule out even the simplest linear parametric model with high‐dimensional regressors with parameter space given by the Euclidean ball with the unit radius:  

$${\mathcal{G}}_N=\{x\mapsto g(x)=x^{\prime }\theta ;\theta \in {\mathbb{R}}^{p_N}:\parallel \theta \parallel \le 1\}.$$

The entropy of this model, as measured by the logarithm of the covering number, grows at the rate $p_N$⁠. Without invoking Donsker conditions, one can still show that terms such as 1.7 vanish as long as the entropy of ${\mathcal{G}}_N$ does not increase with N too rapidly. A fairly general treatment is given by Belloni et al. (2017) who provide a set of conditions under which terms, such as $c^{\ast }$⁠, can vanish making use of the full sample. However, these conditions on the growth of entropy could result in unnecessarily strong restrictions on model complexity, such as very strict requirements on sparsity in the context of lasso estimation, as demonstrated in the IV example mentioned above. Sample splitting allows one to obtain good results under very weak conditions.

Neyman Orthogonality and Moment Conditions

Now we turn to a generalization of the orthogonalization principle above. The first ‘conventional’ estimator ${\widehat{\theta }}_0$ given in 1.3 can be viewed as a solution to estimating equations  

$${\displaystyle \frac{1}{n}\sum _{i\in I}\varphi (W;{\widehat{\theta }}_0,{\widehat{g}}_0)=0,}$$

where φ is a known score function and ${\widehat{g}}_0$ is the estimator of the nuisance parameter g₀. For example, in the partially linear model above, the score function is $\varphi (W;\theta ,g)=(Y-\theta D-g(X))D$⁠. It is easy to see that this score function φ is sensitive to biased estimation of g. Specifically, the Gateaux derivative operator with respect to g does not vanish:⁵ 

$${\partial }_{g}E[\varphi (W;{\theta }_0,g_0)][g-g_0]\ne 0.$$

The proofs of the general results in Section 3 show that this term's vanishing is a key to establishing good behaviour of an estimator for θ₀.

By contrast, the orthogonalized or double/debiased ML estimator ${\check{\theta }}_0$ given in 1.5 solves  

$${\displaystyle \frac{1}{n}\sum _{i\in I}\psi (W;{\check{\theta }}_0,{\widehat{\eta }}_0)=0,}$$

where ${\widehat{\eta }}_0$ is the estimator of the nuisance parameter η₀ and ψ is an orthogonalized or debiased score function that satisfies the property that the Gateaux derivative operator with respect to η vanishes when evaluated at the true parameter values:  

$${\partial }_{\eta }E[\psi (W;{\theta }_0,{\eta }_0)][\eta -{\eta }_0]=0.$$

(1.8)

We refer to property 1.8 as Neyman orthogonality and to ψ as the Neyman orthogonal score function due to the fundamental contributions in Neyman (1959, 1979), where this notion was introduced. Intuitively, the Neyman orthogonality condition means that the moment conditions used to identify θ₀ are locally insensitive to the value of the nuisance parameter, which allows one to plug in noisy estimates of these parameters without strongly violating the moment condition. In the partially linear model 1.1–1.2, the estimator ${\check{\theta }}_0$ uses the score function $\psi (W;\theta ,\eta )=(Y-D\alpha -g(X))(D-m(X))$⁠, with the nuisance parameter being $\eta =(m,g)$⁠. It is easy to see that these score functions ψ are not sensitive to biased estimation of η₀ in the sense that 1.8 holds. The proofs of the general results in Section 3 show that this property and sample splitting are two generic keys that allow us to establish good behaviour of an estimator for θ₀.

Literature Overview

Our paper builds upon two important bodies of research within the semi‐parametric literature. The first is the literature on obtaining $\sqrt{N}$‐consistent and asymptotically normal estimates of low‐dimensional objects in the presence of high‐dimensional or non‐parametric nuisance functions. The second is the literature on the use of sample splitting to relax entropy conditions. We provide links to each of these bodies of literature in turn.

The problem we study is obviously related to the classical semi‐parametric estimation framework, which focuses on obtaining $\sqrt{N}$‐consistent and asymptotically normal estimates for low‐dimensional components with nuisance parameters estimated by conventional non‐parametric estimators, such as kernels or series. See, for example, the work by Levit (1975), Ibragimov and Hasminskii (1981), Bickel (1982), Robinson (1988), Newey (1990, 1994), van der Vaart (1991), Andrews (1994a), Newey et al. (1998, 2004), Robins and Rotnitzky (1995), Linton (1996), Bickel et al. (1998), Chen et al. (2003), van der Laan and Rose (2011) and Ai and Chen (2012). Neyman orthogonality 1.8, introduced by Neyman (1959), plays a key role in optimal testing theory and adaptive estimation, semi‐parametric learning theory and econometrics, and, more recently, targeted learning theory. For example, Andrews (1994a), Newey (1994) and van der Vaart (1998) provide a general set of results on estimation of a low‐dimensional parameter θ₀ in the presence of nuisance parameters η₀. Andrews (1994a) uses Neyman orthogonality (1.8) and Donsker conditions to demonstrate the key equicontinuity condition  

$${\displaystyle \frac{1}{\sqrt{n}}\sum _{i\in I}(\psi (W_i;{\theta }_0,\widehat{\eta })-\int \psi (w;{\theta }_0,\widehat{\eta })\mathit{dP}(w)-\psi (W_i;{\theta }_0,{\eta }_0))\stackrel{p}{\to }0,}$$

which reduces to 1.6 in the PLR model. Newey (1994) gives conditions on estimating equations and nuisance function estimators so that nuisance function estimators do not affect the limiting distribution of parameters of interest, providing a semi‐parametric version of Neyman orthogonality. van der Vaart (1998) discusses use of semi‐parametrically efficient scores to define estimators that solve estimating equations, setting averages of efficient scores to zero. He also uses efficient scores to define k‐step estimators, where a preliminary estimator is used to estimate the efficient score and then updating is done to further improve estimation; see also comments below on the use of sample splitting.

There is also a related targeted maximum likelihood learning approach, introduced in Scharfstein et al. (1999) in the context of treatments effects analysis. This is substantially generalized by van der Laan and Rubin (2006), who use maximum likelihood in a least favourable direction and then perform one‐step or k‐step updates using the estimated scores, in an effort to better estimate the target parameter.⁶ This procedure is like the least favourable direction approach in semi‐parametrics; see, for example, Severini and Wong (1992). The introduction of the likelihood introduces major benefits, such as allowing simple and natural imposition of constraints inherent in the data (e.g. support restrictions when the outcome is binary or censored) and permitting the use of likelihood cross‐validation to choose the nuisance parameter estimator. This data adaptive choice of the nuisance parameter has been dubbed the ‘super learner’ by van der Laan et al. (2007). In subsequent work, van der Laan and Rose (2011) emphasize the use of ML methods to estimate the nuisance parameters for use with the super learner. Much of this work, including recent work such as Luedtke and van der Laan (2016), Toth and van der Laan (2016) and Zheng et al. (2016), focuses on formal results under a Donsker condition, though the use of sample splitting to relax these conditions has also been advocated in the targeted maximum likelihood setting, as discussed below.

The Donsker condition is a powerful classical condition that allows rich structures for fixed function classes $\mathcal{G}$⁠, but unfortunately it is unsuitable for high‐dimensional settings. Examples of function classes where a Donsker condition holds include functions of a single variable that have total variation bounded by 1 and functions $x\mapsto f(x)$ that have $r>dim(x)/2$ uniformly bounded derivatives. As a further example, functions composed from function classes with VC dimensions bounded by p through a fixed number of algebraic and monotonic transforms are Donsker. However, this property will no longer hold if we let $dim(x)$ grow to infinity with the sample size as this increase in dimension would require that the VC dimension also increases with n. More generally, Donsker conditions are easily violated once dimensions become large. A major point of departure of the present work from the classical literature on semi‐parametric estimation is its explicit focus on high‐complexity/entropy cases. One way to analyse the problem of estimation in high‐entropy cases is to see to what degree equicontinuity results continue to hold while allowing moderate growth of the complexity/entropy of ${\mathcal{G}}_N$⁠. Examples of papers taking this approach in approximately sparse settings are Belloni et al. (2016, 2016, 2017), Chernozhukov et al. (2015b), Javanmard and Montanari (2014a), van de Geer et al. (2014) and Zhang and Zhang (2014). In all of these examples, entropy growth must be limited in what can be very restrictive ways. The entropy conditions rule out the contrived overfitting example mentioned above, which does approximate realistic examples, and might otherwise place severe restrictions on the model. For example, in Belloni et al. (2010, 2012), the optimal instrument needs to be sparse of order $s\ll \sqrt{n}$⁠.

A key device that we use to avoid strong entropy conditions is cross‐fitting via sample splitting. Cross‐fitting is a practical, efficient form of data splitting. Importantly, its use here is not simply as a device to make proofs elementary (which it does), but as a practical method to allow us to overcome the overfitting/high‐complexity phenomena that commonly arise in data analysis based on highly adaptive ML methods. Our treatment builds upon the sample‐splitting ideas employed in Belloni et al. (2010, 2012) who considered sample splitting in a high‐dimensional sparse optimal IV model to weaken the sparsity condition mentioned in the previous paragraph to $s\ll n$⁠. This work, in turn, was inspired by Angrist and Krueger (1995). We also build on Ayyagari (2010) and Robins et al. (2013), where ML methods and sample splitting were used in the estimation of a partially linear model of the effects of pollution while controlling for several covariates. We use the term cross‐fitting to characterize our recommended procedure, partly borrowing the jargon from Fan et al. (2012), who employed a slightly different form of sample splitting to estimate the scale parameter in a high‐dimensional sparse regression. Of course, the use of sample splitting to relax entropy conditions has a long history in semi‐parametric estimation problems. For example, Bickel (1982) considered estimating nuisance functions using a vanishing fraction of the sample, and these results were extended to sample splitting into two equal halves and discretization of the parameter space by Schick (1986). Similarly, van der Vaart (1998) uses two‐way sample splitting and discretization of the parameter space to give weak conditions for k‐step estimators using the efficient scores where sample splitting is used to estimate the updates; see also Hubbard et al. (2016). Robins et al. (2008, 2017) use sample splitting in the construction of higher‐order influence function corrections in semi‐parametric estimation. Some recent work in the targeted maximum likelihood literature, e.g. Zheng and van der Laan (2011), also notes the utility of sample splitting in the context of k‐step updating, though this sample splitting approach is different from the cross‐fitting approach we pursue.

Plan of The Paper

We organize the rest of the paper as follows. In Section 2, we formally define Neyman orthogonality and provide a brief discussion that synthesizes various models and frameworks that can be used to produce estimating equations satisfying this key condition. In Section 3, we carefully define DML estimators and develop their general theory. We then illustrate this general theory by applying it to provide theoretical results for using DML to estimate and carry out inference for key parameters in the PLR model, and for using DML to estimate and carry out inference for coefficients on endogenous variables in a partially linear IV model in Section 4. In Section 5, we provide a further illustration of the general theory by applying it to develop theoretical results for DML estimation and inference for average treatment effects (ATEs) and average treatment effects on the treated (ATTEs) under unconfoundedness, and for DML estimation of local average treatment effects (LATEs) in an IV context within the potential outcomes framework; see Imbens and Rubin (2015). Finally, we apply DML in three empirical illustrations in Section 6. In the Appendix, we define additional notation and present proofs.

Notation

The symbols Pr and E denote probability and expectation operators with respect to a generic probability measure that describes the law of the data. If we need to signify the dependence on a probability measure P, we use P as a subscript in ${\mathit{Pr}}_P$ and $E_P$⁠. We use capital letters, such as W, to denote random elements and we use the corresponding lowercase letters, such as w, to denote fixed values that these random elements can take. In what follows, we use ${\parallel ·\parallel }_{P,q}$ to denote the $L^q(P)$ norm; for example, we denote ${\displaystyle {\parallel f\parallel }_{P,q}:={\parallel f(W)\parallel }_{P,q}:={(\int |f(w)|}^q{dP(w))}^{1/q}}$⁠. We use $x^{\prime }$ to denote the transpose of a column vector x. For a differentiable map $x\mapsto f(x)$⁠, mapping ${\mathbb{R}}^d$ to ${\mathbb{R}}^k$⁠, we use ${\partial }_{x^{\prime }}f$ to abbreviate the partial derivatives $(\partial /\partial x^{\prime })f$⁠, and we correspondingly use the expression ${\partial }_{x^{\prime }}f(x_0)$ to mean ${\partial }_{x^{\prime }}f(x){\mid }_{x=x_0}$⁠, etc.

2. Construction of Neyman Orthogonal Score/Moment Functions

Here we formally introduce the model and we discuss several methods for generating orthogonal scores in a wide variety of settings, including the classical Neyman construction. We also use this as an opportunity to synthesize some recent developments in the literature.

2.1. Moment Condition/Estimating Equation Framework

We are interested in the true value θ₀ of the low‐dimensional target parameter $\theta \in \mathrm{\Theta }$⁠, where Θ is a non‐empty measurable subset of ${\mathbb{R}}^{d_{\theta }}$⁠. We assume that θ₀ satisfies the moment conditions  

$$E_P[\psi (W;{\theta }_0,{\eta }_0)]=0,$$

(2.1)

where $\psi ={({\psi }_1,\dots ,{\psi }_{d_{\theta }})}^{\prime }$ is a vector of known score functions, W is a random element taking values in a measurable space $(\mathcal{W},{\mathcal{A}}_{\mathcal{W}})$ with law determined by a probability measure $P\in {\mathcal{P}}_N$⁠, and η₀ is the true value of the nuisance parameter $\eta \in T$⁠, where T is a convex subset of some normed vector space with the norm denoted by ${\parallel ·\parallel }_T$⁠. We assume that the score functions ${\psi }_j:\mathcal{W}\times \mathrm{\Theta }\times T\to \mathbb{R}$ are measurable once we equip Θ and T with their Borel σ‐fields, and we assume that a random sample ${(W_i)}_{i=1}^N$ from the distribution of W is available for estimation and inference.

As discussed in Section 1, we require the Neyman orthogonality condition for the score ψ. To introduce the condition, for $\stackrel{\sim }{T}=\{\eta -{\eta }_0:\eta \in T\}$ we define the pathwise (or the Gateaux) derivative map ${\mathrm{D}}_r:\stackrel{\sim }{T}\to {\mathbb{R}}^{d_{\theta }}$⁠,  

$${\mathrm{D}}_r[\eta -{\eta }_0]:={\partial }_r\{E_P[\psi (W;{\theta }_0,{\eta }_0+r(\eta -{\eta }_0)]\},\eta \in T,$$

for all $r\in [0,1)$⁠, which we assume to exist. For convenience, we also denote  

$${\partial }_{\eta }{E}_P[\psi (W;{\theta }_0,{\eta }_0)][\eta -{\eta }_0]:={\mathrm{D}}_0[\eta -{\eta }_0],\eta \in T.$$

(2.2)

Note that $\psi (W;{\theta }_0,{\eta }_0+r(\eta -{\eta }_0))$ here is well defined because for all $r\in [0,1)$ and $\eta \in T$⁠,  

$${\eta }_0+r(\eta -{\eta }_0)=(1-r){\eta }_0+r\eta \in T,$$

as T is a convex set. In addition, let ${\mathcal{T}}_N\subset T$ be a nuisance realization set such that the estimators ${\widehat{\eta }}_0$ of η₀ specified below take values in this set with high probability. In practice, we typically assume that ${\mathcal{T}}_N$ is a properly shrinking neighbourhood of η₀. Note that ${\mathcal{T}}_N-{\eta }_0$ is the nuisance deviation set, which contains deviations of ${\widehat{\eta }}_0$ from η₀, ${\widehat{\eta }}_0-{\eta }_0$⁠, with high probability. The Neyman orthogonality condition requires that the derivative in 2.2 vanishes for all $\eta \in {\mathcal{T}}_N$⁠.  

Definition 2.1.Neyman Orthogonality

The score $\psi ={({\psi }_1,\dots ,{\psi }_{d_{\theta }})}^{\prime }$ obeys the orthogonality condition at $({\theta }_0,{\eta }_0)$ with respect to the nuisance realization set ${\mathcal{T}}_N\subset T$ if 2.1 holds and the pathwise derivative map ${\mathrm{D}}_r[\eta -{\eta }_0]$ exists for all $r\in [0,1)$ and $\eta \in {\mathcal{T}}_N$ and vanishes at $r=0$⁠; namely,  

$${\partial }_{\eta }{E}_P[\psi (W;{\theta }_0,{\eta }_0)][\eta -{\eta }_0]=0,\text{for}\text{all}\eta \in {\mathcal{T}}_N.$$

(2.3)

We remark here that condition 2.3 holds with ${\mathcal{T}}_N=T$ when η is a finite‐dimensional vector as long as ${\partial }_{\eta }{E}_P[{\psi }_j(W;{\theta }_0,{\eta }_0)]=0$ for all $j=1,\dots ,d_{\theta }$⁠, where ${\partial }_{\eta }{E}_P[{\psi }_j(W;{\theta }_0,{\eta }_0)]$ denotes the vector of partial derivatives of the function $\eta \mapsto E_P[{\psi }_j(W;{\theta }_0,\eta )]$ for $\eta ={\eta }_0$⁠.

Sometimes it will also be helpful to use an approximate Neyman orthogonality condition as opposed to the exact one given in Definition 2.1.  

Definition 2.2.Neyman Near‐Orthogonality

The score $\psi ={({\psi }_1,\dots ,{\psi }_{d_{\theta }})}^{\prime }$ obeys the ${\lambda }_N$ near‐orthogonality condition at $({\theta }_0,{\eta }_0)$ with respect to the nuisance realization set ${\mathcal{T}}_N\subset T$ if 2.1 holds and the pathwise derivative map ${\mathrm{D}}_r[\eta -{\eta }_0]$ exists for all $r\in [0,1)$ and $\eta \in {\mathcal{T}}_N$ and is small at $r=0$⁠; namely,  

$$\parallel {\partial }_{\eta }{E}_P[\psi (W;{\theta }_0,{\eta }_0)][\eta -{\eta }_0]\parallel \le {\lambda }_N,\text{for}\text{all}\eta \in {\mathcal{T}}_N,$$

(2.4)

where ${\{{\lambda }_N\}}_{N\ge 1}$ is a sequence of positive constants such that ${\lambda }_N=o(N^{-1/2})$⁠.

2.2. Construction of Neyman Orthogonal Scores

If we start with a score φ that does not satisfy the orthogonality condition above, we first transform it into a score ψ that does. Here we outline several methods for doing so.

2.2.1. Neyman Orthogonal Scores for Likelihood and Other M‐Estimation Problems With Finite‐Dimensional Nuisance Parameters

First, we describe the construction used by Neyman (1959) to derive his celebrated orthogonal score and $C(\alpha )$‐statistic in a maximum likelihood setting.⁷ Such a construction also underlies the concept of local unbiasedness in the construction of optimal tests in, e.g., Ferguson (1967), and it was extended to non‐likelihood settings by Wooldridge (1991). The discussion of Neyman's construction here draws on Chernozhukov et al. (2015a).

To describe the construction, let $\theta \in \mathrm{\Theta }\subset {\mathbb{R}}^{d_{\theta }}$ and $\beta \in \mathcal{B}\subset {\mathbb{R}}^{d_{\beta }}$⁠, where $\mathcal{B}$ is a convex set, be the target and the nuisance parameters, respectively. Further, suppose that the true parameter values θ₀ and β₀ solve the optimization problem  

$$\underset{\theta \in \mathrm{\Theta },\beta \in \mathcal{B}}{max}{E}_P[\ell (W;\theta ,\beta )],$$

(2.5)

where $\ell (W;\theta ,\beta )$ is a known criterion function. For example, $\ell (W;\theta ,\beta )$ can be the log‐likelihood function associated with observation W. More generally, we refer to $\ell (W;\theta ,\beta )$ as the quasi‐log‐likelihood function. Then, under mild regularity conditions, θ₀ and β₀ satisfy  

$$E_P[{\partial }_{\theta }\ell (W;{\theta }_0,{\beta }_0)]=0,E_P[{\partial }_{\beta }\ell (W;{\theta }_0,{\beta }_0)]=0.$$

(2.6)

Note that the original score function $\varphi (W;\theta ,\beta )={\partial }_{\theta }\ell (W;\theta ,\beta )$ for estimating θ₀ will not generally satisfy the orthogonality condition. Now consider the new score function, which we refer to as the Neyman orthogonal score,  

$$\psi (W;\theta ,\eta )={\partial }_{\theta }\ell (W;\theta ,\beta )-\mu {\partial }_{\beta }\ell (W;\theta ,\beta ),$$

(2.7)

where the nuisance parameter is  

$$\eta ={({\beta }^{\prime },\mathrm{vec}{(\mu )}^{\prime })}^{\prime }\in T=\mathcal{B}\times {\mathbb{R}}^{d_{\theta }{d}_{\beta }}\subset {\mathbb{R}}^p,p=d_{\beta }+d_{\theta }{d}_{\beta },$$

and μ is the $d_{\theta }\times d_{\beta }$ orthogonalization parameter matrix whose true value μ₀ solves the equation  

$$J_{\theta \beta }-\mu J_{\beta \beta }=0$$

(2.8)

for  

$$J=\left(\begin{array}{cc}{J}_{\theta \theta }& J_{\theta \beta }\\ J_{\beta \theta }& J_{\beta \beta }\end{array}\right)={\partial }_{({\theta }^{\prime },{\beta }^{\prime })}{E}_P[{\partial }_{{({\theta }^{\prime },{\beta }^{\prime })}^{\prime }}\ell (W;\theta ,\beta )]{|}_{\theta ={\theta }_0;\beta ={\beta }_0}.$$

The true value of the nuisance parameter η is  

$${\eta }_0={({\beta }_0^{\prime },\mathrm{vec}{({\mu }_0)}^{\prime })}^{\prime };$$

(2.9)

and when $J_{\beta \beta }$ is invertible, 2.8 has the unique solution,  

$${\mu }_0=J_{\theta \beta }{J}_{\beta \beta }^{-1}.$$

(2.10)

The following lemma shows that the score ψ in 2.7 satisfies the Neyman orthogonality condition.  

Lemma 2.1. Neyman Orthogonal Scores for Quasi‐Likelihood Settings

If 2.6 holds, J exists, and $J_{\beta \beta }$ is invertible, then the score ψ in 2.7 is Neyman orthogonal at $({\theta }_0,{\eta }_0)$ with respect to the nuisance realization set ${\mathcal{T}}_N=T$⁠.

 

Remark 2.1. Additional Nuisance Parameters

Note that the orthogonal score ψ in 2.7 has nuisance parameters consisting of the elements of μ in addition to the elements of β, and Lemma 2.1 shows that Neyman orthogonality holds both with respect to β and with respect to μ. We will find that Neyman orthogonal scores in other settings, including infinite‐dimensional ones, have a similar property.

 

Remark 2.2. Efficiency

Note that in this example, μ₀ not only creates the necessary orthogonality but also creates the efficient score for inference on the target parameter θ when the quasi‐log‐likelihood function is the true (possibly conditional) log‐likelihood, as demonstrated by Neyman (1959).

 

Example 2.1. High‐Dimensional Linear Regression

As an application of the construction above, consider the following linear predictive model,  

$$\begin{array}{ccc}& Y=D{\theta }_0+X^{\prime }{\beta }_0+U,& E_P[U{(X^{\prime },D)}^{\prime }]=0,\hfill \end{array}$$

(2.11)

 

$$\begin{array}{ccc}& D=X^{\prime }{\gamma }_0+V,& E_P[VX]=0,\hfill \end{array}$$

(2.12)

where, for simplicity, we assume that θ₀ is a scalar. The first equation here is the main predictive model, and the second equation only plays a role in the construction of the Neyman orthogonal scores. It is well known that θ₀ and β₀ in this model solve the optimization problem 2.5 with  

$$\ell (W;\theta ,\beta )=-\frac{{(Y-D\theta -X^{\prime }\beta )}^2}{2},\theta \in \mathrm{\Theta }=\mathbb{R},\beta \in \mathcal{B}={\mathbb{R}}^{d_{\beta }},$$

where we denote $W={(Y,D,X^{\prime })}^{\prime }$⁠. Hence, equations 2.6 hold with  

$$\partial {\ell }_{\theta }(W;\theta ,\beta )=(Y-D\theta -X^{\prime }\beta )D,\partial {\ell }_{\beta }(W;\theta ,\beta )=(Y-D\theta -X^{\prime }\beta )X,$$

and the matrix J satisfies  

$$J_{\theta \beta }=-E_P[{\mathit{DX}}^{\prime }],J_{\beta \beta }=-E_P[{\mathit{XX}}^{\prime }].$$

The Neyman orthogonal score is then given by  

$$\psi (W;\theta ,\eta )=(Y-D\theta -X^{\prime }\beta )(D-\mu X);\eta ={({\beta }^{\prime },\mathrm{vec}{(\mu )}^{\prime })}^{\prime };$$

 

$$\psi (W;{\theta }_0,{\eta }_0)=U(D-{\mu }_{0}X);{\mu }_0=E_P[{\mathit{DX}}^{\prime }]{(E_P[{XX}^{\prime }])}^{-1}={\gamma }_0^{\prime }.$$

(2.13)

If the vector of covariates X here is high‐dimensional but the vectors of parameters β₀ and γ₀ are approximately sparse, we can use ℓ₁‐penalized least‐squares, ℓ₂‐boosting, or forward selection methods to estimate β₀ and ${\gamma }_0={\mu }_0^{\prime }$⁠, and hence ${\mu }_0={({\beta }_0^{\prime },\mathrm{vec}{({\mu }_0)}^{\prime })}^{\prime }$⁠; see references cited in Section 1.

If $J_{\beta \beta }$ is not invertible, 2.8 typically has multiple solutions. In this case, it is convenient to focus on a minimal norm solution,  

$${\mu }_0=argmin\parallel \mu \parallel \text{such}\text{that}\parallel J_{\theta \beta }-\mu J_{\beta \beta }{\parallel }_q=0$$

for a suitably chosen norm ${\parallel ·\parallel }_q$ on the space of $d_{\theta }\times d_{\beta }$ matrices. With an eye on solving the empirical version of this problem, we can also consider the relaxed version of this problem,  

$${\mu }_0=argmin\parallel \mu \parallel \text{such}\text{that}\parallel J_{\theta \beta }-\mu J_{\beta \beta }{\parallel }_q\le r_N$$

(2.14)

for some $r_N>0$ such that $r_N\to 0$ as $N\to \infty$⁠. This relaxation is also helpful when $J_{\beta \beta }$ is invertible but ill‐conditioned. The following lemma shows that using μ₀ in 2.14 leads to Neyman near‐orthogonal scores. The proof of this lemma can be found in the Appendix.  

Lemma 2.2.Neyman Near‐Orthogonal Scores for Quasi‐Likelihood Settings

If 2.6 holds, J exists, the solution of the optimization problem 2.14 exists, and μ₀ is taken to be this solution, then the score ψ defined in 2.7 is Neyman ${\lambda }_N$ near‐orthogonal at $({\theta }_0,{\eta }_0)$ with respect to the nuisance realization set ${\mathcal{T}}_N=\{\beta \in \mathcal{B}:\parallel \beta -{\beta }_0{\parallel }_q^{\ast }\le {\lambda }_N/r_N\}\times {\mathbb{R}}^{d_{\theta }{d}_{\beta }}$⁠, where the norm ${\parallel ·\parallel }_q^{\ast }$ on ${\mathbb{R}}^{d_{\beta }}$ is defined by ${\parallel \beta \parallel }_q^{\ast }={sup}_A\parallel A\beta \parallel$ with the supremum being taken over all $d_{\theta }\times d_{\beta }$ matrices A such that ${\parallel A\parallel }_q\le 1$⁠.

 

Example 2.1. Continued

In the high‐dimensional linear regression example above, the relaxation 2.14 is helpful when $J_{\beta \beta }=E_P[X{X}^{\prime }]$ is ill‐conditioned. Specifically, if one suspects that $E_P[X{X}^{\prime }]$ is ill‐conditioned, one can define μ₀ as the solution to the following optimization problem:  

$$min\parallel \mu \parallel \text{such}\text{that}\parallel E_P[{\mathit{DX}}^{\prime }]-\mu E_P[{\mathit{XX}}^{\prime }]{\parallel }_{\infty }\le r_N.$$

(2.15)

Lemma 2.2 then shows that using this μ₀ leads to a score ψ that obeys the Neyman near‐orthogonality condition. Alternatively, one can define μ₀ as the solution of the following closely related optimization problem,  

$$\underset{\mu }{min}(\mu E_P[{\mathit{XX}}^{\prime }]{\mu }^{\prime }-\mu E_P[\mathit{DX}]+r_N{\parallel \mu \parallel }_1),$$

whose solution also obeys $\parallel E_P[\mathit{DX}]-\mu E_P[{\mathit{XX}}^{\prime }]{\parallel }_{\infty }\le r_N$⁠, which follows from the first‐order conditions. An empirical version of either problem leads to a Lasso‐type estimator of the regularized solution μ₀; see Javanmard and Montanari (2014a).

 

Remark 2.3. Giving Up Efficiency

Note that the regularized μ₀ in 2.14 creates the necessary near‐orthogonality at the cost of giving up somewhat on the efficiency of the score ψ. At the same time, regularization may generate additional robustness gains as achieving full efficiency by estimating μ₀ in 2.10 may require stronger conditions.

 

Remark 2.4. Concentrating‐Out Approach

The approach for constructing Neyman orthogonal scores described above is closely related to the following concentrating‐out approach, which has been used, for example, in Newey (1994), to show Neyman orthogonality when β is infinite dimensional. For all $\theta \in \mathrm{\Theta }$⁠, let ${\beta }_{\theta }$ be the solution of the following optimization problem:  

$$\underset{\beta \in \mathcal{B}}{max}{E}_P[\ell (W;\theta ,\beta )].$$

Under mild regularity conditions, ${\beta }_{\theta }$ satisfies  

$${\partial }_{\beta }{E}_P[\ell (W;\theta ,{\beta }_{\theta })]=0,\text{for}\text{all}\theta \in \mathrm{\Theta }.$$

(2.16)

Differentiating 2.16 with respect to θ and interchanging the order of differentiation gives  

$$\begin{array}{ccc}\hfill 0& =& {\partial }_{\theta }{\partial }_{\beta }{E}_P[\ell (W;\theta ,{\beta }_{\theta })]={\partial }_{\beta }{\partial }_{\theta }{E}_P[\ell (W;\theta ,{\beta }_{\theta })]\hfill \\ & =& {\partial }_{\beta }{E}_P[{\partial }_{\theta }\ell (W;\theta ,{\beta }_{\theta })+{[{\partial }_{\theta }{\beta }_{\theta }]}^{{}^{\prime }}{\partial }_{\beta }\ell (W;\theta ,{\beta }_{\theta })]\hfill \\ & =& {\partial }_{\beta }{E}_P[\psi (W;\theta ,\beta ,{\partial }_{\theta }{\beta }_{\theta })]{|}_{\beta ={\beta }_{\theta }},\hfill \end{array}$$

where we denote  

$$\psi (W;\theta ,\beta ,{\partial }_{\theta }{\beta }_{\theta }):={\partial }_{\theta }\ell (W;\theta ,\beta )+{[{\partial }_{\theta }{\beta }_{\theta }]}^{\prime }{\partial }_{\beta }\ell (W;\theta ,\beta ).$$

This vector of functions is a score with nuisance parameters $\eta ={({\beta }^{\prime },\text{vec}({\partial }_{\theta }{\beta }_{\theta }))}^{\prime }$⁠. As before, additional nuisance parameters, ${\partial }_{\theta }{\beta }_{\theta }$ in this case, are introduced when the orthogonal score is formed. Evaluating these equations at θ₀ and β₀, it follows from the previous equation that $\psi (W;\theta ,\beta ,{\partial }_{\theta }{\beta }_{\theta })$ is orthogonal with respect to β and from $E_P[{\partial }_{\beta }\ell (W;{\theta }_0,{\beta }_0)]=0$ that we have orthogonality with respect to ${\partial }_{\theta }{\beta }_{\theta }$⁠. Thus, maximizing the expected objective function with respect to the nuisance parameters, plugging that maximum back in, and differentiating with respect to the parameters of interest produces an orthogonal moment condition. See also Section 2.2.3.

2.2.2. Neyman Orthogonal Scores in Generalized Method of Moments Problems

The construction in the previous section gives a Neyman orthogonal score whenever the moment conditions 2.6 hold, and, as discussed in Remark 2.2, the resulting score is efficient as long as $\ell (W;\theta ,\beta )$ is the log‐likelihood function. The question, however, remains about constructing the efficient score when $\ell (W;\theta ,\beta )$ is not necessarily a log‐likelihood function. In this section, we answer this question and describe a generalized method of moments (GMM)‐based method of constructing an efficient and Neyman orthogonal score in this more general case. The discussion here is related to Lee (2005), Bera et al. (2010) and Chernozhukov et al. (2015b).

Because GMM does not require that the moment conditions 2.6 are obtained from the first‐order conditions of the optimization problem 2.5, we use a different notation for the moment conditions. Specifically, we consider parameters $\theta \in \mathrm{\Theta }\subset {\mathbb{R}}^{d_{\theta }}$ and $\beta \in \mathcal{B}\subset {\mathbb{R}}^{d_{\beta }}$⁠, where $\mathcal{B}$ is a convex set, whose true values, θ₀ and β₀, solve the moment conditions  

$$E_P[m(W;{\theta }_0,{\beta }_0)]=0,$$

(2.17)

where $m:\mathcal{W}\times \mathrm{\Theta }\times \mathcal{B}\to {\mathbb{R}}^{d_m}$ is a known vector‐valued function, and $d_m\ge d_{\theta }+d_{\beta }$ is the number of moment conditions. In this case, a Neyman orthogonal score function is  

$$\psi (W;\theta ,\eta )=\mu m(W;\theta ,\beta ),$$

(2.18)

where the nuisance parameter is  

$$\eta ={({\beta }^{\prime },\mathrm{vec}{(\mu )}^{\prime })}^{\prime }\in T=\mathcal{B}\times {\mathbb{R}}^{d_{\theta }{d}_m}\subset {\mathbb{R}}^p,p=d_{\beta }+d_{\theta }{d}_m,$$

and μ is the $d_{\theta }\times d_m$ orthogonalization parameter matrix whose true value is  

$${\mu }_0=(A^{\prime }{\mathrm{\Omega }}^{-1}-A^{\prime }{\mathrm{\Omega }}^{-1}{G}_{\beta }{(G_{\beta }^{\prime }{\mathrm{\Omega }}^{-1}{G}_{\beta })}^{-1}{G}_{\beta }^{\prime }{\mathrm{\Omega }}^{-1}).$$

Here  

$$\begin{array}{ccc}\hfill G_{\gamma }& =& {\partial }_{{\gamma }^{\prime }}{E}_P{[m(W;\theta ,\beta )]|}_{\gamma ={\gamma }_0}\hfill \\ & =& [{\partial }_{{\theta }^{\prime }}{E}_P[m(W;\theta ,\beta )],{\partial }_{{\beta }^{\prime }}{E}_P[m(W;\theta ,\beta )]]{|}_{\gamma ={\gamma }_0}=:[G_{\theta },G_{\beta }],\hfill \end{array}$$

for $\gamma ={({\theta }^{\prime },{\beta }^{\prime })}^{\prime }$ and ${\gamma }_0={({\theta }_0^{\prime },{\beta }_0^{\prime })}^{\prime }$⁠, A is a $d_m\times d_{\theta }$ moment selection matrix, Ω is a $d_m\times d_m$ positive definite weighting matrix, and both A and Ω can be chosen arbitrarily. Note that setting  

$$A=G_{\theta }\text{and}\mathrm{\Omega }={\mathrm{Var}}_P(m(W;{\theta }_0,{\beta }_0)])=E_P[m(W;{\theta }_0,{\beta }_0)m{(W;{\theta }_0,{\beta }_0)}^{\prime }]$$

leads to the efficient score in the sense of yielding an estimator of θ₀ having the smallest variance in the class of GMM estimators Hansen, (1982), and, in fact, to the semi‐parametrically efficient score; see Levit (1975), Nevelson (1977) and Chamberlain (1987). Let ${\eta }_0={({\beta }_0^{\prime },\mathrm{vec}{({\mu }_0)}^{\prime })}^{\prime }$ be the true value of the nuisance parameter $\eta ={({\beta }^{\prime },\mathrm{vec}{(\mu )}^{\prime })}^{\prime }$⁠. The following lemma shows that the score ψ in 2.18 satisfies the Neyman orthogonality condition.  

Lemma 2.3.Neyman Orthogonal Scores for Gmm Settings

If 2.17 holds, $G_{\gamma }$ exists and Ω is invertible, then the score ψ in 2.18 is Neyman orthogonal at $({\theta }_0,{\eta }_0)$ with respect to the nuisance realization set ${\mathcal{T}}_N=T$⁠.

As in the quasi‐likelihood case, we can also consider near‐orthogonal scores. Specifically, note that one of the orthogonality conditions that the score ψ in 2.18 has to satisfy is that ${\mu }_0{G}_{\beta }=0$⁠, which can be rewritten as  

$$A^{\prime }{\mathrm{\Omega }}^{-1/2}(I-L{(L^{\prime }L)}^{-1}{L}^{\prime })L=0,$$

where $L={\mathrm{\Omega }}^{-1/2}{G}_{\beta }$⁠. Here, the part $A^{\prime }{\mathrm{\Omega }}^{-1/2}L{(L^{\prime }L)}^{-1}{L}^{\prime }$ can be expressed as ${\gamma }_0{L}^{\prime }$⁠, where ${\gamma }_0=A^{\prime }{\mathrm{\Omega }}^{-1/2}L{(L^{\prime }L)}^{-1}$ solves the optimization problem  

$${min\parallel \gamma \parallel }_o\text{such}\text{that}{\parallel A^{\prime }{\mathrm{\Omega }}^{-1/2}L-\gamma L^{\prime }L\parallel }_{\infty }=0,$$

for a suitably chosen norm ${\parallel ·\parallel }_o$⁠. When $L^{\prime }L$ is close to being singular, this problem can be relaxed:  

$${min\parallel \gamma \parallel }_o\text{such}\text{that}{\parallel A^{\prime }{\mathrm{\Omega }}^{-1/2}L-\gamma L^{\prime }L\parallel }_{\infty }\le r_N.$$

(2.19)

This relaxation leads to Neyman near‐orthogonal scores.  

Lemma 2.4.Neyman Near‐Orthogonal Scores for Gmm Settings

In the set‐up above, with γ₀ denoting the solution of 2.19, we have for ${\mu }_0:=A^{\prime }{\mathrm{\Omega }}^{-1}-{\gamma }_0{L}^{\prime }{\mathrm{\Omega }}^{-1/2}$ and ${\eta }_0={({\beta }_0^{\prime },\mathrm{vec}{({\mu }_0)}^{\prime })}^{\prime }$ that ψ defined in 2.18 is the Neyman ${\lambda }_N$ near‐orthogonal score at $({\theta }_0,{\eta }_0)$ with respect to the nuisance realization set ${\mathcal{T}}_N=\{\beta \in \mathcal{B}:\parallel \beta -{\beta }_0{\parallel }_1\le {\lambda }_N/r_N\}\times {\mathbb{R}}^{d_{\theta }{d}_m}$⁠.

2.2.3. Neyman Orthogonal Scores for Likelihood and Other M‐Estimation Problems With Infinite‐Dimensional Nuisance Parameters

Here we show that the concentrating‐out approach described in Remark 2.4 for the case of finite‐dimensional nuisance parameters can be extended to the case of infinite‐dimensional nuisance parameters. Let $\ell (W;\theta ,\beta )$ be a known criterion function, where θ and β are the target and the nuisance parameters taking values in Θ and $\mathcal{B}$⁠, respectively, and let us assume that the true values of these parameters, θ₀ and β₀, solve the optimization problem 2.5. The function $\ell (W;\theta ,\beta )$ is analogous to that discussed above but now, instead of assuming that $\mathcal{B}$ is a (convex) subset of a finite‐dimensional space, we assume that $\mathcal{B}$ is some (convex) set of functions, so that β is the functional nuisance parameter. For example, $\ell (W;\theta ,\beta )$ could be a semi‐parametric log‐likelihood where β is the non‐parametric part of the model. More generally, $\ell (W;\theta ,\beta )$ could be some other criterion function such as the negative of a squared residual. Also, let  

$${\beta }_{\theta }=arg\underset{\beta \in \mathcal{B}}{max}{E}_P[\ell (W;\theta ,\beta )]$$

(2.20)

be the concentrated‐out non‐parametric part of the model. Note that ${\beta }_{\theta }$ is a function‐valued function. Now consider the score function  

$$\psi (W;\theta ,\eta )=\frac{d\ell (W;\theta ,\eta (\theta ))}{d\theta },$$

(2.21)

where the nuisance parameter is $\eta :\mathrm{\Theta }\to \mathcal{B}$⁠, and its true value η₀ is given by  

$${\eta }_0(\theta )={\beta }_{\theta },\text{for}\text{all}\theta \in \mathrm{\Theta }.$$

Here, the symbol $d/d\theta$ denotes the full derivative with respect to θ, so that we differentiate with respect to both θ arguments in $\ell (W;\theta ,\eta (\theta ))$⁠. The following lemma shows that the score ψ in 2.21 satisfies the Neyman orthogonality condition.  

Lemma 2.5.Neyman Orthogonal Scores Via Concentrating‐Out Approach

Suppose that 2.5 holds, and let T be a convex set of functions mapping Θ into $\mathcal{B}$ such that ${\eta }_0\in T$⁠. Also, suppose that for each $\eta \in T$⁠, the function $\theta \mapsto \ell (W;\theta ,\eta (\theta ))$ is continuously differentiable almost surely. Then, under mild regularity conditions, the score ψ in 2.21 is Neyman orthogonal at $({\theta }_0,{\eta }_0)$ with respect to the nuisance realization set ${\mathcal{T}}_N=T$⁠.

As an example, consider the partially linear model from Section 1. Let  

$$\ell (W;\theta ,\beta )=-\frac{1}{2}{(Y-D\theta -\beta (X))}^2,$$

and let $\mathcal{B}$ be the set of functions of X with finite mean square. Then  

$$({\theta }_0,{\beta }_0)=arg\underset{\theta \in \mathrm{\Theta },\beta \in \mathcal{B}}{max}{E}_P[\ell (W;\theta ,\beta )]$$

and  

$${\beta }_{\theta }(X)=E_P[Y-D\theta |X],\theta \in \mathrm{\Theta }.$$

Hence, 2.21 gives the following Neyman orthogonal score  

$$\begin{array}{ccc}\hfill \psi (W;\theta ,{\beta }_{\theta })& =& -\frac{1}{2}\frac{d{\{Y-D\theta -E_P[Y-D\theta |X]\}}^2}{d\theta }\hfill \\ & =& (D-E_P[D|X])\times (Y-E_P[Y|X]-(D-E_P[D|X])\theta )\hfill \\ & =& (D-m_0(X))\times (Y-D\theta -g_0(X)),\hfill \end{array}$$

which corresponds to the estimator θ₀ described in Section 1 in 1.5.

It is important to note that the concentrating‐out approach described here gives a Neyman orthogonal score without requiring that $\ell (W;\theta ,\beta )$ is the log‐likelihood function. Except for the technical conditions needed to ensure the existence of derivatives and their interchangeability, the only condition that is required is that θ₀ and β₀ solve the optimization problem 2.5. If $\ell (W;\theta ,\beta )$ is the log‐likelihood function, however, it follows from Newey (1994, p. 1359), that the concentrating‐out approach actually yields the efficient score. An alternative, but closely related, approach to derive the efficient score in the likelihood setting would be to apply Neyman's construction described above for a one‐dimensional least favourable parametric submodel; see Severini and Wong (1992) and Chapter 25 of van der Vaart (1998).

 

Remark 2.5. Generating Orthogonal Scores By Varying

$\mathcal{B}$ When we calculate the concentrated‐out non‐parametric part ${\beta }_{\theta }$⁠, we can use some other set of functions Υ instead of $\mathcal{B}$ on the right‐hand side of 2.20:  

$${\beta }_{\theta }=arg\underset{\beta \in {\rm Y}}{max}{E}_P[\ell (W;\theta ,\beta )].$$

By replacing $\mathcal{B}$ by Υ, we can generate a different Neyman orthogonal score. Of course, this replacement might also change the true value θ₀ of the parameter of interest, which is an important consideration for the selection of Υ. For example, consider the partially linear model and assume that X has two components, X₁ and X₂. Now, consider what would happen if we replaced $\mathcal{B}$⁠, which is the set of functions of X with finite mean square, by the set of functions Υ that is the mean square closure of functions that are additive in X₁ and X₂:  

$${\rm Y}=\overline{\{h(X_1)+h(X_2)\}}.$$

Let ${\overline{E}}_P$ denote the least‐squares projection on Υ. Then, applying the previous calculation with ${\overline{E}}_P$ replacing $E_P$ gives  

$$\psi (W;\theta ,{\beta }_{\theta })=(D-{\overline{E}}_P[D|X])\times (Y-{\overline{E}}_P[Y|X]+(D-{\overline{E}}_P[D|X])\theta ),$$

which provides an orthogonal score based on additive function of X₁ and X₂. Here, it is important to note that the solution to $E_P[\psi (W,\theta ,{\beta }_{\theta })]=0$ will be the true θ₀ only when the true function of X in the partially linear model is additive. More generally, the solution of the moment condition would be the coefficient of D in the least‐squares projection of Y on functions of the form $D\theta +h_1(X_1)+h_1(X_2)$⁠. Note, though, that the corresponding score is orthogonal by virtue of additivity being imposed in the estimation of ${\overline{E}}_P[Y|X]$ and ${\overline{E}}_P[D|X]$⁠.

2.2.4. Neyman Orthogonal Scores for Conditional Moment Restriction Problems With Infinite‐Dimensional Nuisance Parameters

Next we consider the conditional moment restrictions framework studied in Chamberlain (1992). To define the framework, let W, R and Z be random vectors taking values in $\mathcal{W}\subset {\mathbb{R}}^{d_w}$⁠, $\mathcal{R}\subset {\mathbb{R}}^{d_r}$ and $\mathcal{Z}\subset {\mathbb{R}}^{d_z}$⁠, respectively. Assume that Z is a subvector of R and R is a subvector of W, so that $d_z\le d_r\le d_w$⁠. Also, let $\theta \in \mathrm{\Theta }\subset {\mathbb{R}}^{d_{\theta }}$ be a finite‐dimensional parameter whose true value θ₀ is of interest, and let h be a vector‐valued functional nuisance parameter taking values in a convex set of functions $\mathcal{H}$ mapping $\mathcal{Z}$ to ${\mathbb{R}}^{d_h}$⁠, with the true value of h being h₀. The conditional moment restrictions framework assumes that θ₀ and h₀ satisfy the moment conditions  

$$E_P[m(W;{\theta }_0,h_0(Z))\mid R]=0,$$

(2.22)

where $m:\mathcal{W}\times \mathrm{\Theta }\times {\mathbb{R}}^{d_h}\to {\mathbb{R}}^{d_m}$ is a known vector‐valued function. This framework is of interest because it covers a rich variety of models without having to explicitly rely on the likelihood formulation.

To build a Neyman orthogonal score $\psi (W;\theta ,\eta )$ for estimating θ₀, consider the matrix‐valued functional parameter $\mu :\mathcal{R}\to {\mathbb{R}}^{d_{\theta }\times d_m}$ whose true value is given by  

$${\mu }_0(R)=A{(R)}^{\prime }\mathrm{\Omega }{(R)}^{-1}-G(Z)\mathrm{\Gamma }{(R)}^{\prime }\mathrm{\Omega }{(R)}^{-1},$$

(2.23)

where the moment selection matrix‐valued function $A:\mathcal{R}\to {\mathbb{R}}^{d_m\times d_{\theta }}$ and the weighting positive definite matrix‐valued function $\mathrm{\Omega }:\mathcal{R}\to {\mathbb{R}}^{d_m\times d_m}$ can be chosen arbitrarily, and the matrix‐valued functions $\mathrm{\Gamma }:\mathcal{R}\to {\mathbb{R}}^{d_m\times d_{\theta }}$ and $G:\mathcal{Z}\to {\mathbb{R}}^{d_{\theta }\times d_m}$ are given by  

$$\mathrm{\Gamma }(R)={\partial }_{v^{\prime }}{E}_P[m(W;{\theta }_0,v)\mid R]{|}_{v=h_0(Z)}$$

(2.24)

and  

$$G(Z)=E_P[A{(R)}^{\prime }\mathrm{\Omega }{(R)}^{-1}\mathrm{\Gamma }(R)\mid Z]\times {(E_P[\mathrm{\Gamma }{(R)}^{\prime }\mathrm{\Omega }{(R)}^{-1}\mathrm{\Gamma }(R)\mid Z])}^{-1}.$$

(2.25)

Note that μ₀ in 2.23 is well defined even though the right‐hand side of 2.23 contains both R and Z as Z is a subvector of R. Then a Neyman orthogonal score is  

$$\psi (W;\theta ,\eta )=\mu (R)m(W;\theta ,h(Z)),$$

(2.26)

where the nuisance parameter is  

$$\eta =(\mu ,h)\in T={\mathcal{L}}^1(\mathcal{R};{\mathbb{R}}^{d_{\theta }\times d_m})\times \mathcal{H}.$$

Here, ${\mathcal{L}}^1(\mathcal{R};{\mathbb{R}}^{d_{\theta }\times d_m})$ is the vector space of matrix‐valued functions $f:\mathcal{R}\to {\mathbb{R}}^{d_{\theta }\times d_m}$ satisfying $E_P[\parallel f(R)\parallel ]<\infty$⁠. Also, note that even though the matrix‐valued functions A and Ω can be chosen arbitrarily, setting  

$$A(R)={\partial }_{{\theta }^{\prime }}{E}_P[m(W;\theta ,h_0(Z))\mid R]{|}_{\theta ={\theta }_0}$$

(2.27)

and  

$$\mathrm{\Omega }(R)=E_P[m(W;{\theta }_0,h_0(Z))m{(W;{\theta }_0,h_0(Z))}^{\prime }\mid R]$$

(2.28)

leads to an asymptotic variance equal to the semi‐parametric bound of Chamberlain (1992). Let ${\eta }_0=({\mu }_0,h_0)$ be the true value of the nuisance parameter $\eta =(\mu ,h)$⁠. The following lemma shows that the score ψ in 2.26 satisfies the Neyman orthogonality condition.

 

Lemma 2.6.Neyman Orthogonal Scores for Conditional Moment Settings

Suppose that (a) 2.22 holds, (b) the matrices $E_P{[\parallel \mathrm{\Gamma }(R)\parallel }^4]$⁠, $E_P{[\parallel G(Z)\parallel }^4]$⁠, $E_P{[\parallel A(R)\parallel }^2]$ and $E_P{[\parallel \mathrm{\Omega }(R)\parallel }^{-2}]$ are finite, and (c) for all $h\in \mathcal{H}$⁠, there exists a constant $C_h>0$ such that ${Pr}_P(E_P[\parallel m(W;{\theta }_0,h(Z))\parallel \mid R]\le C_h)=1$⁠. Then the score ψ in 2.26 is Neyman orthogonal at $({\theta }_0,{\eta }_0)$ with respect to the nuisance realization set ${\mathcal{T}}_N=T$⁠.

As an application of the conditional moment restrictions framework, let us derive Neyman orthogonal scores in the PLR example using this framework. The PLR model 1.1 is equivalent to  

$$E_P[Y-D{\theta }_0-g_0(X)\mid X,D]=0,$$

which can be written in the form of the conditional moment restrictions framework 2.22 with $W={(Y,D,X^{\prime })}^{\prime }$⁠, $R={(D,X^{\prime })}^{\prime }$⁠, $Z=X$⁠, $h(Z)=g(X)$ and $m(W;\theta ,v)=Y-D\theta -v$⁠. Hence, using 2.27 and 2.28 and denoting $\sigma {(D,X)}^2=E_P[U^2\mid D,X]$ for $U=Y-D{\theta }_0-g_0(X)$⁠, we can take  

$$A(R)=-D,\mathrm{\Omega }(R)=E_P[U^2\mid D,X]=\sigma {(D,X)}^2.$$

With this choice of $A(R)$ and $\mathrm{\Omega }(R)$⁠, we have  

$$\mathrm{\Gamma }(R)=-1,G(Z)=\left(E_P\left[\frac{D}{\sigma {(D,X)}^2}\mid X\right]\right)\times {\left(E_P\left[\frac{1}{\sigma {(D,X)}^2}\mid X\right]\right)}^{-1},$$

and so 2.23 and 2.26 give  

$$\begin{array}{ccc}\hfill \psi (W;\theta ,{\eta }_0)& =& \frac{1}{\sigma {(D,X)}^2}\left(D-\frac{E_P[(D/\sigma {(D,X)}^2)\mid X]}{E_P[(1/\sigma {(D,X)}^2)\mid X]}\right)\hfill \\ & & \times (Y-D\theta -g_0(X)).\hfill \end{array}$$

By construction, the score ψ above is efficient and Neyman orthogonal. Note, however, that using this score would require estimating the heteroscedasticity function $\sigma {(D,X)}^2$⁠, which would require the imposition of some additional smoothness assumptions over this conditional variance function. Instead, if are willing to give up on efficiency to gain some robustness, we can take  

$$A(R)=-D,\mathrm{\Omega }(R)=1;$$

in which case we have  

$$\mathrm{\Gamma }(R)=-1,G(Z)=E_P[D\mid X].$$

Equations 2.23 and 2.26 then give  

$$\begin{array}{ccc}\hfill \psi (W;\theta ,{\eta }_0)& =& (D-E_P[D\mid X])\times (Y-D\theta -g_0(X))\hfill \\ & =& (D-m_0(X))\times (Y-D\theta -g_0(X)).\hfill \end{array}$$

This score ψ is Neyman orthogonal and corresponds to the estimator of θ₀ described in the Introduction in 1.5. Note, however, that this score ψ is efficient only if $\sigma (X,D)$ is a constant.

2.2.5. Neyman Orthogonal Scores and Influence Functions

Neyman orthogonality is a joint property of the score $\psi (W;\theta ,\eta )$⁠, the true parameter value η₀, the parameter set T, and the distribution of W. It is not determined by any particular model for the parameter θ. Nevertheless, it is possible to use semi‐parametric efficiency calculations to construct the orthogonal score from the original score as in Chernozhukov et al. (2016). Specifically, an orthogonal score can be constructed by adding to the original score the influence function adjustment for estimation of the nuisance functions that is analysed in Newey (1994). The resulting orthogonal score will be the influence function of the limit of the average of the original score.

To explain, consider the original score $\varphi (W;\theta ,\beta )$⁠, where β is some function, and let ${\widehat{\beta }}_0$ be a non‐parametric estimator of β₀, the true value of β. Here, β is implicitly allowed to depend on θ, though we suppress that dependence for notational convenience. The corresponding orthogonal score can be formed when there is $\phi (W;\theta ,\eta )$ such that  

$${\displaystyle \int \varphi (w;{\theta }_0,{\widehat{\beta }}_0)\mathit{dP}(w)=\frac{1}{n}\sum _{i=1}^n\phi (W_i;{\theta }_0,{\eta }_0)+o_P(n^{-1/2}),}$$

*

where η is a vector of nuisance functions that includes β, and $\phi (W;\theta ,\eta )$ is an adjustment for the presence of the estimated function ${\widehat{\beta }}_0$ in the original score $\varphi (W;\theta ,\beta )$⁠. The decomposition * typically holds when $\widehat{\beta }$ is either a kernel or a series estimator with a suitably chosen tuning parameter. The Neyman orthogonal score is given by  

$$\psi (W;\theta ,\eta )=\varphi (W;\theta ,\beta )+\phi (W;\theta ,\eta ).$$

(2.29)

Here, $\psi (W;{\theta }_0,{\eta }_0)$ is the influence function of the limit of ${\displaystyle n^{-1}\sum _{i=1}^n\varphi (W_i;{\theta }_0,{\widehat{\beta }}_0)}$⁠, as analysed in Newey (1994), with the restriction $E_P[\psi (W;{\theta }_0,{\eta }_0)]=0$ identifying θ₀.

The form of the adjustment term $\phi (W;\theta ,\eta )$ depends on the estimator ${\widehat{\beta }}_0$ and, of course, on the form of $\varphi (W;\theta ,\beta )$⁠. Such adjustment terms have been derived for various ${\widehat{\beta }}_0$ by Newey (1994). Also, Ichimura and Newey (2015) show how the adjustment term can be computed from the limit of a certain derivative. Any of these results can be applied to a particular starting score $\varphi (W;\theta ,\beta )$ and estimator ${\widehat{\beta }}_0$ to obtain an orthogonal score.

For example, consider again the partially linear model with the original score  

$$\varphi (W;\theta ,\beta )=D(Y-D\theta -g_0(X)).$$

Here, ${\widehat{\beta }}_0={\widehat{g}}_0$ is a non‐parametric regression estimator. From Newey (1994), we know that we obtain the influence function adjustment by taking the conditional expectation of the derivative of the score with respect to $g_0(x)$ (obtaining $-m_0(X)=-E_P[D|X]$⁠) and multiplying the result by the non‐parametric residual to obtain  

$$\phi (W,\theta ,\eta )=-m_0(X)\{Y-D\theta -\beta (X,\theta )\}.$$

The corresponding orthogonal score is then simply  

$$\begin{array}{ccc}\hfill \psi (W;\theta ,\eta )& =& \{D-m_0(X)\}\{Y-D\theta -\beta (X,\theta )\},\hfill \\ \hfill {\beta }_0(X,\theta )& =& E_P[Y-D\theta |X],m_0(X)=E_P[D|X],\hfill \end{array}$$

illustrating that an orthogonal score for the partially linear model can be derived from an influence function adjustment.

Influence functions have been used to estimate functionals of non‐parametric estimators by Hasminskii and Ibragimov (1979) and Bickel and Ritov (1988). Newey et al. (1998, 2004) showed that ${\displaystyle n^{-1/2}\sum _{i=1}^n\psi (W_i;{\theta }_0,{\widehat{\eta }}_0)}$ from 2.29 will have a second‐order remainder in ${\widehat{\eta }}_0$⁠, which is the key asymptotic property of orthogonal scores. Orthogonality of influence functions in semi‐parametric models follows from van der Vaart (1991), as shown for higher‐order counterparts in Robins et al. (2008, 2017). Chernozhukov et al. (2016) point out that, in general, an orthogonal score can be constructed from an original score and non‐parametric estimator ${\widehat{\beta }}_0$ by adding to the original score the adjustment term for estimation of β₀ as described above. This construction provides a way of obtaining an orthogonal score from any initial score $\varphi (W;\theta ,\beta )$ and non‐parametric estimator ${\widehat{\beta }}_0$⁠.

3. Dml: Post‐Regularized Inference Based on Neyman‐Orthogonal Estimating Equations

3.1. Definition of Dml and Its Basic Properties

We assume that we have a sample ${(W_i)}_{i=1}^N$⁠, modelled as independent and identically distributed (i.i.d.) copies of W, whose law is determined by the probability measure P on $\mathcal{W}$⁠. Estimation will be carried out using the finite‐sample analogue of the estimating equation 2.1.

We assume that the true value η₀ of the nuisance parameter η can be estimated by ${\widehat{\eta }}_0$ using a part of the data ${(W_i)}_{i=1}^N$⁠. Different structured assumptions on η₀ allow us to use different machine‐learning tools for estimating η₀, for example:

• (1)

  approximate sparsity for η₀ with respect to some dictionary calls for the use of forward selection, lasso, post‐lasso, ℓ₂‐boosting, or some other sparsity‐based technique;

• (2)

  well‐approximability of η₀ by trees calls for the use of regression trees and random forests;

• (3)

  well‐approximability of η₀ by sparse neural and deep neural nets calls for the use of ℓ₁‐penalized neural and deep neural networks;

• (4)

  well‐approximability of η₀ by at least one model mentioned in (1)–(3) above calls for the use of an ensemble/aggregated method over the estimation methods mentioned in (1)–(3).

There are performance guarantees for most of these ML methods that make it possible to satisfy the conditions stated below. Ensemble and aggregation methods ensure that the performance guarantee is approximately no worse than the performance of the best method.

We assume that N is divisible by K in order to simplify the notation. The following algorithm defines the simple cross‐fitted DML as outlined in the Introduction.  

Definition 3.1.Dml1

(a) Take a K‐fold random partition ${(I_k)}_{k=1}^K$ of observation indices $[N]=\{1,\dots ,N\}$ such that the size of each fold $I_k$ is $n=N/K$⁠. Also, for each $k\in [K]=\{1,\dots ,K\}$⁠, define $I_k^c:=\{1,\dots ,N\}\setminus I_k$⁠. (b) For each $k\in [K]$⁠, construct an ML estimator  

$${\widehat{\eta }}_{0,k}={\widehat{\eta }}_0({(W_i)}_{i\in I_k^c})$$

of η₀, where ${\widehat{\eta }}_{0,k}$ is a random element in T, and where randomness depends only on the subset of data indexed by $I_k^c$⁠. (c) For each $k\in [K]$⁠, construct the estimator ${\check{\theta }}_{0,k}$ as the solution of the following equation:  

$${\mathbb{E}}_{n,k}[\psi (W;{\check{\theta }}_{0,k},{\widehat{\eta }}_{0,k}]=0,$$

(3.1)

where ψ is the Neyman orthogonal score, and $E_{n,k}$ is the empirical expectation over the kth fold of the data; that is, ${\displaystyle E_{n,k}[\psi (W)]=n^{-1}\sum _{i\in I_k}\psi (W_i)}$⁠. If achievement of exact 0 is not possible, define the estimator ${\check{\theta }}_{0,k}$ of θ₀ as an approximate ${\epsilon }_N$‐solution:  

$$\parallel E_{n,k}[\psi (W;{\check{\theta }}_{0,k},{\widehat{\eta }}_{0,k})]\parallel \le \underset{\theta \in \mathrm{\Theta }}{inf}\parallel E_{n,k}[\psi (W;\theta ,{\widehat{\eta }}_{0,k})]\parallel +{\epsilon }_N,{\epsilon }_N=o({\delta }_N{N}^{-1/2}),$$

(3.2)

where ${({\delta }_N)}_{N\ge 1}$ is some sequence of positive constants converging to zero. (4) Aggregate the estimators:  

$${\displaystyle {\stackrel{\sim }{\theta }}_0=\frac{1}{K}\sum _{k=1}^K{\check{\theta }}_{0,k}.}$$

(3.3)

This approach generalizes the 50–50 cross‐fitting method mentioned in the Introduction. We now define a variation of this basic cross‐fitting approach that may behave better in small samples.  

Definition 3.2.Dml2

(a) Take a K‐fold random partition ${(I_k)}_{k=1}^K$ of observation indices $[N]=\{1,\dots ,N\}$ such that the size of each fold $I_k$ is $n=N/K$⁠. Also, for each $k\in [K]=\{1,\dots ,K\}$⁠, define $I_k^c:=\{1,\dots ,N\}\setminus I_k$⁠. (b) For each $k\in [K]$⁠, construct an ML estimator  

$${\widehat{\eta }}_{0,k}={\widehat{\eta }}_0({(W_i)}_{i\in I_k^c})$$

of η₀, where ${\widehat{\eta }}_{0,k}$ is a random element in T, and where randomness depends only on the subset of data indexed by $I_k^c$⁠. (c) Construct the estimator ${\stackrel{\sim }{\theta }}_0$ as the solution to  

$${\displaystyle \frac{1}{K}\sum _{k=1}^K{E}_{n,k}[\psi (W;{\stackrel{\sim }{\theta }}_0,{\widehat{\eta }}_{0,k})]=0,}$$

(3.4)

where ψ is the Neyman orthogonal score, and $E_{n,k}$ is the empirical expectation over the kth fold of the data; that is, ${\displaystyle E_{n,k}[\psi (W)]=n^{-1}\sum _{i\in I_k}\psi (W_i)}$⁠. If achievement of exact 0 is not possible, define the estimator ${\stackrel{\sim }{\theta }}_0$ of θ₀ as an approximate ${\epsilon }_N$‐solution,  

$${\displaystyle ‖\frac{1}{K}\sum _{k=1}^K{E}_{n,k}[\psi (W;{\stackrel{\sim }{\theta }}_0,{\widehat{\eta }}_{0,k})]]‖\le \underset{\theta \in \mathrm{\Theta }}{inf}‖\frac{1}{K}\sum _{k=1}^K{E}_{n,k}[\psi (W;{\theta }_0,{\widehat{\eta }}_{0,k})]]‖+{\epsilon }_N,}$$

(3.5)

for ${\epsilon }_N=o({\delta }_N{N}^{-1/2})$⁠, where ${({\delta }_N)}_{N\ge 1}$ is some sequence of positive constants converging to zero.

 

Remark 3.1. Recommendations

The choice of K has no asymptotic impact under our conditions but, of course, the choice of K may matter in small samples. Intuitively, larger values of K provide more observations in $I_k^c$ from which to estimate the high‐dimensional nuisance functions, which seems to be the more difficult part of the problem. We have found moderate values of K, such as 4 or 5, to work better than $K=2$ in a variety of empirical examples and in simulations. Moreover, we generally recommend DML2 over DML1 though in some problems such as estimation of ATE in the interactive model, which we discuss later, there is no difference between the two approaches. In most other problems, DML2 is better behaved because the pooled empirical Jacobian for the equation in 3.4 exhibits more stable behaviour than the separate empirical Jacobians for the equation in 3.1.

3.2. Moment Condition Models With Linear Scores

We first consider the case of linear scores, where  

$$\psi (w;\theta ,\eta )={\psi }^a(w;\eta )\theta +{\psi }^b(w;\eta ),\text{for}\text{all}w\in \mathcal{W},\theta \in \mathrm{\Theta },\eta \in T.$$

(3.6)

Let $c_0>0$⁠, $c_1>0$⁠, $s>0$ and $q>2$ be some finite constants such that $c_0\le c_1$⁠, and let ${\{{\delta }_N\}}_{N\ge 1}$ and ${\{{\mathrm{\Delta }}_N\}}_{N\ge 1}$ be some sequences of positive constants converging to zero such that ${\delta }_N\ge N^{-1/2}$⁠. Also, let $K\ge 2$ be some fixed integer, and let ${\{{\mathcal{P}}_N\}}_{N\ge 1}$ be some sequence of sets of probability distributions P of W on $\mathcal{W}$⁠.  

Assumption 3.1. Linear Scores With Approximate Neyman Orthogonality

For all $N\ge 3$ and $P\in {\mathcal{P}}_N$⁠, the following conditions hold. (a) The true parameter value θ₀ obeys 2.1. (b) The score ψ is linear in the sense of 3.6. (c) The map $\eta \mapsto E_P[\psi (W;\theta ,\eta )]$ is twice continuously Gateaux‐differentiable on T. (d) The score ψ obeys the Neyman orthogonality or, more generally, the Neyman ${\lambda }_N$ near‐orthogonality condition at $({\theta }_0,{\eta }_0)$ with respect to the nuisance realization set ${\mathcal{T}}_N\subset T$ for  

$${\lambda }_N:=\underset{\eta \in {\mathcal{T}}_N}{sup}‖{\partial }_{\eta }{E}_P\psi (W;{\theta }_0,{\eta }_0)[\eta -{\eta }_0]‖\le {\delta }_N{N}^{-1/2}.$$

(e) The identification condition holds; namely, the singular values of the matrix  

$$J_0:=E_P[{\psi }^a(W;{\eta }_0)]$$

are between c₀ and c₁.

Assumption 3.1 requires scores to be Neyman orthogonal or near‐orthogonal and imposes mild smoothness requirements as well as the canonical identification condition.

 

Assumption 3.2. Score Regularity and Quality of Nuisance Parameter Estimators

For all $N\ge 3$ and $P\in {\mathcal{P}}_N$⁠, the following conditions hold. (a) Given a random subset I of [N] of size $n=N/K$⁠, the nuisance parameter estimator ${\widehat{\eta }}_0={\widehat{\eta }}_0({(W_i)}_{i\in I^c})$ belongs to the realization set ${\mathcal{T}}_N$ with probability at least $1-{\mathrm{\Delta }}_N$⁠, where ${\mathcal{T}}_N$ contains η₀ and is constrained by the next conditions. (b) The moment conditions hold:  

$$\begin{array}{ccc}\hfill m_N& :=& \underset{\eta \in {\mathcal{T}}_N}{sup}(E_P[\parallel \psi (W;{\theta }_0,\eta ){\parallel }^q{])}^{1/q}\le c_1;\hfill \\ \hfill m_N^{\prime }& :=& \underset{\eta \in {\mathcal{T}}_N}{sup}(E_P[\parallel {\psi }^a{(W;\eta )\parallel }^q{])}^{1/q}\le c_1.\hfill \end{array}$$

(c) The following conditions on the statistical rates $r_N$⁠, $r_N^{\prime }$⁠, and ${\lambda }_N^{\prime }$ hold:  

$$\begin{array}{ccc}\hfill r_N& :=& \underset{\eta \in {\mathcal{T}}_N}{sup}\parallel E_P[{\psi }^a(W;\eta )]-E_P[{\psi }^a(W;{\eta }_0)]\parallel \le {\delta }_N,\hfill \\ \hfill r_N^{\prime }& :=& \underset{\eta \in {\mathcal{T}}_N}{sup}(E_P[\parallel \psi (W;{\theta }_0,\eta )-\psi (W;{\theta }_0,{\eta }_0){\parallel }^2{])}^{1/2}\le {\delta }_N,\hfill \\ \hfill {\lambda }_N^{\prime }& :=& \underset{r\in (0,1),\eta \in {\mathcal{T}}_N}{sup}\parallel {\partial }_r^2{E}_P[\psi (W;{\theta }_0,{\eta }_0+r(\eta -{\eta }_0))]\parallel \le {\delta }_N/\sqrt{N}.\hfill \end{array}$$

(d) The variance of the score ψ is non‐degenerate: All eigenvalues of the matrix  

$$E_P[\psi (W;{\theta }_0,{\eta }_0)\psi {(W;{\theta }_0,{\eta }_0)}^{\prime }]$$

are bounded from below by c₀.

Assumptions 3.2(a)–(c) state that the estimator of the nuisance parameter belongs to the realization set ${\mathcal{T}}_N\subset T$⁠, which is a shrinking neighbourhood of η₀, that contracts around η₀ with the rate determined by the ‘statistical’ rates $r_N$⁠, $r_N^{\prime }$ and ${\lambda }_N^{\prime }$⁠. These rates are not given in terms of the norm ${\parallel ·\parallel }_T$ on T, but rather are the intrinsic rates that are most connected to the statistical problem at hand. However, in smooth problems, as discussed below this translates, in the worst cases, to the crude requirement that the nuisance parameters are estimated at the rate $o(N^{-1/4})$⁠.

The conditions in Assumption 3.2 embody refined requirements on the quality of nuisance parameter estimators. In many applications, where $(\theta ,\eta )\mapsto \psi (W;\theta ,\eta )$ is smooth, we can bound  

$$r_N\lesssim {\epsilon }_N,r_N^{\prime }\lesssim {\epsilon }_N,{\lambda }_N^{\prime }\lesssim {\epsilon }_N^2,$$

(3.7)

where ${\epsilon }_N$ is the upper bound on the rate of convergence of ${\widehat{\eta }}_0$ to η₀ with respect to the norm ${\parallel ·\parallel }_T={\parallel ·\parallel }_{P,2}$⁠:  

$$\parallel {\widehat{\eta }}_0{-\eta \parallel }_T\lesssim {\epsilon }_N.$$

Note that ${\mathcal{T}}_N$ can be chosen as the set of η that is within a neighbourhood of size ${\epsilon }_N$ of η₀, possibly with other restrictions, in this case. If only 3.7 holds, Assumption 3.2, particularly ${\lambda }_N^{\prime }=o(N^{-1/2})$⁠, imposes the (crude) rate requirement  

$${\epsilon }_N=o(N^{-1/4}).$$

(3.8)

This rate is achievable for many ML methods under structured assumptions on the nuisance parameters. Among many others, see Bickel et al. (2009), Bühlmann and van de Geer (2011), Belloni et al. (2011, 2012) and Belloni and Chernozhukov (2011, 2013) for ℓ₁‐penalized and related methods in a variety of sparse models, Kozbur (2016) for forward selection in sparse models, Luo and Spindler (2016) for L₂‐boosting in sparse linear models, Wager and Walther (2016) for concentration results for a class of regression trees and random forests, and Chen and White (1999) for a class of neural nets.

However, the presented conditions allow for more refined statements than 3.8. We note that many important structured problems – such as estimation of parameters in PLR models, estimation of parameters in partially linear structural equation models, and estimation of ATEs under unconfoundedness – are such that some cross‐derivatives vanish, allowing more refined requirements than 3.8. This feature allows us to require much finer conditions on the quality of the nuisance parameter estimators than the crude bound 3.8. For example, in many problems  

$${\lambda }_N^{\prime }=0,$$

(3.9)

because the second derivatives vanish,  

$${\partial }_r^2{E}_P[\psi (W;{\theta }_0,{\eta }_0+r(\eta -{\eta }_0))]=0.$$

This occurs in the following important examples:

• (1) the optimal instrument problem (see Belloni et al., 2012);

• (2) the PLR model when $m_0(X)=0$ or is otherwise known (see Section 4);

• (3) the treatment effect examples when the propensity score is known, which includes randomized control trials as an important special case (see Section 5).

If both 3.7 and 3.9 hold, Assumption 3.2, particularly $r_N=o(1)$ and $r_N^{\prime }=o(1)$⁠, imposes the weakest possible rate requirement:  

$${\epsilon }_N=o(1).$$

We note that similar refined rates have appeared in the context of estimation of treatment effects in high‐dimensional settings under sparsity; see Farrell (2015) and Athey et al. (2016) and the related discussion in Remark 5.2. Our refined rate results complement this work by applying to a broad class of estimation contexts, including estimation of ATEs, and to a broad set of ML estimators.

 

Theorem 3.1. Properties of The Dml

Suppose that Assumptions 3.1 and 3.2 hold. In addition, suppose that ${\delta }_N\ge N^{-1/2}$ for all $N\ge 1$⁠. Then the DML1 and DML2 estimators ${\stackrel{\sim }{\theta }}_0$ concentrate in a $1/\sqrt{N}$ neighbourhood of θ₀ and are approximately linear and centred Gaussian,  

$${\displaystyle \sqrt{N}{\sigma }^{-1}({\stackrel{\sim }{\theta }}_0-{\theta }_0)=\frac{1}{\sqrt{N}}\sum _{i=1}^N\overline{\psi }(W_i)+O_P({\rho }_N)⇝N(0,{\mathrm{I}}_d),}$$

(3.10)

uniformly over $P\in {\mathcal{P}}_N$⁠, where the size of the remainder term obeys  

$${\rho }_N:=N^{-1/2}+r_N+r_N^{\prime }+N^{1/2}{\lambda }_N+N^{1/2}{\lambda }_N^{\prime }\lesssim {\delta }_N.$$

(3.11)

Here, $\overline{\psi }(·):=-{\sigma }^{-1}{J}_0^{-1}\psi (·,{\theta }_0,{\eta }_0)$ is the influence function, and the approximate variance is  

$${\sigma }^2:=J_0^{-1}{E}_P[\psi (W;{\theta }_0,{\eta }_0)\psi {(W;{\theta }_0,{\eta }_0)}^{\prime }]{(J_0^{-1})}^{\prime }.$$

The result establishes that the estimator based on the orthogonal scores achieves the root‐N rate of convergence and is approximately normally distributed. It is noteworthy that this convergence result, both the rate of concentration and the distributional approximation, holds uniformly with respect to P varying over an expanding class of probability measures ${\mathcal{P}}_N$⁠. This means that the convergence holds under any sequence of probability distributions ${(P_N)}_{N\ge 1}$ with $P_N\in {\mathcal{P}}_N$ for each N, which in turn implies that the results are robust with respect to perturbations of a given P along such sequences. The same property can be shown to fail for methods not based on orthogonal scores.  

Theorem 3.2. Variance Estimator for Dml

Suppose that Assumptions 3.1 and 3.2 hold. In addition, suppose that ${\delta }_N\ge N^{-[(1-2/q)\wedge 1/2]}$ for all $N\ge 1$⁠. Consider the following estimator of the asymptotic variance matrix of $\sqrt{N}({\stackrel{\sim }{\theta }}_0-{\theta }_0)$⁠:  

$${\displaystyle {\widehat{\sigma }}^2={\widehat{J}}_0^{-1}\frac{1}{K}\sum _{k=1}^K{E}_{n,k}[\psi (W;{\stackrel{\sim }{\theta }}_0,{\widehat{\eta }}_{0,k})\psi {(W;{\stackrel{\sim }{\theta }}_0,{\widehat{\eta }}_{0,k})}^{\prime }]{({\widehat{J}}_0^{-1})}^{\prime },}$$

where  

$${\displaystyle {\widehat{J}}_0=\frac{1}{K}\sum _{k=1}^K{E}_{n,k}[{\psi }^a(W;{\widehat{\eta }}_{0,k})],}$$

and ${\stackrel{\sim }{\theta }}_0$ is either the DML1 or the DML2 estimator. This estimator concentrates around the true variance matrix σ²,  

$${\widehat{\sigma }}^2={\sigma }^2+O_P({\varrho }_N),{\varrho }_N:=N^{-[(1-2/q)\wedge 1/2]}+r_N+r_N^{\prime }\lesssim {\delta }_N.$$

Moreover, ${\widehat{\sigma }}^2$ can replace σ² in the statement of Theorem 3.1 with the size of the remainder term updated as ${\rho }_N=N^{-[(1-2/q)\wedge 1/2]}+r_N+r_N^{\prime }+N^{1/2}{\lambda }_N+N^{1/2}{\lambda }_N^{\prime }$⁠.

Theorems 3.1 and 3.2 can be used for standard construction of confidence regions, which are uniformly valid over a large, interesting class of models:  

Corollary 3.1. Uniformly Valid Confidence Bands

Under the conditions of Theorem 3.2, suppose we are interested in the scalar parameter ${\ell }^{\prime }{\theta }_0$ for some $d_{\theta }\times 1$ vector ℓ. Then the confidence interval  

$$\mathrm{CI}:=({\ell }^{\prime }{\stackrel{\sim }{\theta }}_0\pm {\mathrm{\Phi }}^{-1}(1-\alpha /2)\sqrt{{\ell }^{\prime }{\widehat{\sigma }}^2\ell /N})$$

obeys  

$$\underset{P\in {\mathcal{P}}_N}{sup}|{Pr}_P({\ell }^{\prime }{\theta }_0\in \mathrm{CI})-(1-\alpha )|\to 0.$$

Indeed, the above theorem implies that CI obeys ${\mathit{Pr}}_{P_N}({\ell }^{\prime }{\theta }_0\in \mathrm{CI})\to (1-\alpha )$ under any sequence $\{P_N\}\in {\mathcal{P}}_N$⁠, which implies that these claims hold uniformly in $P\in {\mathcal{P}}_N$⁠. For example, one may choose $\{P_N\}$ such that, for some ${\epsilon }_N\to 0$ 

$$\underset{P\in {\mathcal{P}}_N}{sup}|{\mathit{Pr}}_P({\ell }^{\prime }{\theta }_0\in \mathrm{CI})-(1-\alpha )|\le |{\mathit{Pr}}_{P_N}({\ell }^{\prime }{\theta }_0\in \mathrm{CI})-(1-\alpha )|+{\epsilon }_N\to 0.$$

Next we note that the estimators need not be semi‐parametrically efficient, but under some conditions they can be.

 

Corollary 3.2. Cases With Semi‐Parametric Efficiency

Under the conditions of Theorem 3.1, if the score ψ is efficient for estimating θ₀ at a given $P\in \mathcal{P}\subset {\mathcal{P}}_N$⁠, in the semi‐parametric sense as defined in van der Vaart (1998), then the large sample variance ${\sigma }_0^2$ of ${\stackrel{\sim }{\theta }}_0$ reaches the semi‐parametric efficiency bound at this P relative to the model $\mathcal{P}$⁠.

3.3. Models With Non‐Linear Scores

Let $c_0>0$⁠, $c_1>0$⁠, $a>1$⁠, $v>0$⁠, $s>0$ and $q>2$ be some finite constants, and let ${\{{\delta }_N\}}_{N\ge 1}$⁠, ${\{{\mathrm{\Delta }}_N\}}_{N\ge 1}$ and ${\{{\tau }_N\}}_{N\ge 1}$ be some sequences of positive constants converging to zero. To derive the properties of the DML estimator, we will use the following assumptions.  

Assumption 3.3. Non‐Linear Moment Condition Problem With Approximate Neyman Orthogonality

For all $N\ge 3$ and $P\in {\mathcal{P}}_N$⁠, the following conditions hold. (a) The true parameter value θ₀ obeys 2.1, and Θ contains a ball of radius $c_1{N}^{-1/2}logN$ centred at θ₀. (b) The map $(\theta ,\eta )\mapsto E_P[\psi (W;\theta ,\eta )]$ is twice continuously Gateaux‐differentiable on $\mathrm{\Theta }\times T$⁠. (c) For all $\theta \in \mathrm{\Theta }$⁠, the identification relation  

$$2\parallel E_P[\psi (W;\theta ,{\eta }_0)]\parallel \ge \parallel J_0(\theta -{\theta }_0)\parallel \wedge c_0$$

is satisfied, for the Jacobian matrix  

$$J_0:={\partial }_{{\theta }^{\prime }}\{E_P[\psi (W;\theta ,{\eta }_0)]\}{|}_{\theta ={\theta }_0}$$

having singular values between c₀ and c₁. (d) The score ψ obeys the Neyman orthogonality or, more generally, the Neyman near‐orthogonality with ${\lambda }_N={\delta }_N{N}^{-1/2}$ for the set ${\mathcal{T}}_N\subset T$⁠.

Assumption 3.3 is mild and rather standard in moment condition problems. Assumption 3.3(a) requires θ₀ to be sufficiently separated from the boundary of Θ. Assumption 3.3(b) only requires differentiability of the function $(\theta ,\eta )\mapsto E_P[\psi (W;\theta ,\eta )]$ and does not require differentiability of the function $(\theta ,\eta )\mapsto \psi (W;\theta ,\eta )$⁠. Assumption 3.3(c) implies sufficient identifiability of θ₀. Assumption 3.3(d) is the orthogonality condition that has already been extensively discussed.  

Assumption 3.4. Score Regularity and Requirements on The Quality of Estimation of Nuisance Parameters

Let K be a fixed integer. For all $N\ge 3$ and $P\in {\mathcal{P}}_N$⁠, the following conditions hold. (a) Given a random subset I of $\{1,\dots ,N\}$ of size $n=N/K$⁠, we have that the nuisance parameter estimator ${\widehat{\eta }}_0={\widehat{\eta }}_0({(W_i)}_{i\in I^c})$ belongs to the realization set ${\mathcal{T}}_N$ with probability at least $1-{\mathrm{\Delta }}_N$⁠, where ${\mathcal{T}}_N$ contains η₀ and is constrained by conditions given below. (b) The parameter space Θ is bounded and for each $\eta \in {\mathcal{T}}_N$⁠, the function class ${\mathcal{F}}_{1,\eta }=\{{\psi }_j(·,\theta ,\eta ):j=1,\dots ,d_{\theta },\theta \in \mathrm{\Theta }\}$ is suitably measurable and its uniform covering entropy obeys  

$$\underset{Q}{sup}logN(\epsilon \parallel F_{1,\eta }{\parallel }_{Q,2},{\mathcal{F}}_{1,\eta }{,\parallel ·\parallel }_{Q,2})\le vlog(a/\epsilon ),\text{for}\text{all}0<\epsilon \le 1,$$

(3.12)

where $F_{1,\eta }$ is a measurable envelope for ${\mathcal{F}}_{1,\eta }$ that satisfies $\parallel F_{1,\eta }{\parallel }_{P,q}\le c_1$⁠. (c) The following conditions on the statistical rates $r_N$⁠, $r_N^{\prime }$⁠, and ${\lambda }_N^{\prime }$ hold:  

$$\begin{array}{ccc}\hfill r_N& :=& \underset{\eta \in {\mathcal{T}}_N,\theta \in \mathrm{\Theta }}{sup}\parallel E_P[\psi (W;\theta ,\eta )-E_P[\psi (W;\theta ,{\eta }_0)]\parallel \le {\delta }_N{\tau }_N,\hfill \\ \hfill r_N^{\prime }& :=& \underset{\eta \in {\mathcal{T}}_N,\parallel \theta -{\theta }_0\parallel \le {\tau }_N}{sup}(E_P[\parallel \psi (W;\theta ,\eta )-\psi (W;{\theta }_0,{\eta }_0){\parallel }^2{])}^{1/2}\text{and}{r}_N^{\prime }{log}^{1/2}(1/r_N^{\prime })\le {\delta }_N,\hfill \\ \hfill {\lambda }_N^{\prime }& :=& \underset{r\in (0,1),\eta \in {\mathcal{T}}_N,\parallel \theta -{\theta }_0\parallel \le {\tau }_N}{sup}\parallel {\partial }_r^2{E}_P[\psi (W;{\theta }_0,{\eta }_0+r(\theta -{\theta }_0)+r(\eta -{\eta }_0))]\parallel \le {\delta }_N{N}^{-1/2}.\hfill \end{array}$$

(d) The variance of the score is non‐degenerate. All eigenvalues of the matrix  

$$E_P[\psi (W;{\theta }_0,{\eta }_0)\psi {(W;{\theta }_0,{\eta }_0)}^{\prime }]$$

are bounded from below by c₀.

Assumptions 3.3(a)–(c) state that the estimator of the nuisance parameter belongs to the realization set ${\mathcal{T}}_N\subset T$⁠, which is a shrinking neighbourhood of η₀ that contracts at the statistical rates $r_N$⁠, $r_N^{\prime }$ and ${\lambda }_N^{\prime }$⁠. These rates are not given in terms of the norm ${\parallel ·\parallel }_T$ on T, but rather are intrinsic rates that are connected to the statistical problem at hand. In smooth problems, these conditions translate to the crude requirement that nuisance parameters are estimated at the $o(N^{-1/4})$ rate, as discussed previously in the case with linear scores. However, these conditions can be refined as, for example, when ${\lambda }_N^{\prime }=0$ or when some cross‐derivatives vanish in ${\lambda }_N^{\prime }$⁠; see the linear case in the previous subsection for further discussion. Suitable measurability and pointwise entropy conditions, required in Assumption 3.4(b), are mild regularity conditions that are satisfied in all practical cases. The assumption of a bounded parameter space Θ in Assumption 3.4(b) is embedded in the entropy condition, but we state it separately for clarity. This assumption was not needed in the linear case, and it can be removed in the non‐linear case with the imposition of more complicated Huber‐like regularity conditions. Assumption 3.4(c) is a set of mild growth conditions.  

Remark 3.2. Rate Requirements on Nuisance Parameter Estimators

Similar to the discussion in the linear case, the conditions in Assumption 3.4 are very flexible and embody refined requirements on the quality of the nuisance parameter estimators. The conditions essentially reduce to the previous conditions in the linear case, with the exception of compactness, which is imposed to make the conditions easy to verify in non‐linear cases.

 

Theorem 3.3. Properties of The Dml for Non‐Linear Scores

Suppose that Assumptions 3.3 and 3.4 hold. In addition, suppose that ${\delta }_N\ge N^{-1/2+1/q}logN$ and that $N^{-1/2}logN\le {\tau }_N\le {\delta }_N$ for all $N\ge 1$⁠. Then the DML1 and DML2 estimators ${\stackrel{\sim }{\theta }}_0$ concentrate in a $1/\sqrt{N}$ neighbourhood of θ₀, and are approximately linear and centred Gaussian,  

$${\displaystyle \sqrt{N}{\sigma }^{-1}({\stackrel{\sim }{\theta }}_0-{\theta }_0)=\frac{1}{\sqrt{N}}\sum _{i=1}^N\overline{\psi }(W_i)+O_P({\rho }_N)⇝N(0,\mathrm{I}),}$$

uniformly over $P\in {\mathcal{P}}_N$⁠, where the size of the remainder term obeys  

$${\rho }_N:=N^{-1/2+1/q}logN+r_N^{\prime }{log}^{1/2}(1/r_N^{\prime })+N^{1/2}{\lambda }_N+N^{1/2}{\lambda }_N^{\prime }\lesssim {\delta }_N,$$

$\overline{\psi }(·):=-{\sigma }_0^{-1}{J}_0^{-1}\psi (·,{\theta }_0,{\eta }_0)$ is the influence function, and the approximate variance is  

$${\sigma }^2:=J_0^{-1}{E}_P[\psi (W;{\theta }_0,{\eta }_0)\psi {(W;{\theta }_0,{\eta }_0)}^{\prime }]{(J_0^{-1})}^{\prime }.$$

Moreover, in the statement above, σ² can be replaced by a consistent estimator ${\widehat{\sigma }}^2$⁠, obeying ${\widehat{\sigma }}^2={\sigma }^2+o_P({\varrho }_N)$ uniformly in $P\in {\mathcal{P}}_N$⁠, with the size of the remainder term updated as ${\rho }_N={\rho }_N+{\varrho }_N$⁠. Furthermore, Corollaries 3.1 and 3.2 continue to hold under the conditions of this theorem.

3.4. Finite‐Sample Adjustments To Incorporate Uncertainty Induced By Sample Splitting

The estimation technique developed in this paper relies on subsamples obtained by randomly partitioning the sample: an auxiliary sample for estimating the nuisance functions and a main sample for estimating the parameter of interest. Although the specific sample partition has no impact on estimation results asymptotically, the effect of the particular random split on the estimate can be important in finite samples. To make the results more robust with respect to the partitioning, we propose to repeat the DML estimator S times, obtaining the estimates  

$${\stackrel{\sim }{\theta }}_0^s,s=1,\dots ,S.$$

Features of these estimates may then provide insight into the sensitivity of results to the sample splitting, and we can report results that incorporate features of this set of estimates that should be less driven by any particular sample‐splitting realization.

 

Definition 3.3.Incorporating The Impact of Sample Splitting Using Mean and Median Methods

For point estimation, we define  

$${\displaystyle {\stackrel{\sim }{\theta }}_0^{\text{mean}}=\frac{1}{S}\sum _{s=1}^S{\stackrel{\sim }{\theta }}_0^s\text{or}{\stackrel{\sim }{\theta }}_0^{\text{median}}=\text{median}{\{\stackrel{\sim }{{\theta }_0^s}\}}_{s=1}^S,}$$

where the median operation is applied coordinate‐wise. To quantify and incorporate the variation introduced by sample splitting, we consider variance estimators,  

$${\displaystyle {\widehat{\sigma }}^{2,\text{mean}}=\frac{1}{S}\sum _{s=1}^S\left({\widehat{\sigma }}_s^2+({\widehat{\theta }}_s-{\stackrel{\sim }{\theta }}^{\text{mean}}){({\widehat{\theta }}_s-{\stackrel{\sim }{\theta }}^{\text{mean}})}^{\prime }\right),}$$

(3.13)

and a more robust version,  

$${\widehat{\sigma }}^{2,\text{median}}=\text{median}\{{\widehat{\sigma }}_s^2+{(({\widehat{\theta }}_s-{\stackrel{\sim }{\theta }}^{\text{median}}){({\widehat{\theta }}_s-{\stackrel{\sim }{\theta }}^{\text{median}})}^{\prime }\}}_{s=1}^S,$$

(3.14)

where the median picks out the matrix with median operator norm, which preserve non‐negative definiteness.

We recommend using medians, reporting ${\stackrel{\sim }{\theta }}_0^{\text{median}}$ and ${\widehat{\sigma }}^{2,\mathit{median}}$⁠, as these quantities are more robust to outliers.

 

Corollary 3.3.

If S is fixed, as $N\to \infty$ and maintaining either Assumptions 3.1 and 3.2 or Assumptions 3.3 and 3.4 as appropriate, ${\stackrel{\sim }{\theta }}_0^{\text{mean}}$ and ${\stackrel{\sim }{\theta }}_0^{\text{median}}$ are first‐order equivalent to ${\stackrel{\sim }{\theta }}_0$ and obey the conclusions of Theorems 3.1 and 3.2 or of Theorem 3.3. Moreover, ${\widehat{\sigma }}^{2,\text{median}}$ and ${\widehat{\sigma }}^{2,\text{mean}}$ can replace $\widehat{\sigma }$ in the statement of the appropriate theorems.

It would be interesting to investigate the behaviour under the regime where $S\to \infty$ as $N\to \infty$⁠.

4. Inference in Partially Linear Models

4.1. Inference in Partially Linear Regression Models

Here we revisit the PLR model  

$$\begin{array}{ccc}\hfill Y& =& D{\theta }_0+g_0(X)+U,E_P[U\mid X,D]=0,\hfill \end{array}$$

(4.1)

 

$$\begin{array}{ccc}\hfill D& =& m_0(X)+V,E_P[V\mid X]=0.\hfill \end{array}$$

(4.2)

The parameter of interest is the regression coefficient θ₀. If D is conditionally exogenous (as good as randomly assigned conditional on covariates), then θ₀ measures the average causal/treatment effect of D on potential outcomes.

The first approach to inference on θ₀, which we described in the Introduction, is to employ the DML method using the score function  

$$\psi (W;\theta ,\eta ):=\{Y-D\theta -g(X)\}(D-m(X)),\eta =(g,m),$$

(4.3)

where $W=(Y,D,X)$⁠, and g and m are P‐square‐integrable functions mapping the support of X to $\mathbb{R}$⁠. It is easy to see that θ₀ satisfies the moment condition $E_P\psi (W;{\theta }_0,{\eta }_0)=0$⁠, and also the orthogonality condition ${\partial }_{\eta }{E}_P\psi (W;{\theta }_0,{\eta }_0)[\eta -{\eta }_0]=0$ where ${\eta }_0=(g_0,m_0)$⁠.

A second approach employs the Robinson‐style ‘partialling‐out’ score function  

$$\psi (W;\theta ,\eta ):=\{Y-\ell (X)-\theta (D-m(X))\}(D-m(X)),\eta =(\ell ,m),$$

(4.4)

where $W=(Y,D,X)$ and ℓ and m are P‐square‐integrable functions mapping the support of X to $\mathbb{R}$⁠. This gives an alternative parametrization of the score function above, and using this score is first‐order equivalent to using the previous score. It is easy to see that θ₀ satisfies the moment condition $E_P\psi (W;{\theta }_0,{\eta }_0)=0$⁠, and also the orthogonality condition ${\partial }_{\eta }{E}_P\psi (W;{\theta }_0,{\eta }_0)[\eta -{\eta }_0]=0$⁠, for ${\eta }_0=({\ell }_0,m_0)$⁠, where ${\ell }_0(X)=E_P[Y|X]$⁠.

In the partially linear model, the DML approach complements Belloni et al. (2013, 2014a,b, 2015), Zhang and Zhang (2014), van de Geer et al. (2014) and Javanmard and Montanari (2014b), all of which consider estimation and inference for parameters within the partially linear model using lasso‐type methods without cross‐fitting. By relying upon cross‐fitting, the DML approach allows for the use of a much broader collection of ML methods for estimating the nuisance functions and also allows relaxation of sparsity conditions in the case where lasso or other sparsity‐based estimators are used. Both the DML approach and the approaches taken in the aforementioned papers can be seen as heuristically debiasing the score function $(Y-D\theta -g(X))D$⁠, which does not possess the orthogonality property unless $m_0(X)=0$⁠.

Let ${({\delta }_N)}_{n=1}^{\infty }$ and ${({\mathrm{\Delta }}_N)}_{n=1}^{\infty }$ be sequences of positive constants approaching 0 as before. Also, let c, C and q be fixed strictly positive constants such that $q>4$⁠, and let $K\ge 2$ be a fixed integer. Moreover, for any $\eta =({\ell }_1,{\ell }_2)$⁠, where ℓ₁ is a function ℓ₁ and ℓ₂ are functions mapping the support of X to $\mathbb{R}$⁠, denote ${\parallel \eta \parallel }_{P,q}=\parallel {\ell }_1{\parallel }_{P,q}\vee {\parallel {\ell }_2\parallel }_{P,q}$⁠. For simplicity, assume that $N/K$ is an integer.  

Assumption 4.1. Regularity Conditions for Partially Linear Regression Model

Let $\mathcal{P}$ be the collection of probability laws P for the triple $W=(Y,D,X)$ such that (a) 4.1 and 4.2 hold; (b) ${\parallel Y\parallel }_{P,q}+{\parallel D\parallel }_{P,q}\le C$⁠; (c) ${\parallel UV\parallel }_{P,2}\ge c^2$ and $E_P[V^2]\ge c$⁠; (d) $\parallel E_P[U^2\mid X]{\parallel }_{P,\infty }\le C$ and $\parallel E_P[V^2\mid X]{\parallel }_{P,\infty }\le C$⁠; and (e) given a random subset I of [N] of size $n=N/K$⁠, the nuisance parameter estimator ${\widehat{\eta }}_0={\widehat{\eta }}_0({(W_i)}_{i\in I^c})$ obeys the following conditions for all $n\ge 1$⁠. With P‐probability no less than $1-{\mathrm{\Delta }}_N$⁠, $\parallel {\widehat{\eta }}_0-{\eta }_0{\parallel }_{P,q}\le C$⁠, $\parallel {\widehat{\eta }}_0-{\eta }_0{\parallel }_{P,2}\le {\delta }_N$ and (i) for the score ψ in 4.3, where ${\widehat{\eta }}_0=({\widehat{g}}_0,{\widehat{m}}_0)$⁠, $\parallel {\widehat{m}}_0-m_0{\parallel }_{P,2}\times {\parallel {\widehat{g}}_0-g_0\parallel }_{P,2}\le {\delta }_N{N}^{-1/2}$⁠, and (ii) for the score ψ in 4.4, where ${\widehat{\eta }}_0=({\widehat{\ell }}_0,{\widehat{m}}_0)$⁠, $\parallel {\widehat{m}}_0-m_0{\parallel }_{P,2}\times (\parallel {\widehat{m}}_0-m_0{\parallel }_{P,2}+\parallel {\widehat{\ell }}_0-{\ell }_0{\parallel }_{P,2})\le {\delta }_N{N}^{-1/2}$⁠.

 

Remark 4.1. Rate Conditions for Estimators of Nuisance Parameters

The only non‐primitive condition here is the assumption on the rate of estimating the nuisance parameters. These rates of convergence are available for most often used ML methods and are case‐specific, so we do not restate conditions that are needed to reach these rates.

The following theorem follows as a corollary to the results in Section 3 by verifying Assumptions 3.1 and 3.2, and it will be proven as a special case of Theorem 4.2 below.  

Theorem 4.1. Dml Inference on Regression Coefficients in The Partially Linear Regression Model

Suppose that Assumption 4.1 holds. Then the DML1 and DML2 estimators ${\stackrel{\sim }{\theta }}_0$ constructed in Definitions 3.1 and 3.2 using the score in either 4.3 or 4.4 are first‐order equivalent and obey  

$${\sigma }^{-1}\sqrt{N}({\stackrel{\sim }{\theta }}_0-{\theta }_0)⇝N(0,1),$$

uniformly over $P\in \mathcal{P}$⁠, where ${\sigma }^2={(E_P[V^2])}^{-1}{E}_P[V^2{U}^2]{(E_P[V^2])}^{-1}$⁠. Moreover, the result continues to hold if σ² is replaced by ${\widehat{\sigma }}^2$ defined in Theorem 3.2. Consequently, confidence regions based upon the DML estimators ${\stackrel{\sim }{\theta }}_0$ have uniform asymptotic validity:  

$$\underset{N\to \infty }{lim}\underset{P\in \mathcal{P}}{sup}|{Pr}_P({\theta }_0\in [{\stackrel{\sim }{\theta }}_0\pm {\mathrm{\Phi }}^{-1}(1-\alpha /2)\widehat{\sigma }/\sqrt{N}])-(1-\alpha )|=0.$$

 

Remark 4.2. Asymptotic Efficiency Under Homoscedasticity

Under conditional homoscedasticity, i.e. $E[U^2|Z]=E[U^2]$⁠, the asymptotic variance σ² reduces to $E{[V^2]}^{-1}E[U^2]$⁠, which is the semi‐parametric efficiency bound for θ.

 

Remark 4.3. Tightness of Conditions Under Cross‐Fitting

The conditions in Theorem 4.1 are sharp, though they are simplified for ease of presentation. The sharpness can be understood by examining the case where the regression function g₀ and the propensity function m₀ are sparse with sparsity indices $s^g\ll N$ and $s^m\ll N$⁠. They are estimated by ℓ₁‐penalized estimators ${\widehat{g}}_0$ and ${\widehat{m}}_0$ that have sparsity indices of orders $s^g$ and $s^m$ and converge to g₀ and m₀ at the rates $\sqrt{s^g/N}$ and $\sqrt{s^m/N}$ (ignoring logs). The rate conditions in Assumption 4.1 then require (ignoring logs) that  

$$\sqrt{s^g/N}\sqrt{s^m/N}\ll N^{-1/2}\iff s^g{s}^m\ll N,$$

which is much weaker than the condition  

$${(s^g)}^2+{(s^m)}^2\ll N$$

(ignoring logs) required without sample splitting. For example, if the propensity function m₀ is very sparse (low $s_m$⁠), then the regression function is allowed to be quite dense (high $s_g$⁠), and vice versa. If the propensity function is known (⁠$s^m=0$⁠) or can be estimated at the $N^{-1/2}$ rate, then only consistency for ${\widehat{g}}_0$ is needed. Such comparisons also extend to approximately sparse models.

4.2. Inference in Partially Linear Iv Models

Here we extend the PLR model studied in Section 4.1 to allow for IV identification. Specifically, we consider the model  

$$\begin{array}{ccc}\hfill Y& =& D{\theta }_0+g_0(X)+U,E_P[U\mid X,Z]=0,\hfill \end{array}$$

(4.5)

 

$$\begin{array}{ccc}\hfill Z& =& m_0(X)+V,E_P[V\mid X]=0,\hfill \end{array}$$

(4.6)

where Z is the instrumental variable. As before, the parameter of interest is θ and its true value is θ₀. If $Z=D$⁠, the model 4.5–4.6 coincides with 4.1–4.2 but is otherwise different.

To estimate θ₀ and to perform inference on it, we use the score  

$$\psi (W;\theta ,\eta ):=(Y-D\theta -g(X))(Z-m(X)),\eta =(g,m),$$

(4.7)

where $W=(Y,D,X,Z)$ and g and m are P‐square‐integrable functions mapping the support of X to $\mathbb{R}$⁠. Alternatively, we can use the Robinson‐style score  

$$\psi (W;\theta ,\eta ):=(Y-\ell (X)-\theta (D-r(X)))(Z-m(X)),\eta =(\ell ,m,r),$$

(4.8)

where $W=(Y,D,X,Z)$ and ℓ, m and r are P‐square‐integrable functions mapping the support of X to $\mathbb{R}$⁠. It is straightforward to verify that both scores satisfy the moment condition $E_P[\psi (W;{\theta }_0,{\eta }_0)]=0$ and also the orthogonality condition ${\partial }_{\eta }{E}_P[\psi (W;{\theta }_0,{\eta }_0)][\eta -{\eta }_0]=0$⁠, for ${\eta }_0=(g_0,m_0)$ in the former case and ${\eta }_0=({\ell }_0,m_0,r_0)$ for ℓ₀ and r₀ defined by ${\ell }_0(X)=E_P[Y\mid X]$ and $r_0(X)=E_P[D\mid X]$⁠, respectively, in the latter case.⁸

Note that the score in 4.8 has a minor advantage over the score in 4.7 because all of its nuisance parameters are conditional mean functions, which can be directly estimated by the ML methods. If one prefers to use the score in 4.7, one has to construct an estimator of g₀ first. To do so, one can first obtain a DML estimator of θ₀ based on the score in 4.8, say ${\stackrel{\sim }{\theta }}_0$⁠. Then, using the fact that $g_0(X)=E_P[Y-D{\theta }_0\mid X]$⁠, one can construct an estimator ${\widehat{g}}_0$ by applying an ML method to regress $Y-D{\stackrel{\sim }{\theta }}_0$ on X. Alternatively, one can use assumption‐specific methods to directly estimate g₀, without using the score 4.8 first. For example, if g₀ can be approximated by a sparse linear combination of a large set of transformations of X, one can use the methods of Gautier and Tsybakov (2014) to obtain an estimator of g₀.

Let ${({\delta }_N)}_{n=1}^{\infty }$ and ${({\mathrm{\Delta }}_N)}_{n=1}^{\infty }$ be sequences of positive constants approaching 0 as before. Also, let c, C and q be fixed strictly positive constants such that $q>4$⁠, and let $K\ge 2$ be a fixed integer. Moreover, for any $\eta ={({\ell }_j)}_{j=1}^l$ mapping the support of X to ${\mathbb{R}}^l$⁠, denote ${\parallel \eta \parallel }_{P,q}={max}_{1\le j\le l}{\parallel {\ell }_j\parallel }_{P,q}$⁠. For simplicity, assume that $N/K$ is an integer.  

Assumption 4.2. Regularity Conditions for Partially Linear Iv Model

For all probability laws $P\in \mathcal{P}$ for the quadruple $W=(Y,D,X,Z)$ the following conditions hold: (a) equations 4.5 and 4.6 hold; (b) ${\parallel Y\parallel }_{P,q}+{\parallel D\parallel }_{P,q}+{\parallel Z\parallel }_{P,q}\le C$⁠; (c) ${\parallel UV\parallel }_{P,2}\ge c^2$ and $|E_P[DV]|\ge c$⁠; (d) $\parallel E_P[U^2\mid X]{\parallel }_{P,\infty }\le C$ and $\parallel E_P[V^2\mid X]{\parallel }_{P,\infty }\le C$⁠; and (e) given a random subset I of [N] of size $n=N/K$⁠, the nuisance parameter estimator ${\widehat{\eta }}_0={\widehat{\eta }}_0({(W_i)}_{i\in I^c})$ obeys the following conditions. With P‐probability no less than $1-{\mathrm{\Delta }}_N$⁠, $\parallel {\widehat{\eta }}_0-{\eta }_0{\parallel }_{P,q}\le C$⁠, $\parallel {\widehat{\eta }}_0-{\eta }_0{\parallel }_{P,2}\le {\delta }_N$ and (i) for the score ψ in 4.7, where ${\widehat{\eta }}_0=({\widehat{g}}_0,{\widehat{m}}_0)$⁠, $\parallel {\widehat{m}}_0-m_0{\parallel }_{P,2}\times {\parallel {\widehat{g}}_0-g_0\parallel }_{P,2}\le {\delta }_N{N}^{-1/2}$⁠, and (ii) for the score ψ in 4.8, where ${\widehat{\eta }}_0=({\widehat{\ell }}_0,{\widehat{m}}_0,{\widehat{r}}_0)$⁠, $\parallel {\widehat{m}}_0-m_0{\parallel }_{P,2}\times \left(\parallel {\widehat{r}}_0-r_0{\parallel }_{P,2}+{\parallel {\widehat{\ell }}_0-{\ell }_0\parallel }_{P,2}\right)\le {\delta }_N{N}^{-1/2}$⁠.

The following theorem follows as a corollary to the results in Section 3 by verifying Assumptions 3.1 and 3.2.  

Theorem 4.2. Dml Inference on Regression Coefficients in The Partially Linear Iv Model

Suppose that Assumption 4.2 holds. Then the DML1 and DML2 estimators ${\stackrel{\sim }{\theta }}_0$ constructed in Definitions 3.1 and 3.2 using the score in either 4.7 or 4.8 are first‐order equivalent and obey ${\sigma }^{-1}\sqrt{N}({\stackrel{\sim }{\theta }}_0-{\theta }_0)⇝N(0,1)$ uniformly over $P\in \mathcal{P}$⁠, where ${\sigma }^2={(E_P[DV])}^{-1}{E}_P[V^2{U}^2]{(E_P[DV])}^{-1}$⁠. Moreover, the result continues to hold if σ² is replaced by ${\widehat{\sigma }}^2$ defined in Theorem 3.2. Consequently, confidence regions based upon the DML estimators ${\stackrel{\sim }{\theta }}_0$ have uniform asymptotic validity:  

$$\underset{N\to \infty }{lim}\underset{P\in \mathcal{P}}{sup}|{Pr}_P({\theta }_0\in [{\stackrel{\sim }{\theta }}_0\pm {\mathrm{\Phi }}^{-1}(1-\alpha /2)\widehat{\sigma }/\sqrt{N}])-(1-\alpha )|=0.$$

5. Inference on Treatment Effects in The Interactive Model

5.1. Inference on Ate and Atte

In this section, we specialize the results of Section 3 to estimate treatment effects under the unconfoundedness assumption of Rosenbaum and Rubin (1983). Within this setting, there is a large classical literature focused on low‐dimensional settings that provides methods for adjusting for confounding variables including regression methods, propensity score adjustment methods, matching methods, and ‘doubly‐robust’ combinations of these methods; see, e.g. Robins and Rotnitzky (1995), Hahn (1998), Hirano et al. (2003) and Abadie and Imbens (2006) as well as the textbook overview provided in Imbens and Rubin (2015). In this section, we present results that complement this important classic work as well as the rapidly expanding body of work on estimation under unconfoundedness using ML methods; see, among others, Athey et al. (2016), Belloni et al. (2014a, 2017), Farrell (2015) and Imai and Ratkovic (2013).

We specifically consider estimation of ATEs when treatment effects are fully heterogeneous and the treatment variable is binary, $D\in \{0,1\}$⁠. We consider vectors $(Y,D,X)$ such that  

$$\begin{array}{ccc}\hfill Y& =& g_0(D,X)+U,E_P[U\mid X,D]=0,\hfill \end{array}$$

(5.1)

 

$$\begin{array}{ccc}\hfill D& =& m_0(X)+V,E_P[V\mid X]=0.\hfill \end{array}$$

(5.2)

Because D is not additively separable, this model is more general than the partially linear model for the case of binary D. A common target parameter of interest in this model is the ATE:⁹ 

$${\theta }_0=E_P[g_0(1,X)-g_0(0,X)].$$

Another common target parameter is the ATTE:  

$${\theta }_0=E_P[g_0(1,X)-g_0(0,X)|D=1].$$

The confounding factors X affect the policy variable via the propensity score $m_0(X)$ and the outcome variable via the function $g_0(D,X)$⁠. Both of these functions are unknown and potentially complicated, and we can employ ML methods to learn them.

We proceed to set up moment conditions with scores obeying orthogonality conditions. For estimation of the ATE, we employ  

$$\psi (W;\theta ,\eta ):=(g(1,X)-g(0,X))+\frac{D(Y-g(1,X))}{m(X)}-\frac{(1-D)(Y-g(0,X))}{1-m(X)}-\theta ,$$

(5.3)

where the nuisance parameter $\eta =(g,m)$ consists of P‐square‐integrable functions g and m mapping the support of $(D,X)$ to $\mathbb{R}$ and the support of X to $(\epsilon ,1-\epsilon )$⁠, respectively, for some $\epsilon \in (0,1/2)$⁠. The true value of η is ${\eta }_0=(g_0,m_0)$⁠. This orthogonal moment condition is based on the influence function for the mean for missing data from Robins and Rotnitzky (1995).

For estimation of the ATTE, we use the score  

$$\psi (W;\theta ,\eta )=\frac{D(Y-\overline{g}(X))}{p}-\frac{m(X)(1-D)(Y-\overline{g}(X))}{p(1-m(X))}-\frac{D\theta }{p},$$

(5.4)

where the nuisance parameter $\eta =(\overline{g},m,p)$ consists of P‐square‐integrable functions $\overline{g}$ and m mapping the support of X to $\mathbb{R}$ and to $(\epsilon ,1-\epsilon )$⁠, respectively, and a constant $p\in (\epsilon ,1-\epsilon )$⁠, for some $\epsilon \in (0,1/2)$⁠. The true value of η is ${\eta }_0=({\overline{g}}_0,m_0,p_0)$⁠, where ${\overline{g}}_0(X)=g_0(0,X)$ and $p_0=E_P[D]$⁠. Note that estimating ATTE does not require estimating $g_0(1,X)$⁠. Note also that because p is a constant, it does not affect the DML estimators ${\stackrel{\sim }{\theta }}_0$ based on the score ψ in 5.4, but having p simplifies the formula for the variance of ${\stackrel{\sim }{\theta }}_0$⁠.

Using their respective scores, it can be easily seen that true parameter values θ₀ for ATE and ATTE obey the moment condition $E_P[\psi (W;{\theta }_0,{\eta }_0)]=0$⁠, and also that the orthogonality condition ${\partial }_{\eta }{E}_P[\psi (W;{\theta }_0,{\eta }_0)][\eta -{\eta }_0]=0$ holds.

Let ${({\delta }_N)}_{n=1}^{\infty }$ and ${({\mathrm{\Delta }}_N)}_{n=1}^{\infty }$ be sequences of positive constants approaching 0. Also, let $c,\epsilon ,C$ and q be fixed strictly positive constants such that $q>2$⁠, and let $K\ge 2$ be a fixed integer. Moreover, for any $\eta =({\ell }_1,\dots ,{\ell }_l)$⁠, denote ${\parallel \eta \parallel }_{P,q}={max}_{1\le j\le l}{\parallel {\ell }_j\parallel }_{P,q}$⁠. For simplicity, assume that $N/K$ is an integer.  

Assumption 5.1. Regularity Conditions for Ate and Atte Estimation

For all probability laws $P\in \mathcal{P}$ for the triple $(Y,D,X)$ the following conditions hold: (a) equations 5.1 and 5.2 hold, with $D\in \{0,1\}$⁠; (b) ${\parallel Y\parallel }_{P,q}\le C$⁠; (c) ${Pr}_P(\epsilon \le m_0(X)\le 1-\epsilon )=1$⁠; (d) ${\parallel U\parallel }_{P,2}\ge c$⁠; (e) $\parallel E_P[U^2\mid X]{\parallel }_{P,\infty }\le C$⁠; and (f) given a random subset I of [N] of size $n=N/K$⁠, the nuisance parameter estimator ${\widehat{\eta }}_0={\widehat{\eta }}_0({(W_i)}_{i\in I^c})$ obeys the following conditions. With P‐probability no less than $1-{\mathrm{\Delta }}_N$⁠, $\parallel {\widehat{\eta }}_0-{\eta }_0{\parallel }_{P,q}\le C$⁠, $\parallel {\widehat{\eta }}_0-{\eta }_0{\parallel }_{P,2}\le {\delta }_N,{\parallel {\widehat{m}}_0-1/2\parallel }_{P,\infty }\le 1/2-\epsilon$ and (i) for the score ψ in 5.3, where ${\widehat{\eta }}_0=({\widehat{g}}_0,{\widehat{m}}_0)$ and the target parameter is ATE, $\parallel {\widehat{m}}_0-m_0{\parallel }_{P,2}\times {\parallel {\widehat{g}}_0-g_0\parallel }_{P,2}\le {\delta }_N{N}^{-1/2}$⁠, and (ii) for the score ψ in 5.4, where ${\widehat{\eta }}_0=({\widehat{\overline{g}}}_0,{\widehat{m}}_0,{\widehat{p}}_0)$ and the target parameter is ATTE, $\parallel {\widehat{m}}_0-m_0{\parallel }_{P,2}\times {\parallel {\widehat{\overline{g}}}_0-{\overline{g}}_0\parallel }_{P,2}\le {\delta }_N{N}^{-1/2}$⁠.

 

Remark 5.1.

The only non‐primitive condition here is the assumption on the rate of estimating the nuisance parameters. These rates of convergence are available for the most often used ML methods and are case‐specific, so we do not restate conditions that are needed to reach these rates. The conditions are not the tightest possible, but offer a set of simple conditions under which Theorem 5.1 follows as a special case of the general theorem provided in Section 3. We could obtain more refined conditions by doing customized proofs.

The following theorem follows as a corollary to the results in Section 3 by verifying Assumptions 3.1 and 3.2.

 

Theorem 5.1. Dml Inference on Ate and Atte

Suppose that the target parameter is either ATE, ${\theta }_0=E_P[g_0(1,X)-g_0(0,X)]$⁠, and the score ψ in 5.3 is used, or ATTE, ${\theta }_0=E_P[g_0(1,X)-g_0(0,X)\mid D=1]$⁠, and the score ψ in 5.4 is used. In addition, suppose that Assumption 5.1 holds. Then the DML1 and DML2 estimators ${\stackrel{\sim }{\theta }}_0$⁠, constructed in Definitions 3.1 and 3.2, are first‐order equivalent and obey  

$${\sigma }^{-1}\sqrt{N}({\stackrel{\sim }{\theta }}_0-{\theta }_0)\leftrightsquigarrow N(0,1),$$

(5.5)

uniformly over $P\in \mathcal{P}$⁠, where ${\sigma }^2=E_P[{\psi }^2(W;{\theta }_0,{\eta }_0)]$⁠. Moreover, the result continues to hold if σ² is replaced by ${\widehat{\sigma }}^2$ defined in Theorem 3.2. Consequently, confidence regions based upon the DML estimators ${\stackrel{\sim }{\theta }}_0$ have uniform asymptotic validity:  

$$\underset{N\to \infty }{lim}\underset{P\in \mathcal{P}}{sup}|{Pr}_P({\theta }_0\in [{\stackrel{\sim }{\theta }}_0\pm {\mathrm{\Phi }}^{-1}(1-\alpha /2)\widehat{\sigma }/\sqrt{N}])-(1-\alpha )|=0.$$

The scores ψ in 5.3 and 5.4 are efficient, so both estimators are asymptotically efficient, reaching the semi‐parametric efficiency bound of Hahn (1998).

 

Remark 5.2. Tightness of Conditions

The conditions in Assumption 5.1 are sharp though simplified for ease of presentation. The sharpness can be understood by examining the case where the regression function g₀ and the propensity function m₀ are sparse with sparsity indices $s^g\ll N$ and $s^m\ll N$⁠. They are estimated by ℓ₁‐penalized estimators ${\widehat{g}}_0$ and ${\widehat{m}}_0$ that have sparsity indices of orders $s^g$ and $s^m$ and converge to g₀ and m₀ at the rates $\sqrt{s^g/N}$ and $\sqrt{s^m/N}$ (ignoring logs). Then the rate conditions in Assumption 5.1 require  

$$\sqrt{s^g/N}\sqrt{s^m/N}\ll N^{-1/2}\iff s^g{s}^m\ll N$$

(ignoring logs), which is much weaker than the condition ${(s^g)}^2+{(s^m)}^2\ll N$ (ignoring logs) required without sample splitting. For example, if the propensity score m₀ is very sparse, then the regression function is allowed to be quite dense with $s^g>\sqrt{N}$⁠, and vice versa. If the propensity score is known (⁠$s^m=0$⁠), then only consistency for ${\widehat{g}}_0$ is needed. Such comparisons also extend to approximately sparse models. We note that similar refined rates appeared in Farrell (2015), who considers estimation of treatment effects in a setting where an approximately sparse model holds for both the regression and propensity score functions. In interesting related work, Athey et al. (2016) show that $\sqrt{N}$‐consistent estimation of an ATE is possible under very weak conditions on the propensity score – allowing for the possibility that the propensity score may not be consistently estimated – under strong sparsity of the regression function such that $s_g\ll \sqrt{N}$⁠. Thus, the approach taken in this context and by Athey et al. (2016) are complementary, and either might be preferred, depending on whether or not the regression function can be estimated extremely well based on a sparse method.

5.2. Inference on Local Average Treatment Effects

In this section, we consider estimation of LATEs with a binary treatment variable, $D\in \{0,1\}$⁠, and a binary instrument, $Z\in \{0,1\}$⁠.¹⁰ As before, Y denotes the outcome variable, and X is the vector of covariates.

Consider the functions μ₀, m₀ and p₀, where μ₀ maps the support of $(Z,X)$ to $\mathbb{R}$⁠, and m₀ and p₀, respectively, map the support of $(Z,X)$ and X to $(\epsilon ,1-\epsilon )$ for some $\epsilon \in (0,1/2)$⁠, such that  

$$\begin{array}{ccc}\hfill Y& =& {\mu }_0(Z,X)+U,E_P[U\mid Z,X]=0,\hfill \end{array}$$

(5.6)

 

$$\begin{array}{ccc}\hfill D& =& m_0(Z,X)+V,E_P[V\mid Z,X]=0,\hfill \end{array}$$

(5.7)

 

$$\begin{array}{ccc}\hfill Z& =& p_0(X)+\zeta ,E_P[\zeta \mid X]=0.\hfill \end{array}$$

(5.8)

We are interested in estimating  

$${\theta }_0=\frac{E_P[\mu (1,X)]-E_P[\mu (0,X)]}{E_P[m(1,X)]-E_P[m(0,X)]}.$$

Under the assumptions of Imbens and Angrist (1994) and Frölich (2007), θ₀ is the LATE – the average treatment effect for compliers that are observations that would have $D=1$ if Z were 1 and $D=0$ if Z were 0. To estimate θ₀, we use the score  

$$\begin{array}{ccc}\hfill \psi (W;\theta ,\eta )& :=& \mu (1,X)-\mu (0,X)+\frac{Z(Y-\mu (1,X))}{p(X)}-\frac{(1-Z)(Y-\mu (1,X))}{1-p(X))}\hfill \\ & & -\left(m(1,X)-m(0,X)+\frac{Z(D-m(1,X))}{p(X)}-\frac{(1-Z)(D-m(0,X))}{1-p(X)}\right)\times \theta ,\hfill \end{array}$$

where $W=(Y,D,X,Z)$ and the nuisance parameter $\eta =(\mu ,m,p)$ consists of P‐square‐integrable functions μ, m and p, with μ mapping the support of $(Z,X)$ to $\mathbb{R}$ and m and p, respectively, mapping the support of $(Z,X)$ and X to $(\epsilon ,1-\epsilon )$ for some $\epsilon \in (0,1/2)$⁠. It is easy to verify that this score satisfies the moment condition $E_P[\psi (W;{\theta }_0,{\eta }_0)]=0$ and also the orthogonality condition ${\partial }_{\eta }{E}_P[\psi (W;{\theta }_0,{\eta }_0)][\eta -{\eta }_0]=0$ for ${\eta }_0=({\mu }_0,m_0,p_0)$⁠.

Let ${({\delta }_N)}_{n=1}^{\infty }$ and ${({\mathrm{\Delta }}_N)}_{n=1}^{\infty }$ be sequences of positive constants approaching 0. Also, let c, C and q be fixed strictly positive constants such that $q>4$⁠, and let $K\ge 2$ be a fixed integer. Moreover, for any $\eta =({\ell }_1,{\ell }_2,{\ell }_3)$⁠, where ℓ₁ is a function mapping the support of $(Z,X)$ to $\mathbb{R}$ and ℓ₂ and ℓ₃ are functions respectively mapping the support of $(Z,X)$ and X to $(\epsilon ,1-\epsilon )$ for some $\epsilon \in (0,1/2)$⁠, we denote ${\parallel \eta \parallel }_{P,q}=\parallel {\ell }_1{\parallel }_{P,q}\vee \parallel {\ell }_2{\parallel }_{P,2}\vee {\parallel {\ell }_3\parallel }_{P,q}$⁠. For simplicity, assume that $N/K$ is an integer.  

Assumption 5.2. Regularity Conditions for Late Estimation

For all probability laws $P\in \mathcal{P}$ for the quadruple $W=(Y,D,X,Z)$ the following conditions hold: (a) equations 5.6–5.8 hold, with $D\in \{0,1\}$ and $Z\in \{0,1\}$⁠; (b) ${\parallel Y\parallel }_{P,q}\le C$⁠; (c) ${Pr}_P(\epsilon \le p_0(X)\le 1-\epsilon )=1$⁠; (d) $E_P[m_0(1,X)-m_0(0,X)]\ge c$⁠, (e) $\parallel U-{\theta }_0{V\parallel }_{P,2}\ge c$⁠; (f) $\parallel E_P[U^2\mid X]{\parallel }_{P,\infty }\le C$⁠; and (g) given a random subset I of [N] of size $n=N/K$⁠, the nuisance parameter estimator ${\widehat{\eta }}_0={\widehat{\eta }}_0({(W_i)}_{i\in I^c})$ obeys the following conditions. With P‐probability no less than $1-{\mathrm{\Delta }}_N$⁠, $\parallel {\widehat{\eta }}_0-{\eta }_0{\parallel }_{P,q}\le C$⁠, $\parallel {\widehat{\eta }}_0-{\eta }_0{\parallel }_{P,2}\le {\delta }_N$⁠, $\parallel {\widehat{p}}_0{-1/2\parallel }_{P,\infty }\le 1/2-\epsilon$ and $\parallel {\widehat{p}}_0-p_0{\parallel }_{P,2}\times (\parallel {\widehat{\mu }}_0-{\mu }_0{\parallel }_{P,2}+\parallel {\widehat{m}}_0-m_0{\parallel }_{P,2})\le {\delta }_N{N}^{-1/2}$⁠.

The following theorem follows as a corollary to the results in Section 3 by verifying Assumptions 3.1 and 3.2.

 

Theorem 5.2. Dml Inference on Late

Suppose that Assumption 5.2 holds. Then the DML1 and DML2 estimators ${\stackrel{\sim }{\theta }}_0$ constructed in Definitions 3.1 and 3.2 and based on the score ψ above are first‐order equivalent and obey  

$${\sigma }^{-1}\sqrt{N}({\stackrel{\sim }{\theta }}_0-{\theta }_0)\leftrightsquigarrow N(0,1),$$

(5.9)

uniformly over $P\in \mathcal{P}$⁠, where ${\sigma }^2={(E_P[m(1,X)-m(0,X)])}^{-2}{E}_P[{\psi }^2(W;{\theta }_0,{\eta }_0)]$⁠. Moreover, the result continues to hold if σ² is replaced by ${\widehat{\sigma }}^2$ defined in Theorem 3.2. Consequently, confidence regions based upon the DML estimators ${\stackrel{\sim }{\theta }}_0$ have uniform asymptotic validity:  

$$\underset{N\to \infty }{lim}\underset{P\in \mathcal{P}}{sup}|{Pr}_P({\theta }_0\in [{\stackrel{\sim }{\theta }}_0\pm {\mathrm{\Phi }}^{-1}(1-\alpha /2)\widehat{\sigma }/\sqrt{N}])-(1-\alpha )|=0.$$

6. Empirical Examples

To illustrate the methods developed in the preceding sections, we consider three empirical examples. The first example reexamines the Pennsylvania Reemployment Bonus experiment, which used a randomized control trial to investigate the incentive effect of unemployment insurance. In the second, we use the DML method to estimate the effect of 401(k) eligibility, the treatment variable, and 401(k) participation, a self‐selected decision to receive the treatment that we instrument for with assignment to the treatment state, on accumulated assets. In this example, the treatment variable is not randomly assigned and we aim to eliminate the potential biases due to the lack of random assignment by flexibly controlling for a rich set of variables. In the third, we revisit the IV estimation by Acemoglu et al. (2001) of the effects of institutions on economic growth by estimating a partially linear IV model.

6.1. Effect of Unemployment Insurance Bonus on Unemployment Duration

In this example, we reanalyse the Pennsylvania Reemployment Bonus experiment, which was conducted by the US Department of Labor in the 1980s to test the incentive effects of alternative compensation schemes for unemployment insurance (UI). This experiment has been previously studied by Bilias (2000) and Bilias and Koenker (2002). In these experiments, UI claimants were randomly assigned either to a control group or to one of five treatment groups.¹¹ In the control group, the standard rules of the UI system applied. Individuals in the treatment groups were offered a cash bonus if they found a job within some pre‐specified period of time (qualification period), provided that the job was retained for a specified duration. The treatments differed in the level of the bonus, the length of the qualification period, and whether the bonus was declining over time in the qualification period; see Bilias and Koenker (2002) for further details.

In our empirical example, we focus only on the most generous compensation scheme, treatment 4, and we drop all individuals who received other treatments. In this treatment, the bonus amount is high and the qualification period is long compared to other treatments, and claimants are eligible to enroll in a workshop. Our treatment variable, D, is an indicator variable for being assigned treatment 4, and the outcome variable, Y, is the log of duration of unemployment for the UI claimants. The vector of covariates, X, consists of age group dummies, gender, race, the number of dependents, quarter of the experiment, location within the state, existence of recall expectations and type of occupation.

We report results based on five simple methods for estimating the nuisance functions used in forming the orthogonal estimating equations. We consider three tree‐based methods, labelled ‘Random forest’, ‘Reg. tree’ and ‘Boosting’, one ℓ₁‐penalization based method, labelled ‘Lasso’ and a neural network method, labelled ‘Neural network’. For Reg. tree, we fit a single CART tree to estimate each nuisance function with penalty parameter chosen by tenfold cross‐validation. The results in the Random forest column are obtained by estimating each nuisance function with a random forest that averages over 1000 trees. The results in Boosting are obtained using boosted regression trees with regularization parameters chosen by tenfold cross‐validation. To estimate the nuisance functions using the neural networks, we use two neurons and a decay parameter of 0.02, and we set activation function as logistic for classification problems and as linear for regression problems.¹² Lasso estimates an ℓ₁‐penalized linear regression model using the data‐driven penalty parameter selection rule developed in Belloni et al. (2012). For Lasso, we use a set of 96 potential control variables formed from the raw set of covariates and all second‐order terms (i.e. all squares and first‐order interactions). For the remaining methods, we use the raw set of covariates as features.

We also consider two hybrid methods labelled ‘Ensemble’ and ‘Best’. Ensemble optimally combines four of the ML methods listed above by estimating the nuisance functions as weighted averages of estimates from Lasso, Boosting, Random forest and Neural network. The weights are restricted to sum to one and are chosen so that the weighted average of these methods gives the lowest average mean squared out‐of‐sample prediction error estimated using fivefold cross‐validation. The final column in Table 1 (Best) reports results that combine the methods in a different way. After obtaining estimates from the five simple methods and Ensemble, we select the best methods for estimating each nuisance function based on the average out‐of‐sample prediction performance for the target variable associated with each nuisance function obtained from each of the previously described approaches. As a result, the reported estimate in the last column uses different ML methods to estimate different nuisance functions. Note that if a single method outperformed all the others in terms of prediction accuracy for all nuisance functions, the estimate in the Best column would be identical to the estimate reported under that method.

Table 1.

Estimated effect of cash bonus on unemployment duration

	Lasso	Reg. tree	Random forest	Boosting	Neural network	Ensemble	Best
Panel A: interactive regression model
ATE	−0.081	−0.084	−0.074	−0.079	−0.073	−0.079	−0.078
(twofold)	[0.036]	[0.036]	[0.036]	[0.036]	[0.036]	[0.036]	[0.036]
	(0.036)	(0.036)	(0.036)	(0.036)	(0.036)	(0.036)	(0.036)
ATE	−0.081	−0.085	−0.074	−0.077	−0.073	−0.078	−0.077
(fivefold)	[0.036]	[0.036]	[0.036]	[0.035]	[0.036]	[0.036]	[0.036]
	(0.036)	(0.037)	(0.036)	(0.036)	(0.036)	(0.036)	(0.036)
Panel B: partially linear regression model
ATE	−0.080	−0.084	−0.077	−0.076	−0.074	−0.075	−0.075
(twofold)	[0.036]	[0.036]	[0.035]	[0.035]	[0.035]	[0.035]	[0.035]
	(0.036)	(0.036)	(0.037)	(0.036)	(0.036)	(0.036)	(0.036)
ATE	−0.080	−0.084	−0.077	−0.074	−0.073	−0.075	−0.074
(fivefold)	[0.036]	[0.036]	[0.035]	[0.035]	[0.035]	[0.035]	[0.035]
	(0.036)	(0.037)	(0.036)	(0.035)	(0.036)	(0.035)	(0.035)

Note: Estimated ATE and standard errors from a linear model (Panel B) and heterogeneous effect model (Panel A) based on orthogonal estimating equations. Column labels denote the method used to estimate nuisance functions. Results are based on 100 splits with point estimates calculated the median method. The median standard errors across the splits are reported in brackets and standard errors calculated using the median method to adjust for variation across splits are provided in parentheses. Further details about the methods are provided in the main text.

Table 1 presents DML2 estimates of the ATE on unemployment duration using the median method described in Section 3.4. We report results for the heterogeneous effect model in Panel A and for the partially linear model in Panel B. Because the treatment is randomly assigned, we use the fraction of treated as the estimator of the propensity score in forming the orthogonal estimating equations.¹³ For both the partially linear model and the interactive model, we report estimates obtained using twofold cross‐fitting and fivefold cross‐fitting. All results are based on taking 100 different sample splits. We summarize results across the sample splits using the median method. For comparison, we report two different standard errors: in brackets, we report the median standard errors from across the 100 splits; in parentheses, we report standard errors adjusted for variability across the sample splits using the median method in parentheses.

The estimation results are consistent with the findings of previous studies that have analysed the Pennsylvania Bonus Experiment. The ATE on unemployment duration is negative and significant across all estimation methods at the 5% level regardless of the standard error estimator used. Interestingly, we see that there is no practical difference across the two different standard errors in this example.

6.2. Effect of 401(K) Eligibility and Participation on Net Financial Assets

The key problem in determining the effect of 401(k) eligibility is that working for a firm that offers access to a 401(k) plan is not randomly assigned. To overcome the lack of random assignment, we follow the strategy developed in Poterba et al. (1994a,b). In these papers, the authors use data from the 1991 Survey of Income and Program Participation and they argue that eligibility for enrolling in a 401(k) plan in these data can be taken as exogenous after conditioning on a few observables – of which the most important for their argument is income. The basic idea of their argument is that, at least around the time 401(k) initially became available, people were unlikely to be basing their employment decisions on whether an employer offered a 401(k) but would instead focus on income and other aspects of the job. Following this argument, whether one is eligible for a 401(k) may then be taken as exogenous after appropriately conditioning on income and other control variables related to job choice.

A key component of the argument underlying the exogeneity of 401(k) eligibility is that eligibility can only be taken as exogenous after conditioning on income and other variables related to job choice that might correlate with whether a firm offers a 401(k). Poterba et al. (1994a,b) and many subsequent papers adopt this argument but control only linearly for a small number of terms. One might wonder whether such specifications are able to adequately control for income and other related confounds. At the same time, the power to learn about treatment effects decreases as one allows more flexible models. The principled use of flexible ML tools offers one resolution to this tension. The results presented below thus complement previous results that rely on the assumption that confounding effects can adequately be controlled for by a small number of variables chosen ex ante by the researcher.

In the example in this paper, we use the same data as in Chernozhukov and Hansen (2004). We use net financial assets – defined as the sum of IRA balances, 401(k) balances, checking accounts, US saving bonds, other interest‐earning accounts in banks and other financial institutions, other interest‐earning assets (such as bonds held personally), stocks, and mutual funds less non‐mortgage debt – as the outcome variable, Y, in our analysis. Our treatment variable, D, is an indicator for being eligible to enroll in a 401(k) plan. The vector of raw covariates, X, consists of age, income, family size, years of education, a married indicator, a two‐earner status indicator, a defined benefit pension status indicator, an IRA participation indicator, and a home‐ownership indicator.

In Table 2, we report DML2 estimates of ATE of 401(k) eligibility on net financial assets both in the partially linear model as in 1.1 and allowing for heterogeneous treatment effects using the interactive model outlined in Section 5.1. To reduce the disproportionate impact of extreme propensity score weights in the interactive model, we trim the propensity scores at 0.01 and 0.99. We present two sets of results based on sample splitting as discussed in Section 3 using twofold cross‐fitting and fivefold cross‐fitting. As in the previous section, we consider 100 different sample partitions and summarize the results across different sample splits using the median method. For comparison, we report two different standard errors: in brackets, we report the median standard errors from across the 100 splits; in parentheses, we report standard errors adjusted for variability across the sample splits using the median method. We consider the same methods with the same tuning choices for estimating the nuisance functions as in the previous example, with one exception, and so we do not repeat details for brevity. The one exception is that we implement neural networks with eight neurons and a decay parameter of 0.01 in this example.

Table 2.

Estimated effect of 401(k) eligibility on net financial assets

	Lasso	Reg. tree	Random forest	Boosting	Neural network	Ensemble	Best
Panel A: interactive regression model
ATE	6830	7713	7770	7806	7764	7702	7546
(twofold)	[1282]	[1208]	[1276]	[1159]	[1328]	[1149]	[1360]
	(1530)	(1271)	(1363)	(1202)	(1468)	(1170)	(1533)
ATE	7170	7993	8105	7713	7788	7839	7753
(fivefold)	[1201]	[1198]	[1242]	[1155]	[1238]	[1134]	[1237]
	(1398)	(1236)	(1299)	(1177)	(1293)	(1148)	(1294)
Panel B: partially linear regression model
ATE	7717	8709	9116	8759	8950	9010	9125
(twofold)	[1346]	[1363]	[1302]	[1339]	[1335]	[1309]	[1304]
	(1749)	(1427)	(1377)	(1382)	(1408)	(1344)	(1357)
ATE	8187	8871	9247	9110	9038	9166	9215
(fivefold)	[1298]	[1358]	[1295]	[1314]	[1322]	[1299]	[1294]
	(1558)	(1418)	(1328)	(1328)	(1355)	(1310)	(1312)

Note: Estimated ATE and standard errors from a linear model (Panel B) and heterogeneous effect model (Panel A) based on orthogonal estimating equations. Column labels denote the method used to estimate nuisance functions. Results are based on 100 splits with point estimates calculated the median method. The median standard errors across the splits are reported in brackets and standard errors calculated using the median method to adjust for variation across splits are provided in parentheses. Further details about the methods are provided in the main text.

Turning to the results, it is first worth noting that the estimated ATE of 401(k) eligibility on net financial assets is $19,559 with an estimated standard error of 1413 when no control variables are used. Of course, this number is not a valid estimate of the causal effect of 401(k) eligibility on financial assets if there are neglected confounding variables as suggested by Poterba et al. (1994a,b). When we turn to the estimates that flexibly account for confounding reported in Table 2, we see that they are substantially attenuated relative to this baseline that does not account for confounding, suggesting much smaller causal effects of 401(k) eligibility on financial asset holdings. It is interesting and reassuring that the results obtained from the different flexible methods are broadly consistent with each other. This similarity is consistent with the theory that suggests that results obtained through the use of orthogonal estimating equations and any sensible method of estimating the necessary nuisance functions should be similar. Finally, it is interesting that these results are also broadly consistent with those reported in the original work of Poterba et al. (1994a,b), who used a simple intuitively motivated functional form, suggesting that this intuitive choice was sufficiently flexible to capture much of the confounding variation in this example.

As a further illustration, we also report the LATE in this example where we take the endogenous treatment variable to be participating in a 401(k) plan. Even after controlling for features related to job choice, it seems likely that the actual choice of whether to participate in an offered plan would be endogenous. Of course, we can use eligibility for a 401(k) plan as an instrument for participation in a 401(k) plan under the conditions that were used to justify the exogeneity of eligibility for a 401(k) plan provided above in the discussion of estimation of the ATE of 401(k) eligibility.

We report DML2 results of estimating the LATE of 401(k) participation using 401(k) eligibility as an instrument in Table 3. We employ the procedure outlined in Section 5.2 using the same ML estimators to estimate the quantities used to form the orthogonal estimating equation as we employed to estimate the ATE of 401(k) eligibility outlined previously, so we omit the details for brevity. Looking at the results, we see that the estimated causal effect of 401(k) participation on net financial assets is uniformly positive and statistically significant across all of the considered methods. As when looking at the ATE of 401(k) eligibility, it is reassuring that the results obtained from the different flexible methods are broadly consistent with each other. It is also interesting that the results based on flexible ML methods are broadly consistent with, though somewhat attenuated relative to, those obtained by applying the same specification for controls as used in Poterba et al. (1994a,b) and using a linear IV model, which returns an estimated effect of participation of $13,102 with estimated standard error of (1922). The mild attenuation may suggest that the simple intuitive control specification used in the original baseline specification is too simplistic.

Table 3.

Estimated effect of 401(k) participation on net financial assets

	Lasso	Reg. tree	Random forest	Boosting	Neural network	Ensemble	Best
LATE	8978	11073	11384	11329	11094	11119	10952
(twofold)	[2192]	[1749]	[1832]	[1666]	[1903]	[1653]	[1657]
	(3014)	(1849)	(1993)	(1718)	(2098)	(1689)	(1699)
LATE	8944	11459	11764	11133	11186	11173	11113
(fivefold)	[2259]	[1717]	[1788]	[1661]	[1795]	[1641]	[1645]
	(3307)	(1786)	(1893)	(1710)	(1890)	(1678)	(1675)

Note: Estimated LATE based on orthogonal estimating equations. Column labels denote the method used to estimate nuisance functions. Results are based on 100 splits with point estimates calculated the median method. The median standard errors across the splits are reported in brackets and standard errors calculated using the median method to adjust for variation across splits are provided in parentheses. Further details about the methods are provided in the main text.

Looking at Tables 2 and 3, there are other interesting observations that can provide useful insights into understanding the finite sample properties of the DML estimation method. First, the standard errors of the estimates obtained using fivefold cross‐fitting are lower than those obtained from twofold cross‐fitting for all methods across all cases. This fact suggests that having more observations in the auxiliary sample may be desirable. Specifically, the fivefold cross‐fitting estimates use more observations to learn the nuisance functions than twofold cross‐fitting and thus are likely learn them more precisely. This increase in precision in learning the nuisance functions may then translate into more precisely estimated parameters of interest. While intuitive, we note that this statement does not seem to be generalizable, in that there does not appear to be a general relationship between the number of folds in cross‐fitting and the precision of the estimate of the parameter of interest (see the next example). Secondly, we also see that the standard errors of the Lasso estimates after adjusting for variation due to sample splitting are noticeably larger than the standard errors coming from the other ML methods. We believe that this is due to the fact that the out‐of‐sample prediction errors from a linear model tend to be larger when there is a need to extrapolate. In our framework, if the main sample includes observations that are outside of the range of the observations in the auxiliary sample, the model has to extrapolate to those observations. The fact that the standard errors are lower in fivefold cross‐fitting than in twofold cross‐fitting for the Lasso estimations also supports this hypothesis, because the higher number of observations in the auxiliary sample reduces the degree of extrapolation. We also see that there is a noticeable increase in the standard errors that account for variability due to sample splitting relative to the simple unadjusted standard errors in this case, though these differences do not qualitatively change the results.

6.3. Effect of Institutions on Economic Growth

To demonstrate DML estimation of partially linear structural equation models with instrumental variables, we consider estimation of the effect of institutions on aggregate output following the work of Acemoglu et al. (2001, hereafter AJR). Estimating the effect of institutions on output is complicated by the clear potential for simultaneity between institutions and output. Specifically, better institutions may lead to higher incomes, but higher incomes may also lead to the development of better institutions. To help overcome this simultaneity, AJR use mortality rates for early European settlers as an instrument for institution quality. The validity of this instrument hinges on the arguments that settlers set up better institutions in places where they are more likely to establish long‐term settlements, that where they are likely to settle for the long term is related to settler mortality at the time of initial colonization, and that institutions are highly persistent. The exclusion restriction for the instrumental variable is then motivated by the argument that GDP, while persistent, is unlikely to be strongly influenced by mortality in the previous century, or earlier, except through institutions.

In their paper, AJR note that their IV strategy will be invalidated if other factors are also highly persistent and related to the development of institutions within a country and to the country's GDP. A leading candidate for such a factor, as they discuss, is geography. AJR address this by assuming that the confounding effect of geography is adequately captured by a linear term in distance from the equator and a set of continent dummy variables. Using DML allows us to relax this assumption and to replace it by a weaker assumption that geography can be sufficiently controlled by an unknown function of distance from the equator and continent dummies, which can be learned by ML methods.

We use the same set of 64 country‐level observations as AJR. The data set contains measurements of GDP, settler morality, an index that measures protection against expropriation risk and geographic information. The outcome variable, Y, is the logarithm of GDP per capita and the endogenous explanatory variable, D, is a measure of the strength of individual property rights that is used as a proxy for the strength of institutions. To deal with endogeneity, we use an instrumental variable Z, which is mortality rates for early European settlers. Our raw set of control variables, X, include distance from the equator and dummy variables for Africa, Asia, North America and South America.

We report results from applying DML2 following the procedure outlined in Section 4.2 in Table 4. The considered ML methods and tuning parameters are the same as the previous examples except for the Ensemble method, from which we exclude Neural network, as the small sample size causes stability problems in training the neural network. We use the raw set of covariates and all second‐order terms when doing Lasso estimation, and we simply use the raw set of covariates in the remaining methods. As in the previous examples, we consider 100 different sample splits and report the median estimates of the coefficient and two different standard error estimates. In brackets, we report the median standard errors from across the 100 splits, and we report standard errors adjusted for variability across the sample splits using the median method in parentheses. Finally, we report results from both twofold cross‐fitting and fivefold cross‐fitting as in the other examples.

Table 4.

Estimated effect of institutions on output

	Lasso	Reg. tree	Random forest	Boosting	Neural network	Ensemble	Best
Twofold	0.85	0.81	0.84	0.77	0.94	0.80	0.83
	[0.28]	[0.42]	[0.38]	[0.33]	[0.32]	[0.35]	[0.34]
	(0.22)	(0.29)	(0.30)	(0.27)	(0.28)	(0.30)	(0.29)
Fivefold	0.77	0.95	0.90	0.73	1.00	0.83	0.88
	[0.24]	[0.46]	[0.41]	[0.33]	[0.33]	[0.37]	[0.41]
	(0.17)	(0.45)	(0.40)	(0.27)	(0.30)	(0.34)	(0.39)

Note: Estimated coefficient from a linear IV model based on orthogonal estimating equations. Column labels denote the method used to estimate nuisance functions. Results are based on 100 splits with point estimates calculated the median method. The median standard errors across the splits are reported in brackets and standard errors calculated using the median method to adjust for variation across splits are provided in parentheses. Further details about the methods are provided in the main text.

In this example, we see uniformly large and positive point estimates across all procedures considered, and estimated effects are statistically significant at the 5% level. As in the second example, we see that adjusting for variability across sample splits leads to noticeable increases in estimated standard errors but does not result in qualitatively different conclusions. Interestingly, we see that the estimated standard errors based on fivefold cross‐fitting are larger than those based on twofold cross‐fitting in all procedures except lasso, which differs from the finding in the 401(k) example. Further understanding these differences and the impact of the number of folds on inference for objects of interest seems to be an interesting question for future research. Finally, although the estimated coefficients are somewhat smaller than the baseline estimates reported in AJR – an estimated coefficient of 1.10 with estimated standard error of 0.46 (see AJR, Table 4, Panel A, column 7) – the results are qualitatively similar, indicating a strong and positive effect of institutions on output.

6.4. Comments on Empirical Results

Before closing this section, we want to emphasize some important conclusions that can be drawn from these empirical examples. First, the choice of the ML method used in estimating nuisance functions does not substantively change the conclusion in any of the examples, and we have obtained broadly consistent results regardless of which method we employ. The robustness of the results to the different methods is implied by the theory assuming that all of the employed methods are able to deliver sufficiently high‐quality approximations to the underlying nuisance functions. Secondly, the incorporation of uncertainty due to sample splitting using the median method increases the standard errors relative to a baseline that does not account for this uncertainty, though these differences do not alter the main results in any of the examples. This lack of variation suggests that the parameter estimates are robust to the particular sample split used in the estimation in these examples.

Acknowledgements

We would like to acknowledge research support from the National Science Foundation. We also thank participants of the MIT Stochastics and Statistics seminar, the Kansas Econometrics conference, the Royal Economic Society Annual Conference, the Hannan Lecture at the Australasian Econometric Society meeting, the Econometric Theory lecture at the EC² meetings 2016 in Toulouse, the CORE 50th Anniversary Conference, the Becker–Friedman Institute Conference on Machine Learning and Economics, the INET conferences on Big Data at the University of Southern California in Los Angeles, the World Congress of Probability and Statistics 2016, the Joint Statistical Meetings 2016, the New England Day of Statistics Conference, CEMMAP's Masterclass on Causal Machine Learning, and St Gallen's summer school on Big Data, for many useful comments and questions. We would like to thank Susan Athey, Peter Aronow, Jin Hahn, Guido Imbens, Mark van der Laan and Matt Taddy for constructive comments. We thank Peter Aronow for pointing us to the literature on targeted learning on which we build, along with prior works of Neyman, Bickel, and the many other contributions to semi‐parametric learning theory.

Footnotes

We consider the case where D is a scalar for simplicity. Extension to the case where D is a vector of fixed, finite dimension is accomplished by introducing an equation such as 1.2 for each element of the vector.

For instance, we could use lasso if we believe g₀ is well approximated by a sparse linear combination of pre‐specified functions of X. In other settings, we could, for example, use iterative methods that alternate between random forests, for estimating g₀, and least squares, for estimating θ₀.

In Section 4, we also consider another debiased estimator, based on the partialling‐out approach of Robinson (1988):  

$${\displaystyle {\check{\theta }}_0={\left(\frac{1}{n}\sum _{i\in I}{\widehat{V}}_i{\widehat{V}}_i\right)}^{-1}\frac{1}{n}\sum _{i\in I}{\widehat{V}}_i(Y_i-{\widehat{\ell }}_0(X_i)),{\ell }_0(X)=E[Y|X].}$$

Each of these works differs in terms of detail but can be viewed through the lens of either debiasing or orthogonalization to alleviate the impact of regularization bias on subsequent estimation and inference.

See Section 2 for the definition of the Gateaux derivative operator.

Targeted minimum loss estimation, which shares similar properties, is also discussed in, e.g., van der Laan and Rose (2011) and van der Laan (2015).

The $C(\alpha )$‐statistic, or the orthogonal score statistic, has been explicitly used for testing and estimation in high‐dimensional sparse models in Belloni et al. (2015).

It is interesting to note that the methods for constructing Neyman orthogonal scores described in Section 2 may give scores that are different from those in 4.7 and 4.8. For example, applying the method for conditional moment restriction problems in Section 2.2.4 with $\mathrm{\Omega }(R)=1$ gives the score $\psi (W;\theta ,\eta )=(Y-D\theta -g(X))(r(Z,X)-f(X))$⁠, where the true values of $r(Z,X)$ and $f(X)$ are $r_0(Z,X)=E_P[D\mid Z,X]$ and $f_0(X)=E_P[D\mid X]$⁠, respectively. It may be interesting to compare properties of the DML estimators ${\stackrel{\sim }{\theta }}_0$ based on this score with those based on 4.7 and 4.8 in future work.

Without unconfoundedness/conditional exogeneity, these quantities measure association, and could be referred to as average predictive effect (APE) and average predictive effect for the exposed (APEX). Inferential results for these objects would follow immediately from Theorem 5.1.

Similar results can be provided for the local average treatment effect on the treated (LATT) by adapting the following arguments to use the orthogonal scores for the LATT; see, e.g. Belloni et al. (2017).

There are six treatment groups in the experiments. Following Bilias (2000). we merge groups 4 and 6.

We also experimented with Deep Learning methods from which we obtained similar results for some tuning parameters. However, we ran into stability and computational issues and we have chosen not to report these results in the empirical section.

We also estimated the effects using non‐parametric estimates of the conditional propensity score obtained from the ML procedures given in the column labels. As expected due to randomization, the results are similar to those provided in Table 1 and are not reported for brevity.

Appendix

Proofs of Results

In this appendix, we use C to denote a strictly positive constant that is independent of n and $P\in {\mathcal{P}}_N$⁠. The value of C may change at each appearance. Also, the notation $a_N\lesssim b_N$ means that $a_N\le C{b}_N$ for all n and some C. The notation $a_N\gtrsim b_N$ means that $b_N\lesssim a_N$⁠. Moreover, the notation $a_N=o(1)$ means that there exists a sequence ${(b_N)}_{n\ge 1}$ of positive numbers such that (a) $|a_N|\le b_N$ for all n, (b) $b_N$ is independent of $P\in {\mathcal{P}}_N$ for all n and (c) $b_N\to 0$ as $n\to \infty$⁠. Finally, the notation $a_N=O_P(b_N)$ means that for all $\epsilon >0$⁠, there exists C such that ${Pr}_P(a_N>C{b}_N)\le 1-\epsilon$ for all n. Using this notation allows us to avoid repeating ‘uniformly over $P\in {\mathcal{P}}_N$’ many times in the proofs.

Define the empirical process ${\mathbb{G}}_n(\psi (W))$ as a linear operator acting on measurable functions $\psi :\mathcal{W}\to \mathbb{R}$ such that ${\parallel \psi \parallel }_{P,2}<\infty$ via  

$${\displaystyle {\mathbb{G}}_n(\psi (W)):={\mathbb{G}}_{n,I}(\psi (W)):=\frac{1}{\sqrt{n}}\sum _{i\in I}\psi (W_i)-\int \psi (w)\mathit{dP}(w).}$$

Analogously, we defined the empirical expectation as  

$${\displaystyle E_n[\psi (W)]:=E_{n,I}[\psi (W)]:=\frac{1}{n}\sum _{i\in I}\psi (W_i).}$$

The following lemma is useful particularly in the sample‐splitting contexts.

 

Lemma 6.1.Conditional Convergence Implies Unconditional

Let $\{X_m\}$ and $\{Y_m\}$ be sequences of random vectors. (a) If, for ${\epsilon }_m\to 0$⁠, $\text{Pr}(\parallel X_m\parallel >{\epsilon }_m\mid Y_m){\to }_{\text{Pr}}0$⁠, then $\text{Pr}(\parallel X_m\parallel >{\epsilon }_m)\to 0$⁠. In particular, this occurs if $E[\parallel X_m{\parallel }^q/{\epsilon }_m^q\mid Y_m]{\to }_{\text{Pr}}0$ for some $q\ge 1$⁠, by Markov's inequality. (b) Let $\{A_m\}$ be a sequence of positive constants. If $\parallel X_m\parallel =O_P(A_m)$ conditional on $Y_m$⁠, namely, that for any ${\ell }_m\to \infty$⁠, $Pr(\parallel X_m\parallel >{\ell }_m{A}_m\mid Y_m){\to }_{\text{Pr}}0$⁠, then $\parallel X_m\parallel =O_P(A_m)$ unconditionally, namely, that for any ${\ell }_m\to \infty$⁠, $\text{Pr}(\parallel X_m\parallel >{\ell }_m{A}_m)\to 0$⁠.

 

Proof:

(a) For any $\epsilon >0$$Pr(\parallel X_m\parallel >{\epsilon }_m)\le E[Pr(\parallel X_m\parallel >{\epsilon }_m\mid Y_m)]\to 0$⁠, as the sequence $\{Pr(\parallel X_m\parallel >{\epsilon }_m\mid Y_m)\}$ is uniformly integrable. To show the second part note that $Pr(\parallel X_m\parallel >{\epsilon }_m\mid Y_m)\le E[\parallel X_m{\parallel }^q/{\epsilon }_m^q\mid Y_m]\vee 1\stackrel{p}{\to }0$ by Markov's inequality. (b) This follows from (a). $\square$

Let ${(W_i)}_{i=1}^n$ be a sequence of independent copies of a random element W taking values in a measurable space $(\mathcal{W},{\mathcal{A}}_{\mathcal{W}})$ according to a probability law P. Let $\mathcal{F}$ be a set of suitably measurable functions $f:\mathcal{W}\to \mathbb{R}$⁠, equipped with a measurable envelope $F:\mathcal{W}\to \mathbb{R}$⁠.

 

Lemma 6.2.Maximal Inequality, Chernozhukov EtAl. 2014

Work with the set‐up above. Suppose that $F\ge {sup}_{f\in \mathcal{F}}|f|$ is a measurable envelope for $\mathcal{F}$ with ${\parallel F\parallel }_{P,q}<\infty$ for some $q\ge 2$⁠. Let $M={max}_{i\le n}F(W_i)$ and ${\sigma }^2>0$ be any positive constant such that ${sup}_{f\in \mathcal{F}}{\parallel f\parallel }_{P,2}^2\le {\sigma }^2\le {\parallel F\parallel }_{P,2}^2$⁠. Suppose that there exist constants $a\ge e$ and $v\ge 1$ such that  

$$log\underset{Q}{sup}{N(\epsilon \parallel F\parallel }_{Q,2}{,\mathcal{F},\parallel ·\parallel }_{Q,2})\le vlog(a/\epsilon ),0<\epsilon \le 1.$$

Then  

$$E_P[\parallel {\mathbb{G}}_n{\parallel }_{\mathcal{F}}]\le K\left(\sqrt{v{\sigma }^{2}log\left(\frac{{a\parallel F\parallel }_{P,2}}{\sigma }\right)}+\frac{{v\parallel M\parallel }_{P,2}}{\sqrt{n}}log\left(\frac{{a\parallel F\parallel }_{P,2}}{\sigma }\right)\right),$$

where K is an absolute constant. Moreover, for every $t\ge 1$⁠, with probability $>1-t^{-q/2}$⁠, $\parallel {\mathbb{G}}_n{\parallel }_{\mathcal{F}}\le (1+\alpha )E_P[\parallel {\mathbb{G}}_n{\parallel }_{\mathcal{F}}]+K(q)((\sigma +n^{-1/2}{\parallel M\parallel }_{P,q})\sqrt{t}+{\alpha }^{-1}{n}^{-1/2}{\parallel M\parallel }_{P,2}t),\forall \alpha >0,$ where $K(q)>0$ is a constant depending only on q. In particular, setting $a\ge n$ and $t=logn$⁠, with probability $>1-c{(logn)}^{-1}$⁠,  

$$\parallel {\mathbb{G}}_n{\parallel }_{\mathcal{F}}\le K(q,c)\left(\sigma \sqrt{vlog\left(\frac{{a\parallel F\parallel }_{P,2}}{\sigma }\right)}+\frac{{v\parallel M\parallel }_{P,q}}{\sqrt{n}}log\left(\frac{{a\parallel F\parallel }_{P,2}}{\sigma }\right)\right),$$

A.1

where ${\parallel M\parallel }_{P,q}\le n^{1/q}{\parallel F\parallel }_{P,q}$ and $K(q,c)>0$ is a constant depending only on q and c.

 

Proof of Lemma 2.1:

Because J exists and $J_{\beta \beta }$ is invertible, 2.8 has the unique solution μ₀ given in 2.10, and so we have by 2.6 that $E[\psi (W;{\theta }_0,{\eta }_0)]=0$ for η₀ given in 2.9. Moreover,  

$${\partial }_{{\eta }^{\prime }}{E}_P\psi (W;{\theta }_0,{\eta }_0)=\left([J_{\theta \beta }-{\mu }_0{J}_{\beta \beta }],E[{\partial }_{{\beta }^{\prime }}\ell (W;{\theta }_0,{\beta }_0)]\otimes I_{d_{\theta }\times d_{\theta }}\right)=0,$$

where $I_{d_{\theta }\times d_{\theta }}$ is the $d_{\theta }\times d_{\theta }$ identity matrix and ⊗ is the Kronecker product. Hence, the asserted claim holds by the remark after Definition 2.1. $\square$

 

Proof of Lemma 2.2:

The proof is similar to that of Lemma 2.1, except that now we have to verify 2.4 intead of 2.3. To do so, take any $\beta \in \mathcal{B}$ such that $\parallel \beta -{\beta }_0{\parallel }_q^{\ast }\le {\lambda }_N/r_N$ and any $d_{\theta }\times d_{\beta }$ matrix μ. Denote $\eta ={({\beta }^{\prime },\mathrm{vec}{(\mu )}^{\prime })}^{\prime }$⁠. Then  

$$\begin{array}{ccc}\hfill \parallel {\partial }_{\eta }{E}_P\psi (W,{\theta }_0,{\eta }_0)[\eta -{\eta }_0]\parallel & =& \parallel (J_{\theta \beta }-{\mu }_0{J}_{\beta \beta })(\beta -{\beta }_0)\parallel \hfill \\ & \le & \parallel J_{\theta \beta }-{\mu }_0{J}_{\beta \beta }{\parallel }_q\times {\parallel \beta -{\beta }_0\parallel }_q^{\ast }\le r_n\times ({\lambda }_N/r_N)={\lambda }_N.\hfill \end{array}$$

This completes the proof of the lemma. $\square$

 

Proof of Lemma 2.3:

The proof is similar to that of Lemma 2.1, except that now we have  

$${\partial }_{{\eta }^{\prime }}{E}_P\psi (W,{\theta }_0,{\eta }_0)=[{\mu }_0{G}_{\beta },E_{P}m{(W,{\theta }_0,{\beta }_0)}^{\prime }\otimes I_{d_{\theta }\times d_{\theta }}]=0,$$

where $I_{d_{\theta }\times d_{\theta }}$ is the $d_{\theta }\times d_{\theta }$ identity matrix and ⊗ is the Kronecker product. $\square$

 

Proof of Lemma 2.4:

The proof is similar to that of Lemma 2.2, except that now for any $\beta \in \mathcal{B}$ such that $\parallel \beta -{\beta }_0{\parallel }_1\le {\lambda }_N/r_N$⁠, any $d_{\theta }\times k$ matrix μ, and $\eta ={({\beta }^{\prime },\mathrm{vec}{(\mu )}^{\prime })}^{\prime }$⁠, we have  

$$\begin{array}{ccc}\hfill \parallel {\partial }_{\eta }{E}_P\psi (W,{\theta }_0,{\eta }_0)[\eta -{\eta }_0]\parallel & =& \parallel {\mu }_0{G}_{\beta }(\beta -{\beta }_0)\parallel \hfill \\ & \le & \parallel A^{\prime }{\mathrm{\Omega }}^{-1/2}L-{\gamma }_0{L}^{\prime }{L\parallel }_{\infty }\times {\parallel \beta -{\beta }_0\parallel }_1\hfill \\ & \le & r_n\times ({\lambda }_N/r_N)={\lambda }_N.\hfill \end{array}$$

$\square$

 

Proof of Lemma 2.5:

Take any $\eta \in T$⁠, and consider the function  

$$Q(W;\theta ,r):=\ell (W;\theta ,{\eta }_0(\theta )+r(\eta (\theta )-{\eta }_0(\theta ))),\theta \in \mathrm{\Theta },r\in [0,1].$$

Then  

$$\psi (W;\theta ,{\eta }_0+r(\eta -{\eta }_0))={\partial }_{\theta }Q(W;\theta ,r),$$

and so  

$$\begin{array}{ccc}\hfill {\partial }_r{E}_P[\psi (W;\theta ,{\eta }_0+r(\eta -{\eta }_0))]& =& {\partial }_r{E}_P[{\partial }_{\theta }Q(W;\theta ,r)]\hfill \\ & =& {\partial }_r{\partial }_{\theta }{E}_P[Q(W;\theta ,r)]={\partial }_{\theta }{\partial }_r{E}_P[Q(W;\theta ,r)]\hfill \\ & =& {\partial }_{\theta }{\partial }_r{E}_P[\ell (W;\theta ,{\eta }_0(\theta )+r(\eta (\theta )-{\eta }_0(\theta )))].\hfill \end{array}$$

A.2

Hence,  

$${\partial }_r{E}_P[\psi (W;\theta ,{\eta }_0+r(\eta -{\eta }_0))]{|}_{r=0}=0$$

because  

$${\partial }_r{E}_P[\ell (W;\theta ,{\eta }_0(\theta )+r(\eta (\theta )-{\eta }_0(\theta )))]{|}_{r=0}=0,\text{for}\text{all}\theta \in \mathrm{\Theta },$$

as ${\eta }_0(\theta )={\beta }_{\theta }$ solves the optimization problem  

$$\underset{\beta \in \mathcal{B}}{max}{E}_P[\ell (W;\theta ,\beta )],\text{for}\text{all}\theta \in \mathrm{\Theta }.$$

Here, the regularity conditions are needed to make sure that we can interchange $E_P$ and ${\partial }_{\theta }$ and also ${\partial }_{\theta }$ and ${\partial }_r$ in A.2. $\square$

 

Proof of Lemma 2.6:

First, we demonstrate that ${\mu }_0\in {\mathcal{L}}^1(\mathcal{R};{\mathbb{R}}^{d_{\theta }\times d_m})$⁠. Indeed,  

$$\begin{array}{ccc}\hfill E_P[\parallel {\mu }_0(R)\parallel ]& \le & E_P{[\parallel A(R)}^{\prime }\mathrm{\Omega }{(R)}^{-1}\parallel ]+E_P{[\parallel G(Z)\mathrm{\Gamma }(R)\mathrm{\Omega }(R)}^{-1}\parallel ]\hfill \\ & \le & E_P[\parallel A(R)\parallel \times {\parallel \mathrm{\Omega }(R)\parallel }^{-1}]+E_P[\parallel G(Z)\parallel \times \parallel \mathrm{\Gamma }(R)\parallel \times {\parallel \mathrm{\Omega }(R)\parallel }^{-1}]\hfill \\ & \le & (E_P{[\parallel A(R)\parallel }^2]\times E_P{[\parallel \mathrm{\Omega }(R)\parallel }^{-2}{])}^{1/2}\hfill \\ & & +(E_P{[\parallel G(Z)\parallel }^2\times {\parallel \mathrm{\Gamma }(R)\parallel }^2]\times E_P{[\parallel \mathrm{\Omega }(R)\parallel }^{-2}{])}^{1/2},\hfill \end{array}$$

which is finite by assumptions of the lemma as  

$$E_P{[\parallel G(Z)\parallel }^2\times {\parallel \mathrm{\Gamma }(R)\parallel }^2]\le (E_P{[\parallel G(Z)\parallel }^4{]\times E_P[\mathrm{\Gamma }(R){\parallel }^4])}^{1/2}<\infty .$$

Next, we demonstrate that  

$$E_P[\parallel \psi (W,{\theta }_0,\eta )\parallel ]<\infty \text{for}\text{all}\eta \in T.$$

Indeed, for all $\eta \in T$⁠, there exist $\mu \in {\mathcal{L}}^1(\mathcal{R};{\mathbb{R}}^{d_{\theta }\times d_m})$ and $h\in \mathcal{H}$ such that $\eta =(\mu ,h)$⁠, and so  

$$\begin{array}{ccc}\hfill E_P[\parallel \psi (W,{\theta }_0,\eta )\parallel ]& =& E_P[\parallel \mu (X)m(W,{\theta }_0,h(Z))\parallel ]\hfill \\ & \le & E_P[\parallel \mu (R)\parallel \times \parallel m(W,{\theta }_0,h(Z))\parallel ]\hfill \\ & =& E_P[\parallel \mu (R)\parallel \times E_P[\parallel m(W,{\theta }_0,h(Z))\mid R]]\le C_{h}E[\parallel \mu (R)\parallel ],\hfill \end{array}$$

which is finite by assumptions of the lemma. Further, 2.1 holds because  

$$\begin{array}{ccc}\hfill E_P[\psi (W,{\theta }_0,{\eta }_0)]& =& E_P[{\mu }_0(R)m(W,{\theta }_0,h_0(Z))]\hfill \\ & =& E_P[{\mu }_0(R)E_P[m(W,{\theta }_0,h_0(Z))\mid R]]=0,\hfill \end{array}$$

A.3

where the last equality follows from 2.22.

Finally, we demonstrate that 2.3 holds. To do so, take any $\eta =(\mu ,h)\in {\mathcal{T}}_N=T$⁠. Then  

$$\begin{array}{ccc}& & E_P[\psi (W,{\theta }_0,{\eta }_0+r(\eta -{\eta }_0)]\hfill \\ & & =E_P[({\mu }_0(R)+r(\mu (R)-{\mu }_0(R)))m(W,{\theta }_0,h_0(Z)+r(h(Z)-h_0(Z)))],\hfill \end{array}$$

and so  

$${\partial }_{\eta }{E}_P[\psi (W,{\theta }_0,{\eta }_0)][\eta -{\eta }_0]={\mathcal{I}}_1+{\mathcal{I}}_2,$$

where  

$$\begin{array}{ccc}\hfill {\mathcal{I}}_1& =& E_P[(\mu (R)-{\mu }_0(R))m(W,{\theta }_0,h_0(Z)],\hfill \\ \hfill {\mathcal{I}}_2& =& E_P[{\mu }_0(R){\partial }_{v^{\prime }}m(W,{\theta }_0,v){|}_{v=h_0(Z)}(h(Z)-h_0(Z))].\hfill \end{array}$$

Here, ${\mathcal{I}}_1=0$ by the same argument as that in A.3 and ${\mathcal{I}}_2=0$ because  

$$\begin{array}{ccc}\hfill {\mathcal{I}}_2& =& E_P[{\mu }_0(R)E_P[{\partial }_{v^{\prime }}m(W,{\theta }_0,v){|}_{v=h_0(Z)}\mid X](h(Z)-h_0(Z))]\hfill \\ & =& E_P[{\mu }_0(R)\mathrm{\Gamma }(X)(h(Z)-h_0(Z))]=E_P[E_P[{\mu }_0(R)\mathrm{\Gamma }(R)\mid Z](h(Z)-h_0(Z))]=0\hfill \end{array}$$

because  

$$\begin{array}{ccc}\hfill E_P[{\mu }_0(R)\mathrm{\Gamma }(X)\mid Z]& =& E_P[A{(R)}^{\prime }\mathrm{\Omega }{(R)}^{-1}\mathrm{\Gamma }(R)\mid Z]-E_P[G(Z)\mathrm{\Gamma }{(R)}^{\prime }\mathrm{\Omega }{(R)}^{-1}\mathrm{\Gamma }(R)\mid Z]\hfill \\ & =& E_P[A{(R)}^{\prime }\mathrm{\Omega }{(R)}^{-1}\mathrm{\Gamma }(R)\mid Z]-G(Z)E_P[\mathrm{\Gamma }{(R)}^{\prime }\mathrm{\Omega }{(R)}^{-1}\mathrm{\Gamma }(R)\mid Z]\hfill \\ & =& E_P[A{(R)}^{\prime }\mathrm{\Omega }{(R)}^{-1}\mathrm{\Gamma }(R)\mid Z]-E_P[A{(R)}^{\prime }\mathrm{\Omega }{(R)}^{-1}\mathrm{\Gamma }(R)\mid Z]\hfill \\ & & \times {\left(E_P[\mathrm{\Gamma }{(R)}^{\prime }\mathrm{\Omega }{(R)}^{-1}\mathrm{\Gamma }(R)\mid Z]\right)}^{-1}\times E_P[\mathrm{\Gamma }{(R)}^{\prime }\mathrm{\Omega }{(R)}^{-1}\mathrm{\Gamma }(R)\mid Z]\hfill \\ & =& E_P[A{(R)}^{\prime }\mathrm{\Omega }{(R)}^{-1}\mathrm{\Gamma }(R)\mid Z]-E_P[A{(R)}^{\prime }\mathrm{\Omega }{(R)}^{-1}\mathrm{\Gamma }(R)\mid Z]=0.\hfill \end{array}$$

$\square$

 

Proof of Theorem 3.1 Dml2 Case

: To start with, note that 3.11 follows immediately from the assumptions. Hence, it suffices to show that 3.10 holds uniformly over $P\in {\mathcal{P}}_N$⁠.

Fix any sequence ${\{P_N\}}_{N\ge 1}$ such that $P_N\in {\mathcal{P}}_N$ for all $N\ge 1$⁠. Because this sequence is chosen arbitrarily, to show that 3.10 holds uniformly over $P\in {\mathcal{P}}_N$⁠, it suffices to show that  

$${\displaystyle \sqrt{N}{\sigma }^{-1}({\stackrel{\sim }{\theta }}_0-{\theta }_0)=\frac{1}{\sqrt{N}}\sum _{i=1}^N\overline{\psi }(W_i)+O_{P_N}({\rho }_N)⇝N(0,{\mathrm{I}}_d).}$$

A.4

To do so, we proceed in five steps. Step 1 shows the main argument, and Steps 2–5 present auxiliary calculations. In the proof, it will be convenient to denote by ${\mathcal{E}}_N$ the event that ${\widehat{\eta }}_{0,k}\in {\mathcal{T}}_N$ for all $k\in [K]$⁠. Note that by Assumption 3.2 and the union bound, ${Pr}_{P_N}({\mathcal{E}}_N)\ge 1-K{\mathrm{\Delta }}_n=1-o(1)$ as ${\mathrm{\Delta }}_n=o(1)$⁠.

• Step 1. Denote  

  $$\begin{array}{ccc}\hfill {\widehat{J}}_0& :=& {\displaystyle \frac{1}{K}\sum _{k=1}^K{E}_{n,k}[{\psi }^a(W;{\widehat{\eta }}_{0,k})],R_{N,1}:={\widehat{J}}_0-J_0,}\hfill \\ \hfill R_{N,2}& :=& {\displaystyle \frac{1}{K}\sum _{k=1}^K{E}_{n,k}[\psi (W;{\theta }_0,{\widehat{\eta }}_{0,k})]-\frac{1}{N}\sum _{i=1}^N\psi (W_i;{\theta }_0,{\eta }_0).}\hfill \end{array}$$

  In Steps 2–5 respectively, we will show that  

  $$\begin{array}{ccc}\hfill \parallel R_{N,1}\parallel & =& O_{P_N}(N^{-1/2}+r_N),\hfill \end{array}$$

  A.5

   

  $$\begin{array}{ccc}\hfill \parallel R_{N,2}\parallel & =& O_{P_N}(N^{-1/2}{r}_N^{\prime }+{\lambda }_N+{\lambda }_N^{\prime }),\hfill \end{array}$$

  A.6

   

  $$\begin{array}{ccc}\hfill {\displaystyle \parallel N^{-1/2}\sum _{i=1}^N\psi (W_i;{\theta }_0,{\eta }_0)\parallel }& =& O_{P_N}(1),\hfill \end{array}$$

  A.7

   

  $$\begin{array}{ccc}\hfill \parallel {\sigma }^{-1}\parallel & =& O_{P_N}(1).\hfill \end{array}$$

  A.8

  Because $N^{-1/2}+r_N\le {\rho }_N=o(1)$ and all singular values of J₀ are bounded below from zero by Assumption 3.1, it follows from 189 that with $P_N$‐probability $1-o(1)$⁠, all singular values of ${\widehat{J}}_0$ are bounded below from zero as well. Therefore, with the same $P_N$‐probability,  

  $${\displaystyle {\stackrel{\sim }{\theta }}_0=-{\widehat{J}}_0^{-1}\frac{1}{K}\sum _{k=1}^K{E}_{n,k}[{\psi }^b(W;{\widehat{\eta }}_{0,k})]}$$

  and  

  $$\begin{array}{ccc}\hfill \sqrt{N}({\stackrel{\sim }{\theta }}_0-{\theta }_0)& =& {\displaystyle -\sqrt{N}{\widehat{J}}_0^{-1}\left(\frac{1}{K}\sum _{k=1}^K{E}_{n,k}[{\psi }^b(W;{\widehat{\eta }}_{0,k})]+{\widehat{J}}_0{\theta }_0\right)}\hfill \\ & =& {\displaystyle -\sqrt{N}{\widehat{J}}_0^{-1}\frac{1}{K}\sum _{k=1}^K{E}_{n,k}[\psi (W;{\theta }_0,{\widehat{\eta }}_{0,k})]}\hfill \\ & =& {\displaystyle -{\left(J_0+R_{N,1}\right)}^{-1}\times \left(\frac{1}{\sqrt{N}}\sum _{i=1}^N\psi (W_i;{\theta }_0,{\eta }_0)+\sqrt{N}{R}_{N,2}\right).}\hfill \end{array}$$

  A.9

  In addition, given that  

  $$\begin{array}{ccc}\hfill {(J_0+R_{N,1})}^{-1}-J_0^{-1}& =& {(J_0+R_{N,1})}^{-1}(J_0-(J_0+R_{N,1}))J_0^{-1}\hfill \\ & =& -{(J_0+R_{N,1})}^{-1}{R}_{N,1}{J}_0^{-1},\hfill \end{array}$$

  it follows from 189 that  

  $$\begin{array}{ccc}\hfill \parallel {(J_0+R_{N,1})}^{-1}-J_0^{-1}\parallel & \le & \parallel {(J_0+R_{N,1})}^{-1}\parallel \times \parallel R_{N,1}\parallel \times \parallel J_0^{-1}\parallel \hfill \\ & =& O_{P_N}(1)O_{P_N}(N^{-1/2}+r_N)O_{P_N}(1)=O_{P_N}(N^{-1/2}+r_N).\hfill \end{array}$$

  A.10

  Moreover, because $r_N^{\prime }+\sqrt{N}({\lambda }_N+{\lambda }_N^{\prime })\le {\rho }_N=o(1)$⁠, it follows from 190 and 191 that  

  $$\begin{array}{ccc}\hfill {\displaystyle ‖\frac{1}{\sqrt{N}}\sum _{i=1}^N\psi (W_i;{\theta }_0,{\eta }_0)+\sqrt{N}{R}_{N,2}‖}& \le & {\displaystyle ‖\frac{1}{\sqrt{N}}\sum _{i=1}^N\psi (W_i;{\theta }_0,{\eta }_0)‖+\parallel \sqrt{N}{R}_{N,2}\parallel }\hfill \\ & =& O_{P_N}(1)+o_{P_N}(1)=O_{P_N}(1).\hfill \end{array}$$

  A.11

  Combining 196 and 197 gives  

  $$\begin{array}{ccc}& & {\displaystyle ‖({(J_0+R_{N,1})}^{-1}-J_0^{-1})\times \left(\frac{1}{\sqrt{N}}\sum _{i=1}^N\psi (W_i;{\theta }_0,{\eta }_0)+\sqrt{N}{R}_{N,2}\right)‖}\hfill \\ & & {\displaystyle \le ‖{(J_0+R_{N,1})}^{-1}-J_0^{-1}‖\times ‖\frac{1}{\sqrt{N}}\sum _{i=1}^N\psi (W_i;{\theta }_0,{\eta }_0)+\sqrt{N}{R}_{N,2}‖}\hfill \\ & & =O_{P_N}(N^{-1/2}+r_N).\hfill \end{array}$$

  Now, substituting the last bound into 194 yields  

  $$\begin{array}{ccc}\hfill \sqrt{N}({\stackrel{\sim }{\theta }}_0-{\theta }_0)& =& {\displaystyle -J_0^{-1}\times \left(\frac{1}{\sqrt{N}}\sum _{i=1}^N\psi (W_i;{\theta }_0,{\eta }_0)+\sqrt{N}{R}_{N,2}\right)+O_{P_N}(N^{-1/2}+r_N)}\hfill \\ & =& {\displaystyle -J_0^{-1}\times \frac{1}{\sqrt{N}}\sum _{i=1}^N\psi (W_i;{\theta }_0,{\eta }_0)+O_{P_N}({\rho }_N),}\hfill \end{array}$$

  where in the second line we used 190 and the definition of ${\rho }_N$⁠. Combining this with 192 gives  

  $${\displaystyle \sqrt{N}{\sigma }^{-1}({\stackrel{\sim }{\theta }}_0-{\theta }_0)=\frac{1}{\sqrt{N}}\sum _{i=1}^N\overline{\psi }(W_i)+O_{P_N}({\rho }_N)}$$

  A.12

  by the definition of $\overline{\psi }$ given in the statement of the theorem. In turn, because ${\rho }_N=o(1)$⁠, combining 200 with the Lindeberg–Feller central limit theorem (CLT) and the Cramer–Wold device yields A.4. To complete the proof of the theorem, it remains to establish the bounds 189–192. We do so in the following four steps.
• Step 2. In this step, we establish 189. Because K is a fixed integer, which is independent of N, it suffices to show that for any $k\in [K]$⁠,  

  $$\parallel E_{n,k}[{\psi }^a(W;{\widehat{\eta }}_{0,k})]-E_{P_N}[{\psi }^a(W;{\eta }_0)]\parallel =O_{P_N}(N^{-1/2}+r_N).$$

  A.13

  To do so, fix any $k\in [K]$ and observe that by the triangle inequality,  

  $$\parallel E_{n,k}[{\psi }^a(W;{\widehat{\eta }}_{0,k})]-E_{P_N}[{\psi }^a(W;{\eta }_0)]\parallel \le {\mathcal{I}}_{1,k}+{\mathcal{I}}_{2,k},$$

  A.14

  where  

  $$\begin{array}{ccc}\hfill {\mathcal{I}}_{1,k}& :=& \parallel E_{n,k}[{\psi }^a(W;{\widehat{\eta }}_{0,k})]-E_{P_N}[{\psi }^a(W;{\widehat{\eta }}_{0,k})\mid {(W_i)}_{i\in I_k^c}]\parallel ,\hfill \\ \hfill {\mathcal{I}}_{2,k}& :=& \parallel E_{P_N}[{\psi }^a(W;{\widehat{\eta }}_{0,k})\mid {(W_i)}_{i\in I_k^c}]-E_{P_N}[{\psi }^a(W;{\eta }_0)]\parallel .\hfill \end{array}$$

  To bound ${\mathcal{I}}_{2,k}$⁠, note that on the event ${\mathcal{E}}_N$⁠, which holds with $P_N$‐probability $1-o(1)$⁠,  

  $${\mathcal{I}}_{2,k}\le \underset{\eta \in {\mathcal{T}}_N}{sup}\parallel E_{P_N}[{\psi }^a(W;\eta )]-E_{P_N}[{\psi }^a(W;{\eta }_0)]\parallel =r_N,$$

  and so ${\mathcal{I}}_{2,k}=O_{P_N}(r_N)$⁠. To bound ${\mathcal{I}}_{1,k}$⁠, note that conditional on ${(W_i)}_{i\in I_k^c}$⁠, the estimator ${\widehat{\eta }}_{0,k}$ is non‐stochastic, and so on the event ${\mathcal{E}}_N$⁠,  

  $$\begin{array}{ccc}\hfill E_{P_N}[{\mathcal{I}}_{1,k}^2\mid {(W_i)}_{i\in I_k^c}]& \le & n^{-1}{E}_{P_N}[\parallel {\psi }^a(W;{\widehat{\eta }}_{0,k}){\parallel }^2\mid {(W_i)}_{i\in I_k^c}]\hfill \\ & \le & \underset{\eta \in {\mathcal{T}}_N}{sup}{n}^{-1}{E}_{P_N}[\parallel {\psi }^a(W;\eta ){\parallel }^2]\le c_1^2/n,\hfill \end{array}$$

  where the last inequality holds by Assumption 3.2. Hence, ${\mathcal{I}}_{1,k}=O_{P_N}(N^{-1/2})$ by Lemma 6.1. Combining the bounds ${\mathcal{I}}_{1,k}=O_{P_N}(N^{-1/2})$ and ${\mathcal{I}}_{2,k}=O_{P_N}(r_N)$ with 202 gives 201.
• Step 3. In this step, we establish 190. This is the step where we invoke the Neyman orthogonality (or near‐orthogonality) condition. Again, because K is a fixed integer, which is independent of N, it suffices to show that for any $k\in [K]$⁠,  

  $${\displaystyle E_{n,k}[\psi (W;{\theta }_0,{\widehat{\eta }}_{0,k})]-\frac{1}{n}\sum _{i\in I_k}\psi (W_i;{\theta }_0,{\eta }_0)=O_{P_N}(N^{-1/2}{r}_N^{\prime }+{\lambda }_N+{\lambda }_N^{\prime }).}$$

  A.15

  To do so, fix any $k\in [K]$ and introduce the following additional empirical process notation,  

  $${\displaystyle {\mathbb{G}}_{n,k}[\phi (W)]=\frac{1}{\sqrt{n}}\sum _{i\in I_k}\left(\phi (W_i)-\int \phi (w)d{P}_N\right),}$$

  where ϕ is any $P_N$‐integrable function on $\mathcal{W}$⁠. Then observe that by the triangle inequality,  

  $${\displaystyle ‖E_{n,k}[\psi (W;{\theta }_0,{\widehat{\eta }}_{0,k})]-\frac{1}{n}\sum _{i\in I_k}\psi (W_i;{\theta }_0,{\eta }_0)‖\le \frac{{\mathcal{I}}_{3,k}+{\mathcal{I}}_{4,k}}{\sqrt{n}},}$$

  A.16

  where  

  $$\begin{array}{ccc}\hfill {\mathcal{I}}_{3,k}& :=& \parallel {\mathbb{G}}_{n,k}[\psi (W;{\theta }_0,{\widehat{\eta }}_{0,k})]-{\mathbb{G}}_{n,k}[\psi (W;{\theta }_0,{\eta }_0)]\parallel ,\hfill \\ \hfill {\mathcal{I}}_{4,k}& :=& \sqrt{n}\parallel E_{P_N}[\psi (W;{\theta }_0,{\widehat{\eta }}_{0,k})\mid {(W_i)}_{i\in I_k^c}]-E_{P_N}[\psi (W;{\theta }_0,{\eta }_0)]\parallel .\hfill \end{array}$$

  To bound ${\mathcal{I}}_{3,k}$⁠, note that, as above, conditional on ${(W_i)}_{i\in I_k^c}$⁠, the estimator ${\widehat{\eta }}_{0,k}$ is non‐stochastic, and so on the event ${\mathcal{E}}_N$⁠,  

  $$\begin{array}{ccc}\hfill E_{P_N}[{\mathcal{I}}_{3,k}^2\mid {(W_i)}_{i\in I_k^c}]& =& E_{P_N}[\parallel \psi (W;{\theta }_0,{\widehat{\eta }}_{0,k})-\psi (W;{\theta }_0,{\eta }_0){\parallel }^2\mid {(W_i)}_{i\in I_k^c}]\hfill \\ & \le & \underset{\eta \in {\mathcal{T}}_N}{sup}{E}_{P_N}[\parallel \psi (W;{\theta }_0,\eta )-\psi (W;{\theta }_0,{\eta }_0){\parallel }^2\mid {(W_i)}_{i\in I_k^c}]\hfill \\ & \le & \underset{\eta \in {\mathcal{T}}_N}{sup}{E}_{P_N}[\parallel \psi (W;{\theta }_0,\eta )-\psi (W;{\theta }_0,{\eta }_0){\parallel }^2]={(r_N^{\prime })}^2\hfill \end{array}$$

  by the definition of $r_N^{\prime }$ in Assumption 3.2. Hence, ${\mathcal{I}}_{3,k}=O_{P_N}(r_N^{\prime })$ by Lemma 6.1. To bound ${\mathcal{I}}_{4,k}$⁠, introduce the function  

  $$f_k(r):=E_{P_N}[\psi (W;{\theta }_0,{\eta }_0+r({\widehat{\eta }}_{0,k}-{\eta }_0))\mid {(W_i)}_{i\in I_k^c}]-E_{P_N}[\psi (W;{\theta }_0,{\eta }_0)],r\in [0,1].$$

  Then, by Taylor expansion,  

  $$f_k(1)=f_k(0)+f_k^{\prime }(0)+f_k^{\prime \prime }(\stackrel{\sim }{r})/2,\text{for}\text{some}\stackrel{\sim }{r}\in (0,1).$$

  However, $\parallel f_k(0)\parallel =0$ because  

  $$E_{P_N}[\psi (W;{\theta }_0,{\eta }_0)\mid {(W_i)}_{i\in I_k^c}]=E_{P_N}[\psi (W;{\theta }_0,{\eta }_0)].$$

  In addition, on the event ${\mathcal{E}}_N$⁠, by the Neyman ${\lambda }_N$ near‐orthogonality condition imposed in Assumption 3.1,  

  $$\parallel f_k^{\prime }(0)\parallel =\parallel {\partial }_{\eta }{E}_{P_N}\psi (W;{\theta }_0,{\eta }_0)[{\widehat{\eta }}_{0,k}-{\eta }_0]\parallel \le {\lambda }_N.$$

  Moreover, on the event ${\mathcal{E}}_N$⁠,  

  $$\parallel f_k^{\prime \prime }(\stackrel{\sim }{r})\parallel \le \underset{r\in (0,1)}{sup}\parallel f_k^{\prime \prime }(r)\parallel \le {\lambda }_N^{\prime }$$

  by the definition ${\lambda }_N^{\prime }$ in Assumption 3.2. Hence,  

  $${\mathcal{I}}_{4,k}=\sqrt{n}\parallel f_k(1)\parallel =O_{P_N}(\sqrt{n}({\lambda }_N+{\lambda }_N^{\prime })).$$

  Combining the bounds on ${\mathcal{I}}_{3,k}$ and ${\mathcal{I}}_{4,k}$ with 208 and using the fact that $n^{-1}=O(N^{-1})$ gives 206.
• Step 4. To establish 191, note that  

  $${\displaystyle E_{P_N}\left[{‖\frac{1}{\sqrt{N}}\sum _{i=1}^N\psi (W_i;{\theta }_0,{\eta }_0)‖}^2\right]=E_{P_N}\left[\parallel \psi (W;{\theta }_0,{\eta }_0){\parallel }^2\right]\le c_1^2}$$

  by Assumption 3.2. Combining this with the Markov inequality gives 191.
• Step 5. Here we establish 192. Note that all eigenvalues of the matrix  

  $${\sigma }^2=J_0^{-1}{E}_P[\psi (W;{\theta }_0,{\eta }_0)\psi {(W;{\theta }_0,{\eta }_0)}^{\prime }]{(J_0^{-1})}^{\prime }$$

  are bounded from below by $c_0/c_1^2$ because all singular values of J₀ are bounded from above by c₁ by Assumption 3.1 and all eigenvalues of $E_P[\psi (W;{\theta }_0,{\eta }_0)\psi {(W;{\theta }_0,{\eta }_0)}^{\prime }]$ are bounded from below by c₀ by Assumption 3.2. Hence, given that $\parallel {\sigma }^{-1}\parallel$ is the largest eigenvalue of ${\sigma }^{-1}$⁠, it follows that $\parallel {\sigma }^{-1}\parallel =c_1/\sqrt{c_0}$⁠. This gives 192 and completes the proof of the theorem. $\square$

 

Proof of Theorem 3.1 Dml1 Case

: As in the case of the DML2 version, note that 3.11 follows immediately from the assumptions, and so it suffices to show that 3.10 holds uniformly over $P\in {\mathcal{P}}_N$⁠.

Fix any sequence ${\{P_N\}}_{N\ge 1}$ such that $P_N\in {\mathcal{P}}_N$ for all $N\ge 1$⁠. Because this sequence is chosen arbitrarily, to show that 3.10 holds uniformly over $P\in {\mathcal{P}}_N$⁠, it suffices to show that  

$${\displaystyle \sqrt{N}{\sigma }^{-1}({\stackrel{\sim }{\theta }}_0-{\theta }_0)=\frac{1}{\sqrt{N}}\sum _{i=1}^N\overline{\psi }(W_i)+O_{P_N}({\rho }_N)⇝N(0,{\mathrm{I}}_d).}$$

A.17

To do so, for all $k\in [K]$⁠, denote  

$$\begin{array}{ccc}\hfill {\widehat{J}}_{0,k}& :=& E_{n,k}[{\psi }^a(W;{\widehat{\eta }}_{0,k})],R_{N,1,k}:={\widehat{J}}_{0,k}-J_0,\hfill \\ \hfill R_{N,2,k}& :=& {\displaystyle E_{n,k}[\psi (W;{\theta }_0,{\widehat{\eta }}_{0,k})]-\frac{1}{n}\sum _{i\in I_k}\psi (W_i;{\theta }_0,{\eta }_0).}\hfill \end{array}$$

Because K is a fixed integer, which is independent of n, it follows by the same arguments as those in Steps 2–5 in the Proof of Theorem 3.1 (DML2 case) that  

$$\begin{array}{ccc}\hfill \underset{k\in [K]}{max}\parallel R_{N,1,k}\parallel & =& O_{P_N}(N^{-1/2}+r_N),\hfill \end{array}$$

A.18

 

$$\begin{array}{ccc}\hfill \underset{k\in [K]}{max}\parallel R_{N,2,k}\parallel & =& O_{P_N}(N^{-1/2}{r}_N^{\prime }+{\lambda }_N+{\lambda }_N^{\prime }),\hfill \end{array}$$

A.19

 

$$\begin{array}{ccc}\hfill {max}_{k\in [K]}{\displaystyle \parallel n^{-1/2}\sum _{i\in I_k}\psi (W_i;{\theta }_0,{\eta }_0)\parallel }& =& O_{P_N}(1),\hfill \end{array}$$

A.20

 

$$\begin{array}{ccc}\hfill \parallel {\sigma }^{-1}\parallel & =& O_{P_N}(1).\hfill \end{array}$$

A.21

Because $N^{-1/2}+r_N\le {\rho }_N=o(1)$ and all singular values of J₀ are bounded below from zero by Assumption 3.1, it follows from A.18 that for all $k\in [K]$⁠, with $P_N$‐probability $1-o(1)$⁠, all singular values of ${\widehat{J}}_{0,k}$ are bounded below from zero, and so with the same $P_N$‐probability,  

$${\check{\theta }}_{0,k}=-{\widehat{J}}_{0,k}^{-1}{E}_{n,k}[{\psi }^b(W;{\widehat{\eta }}_{0,k})].$$

Hence, by the same arguments as those in Step 1 in the Proof of Theorem 3.1 (DML2 case), it follows from the bounds A.18–A.21 that for all $k\in [K]$⁠,  

$${\displaystyle \sqrt{n}{\sigma }^{-1}({\check{\theta }}_{0,k}-{\theta }_0)=\frac{1}{\sqrt{n}}\sum _{i\in I_k}\overline{\psi }(W_i)+O_{P_N}({\rho }_N).}$$

Therefore,  

$${\displaystyle \sqrt{N}{\sigma }^{-1}({\stackrel{\sim }{\theta }}_0-{\theta }_0)=\sqrt{N}{\sigma }^{-1}\left(\frac{1}{K}\sum _{k=1}^K{\check{\theta }}_{0,k}-{\theta }_0\right)=\frac{1}{\sqrt{N}}\sum _{i=1}^N\overline{\psi }(W_i)+O_{P_N}({\rho }_N).}$$

A.22

In turn, because ${\rho }_N=o(1)$⁠, combining A.22 with the Lindeberg–Feller CLT and the Cramer–Wold device yields A.17 and completes the proof of the theorem. $\square$

 

Proof of Theorem 3.2:

In this proof, all bounds hold uniformly in $P\in {\mathcal{P}}_N$ for $N\ge 3$⁠, and we do not repeat this qualification throughout. Also, the second asserted claim follows immediately from the first one and Theorem 3.1. Hence, it suffices to prove the first asserted claim.

In the Proof of Theorem 3.1 (DML2 case), we established that $\parallel {\widehat{J}}_0-J_0\parallel =O_P(r_N+N^{-1/2})$⁠. Hence, because $\parallel J_0^{-1}\parallel \le c_0^{-1}$ by Assumption 3.1 and  

$$\parallel E_P[\psi (W;{\theta }_0,{\eta }_0)\psi {(W;{\theta }_0,{\eta }_0)}^{\prime }]\parallel \le E_P[\parallel \psi (W;{\theta }_0,{\eta }_0){\parallel }^2]\le c_1^2$$

by Assumption 3.2, it suffices to show that  

$${\displaystyle ‖\frac{1}{K}\sum _{k=1}^K{\mathbb{E}}_{n,k}[\psi (W;{\stackrel{\sim }{\theta }}_0,{\widehat{\eta }}_{0,k})\psi {(W;{\stackrel{\sim }{\theta }}_0,{\widehat{\eta }}_{0,k})}^{\prime }]-E_P[\psi (W;{\theta }_0,{\eta }_0)\psi {(W;{\theta }_0,{\eta }_0)}^{\prime }]‖=O_P({\varrho }_N).}$$

Moreover, as both K and $d_{\theta }$⁠, the dimension of ψ, are fixed integers, which are independent of N, the last bound will follow if we show that for all $k\in [K]$ and all $j,k\in [d_{\theta }]$⁠,  

$${\mathcal{I}}_{\mathit{kjl}}:=|{\mathbb{E}}_{n,k}[{\psi }_j(W;{\stackrel{\sim }{\theta }}_0,{\widehat{\eta }}_{0,k}){\psi }_l(W;{\stackrel{\sim }{\theta }}_0,{\widehat{\eta }}_{0,k})]-E_P[{\psi }_j(W;{\theta }_0,{\eta }_0){\psi }_l(W;{\theta }_0,{\eta }_0)]|$$

satisfies  

$${\mathcal{I}}_{\mathit{kjl}}=O_P({\varrho }_N).$$

A.23

To do so, observe that by the triangle inequality,  

$${\mathcal{I}}_{\mathit{kjl}}\le {\mathcal{I}}_{\mathit{kjl},1}+{\mathcal{I}}_{\mathit{kjl},2},$$

A.24

where  

$$\begin{array}{ccc}\hfill {\mathcal{I}}_{\mathit{kjl},1}& :=& |E_{n,k}[{\psi }_j(W;{\stackrel{\sim }{\theta }}_0,{\widehat{\eta }}_{0,k}){\psi }_l(W;{\stackrel{\sim }{\theta }}_0,{\widehat{\eta }}_{0,k})]-E_{n,k}[{\psi }_j(W;{\theta }_0,{\eta }_0){\psi }_l(W;{\theta }_0,{\eta }_0)]|,\hfill \\ \hfill {\mathcal{I}}_{\mathit{kjl},2}& :=& |{\mathbb{E}}_{n,k}[{\psi }_j(W;{\theta }_0,{\eta }_0){\psi }_l(W;{\theta }_0,{\eta }_0)]-E_P[{\psi }_j(W;{\theta }_0,{\eta }_0){\psi }_l(W;{\theta }_0,{\eta }_0)]|.\hfill \end{array}$$

We bound ${\mathcal{I}}_{\mathit{kjl},2}$ first. If $q\ge 4$⁠, then  

$$\begin{array}{ccc}\hfill E_P[{\mathcal{I}}_{\mathit{kjl},2}^2]& \le & n^{-1}{E}_P\left[{({\psi }_j(W;{\theta }_0,{\eta }_0){\psi }_l(W;{\theta }_0,{\eta }_0))}^2\right]\hfill \\ & \le & n^{-1}{\left(E_P[{\psi }_j^4(W;{\theta }_0,{\eta }_0)]\times E_P[{\psi }_l^4(W;{\theta }_0,{\eta }_0)]\right)}^{1/2}\hfill \\ & \le & n^{-1}{E}_P[\parallel \psi (W;{\theta }_0,{\eta }_0){\parallel }^4]\le c_1^4,\hfill \end{array}$$

where the second line holds by the Hölder inequality, and the third by Assumption 3.2. Hence, ${\mathcal{I}}_{\mathit{kjl},2}=O_P(N^{-1/2})$⁠. If $q\in (2,4)$⁠, we apply the following von Bahr–Esseen inequality with $p=q/2$⁠: if $X_1,\dots ,X_n$ are independent random variables with mean zero, then for any $p\in [1,2]$⁠,  

$${\displaystyle E\left[{\left|\sum _{i=1}^n{X}_i\right|}^p\right]\le \left(2-\frac{1}{n}\right)\sum _{i=1}^{n}E[|X_i{|}^p];}$$

see DasGupta (2008, p. 650). This gives  

$$\begin{array}{ccc}\hfill E_P[{\mathcal{I}}_{kjl,2}^{q/2}]& \lesssim & n^{-q/2+1}{E}_P[{({\psi }_j(W;{\theta }_0,{\eta }_0){\psi }_l(W;{\theta }_0,{\eta }_0))}^{q/2}]\hfill \\ & \le & n^{-q/2+1}{E}_P[\parallel \psi (W;{\theta }_0,{\eta }_0){\parallel }^q]\lesssim n^{-q/2+1}\hfill \end{array}$$

by Assumption 3.2. Hence, ${\mathcal{I}}_{\mathit{kjl},2}=O_P(N^{2/q-1})$⁠. We conclude that  

$${\mathcal{I}}_{\mathit{kjl},2}=O_P(N^{-[(1-2/q)\wedge (1/2)]}).$$

A.25

Next, we bound ${\mathcal{I}}_{\mathit{kjl},1}$⁠. To do so, observe that for any numbers a, b, $\delta a$ and $\delta b$ such that $|a|\vee |b|\le c$ and $|\delta a|\vee |\delta b|\le r$⁠, we have  

$$|(a+\delta a)(b+\delta b)-ab|\le 2r(c+r).$$

Denoting  

$${\psi }_{\mathit{hi}}:={\psi }_h(W_i;{\theta }_0,{\eta }_0)\text{and}{\widehat{\psi }}_{\mathit{h}\mathit{i}}:={\psi }_h(W_i;{\stackrel{\sim }{\theta }}_0,{\widehat{\eta }}_{0,k}),\text{for}(h,i)\in \{j,l\}\times I_k,$$

and applying the inequality above with $a:={\psi }_{\mathit{j}\mathit{i}}$⁠, $b:={\psi }_{\mathit{l}\mathit{i}}$⁠, $a+\delta a:={\widehat{\psi }}_{\mathit{j}\mathit{i}}$⁠, $b+\delta b:={\widehat{\psi }}_{\mathit{li}}$⁠, $r:=|{\widehat{\psi }}_{\mathit{j}\mathit{i}}-{\psi }_{\mathit{j}\mathit{i}}|\vee |{\widehat{\psi }}_{\mathit{l}\mathit{i}}-{\psi }_{\mathit{l}\mathit{i}}|$ and $c:=|{\psi }_{\mathit{j}\mathit{i}}|\vee |{\psi }_{\mathit{l}\mathit{i}}|$ gives  

$$\begin{array}{ccc}\hfill {\mathcal{I}}_{\mathit{kjl},1}& =& {\displaystyle \left|\frac{1}{n}\sum _{i\in I_k}{\widehat{\psi }}_{\mathit{j}\mathit{i}}{\widehat{\psi }}_{\mathit{l}\mathit{i}}-{\psi }_{\mathit{j}\mathit{i}}{\psi }_{\mathit{l}\mathit{i}}\right|\le \frac{1}{n}\sum _{i\in I_k}|{\widehat{\psi }}_{\mathit{j}\mathit{i}}{\widehat{\psi }}_{\mathit{l}\mathit{i}}-{\psi }_{\mathit{j}\mathit{i}}{\psi }_{\mathit{l}\mathit{i}}|}\hfill \\ & \le & {\displaystyle \frac{2}{n}\sum _{i\in I_k}(|{\widehat{\psi }}_{\mathit{j}\mathit{i}}-{\psi }_{\mathit{j}\mathit{i}}|\vee |{\widehat{\psi }}_{\mathit{l}\mathit{i}}-{\psi }_{\mathit{l}\mathit{i}}|)\times (|{\psi }_{\mathit{j}\mathit{i}}|\vee |{\psi }_{\mathit{l}\mathit{i}}|+|{\widehat{\psi }}_{\mathit{j}\mathit{i}}-{\psi }_{\mathit{j}\mathit{i}}|\vee |{\widehat{\psi }}_{\mathit{l}\mathit{i}}-{\psi }_{\mathit{l}\mathit{i}}|)}\hfill \\ & \le & {\displaystyle {\left(\frac{2}{n}\sum _{i\in I_k}(|{\widehat{\psi }}_{\mathit{j}\mathit{i}}-{\psi }_{\mathit{j}\mathit{i}}{|}^2\vee |{\widehat{\psi }}_{\mathit{l}\mathit{i}}-{\psi }_{\mathit{l}\mathit{i}}{|}^2)\right)}^{1/2}}\hfill \\ & & {\displaystyle \times {\left(\frac{2}{n}\sum _{i\in I_k}(|{\psi }_{\mathit{j}\mathit{i}}|\vee |{\psi }_{\mathit{l}\mathit{i}}|+|{\widehat{\psi }}_{\mathit{j}\mathit{i}}-{\psi }_{\mathit{j}\mathit{i}}|\vee |{\widehat{\psi }}_{\mathit{l}\mathit{i}}-{\psi }_{\mathit{l}\mathit{i}}{|)}^2\right)}^{1/2}.}\hfill \end{array}$$

In addition, the expression in the last line above is bounded by  

$${\displaystyle {\left(\frac{2}{n}\sum _{i\in I_k}|{\psi }_{\mathit{j}\mathit{i}}{|}^2\vee {|{\psi }_{\mathit{l}\mathit{i}}|}^2\right)}^{1/2}+{\left(\frac{2}{n}\sum _{i\in I_k}|{\widehat{\psi }}_{\mathit{j}\mathit{i}}-{\psi }_{\mathit{j}\mathit{i}}{|}^2\vee {|{\widehat{\psi }}_{\mathit{l}\mathit{i}}-{\psi }_{\mathit{l}\mathit{i}}|}^2\right)}^{1/2},}$$

and so  

$${\displaystyle {\mathcal{I}}_{\mathit{kjl},1}^2\lesssim R_N\times \left(\frac{1}{n}\sum _{i\in I_k}{\parallel \psi (W_i;{\theta }_0,{\eta }_0)\parallel }^2+R_N\right),}$$

where  

$${\displaystyle R_N:=\frac{1}{n}\sum _{i\in I_k}{\parallel \psi (W_i;{\stackrel{\sim }{\theta }}_0,{\widehat{\eta }}_{0,k})-\psi (W_i;{\theta }_0,{\eta }_0)\parallel }^2.}$$

Moreover,  

$${\displaystyle \frac{1}{n}\sum _{i\in I_k}{\parallel \psi (W_i;{\theta }_0,{\eta }_0)\parallel }^2=O_P(1),}$$

by the Markov inequality because  

$${\displaystyle E_P\left[\frac{1}{n}\sum _{i\in I_k}{\parallel \psi (W_i;{\theta }_0,{\eta }_0)\parallel }^2\right]=E_P[\parallel \psi (W;{\theta }_0,{\eta }_0{\parallel }^2]\le c_1^2}$$

by Assumption 3.2. It remains to bound $R_N$⁠. We have  

$${\displaystyle R_N\lesssim \frac{1}{n}\sum _{i\in I_k}\parallel {\psi }^a(W_i;{\widehat{\eta }}_{0,k})({\stackrel{\sim }{\theta }}_0-{\theta }_0){\parallel }^2+\frac{1}{n}\sum _{i\in I_k}{\parallel \psi (W_i;{\theta }_0,{\widehat{\eta }}_{0,k})-\psi (W_i;{\theta }_0,{\eta }_0)\parallel }^2.}$$

A.26

The first term on the right‐hand side of A.26 is bounded from above by  

$${\displaystyle \left(\frac{1}{n}\sum _{i\in I_k}{\parallel {\psi }^a(W_i;{\widehat{\eta }}_{0,k})\parallel }^2\right)\times {\parallel {\stackrel{\sim }{\theta }}_0-{\theta }_0\parallel }^2=O_P(1)\times O_P(N^{-1})=O_P(N^{-1}),}$$

and the conditional expectation of the second term given ${(W_i)}_{i\in I_k^c}$ on the event that ${\widehat{\eta }}_{0,k}\in {\mathcal{T}}_N$ is equal to  

$$\begin{array}{ccc}& & E_P[\parallel \psi (W;{\theta }_0,{\widehat{\eta }}_{0,k})-\psi (W;{\theta }_0,{\eta }_0){\parallel }^2\mid {(W_i)}_{i\in I_k^c}]\hfill \\ & & \le \underset{\eta \in {\mathcal{T}}_N}{sup}{E}_P[\parallel \psi (W;{\theta }_0,\eta )-\psi (W;{\theta }_0,{\eta }_0){\parallel }^2\mid {(W_i)}_{i\in I_k^c}]={(r_N^{\prime })}^2.\hfill \end{array}$$

Because the event that ${\widehat{\eta }}_{0,k}\in {\mathcal{T}}_N$ holds with probability $1-{\mathrm{\Delta }}_N=1-o(1)$⁠, it follows that $R_N=O_P(N^{-1}+{(r_N^{\prime })}^2)$⁠, and so  

$${\mathcal{I}}_{\mathit{kjl},1}=O_P(N^{-1/2}+r_N^{\prime }).$$

A.27

Combining the bounds A.25 and A.27 with A.24 gives A.23 and completes the proof of the theorem. $\square$

 

Proof of Theorem 3.3:

We only consider the case of the DML1 estimator and note that the DML2 estimator can be treated similarly.

With the help of Lemma 6.3, which establishes approximate linearity of the subsample DML estimators ${\check{\theta }}_{0,k}$ and is presented below, the proof is the same as that in the linear case presented in the Proof of Theorem 3.1 (DML1 case).

 

Lemma 6.3.Linearization for Subsample Dml in Non‐Linear Problems

Under the conditions of Theorem 3.3, for any $k=1,\dots ,K$⁠, the estimator ${\check{\theta }}_0={\check{\theta }}_{0,k}$ defined by 3.2 obeys  

$${\displaystyle \sqrt{n}{\sigma }_0^{-1}({\check{\theta }}_0-{\theta }_0)=\frac{1}{\sqrt{n}}\sum _{i\in I}\overline{\psi }(W_i)+O_P({\rho }_n^{\prime })}$$

A.28

uniformly over $P\in {\mathcal{P}}_N$⁠, where ${\rho }_n^{\prime }=n^{-1/2}+r_N+r_N^{\prime }+n^{1/2}{\lambda }_N+n^{1/2}{\lambda }_N^{\prime }\lesssim {\delta }_N$ and where $\overline{\psi }(·):=-{\sigma }_0^{-1}{J}_0^{-1}\psi (·,{\theta }_0,{\eta }_0)$⁠.

 

Proof:

Fix any $k=1,\dots ,K$ and any sequence ${\{P_N\}}_{N\ge 1}$ such that $P_N\in {\mathcal{P}}_N$ for all $N\ge 1$⁠. To prove the asserted claim, it suffices to show that the estimator ${\check{\theta }}_0={\check{\theta }}_{0,k}$ satisfies A.28 with P replaced by $P_N$⁠. To do so, we split the proof into four steps. In the proof, we use $E_n$⁠, ${\mathbb{G}}_n$⁠, I and ${\widehat{\eta }}_0$ instead of $E_{n,k}$⁠, ${\mathbb{G}}_{n,k}$⁠, $I_k$ and ${\widehat{\eta }}_{0,k}$⁠, respectively.

• Step 1. (Preliminary rate result) We claim that with $P_N$‐probability $1-o(1)$⁠,  

  $$\parallel {\check{\theta }}_0-{\theta }_0\parallel \le {\tau }_N.$$

  A.29

  To show this claim, note that the definition of ${\check{\theta }}_0$ implies that  

  $$\parallel E_n[\psi (W;{\check{\theta }}_0,{\widehat{\eta }}_0)]\parallel \le \parallel E_n[\psi (W;{\theta }_0,{\widehat{\eta }}_0)]\parallel +{\epsilon }_N,$$

  which in turn implies via the triangle inequality that, with $P_N$‐probability $1-o(1)$⁠,  

  $$\parallel E_{P_N}[\psi (W;\theta ,{\eta }_0)]{|}_{\theta ={\check{\theta }}_0}\parallel \le {\epsilon }_N+2{\mathcal{I}}_1+2{\mathcal{I}}_2,$$

  A.30

  where  

  $$\begin{array}{ccc}\hfill {\mathcal{I}}_1& :=& \underset{\theta \in \mathrm{\Theta },\eta \in {\mathcal{T}}_N}{sup}\parallel E_{P_N}[\psi (W;\theta ,\eta )]-E_{P_N}[\psi (W;\theta ,{\eta }_0)]\parallel ,\hfill \\ \hfill {\mathcal{I}}_2& :=& \underset{\eta \in \{{\eta }_0,{\widehat{\eta }}_0\}}{max}\underset{\theta \in \mathrm{\Theta }}{sup}\parallel E_n[\psi (W;\theta ,\eta )]-E_{P_N}[\psi (W;\theta ,\eta )]\parallel .\hfill \end{array}$$

  Here ${\epsilon }_N=o({\tau }_N)$ because ${\epsilon }_N=o({\delta }_N{N}^{-1/2})$⁠, ${\delta }_N=o(1)$ and ${\tau }_N\ge c_0{N}^{-1/2}logn$⁠. Also, ${\mathcal{I}}_1=r_N\le {\delta }_N{\tau }_N=o({\tau }_N)$ by Assumption 3.4(c). Moreover, applying Lemma 6.2 to the function class ${\mathcal{F}}_{1,\eta }$ for $\eta ={\eta }_0$ and $\eta ={\widehat{\eta }}_0$ defined in Assumption 3.4, conditional on ${(W_i)}_{i\in I^c}$ and $I^c$⁠, so that ${\widehat{\eta }}_0$ is fixed after conditioning, shows that with $P_N$‐probability $1-o(1)$⁠,  

  $$I_2\lesssim N^{-1/2}(1+N^{-1/2+1/q}logn)\lesssim N^{-1/2}=o({\tau }_N).$$

  Hence, it follows from 253 and Assumption 3.3 that with $P_N$‐probability $1-o(1)$⁠,  

  $$\parallel J_0({\check{\theta }}_0-{\theta }_0)\parallel \wedge c_0\le \parallel E_{P_N}[\psi (W;\theta ,{\eta }_0)]{|}_{\theta ={\check{\theta }}_0}\parallel =o({\tau }_N).$$

  A.31

  Combining this bound with the fact that the singular values of J₀ are bounded away from zero, which holds by Assumption 3.3, gives the claim of this step.
• Step 2. (Linearization) Here we prove the claim of the lemma. First, by definition of ${\check{\theta }}_0$⁠, we have  

  $$\sqrt{n}\parallel E_n[\psi (W;{\check{\theta }}_0,{\widehat{\eta }}_0)]\parallel \le \underset{\theta \in \mathrm{\Theta }}{inf}\sqrt{n}\parallel E_n[\psi (W;\theta ,{\widehat{\eta }}_0)]\parallel +{\epsilon }_N\sqrt{n}.$$

  A.32

  Also, it will be shown in Step 4 that  

  $$\begin{array}{ccc}\hfill {\mathcal{I}}_3& :=& \underset{\theta \in \mathrm{\Theta }}{inf}\sqrt{n}\parallel E_n[\psi (W;\theta ,{\widehat{\eta }}_0)]\parallel \hfill \\ & =& O_{P_N}(n^{-1/2+1/q}logn+r_N^{\prime }{log}^{1/2}(1/r_N^{\prime })+{\lambda }_N\sqrt{n}+{\lambda }_N^{\prime }\sqrt{n}).\hfill \end{array}$$

  A.33

  Moreover, for any $\theta \in \mathrm{\Theta }$ and $\eta \in {\mathcal{T}}_N$⁠, we have  

  $$\begin{array}{ccc}\hfill \sqrt{n}{E}_n[\psi (W;\theta ,\eta )]& =& \sqrt{n}{E}_n[\psi (W;{\theta }_0,{\eta }_0)]+{\mathbb{G}}_n[\psi (W;\theta ,\eta )-\psi (W;{\theta }_0,{\eta }_0)]\hfill \\ & & +\sqrt{n}(E_{P_N}[\psi (W;\theta ,\eta )),\hfill \end{array}$$

  A.34

  where we are using the fact that $E_{P_N}[\psi (W;{\theta }_0,{\eta }_0)]=0$⁠. Finally, by Taylor expansion of the function $r\mapsto E_{P_N}[\psi (W;{\theta }_0+r(\theta -{\theta }_0),{\eta }_0+r(\eta -{\eta }_0))]$⁠, which vanishes at $r=0$⁠,  

  $$\begin{array}{ccc}\hfill E_{P_N}[\psi (W;\theta ,\eta )]& =& J_0(\theta -{\theta }_0)+{\partial }_{\eta }{E}_{P_N}\psi (W;{\theta }_0,{\eta }_0)[\eta -{\eta }_0]\hfill \\ & & {\displaystyle +{\int }_0^1{2}^{-1}{\partial }_r^2{E}_{P_N}[W;{\theta }_0+r(\theta -{\theta }_0),{\eta }_0+r(\eta -{\eta }_0)]\mathit{d}\mathit{r}.}\hfill \end{array}$$

  A.35

  Therefore, because $\parallel {\check{\theta }}_0-{\theta }_0\parallel \le {\tau }_N$ and $\eta \in {\mathcal{T}}_N$ with $P_N$‐probability $1-o(1)$⁠, and because by Neyman ${\lambda }_N$‐near orthogonality,  

  $$\parallel {\partial }_{\eta }{E}_{P_N}[\psi (W;{\theta }_0,{\eta }_0)][{\widehat{\eta }}_0-{\eta }_0]\parallel \le {\lambda }_N,$$

  applying 259 with $\theta ={\check{\theta }}_0$ and $\eta ={\widehat{\eta }}_0$⁠, we have with $P_N$‐probability $1-o(1)$⁠,  

  $$\sqrt{n}‖E_n[\psi (W;{\theta }_0,{\eta }_0)]+J_0({\check{\theta }}_0-{\theta }_0)‖\le {\lambda }_N\sqrt{n}+{\epsilon }_N\sqrt{n}+{\mathcal{I}}_3+{\mathcal{I}}_4+{\mathcal{I}}_5,$$

  where by Assumption 3.4,  

  $${\displaystyle {\mathcal{I}}_4:=\sqrt{n}\underset{\parallel \theta -{\theta }_0\parallel \le {\tau }_N,\eta \in {\mathcal{T}}_N}{sup}‖{\int }_0^1{2}^{-1}{\partial }_r^2{E}_{P_N}[W;{\theta }_0+r(\theta -{\theta }_0),{\eta }_0+r(\eta -{\eta }_0)]dr‖\le {\lambda }_N^{\prime }\sqrt{n},}$$

  and by Step 3 below, with $P_N$‐probability $1-o(1)$⁠,  

  $$\begin{array}{ccc}\hfill {\mathcal{I}}_5& :=& \underset{\parallel \theta -{\theta }_0\parallel \le {\tau }_N}{sup}\parallel {\mathbb{G}}_n(\psi (W;\theta ,{\widehat{\eta }}_0)-\psi (W;{\theta }_0,{\eta }_0))\parallel \hfill \\ & \le & r_N^{\prime }{log}^{1/2}(1/r_N^{\prime })+n^{-1/2+1/q}logn.\hfill \end{array}$$

  A.36

  Therefore, because all singular values of J₀ are bounded below from zero by Assumption 3.3(d), it follows that  

  $$\begin{array}{ccc}& & \parallel J_0^{-1}\sqrt{n}{E}_n[\psi (W;{\theta }_0,{\eta }_0)]+\sqrt{n}({\check{\theta }}_0-{\theta }_0)\parallel \hfill \\ & & =O_{P_N}(n^{-1/2+1/q}logn+r_N^{\prime }{log}^{1/2}(1/r_N^{\prime })+({\epsilon }_N+{\lambda }_N+{\lambda }_N^{\prime })\sqrt{n}.\hfill \end{array}$$

  The asserted claim now follows by multiplying both parts of the display by ${\mathrm{\Sigma }}_0^{-1/2}$ (under the norm on the left‐hand side) and noting that singular values of Σ₀ are bounded below from zero by Assumptions 3.3 and 3.4.
• Step 3. Here we derive a bound on ${\mathcal{I}}_5$ in 263. We have  

  $${\mathcal{I}}_5\lesssim \underset{f\in {\mathcal{F}}_2}{sup}|{\mathbb{G}}_n(f)|,{\mathcal{F}}_2=\{{\psi }_j(·,\theta ,{\widehat{\eta }}_0)-{\psi }_j(·,{\theta }_0,{\eta }_0):j=1,\dots ,d_{\theta },\parallel \theta -{\theta }_0\parallel \le {\tau }_n\}.$$

  To bound ${sup}_{f\in {\mathcal{F}}_2}|{\mathbb{G}}_n(f)|$⁠, we apply Lemma 6.2 conditional on ${(W_i)}_{i\in I^c}$ and $I^c$ so that ${\widehat{\eta }}_0$ can be treated as fixed. Observe that with $P_N$‐probability $1-o(1)$⁠, ${sup}_{f\in {\mathcal{F}}_2}{\parallel f\parallel }_{P_N,2}\lesssim r_N^{\prime }$ where we used Assumption 3.4. Thus, an application of Lemma 6.2 to the empirical process $\{{\mathbb{G}}_n(f),f\in {\mathcal{F}}_2\}$ with an envelope $F_2=F_{1,{\widehat{\eta }}_0}+F_{1,{\eta }_0}$ and $\sigma =C{r}_N^{\prime }$ for sufficiently large constant C conditional on ${(W_i)}_{i\in I^c}$ and $I^c$ yields that with $P_N$‐probability $1-o(1)$⁠,  

  $$\begin{array}{c}\hfill \underset{f\in {\mathcal{F}}_2}{sup}|{\mathbb{G}}_n(f)|\lesssim r_N^{\prime }{log}^{1/2}(1/r_N^{\prime })+n^{-1/2+1/q}logn.\end{array}$$

  A.37

  This follows because $\parallel F_2{\parallel }_{P,q}={\parallel F_{1,{\widehat{\eta }}_0}+F_{1,{\eta }_0}\parallel }_{P,q}\le 2C_1$ by Assumption 3.4(b) and the triangle inequality, and  

  $$log\underset{Q}{sup}N(\epsilon \parallel F_2{\parallel }_{Q,2},{\mathcal{F}}_2{,\parallel ·\parallel }_{Q,2})\le 2vlog(2a/\epsilon ),\text{for}\text{all}0<\epsilon \le 1,$$

  because ${\mathcal{F}}_2\subset {\mathcal{F}}_{1,{\widehat{\eta }}_0}-{\mathcal{F}}_{1,{\eta }_0}$ for ${\mathcal{F}}_{1,\eta }$ defined in Assumption 3.4(b), and  

  $$\begin{array}{ccc}& & log\underset{Q}{sup}N(\epsilon \parallel F_{1,{\widehat{\eta }}_0}+F_{1,{\eta }_0}{\parallel }_{Q,2},{\mathcal{F}}_{1,{\widehat{\eta }}_0}-{\mathcal{F}}_{1,{\eta }_0}{,\parallel ·\parallel }_{Q,2})\hfill \\ & & \le log\underset{Q}{sup}N((\epsilon /2)\parallel F_{1,{\widehat{\eta }}_0}{\parallel }_{Q,2},{\mathcal{F}}_{1,{\widehat{\eta }}_0}{,\parallel ·\parallel }_{Q,2})\hfill \\ & & +log\underset{Q}{sup}N((\epsilon /2)\parallel F_{1,{\eta }_0}{\parallel }_{Q,2},{\mathcal{F}}_{1,{\eta }_0}{,\parallel ·\parallel }_{Q,2})\hfill \end{array}$$

  by the proof of Theorem 3 in Andrews (1994b). The claim of this step follows.
• Step 4. Here we derive a bound on ${\mathcal{I}}_3$ in 258. Let ${\overline{\theta }}_0={\theta }_0-J_0^{-1}{E}_n[\psi (W;{\theta }_0,{\eta }_0)]$⁠. Then $\parallel {\overline{\theta }}_0-{\theta }_0\parallel =O_{P_N}(1/\sqrt{n})=o_{P_N}({\tau }_n)$ as $E_{P_N}[\parallel \sqrt{n}{E}_n[\psi (W;{\theta }_0,{\eta }_0)]\parallel ]$ is bounded and the singular values of J₀ are bounded below from zero by Assumption 3.3(d). Therefore, ${\overline{\theta }}_0\in \mathrm{\Theta }$ with $P_N$‐probability $1-o(1)$ by Assumption 3.3(a). Hence, with the same probability,  

  $$\underset{\theta \in \mathrm{\Theta }}{inf}\sqrt{n}\parallel E_n[\psi (W;\theta ,{\widehat{\eta }}_0)]\parallel \le \sqrt{n}\parallel E_n[\psi (W;{\overline{\theta }}_0,{\widehat{\eta }}_0)]\parallel ,$$

  and so it suffices to show that with $P_N$‐probability $1-o(1)$⁠,  

  $$\sqrt{n}\parallel E_n[\psi (W;{\overline{\theta }}_0,{\widehat{\eta }}_0)]\parallel =O(n^{-1/2+1/q}logn+r_N^{\prime }{log}^{1/2}(1/r_N^{\prime })+{\lambda }_N\sqrt{n}+{\lambda }_N^{\prime }\sqrt{n}).$$

  To prove it, substitute $\theta ={\overline{\theta }}_0$ and $\eta ={\widehat{\eta }}_0$ into 259 and use the Taylor expansion in 260. This shows that with $P_N$‐probability $1-o(1)$⁠,  

  $$\begin{array}{ccc}\hfill \sqrt{n}\parallel E_n[\psi (W;{\overline{\theta }}_0,{\widehat{\eta }}_0)]\parallel & \le & \sqrt{n}\parallel E_n[\psi (W;{\theta }_0,{\eta }_0)]+J_0({\overline{\theta }}_0-{\theta }_0)\parallel +{\lambda }_N\sqrt{n}+{\mathcal{I}}_4+{\mathcal{I}}_5\hfill \\ & =& {\lambda }_N\sqrt{n}+{\mathcal{I}}_4+{\mathcal{I}}_5.\hfill \end{array}$$

  Combining this with the bounds on ${\mathcal{I}}_4$ and ${\mathcal{I}}_5$ derived above gives the claim of this step and completes the proof of the lemma. $\square$

 

Proof of Theorems 4.1 and 4.2:

Because Theorem 4.1 is a special case of Theorem 4.2 (with $Z=D$⁠), it suffices to prove the latter. Also, we only consider the DML estimators based on the score 4.7 and note that the estimators based on the score 4.8 can be treated similarly.

Observe that the score ψ in 4.7 is linear in θ:  

$$\begin{array}{ccc}\hfill \psi (W;\theta ,\eta )& =& (Y-D\theta -g(X))(Z-m(X))={\psi }^a(W;\eta )\theta +{\psi }^b(W;\eta );\hfill \\ \hfill {\psi }^a(W;\eta )& =& D(m(X)-Z),{\psi }^b(W;\eta )=(Y-g(X))(Z-m(X)).\hfill \end{array}$$

Therefore, all asserted claims of Theorem 4.2 follow from Theorems 3.1 and 3.2 and Corollary 3.1 as long as we can verify Assumptions 3.1 and 3.2, which we do here. We do so with ${\mathcal{T}}_N$ being the set of all $\eta =(g,m)$ consisting of P‐square‐integrable functions g and m such that  

$$\begin{array}{ccc}\hfill \parallel \eta -{\eta }_0{\parallel }_{P,q}& \le & C,\hfill \\ \hfill \parallel \eta -{\eta }_0{\parallel }_{P,2}& \le & {\delta }_N,\hfill \\ \hfill \parallel m-m_0{\parallel }_{P,2}\times {\parallel g-g_0\parallel }_{P,2}& \le & {\delta }_N{N}^{-1/2}.\hfill \end{array}$$

Also, we replace the constant q and the sequence ${({\delta }_N)}_{N\ge 1}$ in Assumptions 3.1 and 3.2 by $q/2$ and ${({\delta }_N^{\prime })}_{N\ge 1}$ with ${\delta }_N^{\prime }=(C+2\sqrt{C}+2)({\delta }_N\vee N^{-((1-4/q)\wedge (1/2))})$ for all N (recall that we assume that $q>4$⁠, and the analysis in Section 3 only requires that $q>2$⁠; also, ${\delta }_N^{\prime }$ satisfies ${\delta }_N^{\prime }\ge N^{-1/2}$⁠, which is required in Theorems 3.1 and 3.2). We proceed in five steps. All bounds in the proof hold uniformly over $P\in \mathcal{P}$ but we omit this qualifier for brevity).

• Step 1. We first verify Neyman orthogonality. We have that $E_P[\psi (W;{\theta }_0,{\eta }_0)]=0$ by definition of θ₀ of η₀. Also, for any $\eta =(g,m)\in {\mathcal{T}}_N$⁠, the Gateaux derivative in the direction $\eta -{\eta }_0=(g-g_0,m-m_0)$ is given by  

  $$\begin{array}{ccc}\hfill {\partial }_{\eta }{E}_P[\psi (W;{\theta }_0,{\eta }_0)][\eta -{\eta }_0]& =& E_P[(g(X)-g_0(X))(m_0(X)-Z)]\hfill \\ & & +E_P[(m_0(X)-m(X))(Y-D{\theta }_0-g_0(X))]=0,\hfill \end{array}$$

  by the law of iterated expectations, as $V=Z-m_0(X)$ and $U=(Y-D{\theta }_0-g_0(X))$ obey $E_P[V|X]=0$ and $E_P[U|Z,X]=0$⁠. This gives Assumption 3.1(d) with ${\lambda }_N=0$⁠.
• Step 2. Note that  

  $$|J_0|=|E_P[{\psi }^a(W;{\eta }_0)]|=|E_P[D(m_0(X)-Z)]|=|E_P[DV]|\ge c>0$$

  by Assumption 4.2(c). In addition,  

  $$\begin{array}{ccc}\hfill |E_P[{\psi }^a(W;{\eta }_0)]|& =& |E_P[D(m_0(X)-Z)]{|\le \parallel D\parallel }_{P,2}\parallel m_0{(X)\parallel }_{P,2}+{\parallel D\parallel }_{P,2}{\parallel Z\parallel }_{P,2}\hfill \\ & \le & {2\parallel D\parallel }_{P,2}{\parallel Z\parallel }_{P,2}\le {2\parallel D\parallel }_{P,q}{\parallel Z\parallel }_{P,q}\le 2C^2\hfill \end{array}$$

  by the triangle inequality, the Hölder inequality, the Jensen inequality and Assumption 4.2(b). This gives Assumption 3.1(e). Hence, given that Assumptions 3.1(a)–(c) hold trivially, Steps 1 and 2 together show that all conditions of Assumption 3.1 hold.
• Step 3. Note that Assumption 3.2(a) holds by construction of the set ${\mathcal{T}}_N$ and Assumption 4.2(e). Also, note that $\psi (W;{\theta }_0,{\eta }_0)=UV$⁠, and so  

  $$E_P[\psi (W;{\theta }_0,{\eta }_0)\psi {(W;{\theta }_0,{\eta }_0)}^{\prime }]=E_P[U^2{V}^2]\ge c^4>0,$$

  by Assumption 4.2(c), which gives Assumption 3.2(d).
• Step 4. Here we verify Assumption 3.2(b). For any $\eta =(g,m)\in {\mathcal{T}}_N$⁠, we have  

  $$\begin{array}{ccc}\hfill (E_P[\parallel {\psi }^a{(W;\eta )\parallel }^{q/2}{])}^{2/q}& =& \parallel {\psi }^a{(W;\eta )\parallel }_{P,q/2}={\parallel D(m(X)-Z)\parallel }_{P,q/2}\hfill \\ & \le & \parallel D(m(X)-m_0(X)){\parallel }_{P,q/2}+\parallel D{m}_0{(X)\parallel }_{P,q/2}+{\parallel DZ\parallel }_{P,q/2}\hfill \\ & \le & {\parallel D\parallel }_{P,q}\parallel m(X)-m_0{(X)\parallel }_{P,q}+{\parallel D\parallel }_{P,q}\parallel m_0{(X)\parallel }_{P,q}+{\parallel D\parallel }_{P,q}{\parallel Z\parallel }_{P,q}\hfill \\ & \le & {C\parallel D\parallel }_{P,q}+{2\parallel D\parallel }_{P,q}{\parallel Z\parallel }_{P,q}\le 3C^2\hfill \end{array}$$

  by Assumption 4.2(b), which gives the bound on $m_N^{\prime }$ in Assumption 3.2(b). Also, because  

  $$|E_P[(D-r_0(X))(Z-m_0(X))]|=|E_P[DV]|\ge c$$

  by Assumption 4.2(c), it follows that θ₀ satisfies  

  $$\begin{array}{ccc}\hfill |{\theta }_0|& =& \frac{|E_P[(Y-{\ell }_0(X))(Z-m_0(X))]|}{|E_P[(D-r_0(X))(Z-m_0(X))]|}\hfill \\ & \le & c^{-1}{(\parallel Y\parallel }_{P,2}+\parallel {\ell }_0{(X)\parallel }_{P,2}{)(\parallel Z\parallel }_{P,2}+\parallel m_0(X){\parallel }_{P,2})\hfill \\ & \le & 4c^{-1}{\parallel Y\parallel }_{P,2}{\parallel Z\parallel }_{P,2}\le 4C^2/c.\hfill \end{array}$$

  Hence,  

  $$\begin{array}{ccc}& & (E_P[\parallel \psi (W;{\theta }_0,\eta ){\parallel }^{q/2}{])}^{2/q}={\parallel \psi (W;{\theta }_0,\eta )\parallel }_{P,q/2}\hfill \\ & & =\parallel (Y-D{\theta }_0-g(X))(Z-m(X)){\parallel }_{P,q/2}\hfill \\ & & \le {\parallel U(Z-m(X))\parallel }_{P,q/2}+{\parallel (g(X)-g_0(X))(Z-m(X))\parallel }_{P,q/2}\hfill \\ & & \le {\parallel U\parallel }_{P,q}{\parallel Z-m(X)\parallel }_{P,q}+\parallel g(X)-g_0{(X)\parallel }_{P,q}{\parallel Z-m(X)\parallel }_{P,q}\hfill \\ & & \le {(\parallel U\parallel }_{P,q}{+C)\parallel Z-m(X)\parallel }_{P,q}\hfill \\ & & \le (\parallel Y-D{\theta }_0{\parallel }_{P,q}+\parallel g_0(X){\parallel }_{P,q}+C)\hfill \\ & & \times (\parallel Z{\parallel }_{P,q}+\parallel m_0(X){\parallel }_{P,q}+\parallel m(X)-m_0(X){\parallel }_{P,q})\hfill \\ & & \le (2\parallel Y-D{\theta }_0{\parallel }_{P,q}+{C)(2\parallel Z\parallel }_{P,q}+C)\hfill \\ & & \le {(2\parallel Y\parallel }_{P,q}+{2\parallel D\parallel }_{P,q}|{\theta }_0{|+C)(2\parallel Z\parallel }_{P,q}+C)\hfill \\ & & \le 3C(3C+8C^3/c),\hfill \end{array}$$

  where we used the fact that because $g_0(X)=E_P[Y-D{\theta }_0\mid X]$⁠, $\parallel g_0{(X)\parallel }_{P,q}\le {\parallel Y-D{\theta }_0\parallel }_{P,q}$ by the Jensen inequality. This gives the bound on $m_N$ in Assumption 3.2(b). Hence, Assumption 3.2(b) holds.
• Step 5. Finally, we verify Assumption 3.2(c). For any $\eta =(g,m)\in {\mathcal{T}}_N$⁠, we have  

  $$\begin{array}{ccc}\hfill \parallel E_P[{\psi }^a(W;\eta )]-E_P[{\psi }^a(W;{\eta }_0)]\parallel & =& |E_P[{\psi }^a(W;\eta )-{\psi }^a(W;{\eta }_0)]|\hfill \\ & =& |E_P[D(m(X)-m_0(X))]|\hfill \\ & \le & {\parallel D\parallel }_{P,2}{\parallel m(X)-m_0(X)\parallel }_{P,2}\le C{\delta }_N\le {\delta }_N^{\prime },\hfill \end{array}$$

  which gives the bound on $r_N$ in Assumption 3.2(c). Further,  

  $$\begin{array}{ccc}& & (E_P[\parallel \psi (W;{\theta }_0,\eta )-\psi (W;{\theta }_0,{\eta }_0){\parallel }^2{])}^{1/2}={\parallel \psi (W;{\theta }_0,\eta )-\psi (W;{\theta }_0,{\eta }_0)\parallel }_{P,2}\hfill \\ & & =\parallel (U+g_0(X)-g(X))(Z-m(X))-U(Z-m_0(X)){\parallel }_{P,2}\hfill \\ & & \le \parallel U(m(X)-m_0(X)){\parallel }_{P,2}+{\parallel (g(X)-g_0(X))(Z-m(X))\parallel }_{P,2}\hfill \\ & & \le \sqrt{C}\parallel m(X)-m_0{(X)\parallel }_{P,2}+{\parallel V(g(X)-g_0(X))\parallel }_{P,2}\hfill \\ & & +\parallel (g(X)-g_0(X))(m(X)-m_0(X)){\parallel }_{P,2}\hfill \\ & & \le \sqrt{C}\parallel m(X)-m_0{(X)\parallel }_{P,2}+\sqrt{C}{\parallel g(X)-g_0(X)\parallel }_{P,2}\hfill \\ & & +\parallel g(X)-g_0{(X)\parallel }_{P,2}\times {\parallel m(X)-m_0(X)\parallel }_{P,2}\hfill \\ & & \le (2\sqrt{C}+N^{-1/2}){\delta }_N\le (2\sqrt{C}+1){\delta }_N\le {\delta }_N^{\prime },\hfill \end{array}$$

  which gives the bound on $r_N^{\prime }$ in Assumption 3.2(c). Finally, let  

  $$f(r):=E_P[\psi (W;{\theta }_0,{\eta }_0+r(\eta -{\eta }_0)],r\in (0,1).$$

  Then, for any $r\in (0,1)$⁠,  

  $$f(r)=E_P[(U-r(g(X)-g_0(X)))(V-r(m(X)-m_0(X)))],$$

  and so  

  $$\begin{array}{ccc}\hfill \partial f(r)& =& -E_P[(g(X)-g_0(X))(V-r(m(X)-m_0(X)))]\hfill \\ & & -E_P[(U-r(g(X)-g_0(X)))(m(X)-m_0(X))],\hfill \\ \hfill {\partial }^{2}f(r)& =& 2E_P[(g(X)-g_0(X))(m(X)-m_0(X))].\hfill \end{array}$$

  Hence,  

  $$|{\partial }^{2}f(r)|\le 2\parallel g(X)-g_0{(X)\parallel }_{P,2}\times {\parallel m(X)-m_0(X)\parallel }_{P,2}\le 2{\delta }_N{N}^{-1/2}\le {\delta }_N^{\prime }{N}^{-1/2},$$

  which gives the bound on ${\lambda }_N^{\prime }$ in Assumption 3.2(c). Thus, all conditions of Assumptions 3.1 are verified. This completes the proof. $\square$

 

Proof of Theorems 5.1 and 5.2:

The proof of Theorem 5.2 is similar to that of Theorem 5.1 and is therefore omitted. In turn, regarding Theorem 5.1, we show the proof for the case of ATE and note that the proof for the case of ATTE is similar.

Observe that the score ψ in 5.3 is linear in θ:  

$$\begin{array}{ccc}\hfill \psi (W;\theta ,\eta )& =& {\psi }^a(W;\eta )\theta +{\psi }^b(W;\eta ),{\psi }^a(W;\eta )=-1,\hfill \\ \hfill {\psi }^b(W;\eta )& =& (g(1,X)-g(0,X))+\frac{D(Y-g(1,X))}{m(X)}-\frac{(1-D)(Y-g(0,X))}{1-m(X)}.\hfill \end{array}$$

Therefore, all asserted claims of Theorem 5.1 follow from Theorems 3.1 and 3.2 and Corollary 3.1 as long as we can verify Assumptions 3.1 and 3.2, which we do here. We do so with ${\mathcal{T}}_N$ being the set of all $\eta =(g,m)$ consisting of P‐square‐integrable functions g and m such that  

$$\begin{array}{ccc}\hfill \parallel \eta -{\eta }_0{\parallel }_{P,q}& \le & C,\hfill \\ \hfill \parallel \eta -{\eta }_0{\parallel }_{P,2}& \le & {\delta }_N,\hfill \\ \hfill {\parallel m-1/2\parallel }_{P,\infty }& \le & 1/2-\epsilon ;\hfill \\ \hfill \parallel m-m_0{\parallel }_{P,2}\times {\parallel g-g_0\parallel }_{P,2}& \le & {\delta }_N{N}^{-1/2}.\hfill \end{array}$$

Also, we replace the sequence ${({\delta }_N)}_{N\ge 1}$ in Assumptions 3.1 and 3.2 by ${({\delta }_N^{\prime })}_{N\ge 1}$ with ${\delta }_N^{\prime }=C_{\epsilon }({\delta }_N\vee N^{-((1-4/q)\wedge (1/2))})$ for all N, where $C_{\epsilon }$ is a sufficiently large constant that depends only on ε and C (note that ${\delta }_N^{\prime }$ satisfies ${\delta }_N^{\prime }\ge N^{-((1-4/q)\wedge (1/2))}$⁠, which is required in Theorems 3.1 and 3.2). We proceed in five steps. All bounds in the proof hold uniformly over $P\in \mathcal{P}$ but we omit this qualifier for brevity.

• Step 1. We first verify Neyman orthogonality. We have that $E[\psi (W;{\theta }_0,{\eta }_0)]=0$ by definition of θ₀ and η₀. Also, for any $\eta =(g,m)\in {\mathcal{T}}_N$⁠, the Gateaux derivative in the direction $\eta -{\eta }_0=(g-g_0,m-m_0)$ is given by  

  $$\begin{array}{ccc}& & {\partial }_{\eta }{E}_P[\psi (W;{\theta }_0,{\eta }_0)][\eta -{\eta }_0]=E_P[g(1,X)-g_0(1,X)]-E_P[g(0,X)-g_0(0,X)]\hfill \\ & & -E_P\left[\frac{D(g(1,X)-g_0(1,X))}{m_0(X)}\right]+E_P\left[\frac{(1-D)(g(0,X)-g_0(0,X))}{1-m_0(X)}\right]\hfill \\ & & -E_P\left[\frac{D(Y-g_0(1,X))(m(X)-m_0(X))}{m_0^2(X)}\right]\hfill \\ & & -E_P\left[\frac{(1-D)(Y-g_0(0,X))(m(X)-m_0(X))}{{(1-m_0(X))}^2}\right],\hfill \end{array}$$

  which is 0 by the law of iterated expectations, as  

  $$\begin{array}{ccc}\hfill E_P[D\mid X]& =& m_0(X),E_P[1-D\mid X]=1-m_0(X),\hfill \\ \hfill E_P[D(Y-g_0(1,X))\mid X]& =& 0,E_P[(1-D)(Y-g_0(0,X))\mid X]=0.\hfill \end{array}$$

  This gives Assumption 3.1(d) with ${\lambda }_N=0$⁠.

• Step 2. Note that $J_0=-1$⁠, and so Assumption 3.1(e) holds trivially. Hence, given that Assumptions 3.1(a)–(c) hold trivially as well, Steps 1 and 2 together show that all conditions of Assumption 3.1 hold.

• Step 3. Note that Assumption 3.2(a) holds by construction of the set ${\mathcal{T}}_N$ and Assumption 5.1(f). Also,  

  $$\begin{array}{ccc}\hfill E_P[{\psi }^2(W;{\theta }_0,{\eta }_0)]& =& E_P[E_P[{\psi }^2(W;{\theta }_0,{\eta }_0)\mid X]]\hfill \\ & =& E_P\left[E_P[{(g_0(1,X)-g_0(0,X)-{\theta }_0)}^2\mid X]\right.\hfill \\ & & +\left.E_P\left[{\left(\frac{D(Y-g_0(1,X))}{m_0(X)}-\frac{(1-D)(Y-g_0(0,X))}{1-m_0(X)}\right)}^2\mid X\right]\right]\hfill \\ & \ge & E_P\left[{\left(\frac{D(Y-g_0(1,X))}{m_0(X)}-\frac{(1-D)(Y-g_0(0,X))}{1-m_0(X)}\right)}^2\right]\hfill \\ & =& E_P\left[\frac{D^2{(Y-g_0(1,X))}^2}{m_0{(X)}^2}+\frac{{(1-D)}^2{(Y-g_0(0,X))}^2}{{(1-m_0(X))}^2}\right]\hfill \\ & \ge & \frac{1}{{(1-\epsilon )}^2}{E}_P[D^2{(Y-g_0(1,X))}^2+{(1-D)}^2{(Y-g_0(0,X))}^2]\hfill \\ & =& \frac{1}{{(1-\epsilon )}^2}{E}_P[D{U}^2+(1-D)U^2]=\frac{1}{{(1-\epsilon )}^2}{E}_P[U^2]\ge \frac{c^2}{{(1-\epsilon )}^2}.\hfill \end{array}$$

  This gives Assumption 3.2(d).
• Step 4. Here we verify Assumption 3.2(b). We have  

  $$\begin{array}{ccc}\hfill \parallel g_0{(D,X)\parallel }_{P,q}& =& (E_P[|g_0{(D,X)|}^q{])}^{1/q}\hfill \\ & \ge & (E_P[|g_0{(1,X)|}^q{Pr}_P(D=1\mid X)+|g_0{(0,X)|}^q{Pr}_P(D=0\mid X){])}^{1/q}\hfill \\ & \ge & {\epsilon }^{1/q}(E_P[|g_0{(1,X)|}^q]+E_P[|g_0{(0,X)|}^q{])}^{1/q}\hfill \\ & \ge & {\epsilon }^{1/q}(E_P[|g_0{(1,X)|}^q]\vee E_P[|g_0{(0,X)|}^q{])}^{1/q}\hfill \\ & \ge & {\epsilon }^{1/q}(\parallel g_0{(1,X)\parallel }_{P,q}\vee \parallel g_0(0,X){\parallel }_{P,q}),\hfill \end{array}$$

  where in the third line, we used the facts that ${\mathit{Pr}}_P(D=1\mid X)=m_0(X)\ge \epsilon$ and ${Pr}_P(D=0\mid X)=1-m_0(X)\ge \epsilon$⁠. Hence, given that $\parallel g_0{(D,X)\parallel }_{P,q}\le {\parallel Y\parallel }_{P,q}\le C$ by the Jensen inequality and Assumption 5.1(b), it follows that  

  $$\parallel g_0{(1,X)\parallel }_{P,q}\le C/{\epsilon }^{1/q}\text{and}{\parallel g_0(0,X)\parallel }_{P,q}\le C/{\epsilon }^{1/q}.$$

  Similarly, for any $\eta \in (g,m)\in {\mathcal{T}}_N$⁠,  

  $$\parallel g(1,X)-g_0{(1,X)\parallel }_{P,q}\le C/{\epsilon }^{1/q}\text{and}{\parallel g(0,X)-g_0(0,X)\parallel }_{P,q}\le C/{\epsilon }^{1/q}$$

  as $\parallel g(D,X)-g_0{(D,X)\parallel }_{P,q}\le C$⁠. In addition,  

  $$|{\theta }_0|=|E_P[g_0(1,X)-g_0(0,X)]|\le \parallel g_0{(1,X)\parallel }_{P,2}+{\parallel g_0(0,X)\parallel }_{P,2}\le 2C/{\epsilon }^{1/q}.$$

  Therefore, for any $\eta =(g,m)\in {\mathcal{T}}_N$⁠, we have  

  $$\begin{array}{ccc}& & (E_P[|\psi (W;{\theta }_0,\eta ){|}^q{])}^{1/q}={\parallel \psi (W;{\theta }_0,\eta )\parallel }_{P,q}\hfill \\ & & \le (1+{\epsilon }^{-1})\left({\parallel g(1,X)\parallel }_{P,q}+{\parallel g(0,X)\parallel }_{P,q}\right)+{2\parallel Y\parallel }_{P,q}/\epsilon +|{\theta }_0|\hfill \\ & & \le (1+{\epsilon }^{-1})\left(\parallel g(1,X)-g_0{(1,X)\parallel }_{P,q}+{\parallel g(0,X)-g_0(0,X)\parallel }_{P,q}\right)\hfill \\ & & +(1+{\epsilon }^{-1})\left(\parallel g_0{(1,X)\parallel }_{P,q}+{\parallel g_0(0,X)\parallel }_{P,q}\right)+2C/\epsilon +2C/{\epsilon }^{1/q}\hfill \\ & & \le 4C(1+{\epsilon }^{-1})/{\epsilon }^{1/q}+2C/\epsilon +2C/{\epsilon }^{1/q}.\hfill \end{array}$$

  This gives the bound on $m_N$ in Assumption 3.2(b). Also, we have  

  $$(E_P[|{\psi }^a{(W;\eta )|}^q{])}^{1/q}=1.$$

  This gives the bound on $m_N^{\prime }$ in Assumption 3.2(b). Hence, Assumption 3.2(b) holds.
• Step 5. Finally, we verify Assumption 3.2(c). For any $\eta =(g,m)\in {\mathcal{T}}_N$⁠, we have  

  $$\parallel E_P[{\psi }^a(W;\eta )-{\psi }^a(W;{\eta }_0)]\parallel =|1-1|=0\le {\delta }_N^{\prime },$$

  which gives the bound on $r_N$ in Assumption 3.2(c). Further, by the triangle inequality,  

  $$\begin{array}{ccc}\hfill (E_P[\parallel \psi (W;{\theta }_0,\eta )-\psi (W;{\theta }_0,{\eta }_0){\parallel }^2{])}^{1/2}& =& \parallel \psi (W;{\theta }_0,\eta )-\psi (W;{\theta }_0;{\eta }_0){\parallel }_{P,2}\hfill \\ & \le & {\mathcal{I}}_1+{\mathcal{I}}_2+{\mathcal{I}}_3,\hfill \end{array}$$

  where  

  $$\begin{array}{ccc}\hfill {\mathcal{I}}_1& :=& \parallel g(1,X)-g_0{(1,X)\parallel }_{P,2}+{\parallel g(0,X)-g_0(0,X)\parallel }_{P,2},\hfill \\ \hfill {\mathcal{I}}_2& :=& {‖\frac{D(Y-g(1,X))}{m(X)}-\frac{D(Y-g_0(1,X))}{m_0(X)}‖}_{P,2},\hfill \\ \hfill {\mathcal{I}}_3& :=& {‖\frac{(1-D)(Y-g(0,X))}{1-m(X)}-\frac{(1-D)(Y-g_0(0,X))}{1-m_0(X)}‖}_{P,2}.\hfill \end{array}$$

  To bound ${\mathcal{I}}_1$⁠, note that by the same argument as that used in Step 4,  

  $$\parallel g(1,X)-g_0{(1,X)\parallel }_{P,2}\le {\delta }_N/{\epsilon }^{1/2}\text{and}{\parallel g(0,X)-g_0(0,X)\parallel }_{P,2}\le {\delta }_N/{\epsilon }^{1/2},$$

  A.38

  and so ${\mathcal{I}}_1\le 2{\delta }_N/{\epsilon }^{1/2}$⁠. To bound ${\mathcal{I}}_2$⁠, we have  

  $$\begin{array}{ccc}\hfill {\mathcal{I}}_2& \le & {\epsilon }^{-2}{\parallel D{m}_0(X)(Y-g(1,X))-Dm(X)(Y-g_0(1,X))\parallel }_{P,2}\hfill \\ & \le & {\epsilon }^{-2}{\parallel m_0(X)(g_0(1,X)+U-g(1,X))-m(X)U\parallel }_{P,2}\hfill \\ & \le & {\epsilon }^{-2}(\parallel m_0(X)(g(1,X)-g_0(1,X)){\parallel }_{P,2}+\parallel (m(X)-m_0(X))U{\parallel }_{P,2})\hfill \\ & \le & {\epsilon }^{-2}(\parallel g(1,X)-g_0{(1,X)\parallel }_{P,2}+\sqrt{C}\parallel m(X)-m_0(X){\parallel }_{P,2})\hfill \\ & \le & {\epsilon }^{-2}({\epsilon }^{-1/2}+\sqrt{C}){\delta }_N.\hfill \end{array}$$

  Here, the first inequality follows from the bounds $\epsilon \le m_0(X)\le 1-\epsilon$ and $\epsilon \le m(X)\le 1-\epsilon$⁠, the second from the facts that $D\in \{0,1\}$ and for $D=1$⁠, $Y=g_0(1,X)+U$⁠, the third from the triangle inequality, the fourth from the facts that $m_0(X)\le 1$ and $E_P[U^2\mid X]\le C$⁠, and the fifth from 302. Similarly, ${\mathcal{I}}_3\le {\epsilon }^{-2}({\epsilon }^{-1/2}+\sqrt{C}){\delta }_N$⁠. Combining these inequalities shows that  

  $$(E_P[\parallel \psi (W;{\theta }_0,\eta )-\psi (W;{\theta }_0,{\eta }_0){\parallel }^2{])}^{1/2}\le 2({\epsilon }^{-1/2}+{\epsilon }^{-5/2}+\sqrt{C}{\epsilon }^{-2}){\delta }_N\le {\delta }_N^{\prime },$$

  as long as $C_{\epsilon }$ in the definition of ${\delta }_N^{\prime }$ satisfies $C_{\epsilon }\ge 2({\epsilon }^{-1/2}+{\epsilon }^{-5/2}+\sqrt{C}{\epsilon }^{-2})$⁠. This gives the bound on $r_N^{\prime }$ in Assumption 3.2(c). Finally, let  

  $$f(r):=E_P[\psi (W;{\theta }_0,{\eta }_0+r(\eta -{\eta }_0))],r\in (0,1).$$

  Then for any $r\in (0,1)$⁠,  

  $$\begin{array}{ccc}\hfill {\partial }^{2}f(r)& =& E_P\left[\frac{D(g(1,X)-g_0(1,X))(m(X)-m_0(X))}{{(m_0(X)+r(m(X)-m_0(X)))}^2}\right]\hfill \\ & & +E_P\left[\frac{(1-D)(g(0,X)-g_0(0,X))(m(X)-m_0(X))}{{(1-m_0(X)-r(m(X)-m_0(X)))}^2}\right]\hfill \\ & & +E_P\left[\frac{(g(1,X)-g_0(1,X))(m(X)-m_0(X))}{{(m_0(X)+r(m(X)-m_0(X)))}^2}\right]\hfill \\ & & +2E_P\left[\frac{D(Y-g_0(1,X)-r(g(1,X)-g_0(1,X))){(m(X)-m_0(X))}^2}{{(m_0(X)+r(m(X)-m_0(X)))}^3}\right]\hfill \\ & & +E_P\left[\frac{(g(0,X)-g_0(0,X))(m(X)-m_0(X))}{{(1-m_0(X)-r(m(X)-m_0(X)))}^2}\right]\hfill \\ & & -2E_P\left[\frac{(1-D)(Y-g_0(0,X)-r(g(0,X)-g_0(0,X))){(m(X)-m_0(X))}^2}{{(1-m_0(X)-r(m(X)-m_0(X)))}^3}\right],\hfill \end{array}$$

  and so, given that  

  $$\begin{array}{ccc}\hfill D(Y-g_0(1,X))& =& DU,(1-D)(Y-g_0(0,X))=(1-D)U,\hfill \\ \hfill E_P[U\mid D,X]& =& 0,|m(X)-m_0(X)|\le 2,\hfill \end{array}$$

  it follows that for some constant $C_{\epsilon }^{\prime }$ that depends only on ε and C,  

  $$|{\partial }^{2}f(r)|\le C_{\epsilon }^{\prime }\parallel m-m_0{\parallel }_{P,2}\times {\parallel g-g_0\parallel }_{P,2}\le {\delta }_N^{\prime }{N}^{-1/2},$$

  as long as the constant $C_{\epsilon }$ in the definition of ${\delta }_N^{\prime }$ satisfies $C_{\epsilon }\ge C_{\epsilon }^{\prime }$⁠. This gives the bound on ${\lambda }_N^{\prime }$ in Assumption 3.2(c). Thus, all conditions of Assumptions 3.1 are verified. This completes the proof. $\square$