Double Debiased Machine Learning Nonparametric Inference with Continuous Treatments

Double Debiased Machine Learning Nonparametric Inference with Continuous Treatments

We propose a doubly robust inference method for causal effects of continuous treatment variables, under unconfoundedness and with nonparametric or high-dimensional nuisance functions. Our double debiased machine learning (DML) estimators for the average dose-response function (or the average structural function) and the partial effects are asymptotically normal with nonparametric convergence rates. The first-step estimators for the nuisance conditional expectation function and the conditional density can be nonparametric or ML methods. Utilizing a kernel-based doubly robust moment function and cross-fitting, we give high-level conditions under which the nuisance function estimators do not affect the first-order large sample distribution of the DML estimators. We provide sufficient low-level conditions for kernel, series, and deep neural networks. We justify the use of kernel to localize the continuous treatment at a given value by the Gateaux derivative. We implement various ML methods in Monte Carlo simulations and an empirical application on a job training program evaluation. \\ Keywords: Average structural function, cross-fitting, dose-response function, doubly robust. \\ JEL Classification: C14, C21, C55

Introduction

We propose a nonparametric inference method for {\it continuous} treatment or structural effects, under the unconfoundedness assumption and in the presence of nonparametric nuisance functions. We focus on the heterogenous effect with respect to the continuous treatment or policy variable $T$. To identify the causal effects, it is plausible to have a large number of the control variables $X$ that $T$ is randomly assigned conditional on. To achieve doubly robust inference and to employ machine learning (ML) methods, we use a double debiased ML approach that combines a doubly robust moment function and cross-fitting.

We consider a nonparametric and model-free outcome equation $Y = g(T, X, \varepsilon)$. No functional form assumption is imposed on the unobserved disturbances $\varepsilon$, such as restrictions on dimensionality, monotonicity, or separability. The potential outcome is $Y(t) = g(t,X,\varepsilon)$ indexed by the hypothetical treatment value $t$. The object of interest is the {\it average dose-response function} as a function of $t$, defined as the expected value of the potential outcome across observations with the observed and unobserved heterogeneity $(X, \varepsilon)$, i.e.\ $\beta_t \equiv \mathbb{E}[Y(t)] = \int\int g(t, X, \varepsilon) dF_{X\varepsilon}$. It is also known as the {\it average structural function} in nonseparable models in BP03. The well-studied average treatment effect of switching from treatment $t$ to $s$ is $\beta_{s} - \beta_{t}$. We allow the continuous treatment to be multi-dimensional and hence capture the unrestricted heterogenous effects with respect to the multivariate treatment variables. We define the {\it partial (or marginal) effect} of the first component of the continuous variable $T$ at $t = (t_1,... t_{d_T})'$ to be the partial derivative $\theta_t \equiv \partial \beta_t/ \partial t_1$. In program evaluation, the average dose response function $\beta_t$ shows how participants' labor market outcomes vary with the length of exposure to a job training program. In demand analysis when $T$ contains price and income, the average structural function $\beta_t$ can be the Engel curve. The partial effect $\theta_t$ reveals the average price elasticity at given values of price and income. Other examples include the efficacy of political advertisements on campaign contributions in Fong, the effect of nurse staffing on hospital readmissions penalties in Kennedy, etc.

We are among the first to apply the double debiased ML approach to inference on the average structural function $\beta_t$ and the partial effect $\theta_t$ of continuous variables, to our knowledge. They are {\it non-regular nonparametric} objects that cannot be estimated at a root-$n$ convergence rate. We propose a kernel-based {\it double debiased machine learning} (DML) estimator that utilizes a doubly robust moment function and cross-fitting via sample-splitting. The DML estimator uses the moment function

where the conditional expectation function $\gamma(t, x) \equiv \mathbb{E}[Y|T=t, X=x]$ and the conditional density $f_{T|X}(t|x)$ is also known as the generalized propensity score (GPS). A kernel $K_h(T_i-t)$ weights observation $i$ with treatment value around $t$ in a distance of $h$. The number of such observations shrinks as the bandwidth $h$ vanishes with the sample size $n$. A {\it $L$-fold cross-fitting} splits the sample into $L$ subsamples. The nuisance function estimators for $\gamma(t, X_i)$ and $f_{T|X}(t|X_i)$ use observations in the other $(L-1)$ subsamples that do not contain the observation $i$. The DML estimator averages over the subsamples. Then we estimate the partial effect $\theta_t$ by a numerical differentiation.

The doubly robust moment function in equation ((ref)) has appeared in Kallus without asymptotic theory and has been extensively studied in SUZ for Lasso-type estimators. We utilize cross-fitting and provide high-level and low-level conditions that facilitate a variety of nonparametric and ML methods. As each ML method has its strengths and weaknesses depending on the data generating process and applications, flexible employment of various nuisance function estimators is desirable. High-dimensional control variables can be accommodated via the nuisance function estimators; for example, Lasso allows the dimension of $X$ to grow with the sample size. Importantly our inference theory for the proposed DML estimator $\hat \beta_t$ allows one of the nuisance functions to be misspecified. The doubly robust inference is useful, especially when the nuisance functions are estimated under a parametric model or sparsity approximation as Lasso.

We show that the proposed kernel-based DML estimators are asymptotically normal and provide high-level conditions under which the nuisance function estimators for $\gamma(t, X_i)$ and $f_{T|X}(t|X_i)$ do not affect the first-order asymptotic distribution. Specifically the high-level conditions on the convergence rates use a {\it partial $L_2$ norm} that fixes the treatment value at $t$, in contrast to the standard $L_2$ norm that integrates over the joint distribution of $(T,X)$, i.e.\ the root-mean-square rate. We further give low-level conditions for nonparametric kernel, series estimators, and the deep neural networks in Farrell18 that is widely popular in industrial applications. These results on the convergence rates of the nuisance function estimators are new to the literature.

Furthermore, we propose a Simulated DML (SDML) estimator that replaces $\gamma(t, X_i)$ in ((ref)) with $\gamma(U_i, X_i)$ where a simulated variable $U_i$ localizes the realized treatment values around the target value $t$. Introducing such a local variation enables the standard $L_2$ convergence rate of the nuisance function estimator that is available for nonparametric kernel and series estimators, as well as recent ML estimators, such as Lasso in BRT09AS, neural networks in ChenWhite, SH20AS, Farrell18, random forests in SZ20, and empirical $L_2$ rate for boosting in LeoSpindler, as discussed in CCDDHNR (CCDDHNR, hereafter) and CNS21ADML.\footnote{ We are grateful to Whitney Newey for the idea of simulating $U_i$ from the probability density function $f_U(u) = f_{T}\left((u-t)/h\right)/h$. } This result is valuable, as the SDML estimator readily allows for a wider class of ML methods.

In addition, we propose generic ML estimators for the {\it reciprocal of the conditional density} $1/f_{T|X}(t|x)$ when $d_T = 1$ and for $f_{T|X}(t|x)$ when $d_T > 1$ respectively, which may be of independent interest. We also propose a data-driven bandwidth to consistently estimate the optimal bandwidth that minimizes the asymptotic mean squared error.

We aim for a tractable inference procedure that is flexible to employ nonparametric or ML nuisance function estimators and delivers a reliable distributional approximation in practice. Toward that end, the DML method contains two key ingredients: a doubly robust moment function and cross-fitting. The doubly robust moment function reduces sensitivity in estimating $\beta_t$ with respect to nuisance parameters.\footnote{ The double robustness usually refers to consistency of the estimator even if either one of the two nuisance functions is misspecified. The rapidly growing ML literature has utilized the double robustness to reduce regularization and modeling biases in estimating the nuisance functions by ML or nonparametric methods; for example, BCH14RES, Farrell15, BCFH17, Farrell18, CEINR, CCDDHNR, RotheFirpo, and references therein. } Cross-fitting further removes bias induced by overfitting and achieves stochastic equicontinuity without strong entropy conditions.\footnote{ CCDDHNR point out that the commonly used results in empirical process theory, such as Donsker properties, could break down in high-dimensional settings. For example, BCFH17 show how cross-fitting weakens the entropy condition and hence the sparsity assumption on nuisance Lasso estimator. The benefit of cross-fitting is further investigated by WagerAthey for heterogeneous causal effects, NeweyRobins for double cross-fitting, and CJ19ET for cross-fitting bootstrap.} \\[5pt] Related literature:\ Our work builds on the results for semiparametric models in IchimuraNewey22QE, CEINR, and CCDDHNR and extends the literature to nonparametric continuous treatment/structural effects. Note that the doubly robust estimator for a binary/multivalued treatment replaces the kernel $K_h(T_i - t)$ with the indicator function ${\bf 1}\{T_i=t\}$ in equation ((ref)) and has been widely studied, especially in the recent ML literature. We show that the advantageous properties of the DML estimator for the binary treatment carry over to the continuous treatments case. Our DML estimator utilizes the kernel function $K_h(T_i - t)$ for the continuous treatments $T$ of fixed low dimension and averages out the covariates $X$. While our kernel-based estimator appears to be a simple modification of the binary treatment case in practice, we show that one important distinct feature of {\it non-regular nonparametric} parameters is that the Gateaux derivative and the Riesz representer are not unique, depending on how we approximate the continuous treatment distribution approaching a point mass. And the kernel function is a natural choice for localization at $t$. Neyman orthogonality holds as $h \rightarrow 0$ (Neyman). Therefore we provide a foundational justification for the proposed kernel-based DML estimator $\hat\beta_t$, relative to alternative approaches, such as Kennedy and SC; see also VBL. Furthermore, our estimator is doubly robust in the sense that our inference theory is valid, i.e.\ the asymptotic distribution is the same, if either one of the nuisance functions $\mathbb{E}[Y|T,X]$ or $f_{T|X}$ is misspecified, as in Kennedy and Westling. This is a stronger result than the usual doubly robustness on the consistency of the estimator. When one nuisance function is misspecified, we require the other nuisance function to be estimated consistently at a convergence rate faster than $\sqrt{nh^{d_T}}$, which is the convergence rate of $\hat\beta_t$. In contrast, the DML estimators in the semiparametric model of CCDDHNR converge at a regular root-$n$ rate, so their inference theory does not allow this doubly robust property.

There is a small yet rapidly growing literature on employing the DML approach for non-regular nonparametric infinite-dimesional objects. HHLL propose a test for monotonicity of $\beta_t$. CNS21EJ, SC, Fan, and ZimmertLechner study the conditional average binary treatment effect $\mathbb{E}[Y(1) - Y(0)|X_1]$ for a low-dimensional subset $X_1$ of $X$. BonviniKennedy use higher-order influence functions to achieve a faster convergence rate. Despite the advantageous theoretical properties, we note potential drawbacks of the DML approach: Double robustness often requires additional nuisance function estimation, such as the conditional density in the Riesz representer, which could introduce additional variation and implementation complication. Cross-fitting could result in a small effective sample size. For example, the Monte Carlo simulations in Fan show that cross-fitting does not improve the finite-sample performance for Lasso-type estimation.

Our paper adds to the literature on continuous treatment effects estimation. In low-dimensional settings, see Imbens00, HI04, Flores07, and Lee for examples of a class of regression estimators $n^{-1}\sum_{i=1}^n \hat \gamma(t, X_i)$. GW and HHLP study a class of inverse probability weighting estimators. The empirical applications in FFGN12ReStat and KSUZ12 focus on semiparametric results. We extend this literature to a DML framework that enables ML methods for nonparametric inference in practice.

A main contribution of this paper is a formal inference theory for the fully nonparametric causal effects of continuous variables. To uncover the causal effect of the continuous variable $T$ on $Y$, our nonparametric nonseparable model $Y = g(T,X,\varepsilon)$ can be compared to the partially linear model $Y = \theta T + g(X) + \varepsilon$ in Robinson that specifies the homogenous effect by $\theta$ and hence is a semiparametric problem. The important partially linear model has many applications and is one of the leading examples in the recent ML literature, where the nuisance function $g(X)$ can be high-dimensional and estimated by a ML method.\footnote{ See CCDDHNR and references therein. DSLC and OSW extend to more general functional forms. CJN18ET, CJN18JASA, CJM19, and FLM2 propose different approaches. } Another semiparametric parameter of interest is the weighted average of $\beta_t$ or $\theta_t$ over a range of treatment values $t$, such as the average derivative that summarizes certain aggregate effects PSS89ETA and the bound of the average welfare effect in CHN. In contrast, our average structural function $\beta_t$ and the partial effect $\theta_t$ capture the fully nonparametric heterogenous effects of $T$.

The paper proceeds as follows. Section (ref) introduces the framework and estimation procedure. Section (ref) presents the asymptotic theory of the DML estimators and low-level conditions for various nuisance function estimators. Section (ref) introduces the Simulated DML estimator. Section (ref) demonstrates the usefulness of our DML estimator with various ML methods in Monte Carlo simulations and an empirical example on the Job Corps program evaluation. All the proofs are in the Appendix. Additional results, such as uniform inference theory, are in the online supplementary appendix.

Setup and estimation

Let $\{Y_i, T_i, X_i\}_{i=1}^n$ be an i.i.d. sample from $Z = (Y, T^\prime, X^\prime)^\prime \in \mathcal{Z} = \mathcal{Y}\times\mathcal{T}_0\times\mathcal{X} \subseteq \mathcal{R}^{1+d_T+d_X}$ from a population $P$ with a cumulative distribution function (CDF) $F_{Z}(Z)$. Consider a set of treatment values of interest $\mathcal{T}$ to be an interior of the support of $T$.

The commonly used identifying Assumption (ref)(a) based on observational data (also known as unconfoundedness, selection on observables, or ignorability) assumes that conditional on observables, the treatment variable is as good as randomly assigned, or conditionally exogenous.

Define the product kernel as $K_h(T_i-t) \equiv \Pi_{j=1}^{d_T} k((T_{ji} - t_j)/h)/h^{d_T}$, where $T_{ji}$ is the $j^{th}$ component of $T_i$ and the kernel function $k()$ satisfies Assumption (ref). Denote the roughness of $k$ as $R_k \equiv \int_{-\infty}^\infty k(u)^2 du$ and $\kappa \equiv \int_{-\infty}^\infty u^2 k(u) du$.

Assumption (ref) is standard in nonparametric kernel estimation and holds for commonly used kernel functions, such as Epanechnikov and Gaussian. By Assumptions (ref)-(ref) and the same reasoning for the binary treatment, it is straightforward to show the identification,

for $t\in\mathcal{T}$.\footnote{For identification, we only need $\inf_{t\in\mathcal{T}} {\rm ess}\inf_{x\in\mathcal{X}} f_{T|X}(t|x) > 0$ for ((ref)) and additionally $\gamma(t,x)$ and $f_{T|X}(t,x)$ to be continuous in $t$ for ((ref)). Such standard identifying conditions are weaker than Assumption (ref) used for the inference theory.} The expression in equation ((ref)) motivates the class of regression-based (or imputation) estimators, while equation ((ref)) motivates the class of inverse probability weighting estimators; see Section S2.1 in the online supplementary appendix for further discussion. Now we introduce the double debiased machine learning estimator. \\ \\ Estimation procedure:

\setcounter{bean}{0}

The number of folds in cross-fitting $L$ is not random and typically small, such as five or ten in practice; see, e.g.\ Section (ref) or CCDDHNR. When there is no sample splitting ($L=1$), $\hat \gamma_1$ and $\hat f_1$ use all observations in the full sample. Then the DML estimator $\hat\beta_t$ in ((ref)) is the doubly robust estimator considered in Kallus and SUZ.

Our inference theory requires the estimators $\hat\gamma_\ell$ and $\hat f_\ell$ in Step 1 to satisfy Assumption (ref) below. We define the {\it partial $L_2(tX)$ norm} for any $t\in \mathcal{T}$ as $\|\hat\gamma_\ell-\gamma\|_{F_{tX}} \equiv \|\hat\gamma_\ell(T,X)-\gamma(T,X)\|_{F_{tX}} \equiv \big(\int_{\mathcal{X}} (\hat \gamma_\ell(t,x) - \gamma(t,x))^2 f_{TX}(t,x) dx \big)^{1/2}$ and $\|\hat f_\ell -f_{T|X}\|_{F_{tX}} \equiv \|\hat f_\ell(T|X)- f_{T|X}(T|X)\|_{F_{tX}} \equiv \big(\int_{\mathcal{X}} (\hat f_\ell(t|x) - f_{T|X}(t|x))^2 f_{TX}(t, x) dx \big)^{1/2}$, where the joint distribution $F_{TX}(t,X)$ is evaluated at a fixed value of $T$ equal to $t$.

The nuisance function estimators $\hat\gamma_\ell$ and $\hat f_\ell$ converge to some fixed functions $\bar\gamma$ and $\bar f$ respectively in the sense of Assumption (ref)(a)(b). Assumption (ref)(c) allows one of the nuisance functions to be misspecified, or requires at least one of the nuisance function estimators to be consistent. Assumption (ref)(b) and (c) imply that if one nuisance function is misspecified, then the other needs to be estimated consistently at a convergence rate faster than $\sqrt{nh^{d_T}}$. This is the cost of our doubly robust inference. Specifically, if $\bar\gamma \neq \gamma$, then $\|\hat\gamma_\ell - \gamma\|_{F_{tX}} = O_p(1)$. So Assumption (ref)(b) requires $\sqrt{n h^{d_T}} \|\hat f_{\ell} - f_{T|X}\|_{F_{tX}} = o_p(1)$. On the other hand, if $\bar f \neq f_{T|X}$, then we require $\sqrt{n h^{d_T}} \|\hat \gamma_{\ell} - \gamma\|_{F_{tX}} = o_p(1)$. Kennedy show such double robustness using a stronger uniform norm for a local linear estimation.

In Section (ref), we provide sufficient low-level conditions for Assumption (ref) when the nuisance function estimators are kernel estimators, series, and the deep neural networks in Farrell18. The partial $L_2$ convergence rate also appears in Kennedy, the conditional average treatment effect in Fan, and covariate adjustments in regression discontinuity designs in RotheRD.

Conditional density estimation

We propose an estimator for the {\it reciprocal of the generalized propensity score} (GPS), i.e.\ $1/f_{T|X}(t|x)$, when $d_T = 1$. The estimator avoids plugging in a small estimate in the denominator, and the estimate is positive by construction. When $d_T > 1$, we propose an estimator for the GPS. We can use various nonparametric and ML methods designed for the conditional expectation. We provide a root-mean-square convergence rate. In Section (ref), we demonstrate these generic GPS estimators using the deep neural networks in Farrell18 and show how Assumption (ref) can be verified.

The theory of ML methods in estimating the conditional density is less developed compared with estimating the conditional expectation. Alternative estimators for estimating the GPS can be the kernel density estimator, the artificial neural networks in ChenWhite, the Lasso methods in SUZ and BCK19JASA, or the series cross-validated method in Zhang.

It is known that for any CDF $F$, $\frac{d}{du} F^{-1}(u) = \frac{1}{F'(F^{-1}(u))}$ for $u\in(0,1)$. So $\frac{1}{f_{T|X}(t|x)} = \frac{\partial}{\partial u} F^{-1}_{T|X}(u|x)\big|_{u=F_{T|X}(t|x)}$. Inspired by the idea in Koenker94, we estimate $\frac{1}{f_{T|X}(t|x)}$ by a numerical differentiation estimator, labelled as ReGPS,

where $\epsilon = \epsilon_{n}$ is a positive sequence vanishing as $n$ grows and $\hat F_{T|X}(t|x) \pm \epsilon \in (0,1)$. By a standard algebra, the conditional CDF $F_{T|X}(t|x) = \lim_{h_1 \rightarrow 0} \mathbb{E}\left[\Phi\left(\frac{t-T}{h_1}\right)\big|X=x \right]$, where $\Phi$ is the CDF of a standard normal random variable and $h_1 = h_{1n}$ is a bandwidth sequence vanishing as $n$ grows. Let $\hat\mu(W; x)$ be a generic estimator of the conditional expectation $\mathbb{E}\left[W|X= x\right]$ for an outcome variable $W$. Then we estimate $F_{T|X}(t|x)$ by $\hat F_{T|X}(t|x) = \hat \mu\left(\Phi\left(\frac{t-T}{h_1}\right);x\right)$ with a transformed outcome variable of $T$, $W = \Phi\left(\frac{t-T}{h_1}\right)$. The conditional $u$-quantile function $F_{T|X}^{-1}(u|x)$ is estimated by the generalized inverse function $\hat F_{T|X}^{-1}(u|x) = \inf_{t\in\mathcal{T}}\{t: \hat F_{T|X}(t|x) \geq u\}$. When $\hat F_{T|X}(t|x)$ is continuous in $t$, $\hat F^{-1}_{T|X}(u|x)$ is strictly increasing. Then the resulting estimator $\widehat{1/f_{T|X}(t|x)} > 0$.

Denote the standard root-mean-square norm, or the $L_2(X)$ norm, of a random vector $X$ with distribution $F_X$ as $\left\|\hat\mu(W; X) - \mathbb{E}[W|X]\right\|_{F_X} \equiv \Big(\int_\mathcal{X} \big(\hat\mu(W; x) - \mathbb{E}[W|X=x]\big)^2 f_X(x)dx \Big)^{1/2}$ for a random variable $W$.

Assumption (ref) specifies the root-mean-square convergence rate of our conditional density estimator. Thus as long as the root-mean-square convergence rate of a ML method ($R_{1n}$) for the conditional expectation function is available, Assumption (ref) can be satisfied with a suitable bandwidth $h_1$ and $\epsilon$. Then we are able to use such a ML method to estimate the conditional density, as illustrated in Section (ref).

Note that $\hat F^{-1}_{T|X}(u|x)$ estimates the conditional quantile function. The ReGPS estimator inverses the CDF estimate and allows to apply various ML methods for conditional expectation functions. Alternatively we can estimate the conditional quantile function directly by a $\ell_1$-penalized quantile regression; for example, the conditional density function estimation in BCK19JASA.

When $d_T > 1$, we propose a direct estimator for the conditional density function $f_{T|X}(t|x)$ by

labelled as MultiGPS, where the bandwidth $h_1$ is a positive sequence vanishing as $n$ grows, the product kernel $g_{h_1}(T_i-t) \equiv \Pi_{j=1}^{d_T} g((T_{ji} - t_j)/h_1)/h_1^{d_T}$. We can choose $g()$ to be the Gaussian kernel.\footnote{A possible drawback of MultiGPS is that the estimate could be negative or small in finite samples. We may adopt the trimming/flooring approaches to addressing this concern in the literature. For example, following HLL20ER, we can use the estimate $\max\{\hat f_{T|X}(t|X_i), \delta_n\}$ for some positive sequence $\delta_n\rightarrow 0$. }

Asymptotic theory

We present the asymptotically linear representation and asymptotic normality. We provide low-level conditions for estimating the nuisance functions by the deep neural networks in Farrell18 in Section (ref). Conditions for kernel and series estimators are in Section S3.2 in the online supplementary appendix. Section (ref) provides sufficient rate conditions using the standard $L_2$ norm.

Let $\partial_t^\nu g(t, \cdot) \equiv \partial^\nu g(t, \cdot)/\partial t^\nu$ denote the $\nu$th partial derivative of a generic function $g$ with respect to $t$, and $\partial_t \equiv \partial_t^1$.

Note that the second part in the influence function\footnote{ For our non-regular parameters, we borrow the terminology “influence function" in estimating a regular parameter that is $\sqrt{n}$-estimable. An influence function gives the first-order asymptotic effect of a single observation on the estimator. The estimator is asymptotically equivalent to a sample average of the influence function. See Hampel and IchimuraNewey22QE, for example. } in ((ref)) $n^{-1}\sum_{i=1}^n \bar\gamma(t,X_i) - \beta_t = O_p(1/\sqrt{n}) = o_p(1/\sqrt{nh^{d_T}})$ and hence does not contribute to the first-order asymptotic variance $\mathsf{V}_t$. We keep these smaller-order terms to show that the nuisance function estimators have no first-order influence on the asymptotic distribution of $\hat \beta_t$. This is in contrast to the binary treatment case where $K_h(T_i-t)$ is replaced by ${\bf 1}\{T_i - t\}$ in $\hat \beta_t$, so $\hat \beta_t$ converges at a root-$n$ rate. Then the second part in ((ref)) is of first-order for a binary treatment, resulting in the well-studied efficient influence function in estimating the binary treatment effect in Hahn98ETA.

Theorem (ref) is fundamental for inference, such as constructing confidence intervals and the optimal bandwidth $h$ that minimizes the asymptotic mean squared error. We propose an estimator for the leading bias $\mathsf{B}_t$, inspired by the idea in PS96. Let the notation $\hat\beta_t = \hat\beta_{t,b}$ be explicit on the bandwidth $b$ and \[ \hat{\mathsf{B}}_t \equiv \frac{\hat\beta_{t,b} - \hat\beta_{t,ab}}{b^2 (1-a^2)} \] with a pre-specified fixed scaling parameter $a \in (0, 1)$. Theorem (ref) below shows the consistency of $\hat{\mathsf{B}}_t$ under Assumption (ref)(d).

We can estimate the asymptotic variance $\mathsf{V}_t$ by the sample variance of the estimated influence function $\hat{\mathsf{V}}_t \equiv h^{d_T} n^{-1} \sum_{\ell=1}^L\sum_{i \in I_\ell} \hat \psi_{i\ell}^2$, where $\hat \psi_{i\ell} \equiv K_h(T_i-t) (Y_i - \hat \gamma_\ell(t,X_i) )/\hat f_\ell(t|X_i) + \hat \gamma_\ell(t,X_i) - \hat\beta_t$. Then we propose a data-driven bandwidth $\hat h_t \equiv \big(d_T \hat{\mathsf{V}}_t\big/\big(4\hat{\mathsf{B}}_t^2\big)\big)^{1/(d_T+4)} n^{-1/(d_T+4)}$ to consistently estimate the optimal bandwidth that minimizes the asymptotic mean squared error (AMSE) given in Theorem (ref).

Assumption (ref)(a)-(c) are for the consistency of $\hat{\mathsf{V}}_t$. The condition (a) strengthens Assumption (ref)(a), and (b) is mild boundedness conditions that are implied when $\hat\gamma_\ell$ and $\hat f_\ell$ are bounded uniformly. In practice, we may use different Step 1 estimators in $\hat{\mathsf{V}}_t$ and $\hat\beta_t$ due to different high-level conditions.

A common approach is to choose an undersmoothing bandwidth $h$ smaller than $h^\ast_t$ such that the bias is first-order asymptotically negligible, i.e.\ $h^2\sqrt{nh^{d_T}} \rightarrow 0$. Then we can construct the usual $(1-\alpha)\times 100\%$ point-wise confidence interval $\Big[\hat \beta_t \pm \Phi^{-1}(1-\alpha/2) \sqrt{\hat{\mathsf{V}}_t/(nh^{d_T})} \Big]$, where $\Phi$ is the CDF of $\mathcal{N}(0,1)$. Alternatively, we may consider a further bias correction to allow for a wider range of bandwidth choice so that we may implement $\hat h_t$ in practice. Specifically we may use the above bias estimator $\hat{\mathsf{B}}_t$ and account for its variation in the asymptotic theory of the bias-corrected estimator $\hat\beta_t - h^2\hat{\mathsf{B}}_t$. CCF show that the AMSE optimal bandwidth of the original estimator is feasible in different contexts. Westling develop such robust bias-corrected inference for the local linear estimator in Kennedy.

Next we present the asymptotic theory for $\hat \theta_t$.

The conditions (a) and (b) in Theorem (ref) strengthen Assumption (ref) for $\beta_t$ and imply that $\eta$ cannot be too small and depends on the precision of the nuisance function estimators. The bias $\mathsf{B}_{\theta_t}^{mis}$ is due to misspecifying one of the nuisance functions and is zero when both nuisance functions are correctly specified or estimated by nonparametric methods.

We propose a data-driven bandwidth $\hat h_{\theta_t} \equiv \big((d_T+2) \hat{\mathsf{V}}_t^\theta\big/\big(4(\hat{\mathsf{B}}_t^\theta)^2\big)\big)^{1/(d_T+6)} n^{-1/(d_T+6)}$ to consistently estimate the optimal bandwidth that minimizes the AMSE given in Theorem (ref). Following the same procedure for $\hat\beta_t$, let $\hat{\mathsf{V}}_t^\theta \equiv \hat{\mathsf{V}}_t \int k'(u)^2du R_k^{-1}$ and $\hat{\mathsf{B}}^\theta_t \equiv \big(\hat\theta_{t,b} - \hat\theta_{t,ab}\big)/\big(b^2 (1-a^2)\big)$. In practice, we may choose $\eta = h n^{-a}$, where $a > 1/(d_T + 6)$, such that the conditions in Theorem (ref) are satisfied.

Nuisance function estimators

We show that the high-level conditions on the convergence rates in Assumption (ref) are attainable by the nonparametric and ML methods: kernel, series, and the deep neural networks in Farrell18, where the dimension of the control variables $d_X$ is fixed. Lasso methods have been extensively studied in SUZ, SU_PRTE, and SUZ_UQR, where $d_X$ can grow with $n$. These ML methods require different low-level conditions, such as dimensionality, smoothness, and tuning parameters.

Consider the conditional density estimator MultiGPS $\hat f_{T|X}$ given in Section (ref) for example. By Lemma (ref), Assumption (ref)(b) requires

Therefore, we need to obtain the partial $L_2(tX)$ convergence rate $\|\hat\gamma-\gamma\|_{F_{tX}}$ and the standard $L_2(X)$ convergence rate $R_{1n}$.

We seek theoretical results, such as the above rate conditions, for insights on selection of the tuning parameters in practice that is challenging and under-developed in the ML literature. We may use the optimal choices for the nuisance function estimators, as they do not affect the first-order asymptotics. A common method is cross-validation. We may choose the optimal rates for the bandwidths for $\hat\gamma$ and $\hat\mu(h_1^{d_T}g_{h_1}(T-t); x)$ that respectively minimize $\|\hat\gamma-\gamma\|_{F_{tX}}$ and $\big\| \hat\mu\big(h_1^{d_T}g_{h_1}(T-t); X\big) - {\mathbb{E}}[h_1^{d_T}g_{h_1}(T-t)|X] \big\|_{F_X}$, which might be available in the literature. Similarly we can derive the optimal $h_1^\ast \propto R_{1n}^{1/(2+d_T)}$.

Next we propose a deep MLP-ReLU network kernel estimator for $\gamma(t, x)$ (labelled as Kernel NN) and derive its $L_2(tX)$ convergence rate. We illustrate the low-level conditions of conventional kernel and series estimators in Section S3 in the online supplementary appendix, as the calculations are rather standard. These results on the $L_2(tX)$ convergence rates for neural networks and series are new and non-trivial extensions of existing results in the literature.

We consider the deep neural networks in Farrell18 (FLM, hereafter) that use the fully connected feedforward neural networks (multilayer perceptron, or MLP) and the nonsmooth rectified linear units (ReLU) activation function. We propose a deep MLP-ReLU network kernel estimator for $\gamma(t, x)$. The proposed estimator serves the purpose to conveniently apply the $L_2(TX)$ convergence rate given in FLM to obtain the $L_2(tX)$ convergence rate. So we can deliver valid asymptotic inference for $\beta_t$ and $\theta_t$ following deep learning. In this section, we closely follow the notations in FLM for easy reference, by slightly abusing our notations.

We consider a kernel-weighted loss function for any $t\in\mathcal{T}$,

where a product kernel $\mathsf{K}_b(T-t) \equiv \Pi_{j=1}^{d_T} \mathsf{k}((T_j - t_j)/b)/b^{d_T}$ with a kernel function $\mathsf{k}$ and a positive sequence of bandwidth $b = b_n$ vanishing as $n$ grows. We define the {\it deep MLP-ReLU network kernel estimator} for any $t\in\mathcal{T}$ as

where $\mathcal{F}_{MLP}$ is the MLP class, $M$ is an absolute constant, and $\theta$ depending on $t$ collects the weights and constants over all nodes. We refer the details of the MLP-ReLU network estimators to FLM. Then we obtain $\hat \gamma(t, x) = \hat f_{tb}(x)$.

Denote the derivative of a function $f(x)$ as $\mathtt{D}^{\alpha}_x f(x) = \frac{\partial^{|\bm{\alpha}|} f(x)}{\partial x_1^{\alpha_1}\cdots \partial x_{d_X}^{\alpha_{d_X}}}$, where $\bm{\alpha} = (\alpha_1,...,\alpha_{d_X})$ and $|\bm{\alpha}| = \alpha_1 +...+\alpha_{d_X}$.

Assumption (ref)(a) is due to the kernel weight $\mathsf{K}_b(T-t)$ in the loss function. Assumption (ref)(b)-(d) collect assumptions for applying Theorem 1 in FLM. Detailed discussion on these assumptions is referred to FLM. As discussed in FLM, it is standard in nonparametric analysis to assume the true function to be estimated is bounded. The choice of $M$ may be arbitrarily large and is simply a formalization of the requirement that the optimizer is not allowed to diverge on the function level in the sup-norm sense. For practical implementation, we do not impose such bound. But there is a practice for rescaling the output variable which is generally dependent on the activation function being used, e.g.\ the domain of the activation function. A common practice for variable transformation is standardization (subtracting mean and dividing by standard deviation) or scaling to a specific range by an affine transformation (generally chosen between 0 and 1).

We can apply deep neural networks to the GPS estimation proposed in Section (ref). The ReGPS estimator can use $\hat\mu\left(\Phi((t-T)/h_1); x\right) = \hat f_{MLP-ReGPS}(x)$ the MLP estimator in FLM with the unweighted loss function:

The MultiGPS estimator can use $\hat\mu\left(h_1^{d_T} g_{h_1}(T-t); x\right) = \hat f_{MLP-MultiGPS}(x)$:

We are ready to show that Assumption (ref) is attainable by the MLP-ReLU network estimators. Take the MultiGPS estimator in ((ref)) for example, with $d_T=1$ for simplicity. Assumption (ref)(b) is $nh \Big( (nb^{2})^{-\frac{r}{r+d_X}}\log^8n +\log\log n/(nb) + b^2 \Big) \Big( h_1^{-2}\big( n^{-\frac{r}{r+d_X}}$\\ $\log^8n + \log\log n/n\big) + h_1^4 \Big) \rightarrow 0$ by Theorem (ref), Lemmas (ref) and (ref). Assumption (ref)(a) is implied by $nb^2\rightarrow\infty$ and $n^{r/(r+d_X)}h_1^2\rightarrow\infty$. When $h = h_1 = b$, Assumption (ref)(a) holds by letting $n^{r/(r+d_X)}h^2 \rightarrow\infty$, and (ii) holds by letting smoothness $r > {d_X}$. FLM discuss the same condition $r > {d_X}$ for the average treatment effect of a binary treatment variable. Similarly we note that this condition is not minimal but is sufficient to justify the practical use of the MLP-ReLU network estimators for valid inference on the average structural function and the partial effect of continuous treatments by our approach.

Standard $L_2$ norm rate

We provide sufficient rate conditions using the standard $L_2(TX)$ norm in Assumption (ref) to replace the partial $L_2(tX)$ norm in Assumption (ref). Let the $L_2(TX)$ norm be $\|\hat\gamma_\ell-\gamma\|_{F_{TX}} \equiv \big(\int_\mathcal{T}\int_\mathcal{X} \big(\hat\gamma_\ell(t,x) - \gamma(t,x) \big)^2 f_{TX}(t,x)dx dt\big)^{1/2}$. We do not need to modify the rate condition on $\hat f_{\ell}$, as it equivalently uses the standard $L_2(X)$ norm for a given $t$. The cost of the more commonly used $L_2(TX)$ norm is losing the doubly robust inference and stricter regularity conditions.

Note that the rate condition Assumption (ref)(b) is the same rate condition for the regular semiparametric models in CCDDHNR and CNS21ADML, e.g.\ $\hat f_{T|X}(t|x)$ is replaced with the propensity score $P(T=t|X=x)$ for a discrete treatment. We discuss the cost and implications of Assumption (ref) based on the standard $L_2$ norm. First under Assumption (ref), the inference theory cannot allow $\gamma$ to be misspecified, while $f_{T|X}$ can be misspecified. Second, the convergence rates are faster than those required in Assumption (ref). To learn intuition on the rate condition (a), we utilize the bounded kernel in the DML estimator, resulting in a penalty $h^{-d_T/2}$ in the loose bound $h^{-d_T/2} \|\hat \gamma_\ell - \gamma\|_{F_{TX}}$. Third, the bandwidth $h$ choice is more restrictive. Assumption (ref)(c) implies undersmoothing, i.e.\ the leading bias of $\hat\beta_t$ is first-order ignorable by $h^2\sqrt{nh^{d_T}} \rightarrow 0$.

For a specific example, consider $\hat\gamma$ to be FLM's estimator for the conditional expectation function $\mathbb{E}[Y|T,X]$, i.e.\ using the unweighted loss function $(Y-f(T,X))^2/2$. FLM provide the corresponding $L_2(TX)$ convergence rate $\|\hat \gamma_\ell -\gamma\|_{F_{TX}}^2 = O_p\left( n^{-\frac{r}{r+d_X+d_T}}\log^8n +\log\log n/n\right)$ with the smoothness $r$ defined in Assumption (ref)(c). We need $r > d_T + d_X$ that is stronger than $r > d_X$, as discussed in Section (ref). Moreover comparing the rate conditions in Assumption (ref)(b) and Assumption (ref)(b) for $d_T=1$, we can show that using the partial $L_2(tX)$ norm results in a tighter bound by $\sqrt{nh}\|\hat\gamma_\ell - \gamma\|_{F_{tX}} = o_p(\sqrt{n}\|\hat\gamma_\ell - \gamma\|_{F_{TX}})$ with $r = d_x + 1$, $b=h \propto n^{-a}$ and $a < 1/2$.

Simulated DML estimator

We introduce Simulated DML estimator $\check\beta_t$ that enables the high-level rate conditions based on the standard $L_2(TX)$ norm, rather than the partial $L_2(tX)$ norm, and also permits the doubly robust inference. Following the estimation procedure given in Section (ref), Step 1 computes the nuisance function estimators $\hat \gamma_\ell$ and $\hat f_\ell$. In Step 2, let $U_i \equiv \mathsf{T}_i h_0 + t$ where $\{\mathsf{T}_i\}_{i\in I_\ell}$ are i.i.d. draws from $\{T_i\}_{i\in I_\ell}$ with replacement, for $i \in I_\ell$, $\ell\in\{1,..., L\}$, and $h_0$ is a positive sequence converging to zero as $n\rightarrow\infty$. The Simulated DML (SDML) estimator is defined as

The corresponding partial effect estimator $\check\theta_t \equiv (\check\beta_{t^+} -\check\beta_{t^-})/\eta$ as in Step 3.

To get intuition, the SDML estimator $\check\beta_t$ uses $\hat\gamma_\ell(U_i, X_i)$ rather than $\hat\gamma_\ell(t, X_i)$ as in $\hat\beta_t$ that fixes the treatment value at the target $t$. The simulated $U_i$ localizes the realized treatment values around $t$. Introducing such local variation enables the standard $L_2$ rate of $\hat\gamma$. Specifically $\mathsf{T}_i$ defined above follows the empirical distribution function of $\{T_i\}_{i\in I_\ell}$, denoted as $\hat F_{T\ell}$. Therefore $U_i$ follows a CDF conditional on the sample $\{Z_i=(Y_i, X_i, T_i)\}_{i=1}^n$, $P(U_i \leq u|\{Z_i\}_{i=1}^n) = P(\mathsf{T}_i h_0 + t \leq u|\{Z_i\}_{i=1}^n) = \hat F_{T\ell}((u-t)/h_0)$.

Assumption (ref) gives the high-level conditions on the nuisance function estimators.

Assumption (ref)(a) strengthens the consistency condition in Assumption (ref)(a) with a penalty $h_0^{-d_T/2}$. Interestingly when $h_0 = h$, the rate condition (b) $\sqrt{n}\|\hat\gamma_\ell-\gamma\|_{F_{TX}} \|\hat f_\ell - f_{T|X}\|_{F_{tX}} = o_p(1)$, which is the same rate condition for the regular semiparametric models in CCDDHNR and CNS21ADML. Assumption (ref)(c) is a boundedness condition that is implied if $\hat\gamma$ is uniformly bounded; for example, deep neural networks in FLM for a uniformly bounded $Y$. Specifically, FLM provide the corresponding $L_2(TX)$ convergence rate $\|\hat \gamma_\ell -\gamma\|_{F_{TX}}$ and illustrate its usefulness for semiparametric inference on the average treatment effect of a binary treatment. We could use this rate to verify the high-level conditions in Assumption (ref) for the SDML estimator or Assumption (ref) for the DML estimator as discussed in Section (ref).

We can estimate the AMSE optimal bandwidth $h_t^\ast$ given in Theorem (ref) by the SDML approach. We can estimate the leading bias $\mathsf{B}_t$ by $\check{\mathsf{B}}_t \equiv \frac{\check\beta_{t,b} - \check\beta_{t,ub}}{b^2 (1-a^2)}$ with a pre-specified positive scaling parameter $a\in(0,1)$. We can estimate the asymptotic variance $\mathsf{V}_t$ by $\check{\mathsf{V}}_t \equiv h^{d_T} n^{-1} \sum_{\ell=1}^L\sum_{i \in I_\ell} \check \psi_{i\ell}^2$, where $\check \psi_{i\ell} \equiv K_h(T_i-t) (Y_i - \hat \gamma_\ell(U_i,X_i) )/\hat f_\ell(t|X_i) + \hat \gamma_\ell(U_i,X_i) - \check\beta_t$. Then a data-driven bandwidth $\check h_t \equiv \big(d_T \check{\mathsf{V}}_t\big/\big(4\check{\mathsf{B}}_t^2\big)\big)^{1/(d_T+4)} n^{-1/(d_T+4)}$. Theorem (ref) below shows the consistency of $\check{\mathsf{B}}_t$, $\check{\mathsf{V}}_t$, and $\check h_t$, under Assumption (ref) that is modified from Assumption (ref).

We can similarly obatin the asymptotic theory for $\check\theta_t$.

Numerical examples

This section provides numerical examples of Monte Carlo simulations and an empirical illustration. The estimation procedure of the proposed DML estimator is described in Section (ref). To estimate the first-step conditional expectation function $\gamma(t,x) = \mathbb{E}[Y|T=t, X=x]$ and the conditional density $f_{T|X}$ by MultiGPS as described in Section (ref), we employ three methods: Lasso, the deep neural networks (NN) based on Farrell18, and the Kernel NN proposed in Section (ref). We implement our DML estimator with these algorithms respectively. Note that Lasso assumes certain sparsity specifications and NN is a nonparametric method. Our doubly robust inference theory allows one of the first-step functions to be misspecified. Software we develop is available at \url{https://github.com/KColangelo/Double-ML-Continuous-Treatment}. Section S1 in the online supplementary appendix provides the implementation details and additional results of ReGPS.

Simulation study

We consider the data-generating process: $\nu \sim \mathcal{N}(0,1)$, $\varepsilon \sim \mathcal{N}(0,1)$,

where $\theta_j = 1/j^2$, $diag(\Sigma) = 1$, the $(i,j)$-entry $\Sigma_{ij} = 0.5$ for $|i-j|=1$ and $\Sigma_{ij} = 0$ for $|i-j|>1$ for $i,j=1,...,100$, and $\Phi$ is the CDF of $\mathcal{N}(0,1)$. Thus this is a nonseparable model and the potential outcome $Y(t) = 1.2t + 1.2 X'\theta + t^2 + tX_1 + \varepsilon*\sqrt{0.5+\Phi(X_1)}$. The parameters of interest are the average dose response function and the partial effect at $t = 0$, i.e.\ $\beta_0 = \mathbb{E}[Y(0)] = 0$ and $\theta_0 = \partial \mathbb{E}[Y(t)]/\partial t|_{t=0} = 1.2$.

We compare estimations with cross-fitting and without cross-fitting, and with a range of bandwidths to demonstrate robustness to bandwidth choice. We consider sample size $n \in \{1000,10000\}$ and the number of subsamples used for cross-fitting $L \in\{1, 5\}$. We use the second-order Epanechnikov kernel with bandwidth $h$. For the MultiGPS estimator described in Section (ref), we choose bandwidth $h_1=h$. Let the bandwidth $h=c \sigma_T n^{-0.2}$ for a constant $c\in\{0.75, 1.0,1.25, 1.5\}$ and the standard deviation $\sigma_T$ of $T$. We computed the AMSE-optimal bandwidth $h^\ast_0$ given in Theorem (ref) that has the corresponding $c^\ast = 1.45$. Thus using some undersmoothing bandwidth with $c < c^\ast$, the 95% confidence interval $\big[\hat \beta_t \pm 1.96 s.e.\big]$ is asymptotically valid, where the standard error ($s.e.$) is computed using the sample analogue of the estimated influence function, as described in Section (ref).

Table (ref) reports the results based on 1,000 Monte Carlo replications. Under no cross-fitting ($L=1$), the confidence intervals generally have lower coverage rates and the bias is larger than under cross-fitting. The estimators using NN and Kernel NN perform well in the case of fivefold cross-fitting, with coverage rates (Cov.) near the nominal 95%, especially in a smaller sample size $n=1,000$. Intuitively Kernel NN estimates the conditional expectation function $\gamma(t,x)$ locally at $t$, resulting in a smaller effective sample size $nb^{d_T}$ than the full sample size $n$ used by NN. Therefore NN may outperform Kernel NN in small samples. The coverage rate and bias are improved the most for NN and Kernel NN with cross-fitting, but only marginally for Lasso. Cross-fitting should improve our estimation in the case that the machine learning algorithm is over-fitting. Given that cross-fitting does not improve Lasso, it might suggest that Lasso does not have a severe over-fitting problem for this data-generating process.

These methods seem robust to bandwidth choice under cross-fitting. Overall these results demonstrate consistency with the theoretical results of this paper, confirming the usefulness of cross-fitting for ML methods.

{

}

Empirical illustration

We illustrate our method by re-analyzing the Job Corps program in the United States, which was conducted in the mid-1990s. The Job Corps program is the largest publicly funded job training program, which targets disadvantaged youth. The participants are exposed to different numbers of actual hours of academic and vocational training. The participants' labor market outcomes may differ if they accumulate different amounts of human capital acquired through different lengths of exposure. We estimate the average dose response functions to investigate the relationship between employment and the length of exposure to academic and vocational training. As our analysis builds on FFGN12ReStat, HHLP, and Lee, we refer the readers to the reference therein for further details of Job Corps.

We use the same dataset in HHLP. We consider the outcome variable ($Y$) to be the proportion of weeks employed in the second year following the program assignment. The continuous treatment variable ($T$) is the total hours spent in academic and vocational training in the first year. We follow the literature to assume the conditional independence Assumption (ref)(a), meaning that selection into different levels of the treatment is random, conditional on a rich set of observed covariates, denoted by $X$. The identifying Assumption (ref) is indirectly assessed in FFGN12ReStat. Our sample consists of 4,024 individuals who completed at least 40 hours (one week) of academic and vocational training. There are 40 covariates measured at the baseline survey. In the online supplementary appendix, Figure S1 shows the distribution of $T$ by a histogram, and Table S2 provides brief descriptive statistics. \\[5pt] Implementation details:\ We estimate the average dose response function $\beta_t = \mathbb{E}[Y(t)]$ and partial effect $\theta_t = \partial \mathbb{E}[Y(t)]/\partial t$ by the proposed DML estimator with fivefold cross-fitting. We implement three DML estimators: Lasso, the generalized random forests in ATW19AS, the neural networks (NN) based on Farrell18, and the Kernel NN proposed in Section (ref). The parameters for these methods are selected as described in Section S1 in the online supplementary appendix.

We use the second-order Epanechnikov kernel with bandwidth $h$. For the MultiGPS estimator, we use the Gaussian kernel with bandwidth $h_1=h$. We compute the optimal bandwidth that minimizes an asymptotic integrated MSE. For practical implementation, consider a weight function $w(t) = {\bf 1}\{ t \in [\underline{t}, \bar t]\}/(\bar t - \underline{t})$ that is the density of $Uniform[\underline{t}, \bar t]$ on a subset of the support of $T$. The bandwidth that minimizes the asymptotic integrated MSE $ \int_{\mathcal{T}} \big( \mathsf{V}_t/(nh^{d_T}) + h^4 \mathsf{B}_t^2\big) w(t) dt$ for an integrable weight function $w(t): \mathcal{T} \mapsto \mathcal{R}$ is $h^\ast_w = \big(d_T\mathsf{V}_w\big/\big(4\mathsf{B}_w\big)\big)^{1/(d_T+4)} n^{-1/(d_T+4)}$, where $\mathsf{V}_w \equiv \int_{\mathcal{T}} \mathsf{V}_t w(t) dt$ and $\mathsf{B}_w \equiv \int_{\mathcal{T}} \mathsf{B}_t^2 w(t) dt$, following Theorem (ref). Set $m$ equally spaced grid points over $[\underline{t}, \bar t]$: $\big\{\underline{t} = t_1, t_2,..., t_m = \bar t\big\}$. Following the approach given in Section (ref), we estimate ${\mathsf{V}}_{t_j}$ with $h=3\hat\sigma_Tn^{-0.2} = 548.52$ and ${\mathsf{B}}_{t_j}$ with $b = 2h$ and $a=0.5$, for $j=1,...,m$. A plug-in estimator $\hat h^\ast_w = \big(\hat{\mathsf{V}}_w\big/\big(4\hat{\mathsf{B}}_w\big)\big)^{1/5} n^{-1/5}$, where $\hat{\mathsf{V}}_w = m^{-1}\sum_{j=1}^m \hat{\mathsf{V}}_{t_j}$ and $\hat{\mathsf{B}}_w = m^{-1}\sum_{j=1}^m \hat{\mathsf{B}}_{t_j}^2$. We use $[\underline{t}, \bar t] = [160, 1840]$ and $t_j - t_{j-1} = 40$ in this empirical application. We then obtain under-smoothing bandwidths $0.8 \hat h^\ast_w$ that are 213.45 for Lasso, 224.36 for the generalized random forest, 223 for NN, and 225.11 for Kernel NN.

Results:\ Figure (ref) presents the estimated average dose response function $\beta_t$ along with 95% point-wise confidence intervals. The results for the three ML nuisance function estimators have similar patterns. The estimates suggest an inverted-U relationship between the employment and the length of participation. DNN estimates appear to be the most erratic, possibly due to the smaller bandwidth compared with other estimators.

Figure (ref) reports the partial effect estimates $\hat\theta_t$ with step size $\eta =160$ (one month). Across all procedures, we see positive partial effects when hours of training are less than around 500 (three months) and negative partial effect around 1,500 hours (9 months). Taking the estimates by NN for example, $\hat \beta_{400} = 46.83$ with standard error $s.e. = 1.38$ and $\hat \theta_{400} = 0.0217$ with $s.e. = 0.0123$ computed based on the result of Theorem (ref). This estimate implies that increasing the training from two months to three months increases the average proportion of weeks employed in the second year by $3.47\%$ (nearly two weeks) with $s.e. = 1.962\%$.

Lee09RES finds that the program had a negative impact on employment propensities in the short term (104 weeks since random assignments) and a positive effect in the long term (104-208 weeks). Lee09RES considers a binary treatment variable of being in the program or not, with the outcome variable ${\bf 1}\{Y \geq 0\}$ in our notations. We focus on the employment proportion in the second year following the program assignment (52-104 weeks) and estimate the heterogenous effects of the total hours spent in academic and vocational training in the first year.

The empirical practice has focused on semiparametric estimation; see FFGN12ReStat, HHLP, Lee, for example. The semiparametric methods are subject to the risk of misspecification. Our DML estimator provides a feasible approach to implementing a fully nonparametric inference in practice.

Conclusion and outlook

For a future extension, our DML estimator serves as the preliminary element for policy learning and optimization with a continuous decision, following Manski04, HiranoPorter, KT18ECMA, Kallus, DSLC, AW_PL, Farrell18, among others.

Another extension is robustness against multiway clustering, where the conventional cross-fitting does not ensure the independence between observations in $I_\ell$ from $I_\ell^c$. We may adopt the $K^2$-fold multiway cross-fitting proposed by Chiang that focus on regular DML estimators as in CCDDHNR. Since the form of estimators and the proofs of asymptotic theories for our continuous treatment case are similar to those studied in Chiang, we expect that their proposed algorithm works for $\hat \beta_t$; a formal extension is out of the scope of this paper.

When unconfoundedness is violated, we can use the control function approach in triangular simultaneous equations models by including in the covariates some estimated control variables using instrumental variables. For example, Lee09RES studies the issue of sample selection for the wage effects of the Job Corps program. To extend our empirical application to the wage effect of the length of exposure to the program, we may follow Lee09RES to estimate bounds on the wage effect of the continuous treatment using the excess number of individuals who are induced to be selected. A closer approach to our estimator is DNV03, who show that a nonparametric control function method accounts for both selection and endogeneity. IN09ETA show that the conditional independence assumption holds when the covariates $X$ include the additional control variable $V = F_{T|Z}(T|Z)$, the conditional distribution function of the endogenous variable given the instrumental variables $Z$. The influence function that accounts for estimating the control variables as generated regressors has derived in Corollary 2 in Lee15. Lee15 shows that the adjustment terms for the estimated control variables are of smaller order in the influence function of the final estimator, but it may be important to include them to achieve local robustness. This is a distinct feature of the average structural function of continuous treatments, as discussed in Section (ref). Using such an influence function to construct the corresponding DML estimator is left for future research.

Acknowledgements

This work is not related to Kyle Colangolo's position at Amazon, and was completed outside of the regular duties of the position. We are grateful to Max Farrell, Whitney Newey, Takuya Ura, Ted Westling, and Yichong Zhang for valuable discussion. We thank Parush Arora for assistance.