Debiased Machine Learning when Nuisance Parameters Appear in Indicator Functions

Debiased Machine Learning when Nuisance Parameters Appear in Indicator Functions

Keywords: Debiased Machine Learning, High-dimensional Regression, Non-differentiability, Smoothing

Introduction

This paper studies debiased machine learning (DML) when nuisance parameters appear in indicator functions. An important example is the maximized average welfare gain under optimal treatment assignment rules. Provided that unconfoundedness assumption holds, the conditional average treatment effect (CATE) function is identified. The parameter of interest is represented by the expectation of a moment function where the moment function consists of indicator functions. For asymptotically valid inference for a parameter of interest, the current literature on debiased machine learning relies on Gateaux differentiability of the functions inside moment conditions. However, Gateaux differentiability does not hold when nuisance parameters appear in indicator functions, which makes the development of valid inference procedures an open problem.

Let $W=\left(Y,X^{'}\right)^{'}$ denote an observation where $Y$ is an outcome variable with a finite second moment and $X$ is a high-dimensional vector of covariates. Let \[ \boldsymbol{\gamma}_{0}\left(x\right)\equiv\mathbb{E}\left[Y\mid X=x\right] \] be the conditional expectation of $Y$ given $X\in\mathcal{X}.$ Let $\boldsymbol{\gamma}:\mathcal{X}\rightarrow\mathbb{R}$ be a function of $X.$ Define $m\left(w,\boldsymbol{\gamma}\right)$ as a function of the function $\boldsymbol{\gamma}$ (i.e. a functional of $\boldsymbol{\gamma}$), which depends on an observation $w.$ The parameter of interest $\theta_{0}$ has the following expression: \[ \theta_{0}=\mathbb{E}\left[m\left(W,\boldsymbol{\gamma}_{0}\right)\right]. \]

Chernozhukov.et.al.2022a shows an asymptotic distribution theory for DML when $m\left(w,\boldsymbol{\gamma}\right)$ is linear or nonlinear in $\boldsymbol{\gamma}.$ When $m\left(W,\boldsymbol{\gamma}\right)$ is nonlinear in $\boldsymbol{\gamma},$ Chernozhukov.et.al.2022a linearizes it and extends results for the linear case to the linearized function. The key assumption is that the remainder in the linearization, employing Gateaux differentiability in a neighborhood of the true parameter, is bounded by a constant. This assumption is crucial in showing asymptotic normality of the DML estimator in the nonlinear case as it renders the remainder term negligible. On the other hand, we focus on problems where $\boldsymbol{\gamma}$ appears in indicator functions, and hence $m\left(w,\boldsymbol{\gamma}\right)$ is not Gateaux differentiable in $\boldsymbol{\gamma}.$ This motivates us to propose an alternative approach in which a smoothing function is used to smooth the indicator function.

DML has been widely studied in econometrics literature. Chernozhukov.et.al.2017 and Chernozhukov.et.al.2018 propose a general DML approach for valid inference in the context of estimating causal and structural effects. Semenova.and.Chernozhukov2020 studies DML estimation of the best linear predictor (approximation) for structural functions including conditional average structural and treatment effects. Recently, Chernozhukov.et.al.2022a proposes automatic DML for linear and nonlinear functions of regression equations. They provide the average policy effect, weighted average derivative, average treatment effect and the average equivalent variation bound as examples of linear functions of regression equations. As an example of a nonlinear function, they discuss the causal mediation analysis of Imai.et.al.2010. Chernozhukov.et.al.2022b derives the convergence rate and asymptotic results for linear functionals of regression equations. Chernozhukov.et.al.2022c proposes a general method to construct locally robust moment functions for generalized method of moments estimation. The two key features of DML are orthogonal moments functions and cross-fitting. (Neyman) Orthogonal moment functions are used to mitigate regularization and/or model selection bias, and to avoid using plug-in estimators. When employing regularized machine learners in causal or structural estimation, squared bias may decrease at a slower rate than variance. As a result, confidence interval coverage can be poor and estimators will not be $\sqrt{n}$-consistent (Chernozhukov.et.al.2017, Chernozhukov.et.al.2018, and Chernozhukov.et.al.2022c). Combining orthogonal moment functions with cross-fitting makes inference available when regression estimators are high-dimensional.

This paper contributes to the expanding DML literature by considering estimation and inference for non-differentiable functions. We propose smoothing the indicator functions, and develop an asymptotic distribution theory for a class of models which involves indicator functions. We introduce a sigmoid function to smooth the indicator function, and a smoothing parameter which controls the smoothness of the sigmoid function. Asymptotically, the proposed estimator exhibits a trade-off between bias and variance. Its asymptotic behavior depends on two terms. One is a random component which is related to the sampling distribution of the proposed estimator. The other is a nonrandom term which represents the error introduced by approximating the indicator function with a sigmoid function. Smoothing affects the two components in different ways. The random component characterizes the variance of the estimator, and it shrinks as we smooth the indicator function. On the other hand, the nonrandom term represents the bias of the estimator, and blows up as we make the indicator function smoother. The bias order depends on the distribution of the CATE function, and we control its magnitude by imposing a margin assumption (Kitagawa.and.Tetenov2018). In light of the trade-off between bias and variance, we study an optimal choice for the smoothing parameter by minimizing an upper bound of the asymptotic mean squared error. Armstrong.and.Kolesar2020 proposes a method of constructing bias-aware confidence intervals. We construct a feasible version of this confidence interval in our setting. In addition, we derive theoretical results when the margin assumption does not hold.

We conduct Monte Carlo simulations to verify our theoretical results, and an empirical analysis to illustrate implementation of the methods. The simulation results show the asymptotic normality of our estimator and the applicability of the standard inference when a smoothing parameter is chosen optimally. In our empirical analysis, we use experimental data from the National Job Training Partnership Act (JTPA) Study (Bloom.et.al.1997). When the smoothing parameter is chosen optimally, we find that the estimate of maximized welfare gain is similar to previous studies.

Non-differentiability often arises in causal inference problems involving treatment assignment rules. D'Adamo.2022 studies the estimation of optimal treatment rules under partial identification. Christensen.et.al.2023 estimates treatment rules under directional differentiability with respect to a finite-dimensional nuisance parameter. Many papers propose various approaches to handle non-differentiability in specific settings. Horowitz1992 uses the sigmoid function to smooth the indicator function in analyzing the binary response model. The parameters of interest in Horowitz1992 are coefficients of the single index model, while our parameter of interest is the value of a criterion function. Zhou.et.al2017 replaces the 0-1 loss with the smoothed ramp loss in the framework of residual weighted learning to estimate individualized treatment rules. They construct the smoothed ramp loss by replacing the sharp cutoff on the interval $\left[-1,1\right]$ with a quadratic smoothing function. In our approach, we employ a smoothing parameter that depends on the sample size to smooth the indicator function using the sigmoid function. On the other hand, Chen.et.al.2003 studies a class of semiparametric optimization estimators when criterion functions are not smooth, and does not introduce a smoothing function when deriving the asymptotic distribution. Hirano.and.Porter.2012 shows that if the target object is non-differentiable in the parameters of the data distribution, there exist no estimator sequences that are locally asymptotically unbiased. Non-differentiability also arises in the nonparametric IV quantile regression through the non-smooth generalized residual functions. In response, Chen.and.Pouzo.2012 proposes a class of penalized sieve minimum distance estimators. Levis.et.al.2023 studies the covariate-assisted version of the Balke and Pearl bounds (Balke.and.Pearl.1997), which are characterized as non-smooth functionals (specifically, a max function). They smooth the max function using the log-sum-exp (LSE) function and provide an estimator based on the nonparametric efficient influence function for the smoothed functional. In contrast, we smooth the indicator function using a sigmoid function.

Standard bootstrap consistency fails in the presence of non-differentiability, which has led researchers to propose alternative bootstrap methods. Andrews2000 shows inconsistency of the bootstrap if the parameter lies on the boundary of a parameter space defined by linear or nonlinear inequality constraints. Fang.and.Santos.2018 studies inference for (Hadamard) directionally differentiable functions, and Hong.and.Li.2018 proposes a numerical derivative based Delta method to show consistent inference for functions of parameters that are only directionally differentiable. Recently, Kitagawa.et.al2020 characterizes the asymptotic behavior of the posterior distribution of functions which are locally Lipschitz continuous but possibly non-differentiable.

Various works study inference for welfare under optimal treatment assignment rules. Chen.et.al.2023 proposes similar approaches to ours, wherein they smooth the arg maximum operator using the soft-maximum operator. However, our method differs in several aspects. First, unlike their estimator, we propose a DML estimator where the orthogonal moment function involves a Riesz representer for the expectation of the derivative of the smoothing function. As Chernozhukov.et.al.2022a points out, this type of orthogonal moment function, consisting of a Riesz representer, can be understood as the efficient influence function. Second, our approach optimizes the mean squared error criterion in the smoothing parameter and offers a choice of the smoothing parameter in practice. Third, we construct a feasible version of bias-aware confidence intervals, while Chen.et.al.2023 eliminates bias by undersmoothing. Luedtke.and.van.der.Laan2016 studies inference for the mean outcome under optimal treatment rules by developing a regular and asymptotically linear estimator. Semenova2023a and Semenova2024 study estimation and inference for objects involving maximum or minimum of nuisance functions. These works consider plugging in the machine learning estimates of the nuisance functions into non-differentiable functions without any smoothing when the margin assumption is imposed. Our analysis shows that an optimal choice of the smoothing parameter is generally an interior value. In particular, when the margin assumption fails, the performance of plug-in based methods is not guaranteed. In this case, we derive more conservative, bias-aware confidence intervals using an optimal smoothing parameter chosen specifically for the setting without the margin assumption. This provides policy makers with a conservative inference strategy under weaker assumptions, offering an alternative when plug-in based methods are not applicable. Kitagawa.and.Tetenov2018 proposes a resampling based inference procedure for optimized welfare with potentially conservative coverage. Andrews.et.al2023 studies estimators and confidence intervals for the welfare at an estimated policy that controls the winner’s curse.

In addition, the average welfare gain from the unrestricted optimal treatment plays a critical benchmarking role in policy learning. Manski2004 evaluates the performance of statistical treatment rules in terms of their maximum regret and provides finite‐sample regret bounds for conditional empirical success (CES) rules. Kitagawa.and.Tetenov2018 assesses the properties of estimated treatment rules by their average welfare regret relative to the maximum feasible welfare gain by using the empirical welfare maximization method. Athey.and.Wager2021 studies policy learning using observational data with the doubly-robust approach from Chernozhukov.et.al.2022c. While much of the literature focuses on identifying the optimal treatment rule within a restricted class, the unrestricted optimal policy gain quantifies the maximum achievable benefit if treatments were allocated perfectly according to true conditional average treatment effects. Chernozhukov.et.al.2024a views it as a measure of the heterogeneity of treatment effects which quantifies the potential improvement over the average effect achievable through optimally tailored treatment assignments. Our estimator provides a consistent measure of the average welfare gain under the unrestricted optimal treatment rule, suggesting that our estimated benchmark can then be used to compare with the average welfare achieved by restricted treatment rules. The difference between the two may serve as an indication of the welfare loss incurred when policies are restricted to a particular class, thereby quantifying the potential benefit of allowing for more tailored treatment assignments.

The rest of the paper is organized as follows. Section (ref) introduces maximized average welfare gain under optimal treatment rules as an example where the parameter of interest is a non-differentiable function of regression equations. Section (ref) presents the estimation method and inference. Section (ref) gives simulation results. Section (ref) presents an empirical example. Section (ref) concludes the paper.

Non-differentiable Effects

In some cases the parameter of interest is the expectation of a non-differentiable function of regression equations. An important example is maximized average welfare gain under optimal treatment assignment rules. Consider the following potential outcomes framework. Let $D$ be a binary treatment status indicator and $Y\left(D\right)$ be a potential outcome. $Y\left(1\right)$ denotes the potential outcome upon receipt of the treatment, and $Y\left(0\right)$ represents the potential outcome without receipt of the treatment. The observed outcome $Y$ is written as \[ Y=Y\left(D\right)=DY\left(1\right)+\left(1-D\right)Y\left(0\right) \] Let $W=\left(Y,X^{'}\right)^{'}$ denotes an observation, with $X=\left(D,Z^{'}\right)^{'}$ where $Z$ is a high-dimensional vector of covariates. High-dimensional covariates are often considered in recent causal inference literature including Semenova.and.Chernozhukov2020. Let $\delta\left(Z\right)\in\left\{ 0,1\right\} $ be a treatment assignment function, where $\delta\left(Z\right)=1$ if treatment is assigned and $\delta\left(Z\right)=0$ if not. The maximized average welfare with respect to $\delta\left(Z\right)$ is expressed as follows. \[ \mathbb{E}\left[Y\left(1\right)\delta\left(Z\right)+Y\left(0\right)\left(1-\delta\left(Z\right)\right)\right]. \] Under the unconfoundedness assumption $\left(Y\left(1\right),Y\left(0\right)\right)\perp D\mid Z$ and the overlap condition, we can identify the CATE function \[ \tau\left(Z\right)=\mathbb{E}\left[Y\left(1\right)-Y\left(0\right)\mid Z\right]. \] If $\tau\left(Z\right)$ is known, the optimal treatment assignment rule $\delta^{*}\left(Z\right)$ is \[ \delta^{*}\left(Z\right)=\mathds{1}\left\{ \tau\left(Z\right)>0\right\} . \] The welfare gain relative to the no-one treated policy is

where the second equality holds by the law of the iterated expectations. The parameter of interest can be expressed as $\theta_{0}=\mathbb{E}\left[m\left(W,\boldsymbol{\gamma}_{0}\right)\right]$ with

Hirano.and.Porter.2012 shows impossibility results for the estimation of non-differentiable functionals of the data distribution. In particular, when the target object is non-differentiable in the parameters of the data distribution, there exist no sequence of estimators that achieves local asymptotic unbiasedness. Even though $m\left(W,\boldsymbol{\gamma}\right)$ is non-differentiable in $\boldsymbol{\gamma},$ $\mathbb{E}\left[m\left(W,\boldsymbol{\gamma}\right)\right]$ is not necessarily non-differentiable. For example, if $\tau\left(Z\right)$ follows a normal distribution with mean $\mu$ and variance $\sigma^{2},$ we have

where $\Phi\left(\cdot\right)$ and $\phi\left(\cdot\right)$ are, respectively, the cdf and the probability density function (pdf) of the standard normal distribution. The target object is thus differentiable in the parameters of the data distribution. In general, if $\tau\left(Z\right)$ has a continuous density function (i.e., when the margin assumption holds), then $\theta_{0}$ will be differentiable in the parameters of the data distribution as $\mathbb{E}\left[\tau\left(Z\right)\mathds{1}\left\{ \tau\left(Z\right)>0\right\} \right]$ is proportional to the truncated mean of $\tau\left(Z\right).$ On the other hand, the target parameter can be non-differentiable when the margin assumption fails to hold. We will see in Section (ref) that valid inference depends on how $\tau\left(Z\right)$ behaves in the neighborhood $\tau\left(Z\right)=0.$

Since the $m\left(W,\boldsymbol{\gamma}\right)$ is non-differentiable in $\boldsymbol{\gamma},$ the results of Chernozhukov.et.al.2022a cannot be directly applied. To be specific, $m\left(W,\boldsymbol{\gamma}\right)$ is not Gateaux differentiable at $\boldsymbol{\gamma}=\left(c,c\right)$ for $c\in\mathbb{R}.$ To see this, consider the Gateaux differential $dm\left(\left(c,c\right);\psi\right)$ of $m$ at $\left(c,c\right)$ in the direction $\psi=\left(\psi_{1},\psi_{2}\right)$ as follows. \[ dm\left(\left(c,c\right);\psi\right)\equiv\lim_{\delta\rightarrow0}\frac{m\left(\left(c,c\right)+\delta\psi\right)-m\left(\left(c,c\right)\right)}{\delta}. \] If the limit exists for all directions $\psi,$ then $m$ is Gateaux differentiable at $\left(c,c\right).$ However, it is clear that the limit does not exist. Notice that for $\delta\neq0,$

Then,

and

The left and right limits are not the same, and hence $m\left(W,\boldsymbol{\gamma}\right)$ is not Gateaux differentiable at $\boldsymbol{\gamma}=\left(c,c\right)$ for $c\in\mathbb{R}.$

When $m\left(W,\boldsymbol{\gamma}\right)$ is nonlinear in $\boldsymbol{\gamma},$ Chernozhukov.et.al.2022a linearizes the function and extends results for the linear case to the linearized function. The key assumption is that the remainder in the linearization, employing Gateaux differentiability in the neighborhood of the true parameter, is bounded by a constant. This assumption is crucial in deriving asymptotic normality of the DML estimator in the nonlinear case as it renders the remainder term negligible.

We introduce the sigmoid function to smooth the indicator function\footnote{The existence of Gateaux derivative of a functional is guaranteed as long as we stay in $C^{1}$ class. To be specific, consider a functional $\mathscr{F}:V\rightarrow\mathbb{R}$ is given as $\mathscr{F}\left(u\right)=\int_{\Omega}F\left(x,u\left(x\right),Du\left(x\right)\right)dx,\;\Omega\subseteq\mathbb{R}^{n},$ with functions $u:\Omega\rightarrow\mathbb{R}^{n}$ contained in some open subset $V$ of a function space $U\subseteq C^{1}\left(\Omega,\mathbb{R}^{n}\right).$ Under this specification, first variation (functional differential) of $\mathscr{F}$ exists and coincides with the Gateaux derivative of $\mathscr{F}.$ Thus, whenever $u$ and the integrand $F$ are of class $C^{1},$ the Gateaux derivative exists. In our setting, since the sigmoid function is of class $C^{1},$ the existence of Gateaux derivative of $m\left(W,\boldsymbol{\gamma}\right)$ is guaranteed.}. This smoothing function is characterized by a smoothing parameter which depends on the sample size. A notable feature is that the sigmoid function can be interpreted as the cumulative distribution function (cdf) of the logistic distribution, with the smoothing parameter scaling the distribution. Therefore, it is relatively convenient to derive an analytic expression for the approximation error introduced by smoothing using the analytic properties of the logistic distribution. This explicit analytic form facilitates the selection of an optimal smoothing parameter to control the trade-off between bias and variance of our estimator. Proper application of smoothing transforms a non-differentiable function into a differentiable one, while preserving the overall shape of the original function.

Despite the differentiability of the sigmoid function, the presence of the smoothing parameter means the results from Chernozhukov.et.al.2022a do not directly carry over. Chernozhukov.et.al.2022a assumes that the remainder term from linearization is bounded by a constant. When linearizing the sigmoid function, the remainder term depends on the smoothing parameter. This dependence causes the bound of the remainder term to increase as the sample size grows, thereby precluding straightforward application of the results of Chernozhukov.et.al.2022a.

Estimation and Inference

The sigmoid function is defined as $f\left(t\right)=\frac{1}{1+\exp\left(-s_{n}t\right)}$ where $s_{n}>0$ can be interpreted as a smoothing parameter which depends on the sample size. As shown in Figure (ref), the sigmoid function approaches the indicator function as $s_{n}\rightarrow\infty.$

Let \[ m_{\mathrm{sig}}\left(W,\boldsymbol{\gamma}\right)=\frac{\gamma_{1}\left(X\right)-\gamma_{2}\left(X\right)}{1+\exp\left(-s_{n}\left(\gamma_{1}\left(X\right)-\gamma_{2}\left(X\right)\right)\right)} \] be the smoothing function of \[ m\left(W,\boldsymbol{\gamma}\right)=\left[\gamma_{1}\left(X\right)-\gamma_{2}\left(X\right)\right]\mathds{1}\left\{ \gamma_{1}\left(X\right)-\gamma_{2}\left(X\right)>0\right\} . \] Figure (ref) shows that $m_{\mathrm{sig}}\left(W,\boldsymbol{\gamma}\right)$ approaches $m\left(W,\boldsymbol{\gamma}\right)$ as $s_{n}\rightarrow\infty.$

Denote

$\overline{\theta}$ can be viewed as the true parameter of interest, and $\overline{\theta}_{\mathrm{sig}}$ as a pseudo-true parameter.

Estimators

As in Chernozhukov.et.al.2022a, we can construct a DML estimator $\hat{\theta}_{\mathrm{sig}}.$ Consider the following decomposition for any $n\in\mathbb{N}.$

Equation (ref) shows how the asymptotic distribution behaves when $\overline{\theta}$ is estimated by $\hat{\theta}_{\mathrm{sig}}.$ The first term on the right hand side $\sqrt{\frac{n}{s_{n}^{2}}}\left(\hat{\theta}_{\mathrm{sig}}-\overline{\theta}_{\mathrm{sig}}\right)$ is a random term which generates the asymptotic distribution of the estimator $\hat{\theta}_{\mathrm{sig}}$ around a pseudo-true parameter $\overline{\theta}_{\mathrm{sig}},$ and the second term $\sqrt{\frac{n}{s_{n}^{2}}}\left(\overline{\theta}_{\mathrm{sig}}-\overline{\theta}\right)$ is a nonrandom term which accounts for the error introduced by approximating the indicator function with the sigmoid function. As $\hat{\theta}_{\mathrm{sig}}$ and $\overline{\theta}_{\mathrm{sig}}$ involve $s_{n},$ both two terms are affected by the smoothing parameter $s_{n}.$ Another feature of the asymptotic behavior is that we multiply $\sqrt{\frac{n}{s_{n}^{2}}}$ instead of $\sqrt{n}.$ Kernel density estimation has a similar expression when bandwidth selection is involved. The results of Theorem (ref), presented later in this section, make inference feasible when $s_{n}$ is chosen optimally.

Following equation (ref), the DML estimator $\hat{\theta}_{\mathrm{sig}}$ is constructed as follows\footnote{In appendix, we briefly review DML where the parameter of interest depends linearly on a conditional expectation or nonlinearly on multiple conditional expectations}:

where the data $W_{i},$ $i=1,\cdots,n,$ are i.i.d., $I_{\ell},$ $\ell=1,\cdots L,$ is a partition of the observation index set $\left\{ 1,\cdots,n\right\} $ into $L$ distinct subsets of roughly equal size, and $\hat{\boldsymbol{\gamma}}_{\ell}=\left(\hat{\gamma}_{1\ell}\left(X_{1i}\right),\hat{\gamma}_{2\ell}\left(X_{2i}\right)\right)^{'}$ is the vector of regressions constructed by the observations not in $I_{\ell}.$ The estimator $\hat{\alpha}_{k\ell}\left(X_{ki}\right)$ of the Riesz representer specific to each regression is also constructed by the observations not in $I_{\ell}.$ Each $\hat{\alpha}_{k\ell}\left(X_{ki}\right)$ is obtained as follows. For each $k,$ denote $b_{k}\left(x_{k}\right)=\left(b_{k1}\left(x_{k}\right),\cdots,b_{kp}\left(x_{k}\right)\right)^{'}$ as a $p\times1$ dictionary vector specific to the $k$th regression $\gamma_{k}\left(x_{k}\right),$ and let $\hat{\boldsymbol{\gamma}}_{\ell,\ell^{'}}$ be the vector of regressions constructed by all observations not in either $I_{\ell}$ or $I_{\ell^{'}}.$ Also, let $\eta$ be a scalar, and $e_{k}$ be the $k$th column of the $2\times2$ identity matrix. Then, as in the equation (5.2) of Chernozhukov.et.al.2022a,

where $n_{\ell}$ is the number of observations in $I_{\ell},$ $b_{kj}$ is the $j$th element of the dictionary $b_{k}\left(x_{k}\right)$ as a function of $x_{k},$ and $r_{k}$ is the penalty size which must be chosen to be larger than $\sqrt{\ln\left(p\right)/n}.$

In the context of CATE estimation, this estimator can be categorized as a regularized T-learner, taking the difference between two conditional expectations and incorporating an additional debiasing correction term. Researchers may also consider alternative approaches. For example, Athey.and.Wager2019 treats the CATE function as a nuisance parameter and subsequently debiases it. It is generally difficult to argue that one method uniformly dominates another. Künzel.el.al.2019 provides comparisons of multiple learners and discusses the advantages of each method.

Our proposed estimator is an automatic DML for nonlinear functionals following Chernozhukov.et.al.2022a. Its advantage over generic DML (e.g, Chernozhukov.et.al.2018) is that it does not require prior knowledge of the explicit form of the correction term, that is, the Riesz representer. Moreover, even when a closed form is available, the generic DML approach, which first estimates the nuisance parameter such as the propensity score and then applies its analytical functional form, may not be optimal because of structural issues. For instance, to avoid numerical instability, Klosin2021 estimates continuous treatment effects using automatic DML, where the correction term is given as the multiplicative inverse of the (generalized) propensity score. In our setting, the form of the correction term depends on the smoothing parameter. This dependence makes the generic approach highly sensitive to the tuning of the smoothing parameter and may amplify estimation errors. In contrast, the automatic DML approach is designed to balance the trade-off between bias and variance associated with the smoothing parameter in an optimal manner, resulting in a more stable estimate of the correction term.

Our analysis is based on the automatic DML using a sparse linear approximation of the Riesz representer as in Chernozhukov.et.al.2022a. One may consider an adversarial approach to estimate the Riesz representer (Hirshberg.and.Wager.2021 and Chernozhukov.et.al.2024b) within a broader functional class, but adversarial learning methods incur additional computational burdens. By introducing an approximate sparse specification of the Riesz representer, we can avoid these computational challenges by controlling the mean square approximation error and using Lasso. This approach also allows the identity of the important elements in the dictionary $b$ to remain unknown while still achieving the sparse approximation rate. Such a property is particularly useful for a policy maker who does not have prior knowledge of which elements of the dictionary are most relevant to the parameter of interest.

In contrast to the framework in Chernozhukov.et.al.2024a, which considers a setting where the policy maker pre-selects a low-dimensional subset of covariates $Z_{\mathrm{sub}}$ (e.g., income level) from a high-dimensional set $Z$ and then identifies the CATE as

with $\varphi_{0}\left(Y,Z,D\right)$ being the augmented inverse propensity weighted score (Robins et al., 1994), and defines the value function as

so that the maximized average welfare gain is given by

our framework differs in an important way. In our approach, the maximized average welfare gain is defined as $\mathbb{E}\left[\mathds{1}\left(\tau\left(Z\right)\geq0\right)\tau\left(Z\right)\right]$ over the entire high-dimensional set $Z$, and the policy maker is not required to specify in advance which covariates are most relevant. The automatic DML procedure detects the relevant factors through estimation, thereby providing an alternative and flexible benchmark for policy evaluation.

Theoretical Results

We impose the following regularity conditions (Chernozhukov.et.al.2022a). For a matrix $A,$ define the norm $\left\Vert A\right\Vert _{1}=\sum_{i,j}\left|a_{ij}\right|.$ For a $p\times1$ vector $\rho,$ let $\rho_{J}$ be a $J\times1$ subvector of $\rho,$ and $\rho_{J^{c}}$ be the vector consisting of components of $\rho$ that are not in $\rho_{J}.$

This assumption restricts the convergence rate of each $\hat{\gamma}_{k}.$ This is based on the results of Newey1994, which shows that estimators which rely nonlinearly on unknown functions need to converge faster than $n^{-\frac{1}{4}}$ in terms of the norm.

This assumption is a sparse eigenvalue condition, which is generally assumed in Lasso literature (Bickel.et.al.2009, Rudelson.and.Zhou2013, and Belloni.and.Chernozhukov2013).

This assumption characterizes the Lasso penalty size $r_{k},$ and restricts the growth rate of $p$ to be slower than some power of $n.$

This assumption imposes the finite moment conditions.

This assumption restricts the convergence rate of the smoothing parameter $s_{n}.$

The next proposition characterizes the asymptotic distribution of the estimator $\hat{\theta}_{\mathrm{sig}}$ around the pseudo-true parameter $\overline{\theta}_{\mathrm{sig}}.$ We multiply by $\sqrt{\frac{n}{s_{n}^{2}}}$ instead of $\sqrt{n}$ due to the dependence of $\mathrm{Var}\left(\psi_{\mathrm{sig}}\left(w\right)\right)$ on $s_{n}.$ The asymptotic variance $V$ depends on the variance of the CATE function $\overline{\tau}\left(X\right).$

The following proposition characterizes $\overline{\theta}_{\mathrm{sig}}-\overline{\theta},$ which accounts for the approximation bias of our estimator.

Proposition (ref) shows that the behavior of $\overline{\theta}_{\mathrm{sig}}-\overline{\theta}$ depends on the distribution of $\overline{\tau}$ around the cutoff point. This is because \[ \int_{0}^{u}\overline{\tau}f_{\overline{\tau}}\left(\overline{\tau}\right)d\overline{\tau}=\mathrm{Pr}\left(0<\overline{\tau}<u\right)\mathbb{E}\left[\overline{\tau}\mid0<\overline{\tau}<u\right]. \] That is, the convergence rate of $\overline{\theta}_{\mathrm{sig}}-\overline{\theta}$ can depend on the distribution of $\overline{\tau}.$ Example (ref) shows a case where $\overline{\theta}_{\mathrm{sig}}-\overline{\theta}$ has an analytic expression if $\overline{\tau}$ follows a logistic distribution with scale parameter $\frac{1}{\lambda}$ and $\lambda=1.$

$\quad$

$\quad$

Proposition (ref) is difficult to justify in practice as the convergence rate of $\overline{\theta}_{\mathrm{sig}}-\overline{\theta}$ cannot be determined without knowledge of the distribution of $\overline{\tau}.$ Even if the distribution of $\overline{\tau}$ were known, it would still be unclear how fast $\overline{\theta}_{\mathrm{sig}}-\overline{\theta}$ converges to zero. As seen in Example (ref) knowledge of the distribution of $\overline{\tau}$ does not necessarily pin down the convergence rate of $\overline{\theta}_{\mathrm{sig}}-\overline{\theta}.$ Instead, researchers may be interested in the worst-case: the upper bound of $\left|\overline{\theta}_{\mathrm{sig}}-\overline{\theta}\right|.$ To characterize the bounds of $\left|\overline{\theta}_{\mathrm{sig}}-\overline{\theta}\right|,$ we impose an additional assumption, known as the margin assumption.

This assumption explains the behavior of the distribution $\overline{\tau}$ in the neighborhood of $\overline{\tau}=0.$ Kitagawa.and.Tetenov2018 considers the margin assumption in the context of the empirical welfare maximization to improve the convergence rate of welfare loss. Example 2.4 of Kitagawa.and.Tetenov2018 notes that, when the pdf of $\overline{\tau}\left(X\right)$ is bounded from above by $p_{\overline{\tau}}<\infty,$ the upper bound of the margin assumption is satisfied with $\alpha_{4}=1$ and $c_{4}=p_{\overline{\tau}}.$ This choice of $\alpha_{4}$ and $c_{4}$ can be considered as a benchmark. In practice, researchers need to specify or estimate $c_{4}.$ This implementation is explained at the end of this section. We impose a lower bound in the margin assumption in order to characterize the order of bias. The next proposition shows the bounds of $\left|\overline{\theta}_{\mathrm{sig}}-\overline{\theta}\right|$ provided the margin assumption holds.

Proposition (ref) provides an upper and lower bound of $\left|\overline{\theta}_{\mathrm{sig}}-\overline{\theta}\right|.$ The order of the negative bias is $\left(\frac{1}{s_{n}}\right)^{\alpha_{4}+1}.$ Given the bounds of $\left|\overline{\theta}_{\mathrm{sig}}-\overline{\theta}\right|,$ a bias-aware confidence interval can be constructed for $\overline{\theta}.$ Armstrong.and.Kolesar2020 proposes a method of constructing confidence intervals that take into account bias. Following this approach, a confidence interval can be constructed as

where $\mathrm{se}\left(\hat{\theta}_{\mathrm{sig}}\right)$ denotes the standard error, $\widehat{\overline{\mathrm{bias}}}\left(\hat{\theta}_{\mathrm{sig}}\right)$ stands for an estimate of the absolute value of the worst-case bias, which we write as $\overline{\mathrm{bias}}\left(\hat{\theta}_{\mathrm{sig}}\right),$ and $\mathrm{cv}_{1-\alpha}\left(A\right)$ is the $1-\alpha$ quantile of the folded normal distribution, $\left|\mathcal{N}\left(A,1\right)\right|.$ As Armstrong.and.Kolesar2020 points out, this confidence interval has a critical value $\mathrm{cv}_{1-\alpha}\left(\frac{\widehat{\overline{\mathrm{bias}}}\left(\hat{\theta}_{\mathrm{sig}}\right)}{\mathrm{se}\left(\hat{\theta}_{\mathrm{sig}}\right)}\right),$ which is larger than the usual normal quantile $z_{1-\frac{\alpha}{2}}.$ Correct coverage of this confidence interval can be derived from Theorem 2.2 of Armstrong.and.Kolesar2020. For notational convenience, let $\mathrm{sd}\left(\hat{\theta}_{\mathrm{sig}}\right)$ denote the standard deviation of $\hat{\theta}_{\mathrm{sig}}.$

To implement this confidence interval, it is necessary to estimate the worst-case bias and standard deviation of $\hat{\theta}_{\mathrm{sig}}.$ In Proposition (ref), the upper bound of $\overline{\theta}_{\mathrm{sig}}-\overline{\theta}$ is expressed as \[ c_{4}\left(\frac{1}{s_{n}}\right)^{\alpha_{4}+1}\pi^{\alpha_{4}+1}\left(2^{\alpha_{4}+1}-2\right)\left|B_{\alpha_{4}+1}\right|. \] The constants $c_{4}$ and $\alpha_{4}$ must be specified or estimated. The (asymptotic) standard deviation involves $\mathrm{Var}\left(\overline{\tau}\left(X\right)^{2}\right).$ As these constants are also utilized in selecting the optimal smoothing parameter, we discuss how they can be estimated after presenting the optimal smoothing parameter in the next theorem.

Theorem (ref) shows the asymptotic distribution in the worst-case scenario when the optimal smoothing parameter is chosen to balance out the trade-off between bias and variance. The asymptotic distribution exhibits negative bias as $\overline{\theta}_{\mathrm{sig}}-\overline{\theta}$ is negative. A notable feature is that, when the smoothing parameter is chosen optimally, the bias consists of constants $c_{2,\mathrm{opt}},$ $c_{4},$ and $\alpha_{4}.$ The constants $c_{4}$ and $\alpha_{4}$ comes from the upper bound of the margin assumption, and can be estimated by checking the margin assumption. As discussed earlier, when the pdf of $\overline{\tau}\left(X\right)$ is bounded from above by $p_{\overline{\tau}}<\infty,$ the upper bound of the margin assumption is satisfied with $\alpha_{4}=1$ and $c_{4}=p_{\overline{\tau}}.$ The constant $c_{2,\mathrm{opt}}$ can be viewed as a tuning parameter, and it involves $c_{4},$ $\alpha_{4},$ and $\mathrm{Var}\left(\overline{\tau}\left(X\right)^{2}\right).$ In practice, $c_{4},$ $\alpha_{4},$ and $\mathrm{Var}\left(\overline{\tau}\left(X\right)^{2}\right)$ must be specified or estimated in order to choose the tuning parameter $c_{2,\mathrm{opt}}.$

Tuning Parameter Selection

As discussed in the previous subsection, researchers need to select tuning parameters. We provide a practical way to implement our procedure.

$\;$

$\;$

Although Silverman's rule of thumb may yield inaccurate results when the true distribution significantly deviates from normality, it is straightforward to implement, and remains widely used in practice. With tuning parameters chosen based on Silverman's rule of thumb, $c_{2,\mathrm{opt}}$ is \[ c_{2,\mathrm{opt}}=\left\{ \frac{2\left[\hat{p}_{\overline{\tau}}\pi^{2}\left(2^{2}-2\right)\left|B_{2}\right|\right]^{2}}{\frac{1}{16}\widehat{\mathrm{Var}\left(\overline{\tau}\left(X\right)^{2}\right)}}\right\} ^{\frac{1}{6}} \] where $\hat{p}_{\overline{\tau}}$ and $\widehat{\mathrm{Var}\left(\overline{\tau}\left(X\right)^{2}\right)}$ are defined above. Thus, \[ s_{n}^{*}=c_{2,\mathrm{opt}}n^{\frac{1}{6}}. \] With tuning parameters chosen by the rule of thumb, the confidence interval in equation (ref) can be calculated using

Without Margin Assumption

One may consider how to choose the smoothing parameter when the margin assumption does not hold (Levis.et.al.2023). In this case, by slightly adjusting the proof of Proposition (ref), we can show that the upper bound for $\overline{\theta}_{\mathrm{sig}}-\overline{\theta}$ is characterized as \[ \frac{1}{s_{n}}2\log2, \] which is of order $\frac{1}{s_{n}}$. This rate is slower than that obtained under the margin assumption, where the bound is of order $\left(\frac{1}{s_{n}}\right)^{1+\alpha_{4}}$. Therefore, an optimal smoothing parameter in the absence of the margin condition can be chosen as \[ s_{n}^{*\;\mathrm{no}\;\mathrm{margin}}=c_{2,\mathrm{opt}}^{\mathrm{no}\;\mathrm{margin}}n^{\frac{1}{4}}, \] with \[ c_{2,\mathrm{opt}}^{\mathrm{no}\;\mathrm{margin}}=\left(\frac{\left(2\log2\right)^{2}}{\frac{1}{16}\mathrm{Var}\left(\overline{\tau}\left(X\right)^{2}\right)}\right)^{\frac{1}{4}}. \] The asymptotic distribution in the absence of the margin assumption is given by \[ \sqrt{\frac{n}{s_{n}^{2}}}\left(\hat{\theta}_{\mathrm{sig}}-\overline{\theta}\right)\overset{d}{\rightarrow}\mathcal{N}\left(-c_{3}^{\mathrm{no}\;\mathrm{margin}},\frac{1}{16}\mathrm{Var}\left(\tau\left(X\right)^{2}\right)\right) \] where \[ c_{3}^{\mathrm{no}\;\mathrm{margin}}\leq\frac{2\log2}{\left(c_{2,\mathrm{opt}}^{\mathrm{no}\;\mathrm{margin}}\right)^{2}} \] The bias-aware confidence interval in equation (ref) can also be constructed by using the optimal smoothing parameter in the absence of the margin assumption:

Thus, our smoothing methods can provide a conservative inference strategy under weaker assumptions, offering an alternative when plug-in based methods are not applicable.

Miscellaneous Estimands

We note that the construction of our DML estimator suggests alternative approaches for estimating some interesting estimands. Two examples are provided below.

Probability that CATE is positive

The proportion of individuals with a positive CATE is of interest to policy makers, as it represents the fraction of treated individuals when the optimal policy is implemented in the standard binary treatment assignment setting. Kitagawa.and.Tetenov2018 reports the share of the population to be treated in Table 1 of their paper. This quantity can be computed using our DML estimator. Recall from equation (ref) that the DML estimator is defined as \[ \hat{\theta}_{\mathrm{sig}}=\frac{1}{n}\sum_{\ell=1}^{L}\sum_{i\in I_{\ell}}\hat{\psi}_{i\ell} \] where \[ \hat{\psi}_{i\ell}\equiv m_{\mathrm{sig}}\left(W_{i},\hat{\boldsymbol{\gamma}}_{\ell}\right)+\sum_{k=1}^{2}\hat{\alpha}_{k\ell}\left(X_{ki}\right)\left[Y_{ki}-\hat{\gamma}_{k\ell}\left(X_{ki}\right)\right]. \] The term $\hat{\psi}_{i\ell}$ can be interpreted as an estimate of the debiased outcome for \[ m_{\mathrm{sig}}\left(W_{i},\boldsymbol{\gamma}\right)=\frac{\tau\left(X_{i}\right)}{1+\exp\left(-s_{n}\tau\left(X_{i}\right)\right)} \] Since the sign of $m_{\mathrm{sig}}\left(W_{i},\boldsymbol{\gamma}\right)$ is consistent with that of $\tau\left(X_{i}\right)$, the proportion of positive CATE values can be computed as the fraction of positive $\hat{\psi}_{i\ell}$. We also report this value in our empirical analysis.

Half of ATE

Throughout the paper, we let $s_{n}\rightarrow\infty$ in the estimator to derive asymptotic results. It is also interesting to examine how the estimator is constructed when $s_{n}\rightarrow0$. For the vector of covariates $Z$ and the binary treatment status indicator $D$ with $X=\left(D,Z^{'}\right)^{'}$, consider an appropriate dictionary $b\left(x\right)=b\left(d,z\right)$. First, note that \[ \lim_{s_{n}\rightarrow0}m_{\mathrm{sig}}\left(W,\boldsymbol{\gamma}\right)=\frac{1}{2}\left[\gamma_{1}\left(X\right)-\gamma_{2}\left(X\right)\right], \] so that taking the expectation yields an estimand equal to half of the ATE. Next, let us examine how the Riesz representer changes. For notational convenience, we suppress the cross-validation notation. Recall from equation (ref) that when $s_{n}\rightarrow0$, the components $\hat{M}_{1j}$ and $\hat{M}_{2j}$ are given by

Thus, we have

If we define \[ m_{\mathrm{ATE},\mathrm{half}}\left(W,\boldsymbol{\gamma}\right)\equiv\frac{1}{2}\left[\gamma_{1}\left(X\right)-\gamma_{2}\left(X\right)\right] \] and compute its Gateaux derivative with respect to the dictionary, we obtain the equivalents results:

This is because the moment function becomes linear when $s_{n}\rightarrow0$. In other words, the target moment function coincides with its own Gateaux derivative. This observation shows that in the limit $s_{n}\rightarrow0$, the DML estimator targets half of the ATE. In this linear case one can then rely on the automatic DML construction for linear functionals as described in Chernozhukov.et.al.2022a.

Alternative Smoothing Function

Our target parameter is defined as the expectation of the moment function \[ m\left(W,\boldsymbol{\gamma}\right)=\left(\gamma_{1}\left(X\right)-\gamma_{2}\left(X\right)\right)\mathds{1}\left\{ \gamma_{1}\left(X\right)-\gamma_{2}\left(X\right)>0\right\} \] where the moment function is equivalent to $\max\left\{ \gamma_{1}\left(X\right)-\gamma_{2}\left(X\right),0\right\} $. In our approach, we smooth the indicator function by using a sigmoid function, thereby obtaining the smoothed moment function \[ m_{\mathrm{sig}}\left(W,\boldsymbol{\gamma}\right)\equiv\frac{\gamma_{1}\left(X\right)-\gamma_{2}\left(X\right)}{1+\exp\left(-s_{n}\left(\gamma_{1}\left(X\right)-\gamma_{2}\left(X\right)\right)\right)}. \] Alternatively, one may smooth the maximum directly via the log-sum-exp (LSE) function (Levis.et.al.2023). In that case the smoothing function is defined as \[ m_{\mathrm{LSE}}\left(W,\boldsymbol{\gamma}\right)\equiv\frac{1}{h_{n}}\log\left(1+\exp\left(h_{n}\left(\gamma_{1}\left(X\right)-\gamma_{2}\left(X\right)\right)\right)\right), \] where $h_{n}$ plays the same role as the smoothing parameter in our approach. Notably, the Gateaux derivative of $m_{\mathrm{LSE}}\left(W,\boldsymbol{\gamma}\right)$ in the direction of the true treatment effect difference is exactly $m_{\mathrm{sig}}\left(W,\boldsymbol{\gamma}\right)$. Therefore, when assessing the approximation error introduced by smoothing, both the sigmoid‐based and LSE‐based approaches are fundamentally linked to the logistic distribution and exhibit equivalent theoretical properties.

On the other hand, one can observe that \[ m_{\mathrm{sig}}\left(W,\boldsymbol{\gamma}\right)\leq m\left(W,\boldsymbol{\gamma}\right)\leq m_{\mathrm{LSE}}\left(W,\boldsymbol{\gamma}\right), \] with the equalities holding at the cutoff point. As a result, estimates based on the sigmoid smoothing are likely to be smaller than those based on the LSE smoothing. Which smoothing method to adopt can ultimately depend on the policy maker's preference. For example, if one wishes to avoid overestimating the welfare gain, a conservative policy maker may choose the sigmoid function. Furthermore, a useful by-product of the sigmoid‐based approach is that it enables the computation of the proportion of individuals with a positive conditional average treatment effect by leveraging sign consistency, as discussed in the previous subsection.

Simulation

We provide simulation results for a process where $\overline{\tau}\left(X\right)$ follows a logistic distribution with mean 0 and variance 1. The data generating process is as follows. Consider covariates $X=\left(X_{1},\cdots,X_{\frac{p_{0}}{2}},X_{\frac{p_{0}}{2}+1},\cdots,X_{p_{0}},X_{p_{0}+1},\cdots,X_{p}\right)$ where each $X_{j}$ is an i.i.d. exponential random variable with rate parameter $\lambda_{j}=\frac{2}{p_{0}}$ for $j=1,\cdots,p_{0}.$ Here, $p_{0}$ controls sparsity and is set as 6. It can be easily verified\footnote{$\mathrm{Pr}\left(\min\left\{ X_{1},\cdots,X_{\frac{p_{0}}{2}}\right\} \geq x\right)=\mathrm{Pr}\left(X_{1}\geq x,\cdots,X_{\frac{p_{0}}{2}}\geq x\right)=\mathrm{Pr}\left(X_{1}\geq x\right)\times\cdots\times\mathrm{Pr}\left(X_{\frac{p_{0}}{2}}\geq x\right).$ Since each $X_{j}$ is i.i.d. exponential random variable with the rate parameter $\lambda_{j}=\frac{2}{p_{0}},$ we obtain $\mathrm{Pr}\left(\min\left\{ X_{1},\cdots,X_{\frac{p_{0}}{2}}\right\} \geq x\right)=\exp\left(-\left(\frac{2}{p_{0}}\times\frac{p_{0}}{2}\right)x\right)=\exp\left(-x\right).$ This immediately implies $\min\left\{ X_{1},\cdots,X_{\frac{p_{0}}{2}}\right\} \sim\mathrm{Exp}\left(1\right).$} that

Potential outcomes are set to

where $\epsilon_{1}\sim\mathscr{N}\left(0,0.1^{2}\right)$ and $\epsilon_{2}\sim\mathscr{N}\left(0,0.1^{2}\right)$ are independent of $X.$ From properties of the exponential and logistic distributions\footnote{A quick computation shows $\mathrm{Pr}\left(\frac{S_{1}}{S_{2}}\leq x\right)=\frac{x}{x+1}.$ Note that the log function is strictly increasing, and its inverse function is the exponential function. Thus, $\mathrm{Pr}\left(\ln\left(\frac{S_{1}}{S_{2}}\right)\leq x\right)=\frac{e^{x}}{e^{x}+1},$ which is the cdf of the logistic distribution.}, we have \[ \ln\left(\frac{S_{1}}{S_{2}}\right)\sim\mathrm{Logistic}\left(0,1\right) \] when $S_{1}$ and $S_{2}$ are i.i.d. exponential random variables with rate parameter 1. The CATE function $\overline{\tau}\left(X\right)$ is then

The true parameter is

where a detailed derivation is included in Appendix.

When running the simulation (and also analyzing empirical data), there are three major tuning parameters and a dictionary which must be selected. For the choice of dictionary, we consider four specifications as follows:

Specification (1): Includes an intercept, six base covariates, and squared terms for the base covariates. The dimension of the dictionary is 13.

Specification (2): Extends Specification (1) by adding all first-order interaction terms and cubic terms for the base covariates. The dimension of the dictionary is 34.

Specification (3): Extends Specification (2) by adding the fourth- and fifth- and sixth-order terms for the base covariates, and six normal random error terms. The dimension of the dictionary is 58.

The sample size is $n=2,000$ and the iteration number is $1,000$ for all specifications. The first tuning parameter is the penalty degree for estimating the conditional expectation $\boldsymbol{\gamma}\left(X\right).$ When the conditional expectation is estimated by Lasso, Chernozhukov.et.al.2022a provides theoretical justification for choosing the penalty degree that results in the fastest possible mean square convergence rate, which produces the optimal trade-off between bias and variance. We choose the penalty parameter as $\sqrt{\frac{\ln\left(p+1\right)}{n}}$ where $p+1$ is the dimension of the dictionary and $n$ is the sample size. The second tuning parameter is $r_{k}$ in equation (ref) for estimating $\hat{\rho}_{k\ell}.$ Chernozhukov.et.al.2022a argues that this parameter must be larger than $\sqrt{\frac{\ln\left(p+1\right)}{n}}$ when $m\left(w,\boldsymbol{\gamma}\right)$ is not linear on $\boldsymbol{\gamma}.$ They propose choosing $r_{k}$ to be proportional to $n^{-\frac{1}{4}}$ and we set the $r_{k}$ as $n^{-\frac{1}{4}}.$ The third tuning parameter is the optimal smoothing parameter $s_{n}^{*}=c_{2}n^{\frac{1}{2\left(\alpha_{4}+2\right)}}.$ In this example, the pdf of the CATE function is bounded from above by $p_{\overline{\tau}}<\infty.$ Hence, the margin assumption is satisfied with $\alpha_{4}=1.$ $c_{2}$ is chosen as \[ c_{2,\mathrm{opt}}=\left\{ \frac{2\left[p_{\overline{\tau}}\pi^{2}\left(2^{2}-2\right)\left|B_{2}\right|\right]^{2}}{\frac{1}{16}\mathrm{Var}\left(\overline{\tau}\left(X\right)^{2}\right)}\right\} ^{\frac{1}{6}} \] where $p_{\overline{\tau}}=0.25$ and $\mathrm{Var}\left(\overline{\tau}\left(X\right)^{2}\right)=\frac{16}{45}\pi^{4}$ when $\overline{\tau}\left(X\right)\sim\mathrm{Logistic}\left(0,1\right).$

Figure (ref) shows the sampling distribution of three estimators. The red dashed line represents the true parameter $\ln2.$ The first estimator is the DML estimator $\hat{\theta}_{\mathrm{sig}}$ with the optimal tuning parameter $c_{2,\mathrm{opt}}.$ The second estimator is a naive estimator $\hat{\theta}_{\mathrm{naive}}$ defined as

Notice that $\hat{\theta}_{\mathrm{naive}}$ is a sample analogue estimator of $\overline{\theta}$ with neither debiasing nor cross-fitting. As discussed in Section (ref), $\hat{\theta}_{\mathrm{naive}}$ may exhibit large biases when $\widehat{\tau\left(X\right)}$ entails regularization and/or model selection. (Chernozhukov.et.al.2017, Chernozhukov.et.al.2018, and Chernozhukov.et.al.2022c). On the other hand, the DML estimator $\hat{\theta}_{\mathrm{sig}}$ involves negative bias which can be controlled along with variance. The third estimator is the maximum bias DML (MB-DML) estimator, $\hat{\theta}_{\mathrm{sig}}+\hat{c}_{3,\mathrm{max}},$ where $\hat{c}_{3,\mathrm{max}}$ is the estimate of the worst-case bias $c_{3}.$ As expected, the DML estimator $\hat{\theta}_{\mathrm{sig}}$ shows negative bias, and the naive estimator $\hat{\theta}_{\mathrm{naive}}$ produces large bias. For the third estimator $\hat{\theta}_{\mathrm{sig}}+\hat{c}_{3,\mathrm{max}},$ with an estimate of maximal bias plugged in, the center of the distribution for the MB-DML estimator $\hat{\theta}_{\mathrm{sig}}+\hat{c}_{3}$ is above the true parameter. This is consistent with the bias bound presented in Theorem (ref) being the worst-case. By adding an estimate of this worst-case bound, we over adjust when the true bias is less than the worst-case.

Table (ref) shows the Monte Carlo bias, standard error (SE), and root-mean-square error (RMSE), as well as the coverage rate in Specification (1). The confidence level is 0.95, and the coverage rate of the DML estimator $\hat{\theta}_{\mathrm{sig}}$ is around 95%. The bias-aware confidence interval uses a larger critical value in order to take into account the bias.

Figures (ref) and (ref) present quantile-quantile plots (Q-Q plots) for the naive estimator and the DML estimator in Specification (1). Both Q-Q plots show relatively $45^{\circ}$ straight lines. However, the naive estimator is not valid for inference because of its large bias, as the literature has consistently pointed out.

Table (ref) presents the results for various dictionary specifications. As expected, the standard error increases when irrelevant terms are included. This suggests that the efficiency of the estimator can be improved when a policymaker has some knowledge of which factors are important for the target parameter.

Finally, we present similar results based on the LSE-smoothing method. Table (ref) displays the outcomes using the LSE-based smoothing function introduced in Section (ref). Our findings indicate that both methods exhibit equivalent performance, although the LSE-based approach tends to produce higher estimates than the sigmoid-based approach, as discussed in Section (ref).

Empirical Analysis

We apply our method to experimental data from the National Job Training Partnership Act (JTPA) Study, and predominantly follow the empirical strategies of Kitagawa.and.Tetenov2018. The sample consists of 9,223 observations. There are the outcome variable (income), and a binary treatment (assignment to a job training program). Also, there are 5 base covariates: age, education, black indicator, Hispanic indicator, and earnings in the year prior to the assignment (pre-earnings). Kitagawa.and.Tetenov2018 only uses two covariates: education and pre-earnings in the context of the valid empirical welfare maximization (EWM) method. Our target parameter can be viewed as average welfare gain under the optimal treatment assignment rules, and we use more covariates to exploit an appealing feature of our method. For the choice of dictionary, we consider four specifications as follows:

Specification (1): Includes an intercept, five base covariates, squared terms for age, education, and pre-earnings, as well as first-order interaction terms for all base covariates. The dimension of the dictionary is 19.

Specification (2): Includes an intercept, five base covariates, and squared and cubic terms for age, education, and pre-earnings. The dimension of the dictionary is 12.

Specification (3): Extends Specification (2) by adding quadratic terms for age and education. The dimension of the dictionary is 14.

Specification (4): Extends Specification (3) by adding all first-order interaction terms and the fifth- and sixth-order terms for age and education. The dimension of the dictionary is 28.

The covariates are standardized. Tuning parameters are selected by rule-of-thumb as described in Section (ref).

Table (ref) summarizes the estimation results. The confidence interval widens as we include higher-order terms. In Kitagawa.and.Tetenov2018, the corresponding estimate is \$1,340 with 95% CI (\$441, \$2,239) for the EWM quadrant rule, \$1,364 with 95% CI (\$398, \$2,330) for the EWM linear rule, and \$1,489 with 95% CI (\$374, \$2,603) for the EWM linear rule with squared and cubic terms for education. Our confidence intervals broadly align with these values. Additionally, Kitagawa.and.Tetenov2018 reports that the share of the population to be treated ranges between 0.88 and 0.96, depending on their EWM rules, which is also consistent with our results.

Conclusion

This paper focuses on debiased machine learning when nuisance parameters appear in indicator functions and there is a high-dimensional vector of covariates. We propose a DML estimator where the indicator function is smoothed. The asymptotic distribution theory demonstrates that an optimal choice of the smoothing parameter enables standard inference by balancing the trade-off between squared bias and variance. Simulations and empirical exercise corroborate these results.

There are several ways in which the proposed procedure could be developed further. The effectiveness of the proposed procedure relies significantly on the nature of non-differentiable and smoothing functions. The class of non-differentiable functions is large, and formulating a general theory for DML for non-differentiable functions is not straightforward. In addition, it may be possible to construct a tighter confidence. Finally, a formal coverage guarantee for a feasible procedure with estimated bias and variances has yet to be established.