Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
82,697 characters · 18 sections · 8 citation commands
Inference on Welfare and Value Functionals under Optimal Treatment Assignment
\global\long\global\long\global\long\global\long\global\long\global\long\global\long\global\long\global\long\global\long\global\long\global\long\global\long\global\long\global\long\global\long\global\long\global\long\global\long\global\long\global\long\global\long\global\long\global\long\global\long\global\long\global\long\global\long\global\long\global\long\global\long\global\long\global\long\global\long\global\long\global\long\global\long\global\long\def\mc#1{\mathscr{#1}} \global\long\global\long\global\long\global\long\global\long\global\long\global\long\global\long\global\long\global\long\global\long\global\long\global\long\global\long\global\long\def\abs#1{\left|#1\right|} \global\long\def\norm#1{\left\Vert #1\right\Vert } \global\long\def\rest#1{\left.#1\right|} \global\long\def\bracket#1#2{\left\langle #1\middle\vert#2\right\rangle } \global\long\def\sandvich#1#2#3{\left\langle #1\middle\vert#2\middle\vert#3\right\rangle } \global\long\def\turd#1{\frac{#1}{3}} \global\long\global\long\def\sand#1{\left\lceil #1\right\vert } \global\long\def\wich#1{\left\vert #1\right\rfloor } \global\long\def\sandwich#1#2#3{\left\lceil #1\middle\vert#2\middle\vert#3\right\rfloor } \global\long\def\abs#1{\left|#1\right|} \global\long\def\norm#1{\left\Vert #1\right\Vert } \global\long\def\rest#1{\left.#1\right|} \global\long\def\inprod#1{\left\langle #1\right\rangle } \global\long\def\ol#1{\overline{#1}} \global\long\def\ul#1{#1} \global\long\def\td#1{\tilde{#1}} \global\long\def\bs#1{\boldsymbol{#1}} \global\long\global\long\global\long\global\long\global\long{6pt} {6pt}
In this paper, we study the estimation and inference on welfare and value functionals of a given treatment under optimal (“first-best”) treatment assignment.
Let $D_{i}\in\left\{ 0,1\right\} $ denote a certain binary treatment for subject $i$, $\left(Y_{i}\left(0\right),Y_{i}\left(1\right)\right)$ denote the potential outcomes of interest, and $Y_{i}:=Y_{i}\left(D_{i}\right)\in \mathbb{R}$ denote the observed outcome. Let $X_{i}\in \mathbb{R}^d$ denote subject $i$'s observable characteristics of $i$, which is distributed with density $f_0$. We suppose that researchers have access to a random sample of training data $\left\{(D_{i},Y_{i},X_{i})\right\}_{i=1}^{n}$.
Under the standard conditional unconfoundedness assumption $\rest{\left(Y_{i}\left(0\right),Y_{i}\left(1\right)\right)\perp D_{i}}X_{i}$ and the overlap condition $p_0\left(x\right):=\mathbb{E}\left[\rest{D_{i}}X_{i}=x\right]\in\left(0,1\right)$, the conditional average treatment effect (CATE) defined by \[ \text{CATE}\left(x\right):=\mathbb{E}\left[\rest{Y_{i}\left(1\right)-Y_{i}\left(0\right)}X_{i}=x\right] \] is identified from data by
where $h_0:\mathbb{R}^d \mapsto \mathbb{R}$, and
is the nonparametric regression function of the outcome $Y_{i}$ on $X_{i}$ and $D_{i}$. We maintain the conditional unconfoundedness assumption and the overlap condition, and will thereafter simply refer to $h_{0}$ as the CATE function. We further assume that $h_{0}$ belongs to a Holder class of functions with smoothness $s>1$.
We consider a standard scenario where policymakers can assign treatments based on covariates, and focus on the following two core types of welfare and value parameters. The first type is the maximized welfare of the target population under optimal treatment assignment:
where $f$ is the marginal density of $x$ in the target population and $\left[t\right]_{+}:=\max\left(t,0\right)$ is the rectified linear unit (ReLU) function. $W\left(h_{0}\right)$ averages the CATE over the population under the “first-best” treatment assignment rule: “treat type $x$ if and only if $\text{CATE}\left(x\right)\geq0$,” and is thus often referred to as the welfare under optimal treatment assignment. We will thereafter refer to $W\left(h_{0}\right)$ as the welfare functional.
Alternatively, one may also be interested in evaluating the average of a value other than CATE over the population under optimal treatment assignment:
where $v_0:\mathbb{R}^d \mapsto \mathbb{R}$ is a user-defined function that may be a utility function, a cost function, or any economically meaningful function of the observed covariate $x$. For example, setting $v_0\left(x\right)=a'x$ endows $V\left(h_{0}\right)$ with the interpretation as certain aggregate characteristics of the treated population under optimal treatment assignment, and setting $v_0\equiv1$ implies that $V\left(h_{0}\right)\equiv\mathbb{P}_{f}\left(h_{0}\left(X_{i}\right)\geq0\right)$ becomes the share of the target population to be treated (i.e. with nonnegative CATE).
In the formulation of $W\left(h_{0}\right)$ and $V\left(h_{0}\right)$ above, we take the density $f\left(x\right)$ to be known. This is in itself relevant in settings where the covariate density of the target population is configured or known/estimated from other sources than the training sample used to estimated CATE. For example, CATE may be estimated from a smaller pilot program, while the policymakers are contemplating to implement the policy on a statewide or nation-wide basis with a much larger population. That said, in this paper we also consider an important case where $f$ is unknown and set to be the covariate density in the underlying population of the training sample. In this case, $f$ does not need to be estimated, as the integral with respect to $f$ can be naturally approximated via sample average in the training sample. This corresponds more closely to the “empirical welfare” as considered in \citet*{kitagawa2018should}. We provide results for this setting as well.
In this paper, we establish inference results for the welfare and value functionals $W\left(h_{0}\right)$ and $V\left(h_{0}\right)$ with nonparametric estimated CATE $h_0$. The results can be summarized informally as follows.
For the welfare functional $W\left(h_{0}\right)$, we establish semiparametric plug-in estimators of the welfare functional are asymptotically normal at the parametric $\sqrt{n}$ rate and derive closed-form asymptotic variance formulas, along with consistent asymptotic variance estimators. The key insight of the $\sqrt{n}$ rate is geometric: by the definition of the welfare functional, the integrand $h_{0}$$\left(x\right)$ vanishes on the boundary of the integration $\left\{ x:h_{0}\left(x\right)=0\right\} $, neutralizing the non-smoothness of the indicator function $\mathbf{\mathbbm1}\left\{ h_{0}\left(x\right)\geq0\right\} $.
In contrast, for the value functional $V\left(h_{0}\right)$ with general weight $v_{0}$ that does not vanish on the boundary $\left\{ x:h_{0}\left(x\right)=0\right\} $, we show that the rate of convergence is slower than $\sqrt{n}$, and is instead given by the 1-dimensional nonparametric regression rate $n^{-\frac{s}{2s+1}}$ under appropriate conditions. We establish asymptotic normality of the semiparametric plug-in estimator under this irregular convergence rate, and provide a consistent variance estimator based on the sieve Riesz representer. In particular, the consistent variance estimator features a Hausdorff integral on the boundary submanifold $\left\{ x:h_{0}\left(x\right)=0\right\} $, for which we provide a numerical integration and differentiation procedure for the computation of submanifold integrals.
We conduct an array of Monte Carlo experiments to document the good finite-sample accuracy of our theoretical inferential results. We show that the proposed standard error estimators perform well in finite sample, and the corresponding confidence intervals based on the asymptotic normality result and the standard error estimators have coverage probabilities close to their nominal levels. These findings hold not only for the $\sqrt{n}$-estimable welfare functional, but also for the value functional, which is estimated at slower-than-$\sqrt{n}$ rate with standard errors computed through numerical integration and differentiation.
We also apply our results to empirical data from the Job Training Partnership Act (JTPA) data set. Following \citet*{kitagawa2018should}, we take 30-month post-program earning as the outcome variable and consider two covariates: pre-program earning and education. We then provide empirical estimates and confidence intervals for two parameters: the welfare under first-best treatment assignment, which is $\sqrt{n}$ estimable, and the share of population to be treated under first-best treatment assignment, which is not $\sqrt{n}$-estimable. As in \citet*{kitagawa2018should}, we also consider two different scenarios: one with the cost of the treatment incorporated, and one without. These parameters have also been estimated in \citet*{kitagawa2018should} under the label of “nonparametric plug-in rule” using kernel first-stages, but \citet*{kitagawa2018should} only provide point estimates with no confidence intervals for them. We use sieve (B-spline) first-stage nonparametric estimators and find similar results to those in \citet*{kitagawa2018should}, and further provide informative confidence intervals for both the welfare and the share parameters.
This work is a companion paper to \citet*{chen2025semiparametric}, and directly applies the general theoretical results there to handle differentiation with respect to region of integration and semiparametric estimation of integrals over submanifolds. Specifically, this work focuses on the important context of treatment assignment problems, and deal with two features that are not discussed in \citet*{chen2025semiparametric}. First, CATE is defined as the difference of two nonparametric regression functions between the treated and untreated subpopulation (or the difference of two point evaluations of a nonparametric function with treatment status defined as an argument too), and is often estimated as the difference of two first-stage nonparametric estimators. Second, in treatment assignment problems, researchers may face two different scenario in terms of the distribution of the covariates: sometimes, this distribution can be treated as known and may be different from the covariate distribution in the experimental/observational population from which CATE is estimated; in other times, one may want to treat the covariate distribution as unknown and the same as the experimental/observational population, and uses sample averages to automatically incorporate the covariate distribution. In this paper, we takes into account these special structures of the problem, and establish the inference results by providing lower-level sufficient conditions to the general theory in \citet*{chen2025semiparametric}.\footnote{Two concurrent papers by cattaneo2025dist,cattaneo2025loc also feature submanifold integrals, but focuses on 1-dimensional cases that arise from the specific context of boundary discontinuity designs. Our current paper also differs significantly from cattaneo2025dist,cattaneo2025loc, who focus on boundary discontinuity designs and boundary treatment effects, which are very different objects from the welfare and value functionals considered here under first-best treatment assignments. Consequently, the submanifolds (boundaries) in their settings are given or known based on the locations or a distance function, while in our current project the boundary submanifold is defined by the unknown and estimated CATE function.}
Our paper makes new contributions to the literature on estimation and inference on a general value functional of a policy under optimal treatment assignment. To the best of our knowleadge, all the existing work on limiting distributions of functionals of a policy under optimal treatment assignments considered $\sqrt{n}$-normality only. Our paper is the first to establish slower-than-root-n limiting distribution and inference results for irregular value functionals of a policy under optimal treatment assignment. Previously, \citet*{bhattacharya2012inferring} establishes $\sqrt{n}$-normality of the optimal welfare value under budget constraint, using Nadaraya-Waston kernel estimator for first-stage estimation of CATE. The important work of \citet*{kitagawa2018should} focuses on a related but slightly different topic: empirical welfare maximization within a constrained class of policy rules. That said, \citet*{kitagawa2018should} also considers and reports empirical estimates on the nonparametric plug-in rule, which can be interpreted as a plug-in estimator of the first-best welfare, though there were no theoretical results or confidence intervals for this estimate. The recent papers by \citet*{park2024debiased} and \citet*{whitehouse2025inference} use soft-max functions to smooth the max function, and apply the debiased machine learning approach to establish asymptotic normality and provide inferential results, and \citet*{whitehouse2025inference} establishes $\sqrt{n}$ asymptotic normality of their soft-max welfare functional. \citet*{feng2024statistical} analyzes binary treatment assignment under constraints and the asymptotic property of the welfare of the policy under optimal cutoff choice, and establishes root-$n$ asymptotic normality of the functionals. Our paper complements these existing work in several ways. First, we clarify that the $\sqrt{n}$-estimability of the welfare functional is due to the fact that the integrand (CATE) by construction vanishes on the boundary of the optimally treated population $\left\{ x:\text{CATE}\left(x\right)=0\right\} $, which is specific to the welfare functional but generally not satisfied for other types of value functionals. Second, we demonstrate that the $\sqrt{n}$-estimability of the welfare functional can be attained without the use of smoothing/soft-max functions. Third and most importantly, our results extend well beyond the welfare functional, and cover generic value functionals that may be slower than $\sqrt{n}$-estimable.
The rest of the paper is organized as follows. Section (ref) lays out the main model setup, and provides a conceptual explanation of why the welfare functional $W\left(h_{0}\right)$ can be $\sqrt{n}$-estimable while the value functional $V\left(h_{0}\right)$ is not in general. Section (ref) then establishes the inference results for the welfare functional $W\left(h_{0}\right)$, while Section (ref) provides corresponding results for the value function $V\left(h_{0}\right)$. We report numerical results from Monte Carlo simulations in Section (ref), and conduct an empirical illustration in Section (ref). Proofs of theoretical results in the main text are available in Appendix (ref).
We first introduce the standard treatment effect model, the general welfare functional of interest and the maintained assumptions in Subsection (ref)
We first state the maintained assumption in the paper. Let $\hat{h}\left(x\right):=\hat{\mu}\left(x,1\right)-\hat{\mu}\left(x,0\right)$ be a nonparametric estimator of the CATE function, in which $\hat{\mu}\left(x,d\right)$ is a first-stage estimator of the nonparametric regression model
We impose the following basic assumptions in this paper.
We first provide a heuristic overview of our theoretical analysis, and explain the key intuition why the welfare functional $W\left(h_{0}\right)$ may be $\sqrt{n}$-estimable (i.e. $W$ is a regular functionial) while $V\left(h_{0}\right)$ is generally not (i.e., $V$ is a irregular functional).
We first introduce a more general value functional $\Phi \left(h_{0}\right)$ that nests $W\left(h_{0}\right)$ and $V\left(h_{0}\right)$ as special cases:
where $\phi:\mathbb{R}\times \mathbb{R}^d \to\mathbb{R}$ is a known measurable mapping. We note that
By the standard semiparametric theory on the estimation of functionals of nonparametric regression functions, the asymptotic property of the semiparametric plug-in estimator $\Phi\left(\hat{h}\right)$ can be analyzed via the functional derivative of $\Phi\text{\ensuremath{\left(h\right)}}$ w.r.t. $h_0$ in the direction of $h-h_{0}$, i.e., writing $h_{t}:=h_{0}+t\left(h-h_{0}\right)$,
The presence of the ReLU/max function $\left[t\right]_{+}$ in $D_{h}W\left(h_{0}\right)$ and the indicator function $\mathbf{\mathbbm1}\left\{ t\geq0\right\} $ induces a point of nonsmoothness at $t=0$, where $\left[t\right]_{+}$ is nondifferentiable and $\mathbf{\mathbbm1}\left\{ t\geq0\right\} $ is discontinuous. This complicates the calculation of the functional derivatives in (ref) and (ref), though to different degrees.
We now provide an overview of the key difference between the welfare and value functions from the perspective of the generalized Leibniz rule, which has the following generic form regarding the total time derivative of integrals with changing integrand and changing region of integration:
The generalized Leibniz rule,\footnote{See, for example, Theorem 4.2 of delfour2001shapes.} states that, under mild regularity conditions,
where term (I) captures the effect of the change in the integrand $G_{t}\left(x\right)$ with the region of integration $\Omega_{t}$ held fixed, while term (II) captures the effect of the change in the region of integration $\Omega_{t}$ with the integrand $G_{t}\left(x\right)$ held fixed. The somewhat “nonstandard” term (II) warrants some more explanations: $\partial\Omega_{t}$ denotes the boundary of $\Omega_{t}$, ${\bf n}_{t}\left(x\right)$ is the outward-pointing unit normal vector, ${\bf v}_{t}\left(x\right)$ is the velocity vector associated with the time movement in the $\partial\Omega_{t}$, and $S_{\partial\Omega_{t}}\left(x\right)$ denotes the surface measure on the boundary $\partial\Omega_{t}$. Note that, when $x$ is one-dimensional and $\Omega_{t}=\left[a_{t},b_{t}\right]$, (ref) specializes to the standard Leibniz rule: \[ \frac{d}{dt}\int_{a_{t}}^{b_{t}}G_{t}\left(x\right)dx=\int_{a_{t}}^{b_{t}}\frac{\partial}{\partial t}G_{t}\left(x\right)dx+G_{t}\left(b_{t}\right)\frac{d}{dt}b_{t}-G_{t}\left(a_{t}\right)\frac{d}{dt}a_{t}. \]
We observe that $D_{h}\Phi\left(h_{0}\right)$, $D_{h}W\left(h_{0}\right)$ and $D_{h}V\left(h_{0}\right)$ are all of the form (ref), with the same parametrized region of integration \[ \Omega_{t}:=\left\{ x\in \mathbb{R}^d:h_{t}\left(x\right)\geq0\right\}~,~~\Omega_{0}:=\left\{ x\in \mathbb{R}^d:h_{0}\left(x\right)\geq0\right\} \] and $\partial\Omega_{0}= \left\{ x\in \mathbb{R}^d:h_{0}\left(x\right)=0\right\} $ under mild regularity conditions on $h_0$.
Applying the generalized Leibniz rule (ref) to the functional $\Phi(h)$, we obtain:
where the first term (I) is a full-dimensional Lebesgue integral (in $ \mathbb{R}^d$), while the second term (II) is a lower-dimensional boundary integral. In particular, when $\phi\left(h_{0}(x),x\right) f(x)$ does not vanish on $\partial\Omega_{0}=\left\{ x\in \mathbb{R}^d:h_{0}\left(x\right)=0\right\}$, the second term (II) cannot be ignored, despite $\partial\Omega_{0}$ has Lebesgue measure zero in $ \mathbb{R}^d$. In fact, under Assumption (ref)(c) (see, e.g., \citet*{chen2025semiparametric}), the second term (II) of (ref) can be expressed as
where ${\cal H}^{d-1}$ denotes the $(d-1)$ dimensional Hausdorff measure (see \citet*{chen2025semiparametric}), which coincides with the $(d-1)$ dimensional Lebesgue measure in $\mathbb{R}^{d-1}$. Thus the boundary integral term (II) of (ref) is a lower-dimensional integral functional that only extracts information about $h_{0}$ on a Lebesgue measure-$0$ set (in $\mathbb{R}^d$).
For the welfare function $\Phi(h_0)=W(h_0)$, plugging $\phi(h_0(x),x) = h_0(x)$ into (ref) yields
where the term (II) vanishes since $\phi(h_0(x),x) = h_0(x) = 0$ on the boundary $\partial\Omega_{0}$, and the term (I) is a full-dimensional Lebesgue integral functional of $h-h_{0}$.
For the value function $\Phi(h_0)=V(h_0)$, plugging $\phi(h_0(x),x) = v_0(x)$ with $\partial_1 \phi\left(h_{0}(x),x\right)=0$ into (ref) and (ref) yields
which is not $0$ as long as $v_0(x)f(x)$ does not vanish on the boundary $\partial\Omega_{0}= \left\{ x\in \mathbb{R}^d:h_{0}\left(x\right)=0\right\}$. For example, setting $v_0\left(x\right)\equiv1$ yields $V\left(h_{0}\right)=P_{f}\left(h_{0}\left(X_{i}\right)\geq0\right)$, the share of population with nonnegative CATE, and the boundary integral term (II) does not vanish. In fact, as long as $v_0(x)f(x)\neq 0$ on the boundary $\partial\Omega_{0}= \left\{ x\in \mathbb{R}^d:h_{0}\left(x\right)=0\right\}$, $D_{h}V\left(h_{0}\right)$ is a non-zero $(d-1)$ dimensional integral functional that only extracts information about $h_{0}$ on a Lebesgue measure-$0$ set (in $\mathbb{R}^d$), akin to a point evaluation of a nonparametric estimation.
Let $L^2(f)$ denote the Hilbert space of square integrable (against $f$) functions with the inner product $\left\langle g,h\right\rangle _{2,f}:=\int g\left(x\right)h\left(x\right)f\left(x\right)dx$. For the CATE function $h_0 \in L^2(f)$, it is well-known that a linear functional $L\left[h-h_{0}\right]$ is bounded (or equivalently, continuous) if and only if \[ \sup_{\nu \neq 0, \nu \in L^2 (f)}\frac{|L\left[\nu (\cdot)\right]|^2}{E_f[|\nu (X)|^2]}<\infty \] which is a necessary and sufficient condition for the existence of a Riesz representer $\nu^*\in L^2(f)$ such that \[ L[\nu]=\left\langle \nu^*,\nu\right\rangle _{2,f}~~~\text{for all}~\nu \in L^2(f) \] This in turn is a necessary condition for any plug-in estimator of the linear functional $L\left[\hat{h}-h_{0}\right]=\left\langle \nu^*,\hat{h}-h_{0}\right\rangle _{2,f}$ to converge to zero at a root-$n$ rate.
For the welfare functional, its linear directional derivative functional $L[\nu]=D_{h}W\left(h_{0}\right)\left[\nu\right]$ given in (ref), we immediately see that $\nu^{*}\left(x\right):=\mathbf{\mathbbm1}\left\{ h_{0}\left(x\right)\geq0\right\} $ is the Riesz representer of the linear functional $D_{h}W\left(h_{0}\right)\left[\nu\right]$ in the Hilbert space $L^2(f)$, and that this Riesz representer has bounded norm \[ \norm{\nu^{*}}^{2}:=\int\mathbf{\mathbbm1}^{2}\left\{ h_{0}\left(x\right)\geq0\right\} f\left(x\right)dx\leq1, \] and thus, by well-known results in, say, \citet*{chen2014sieveIrregular}, \citet*{chen2014sieve} and chenpouzo2015sieve, the linear functional $D_{h}W\left(h_{0}\right)\left[\nu\right]$ is a regular (i.e., $\sqrt{n}$-estimable) functional under appropriate conditions.
In contrast, for the general value functional, the linear functional corresponding to its directional derivative $D_{h}V\left(h_{0}\right)\left[\nu \right]$ given in (ref) does not have a well-defined Riesz representer in the Hilbert space $L^2(f)$. It is well-known that, according to Lemma 3.3 of chenpouzo2015sieve, Consequently, the functional $V$ becomes an irregular functional that cannot be estimated at $\sqrt{n}$ rate.
The above provides an intuitive explanation of why the welfare functional $W\left(h_{0}\right)$ is very special relative to general types of value functionals $V\left(h_{0}\right)$ or $\Phi\left(h_{0}\right)$, and clarifies why $W\left(h_{0}\right)$ could be $\sqrt{n}$-estimable while others generally cannot. In subsequent sections, we provide formal conditions and theorems that establish the $\sqrt{n}$-normality of plug-in estimators of the welfare functional $W\left(h_{0}\right)$, as well as the slower-than-$\sqrt{n}$ asymptotic normality for the value functional $V\left(h_{0}\right)$.
In this section, we focus on the welfare functional defined in (ref), which can also be equivalently written as a functional of $\mu_{0}$ as follows:
We write out the two equivalent definitions of the functionals based on $h_{0}$ and $\mu_{0}$, since each representation has its own merit. The representations $W\left(h_{0}\right)$ based on the CATE function $h_{0}$ is clearer in terms of its interpretation: the condition $h_{0}\left(x\right)\geq0$ is a direct optimal treatment assignment, and this representation has been adopted in previous work such as \citet*{kitagawa2018should}. On the other hand, the representations $\ol W\left(\mu_{0}\right)$ based on $\mu_{0}$ is clearer in terms of the underlying nonparametric regression function, which is notationally easier to work with in our subsequent semiparametric asymptotic analysis.
We consider two different setups for the estimation of $W\left(h_{0}\right)$, depending on how we treat the marginal density $f\left(x\right)$.
In the first setup, we treat $f$ as known and the functional $W\left(\cdot\right)$ as a known transformation of $h_{0}$. This is relevant in cases where the covariate density $f\left(x\right)$ of the target population is either configured or known/estimated from other sources than the training sample used to estimated CATE. In the formulation of $W\left(h_{0}\right)$ and $V\left(h_{0}\right)$ above, we take the density $f\left(x\right)$ to be known. For example, CATE may be estimated from a smaller pilot program (with sample size $n$), while the policymakers are contemplating to implement the policy on a statewide or nation-wide basis with a much larger target population, whose covariate density may either be known at the population level or estimated from an alternative data source with much larger sample size $N>>n$ (so that the sampling uncertainty in the estimation of $f$ becomes negligible relative to that in the estimation of $h_{0}$).
In the second setup, we take $f$ to be unknown and set it to $f_{0}$, the covariate density in the underlying population of the training sample. In this case, $f$ does not need to be explicitly estimated, but the integral in $W\left(\cdot\right)$ with respect to $f$ can be naturally approximated via sample average in the training sample. This corresponds more closely to the “empirical welfare” as considered in \citet*{kitagawa2018should}. We also provide results for this setting as well.
We start with the first setup, where the covariate density $f$ is taken to be known and $W\left(\cdot\right)$ is treated as a deterministic known integral of $h_{0}$ w.r.t. $f$. In this case, we can define a simple nonparametric plug-in estimators of $W\left(h_{0}\right)$ as \[ W\left(\hat{h}\right)\equiv\ol W\left(\hat{\mu}\right)=\int\left[\hat{h}\left(x\right)\right]_{+}f\left(x\right)dx. \]
We first state some key assumptions. Let $\hat{h}\left(x\right):=\hat{\mu}\left(x,1\right)-\hat{\mu}\left(x,0\right)$ be a nonparametric estimator of the CATE function, in which $\hat{\mu}\left(x,d\right)$ is a first-stage estimator of the nonparametric regression model
Theorem (ref) suggests the following natural estimator for the asymptotic variance $\sigma_{W}^{2}$:
where $\hat{u}_{i}:=Y_{i}-\hat{h}\left(X_{i}\right)$. This requires nonparametric estimation of propensity score function $p(x)$, as well as the knowledge (or nonparametric estimation) of the density $f_0$ or the Radon-Nikodym derivative $\lambda(x)$.
Alternatively and more preferably, we can use the sieve-based asymptotic variance estimator, which does not require estimation or knowledge of $p(x)$ and $\lambda(x)$. To do so, we first clarify some subtlety in the definition of the sieve in the first-stage nonparametric estimation of CATE as $\hat{\mu}(x,1)-\hat{\mu}(x,0)$, where $\hat{\mu}$ is a linear sieve estimator of $\mu_{0}\left(x,d\right)$ under random design on $\left(X_{i},D_{i}\right)$ with $D_{i}$ being binary.
In practice, $\hat{\mu}\left(x,1\right)$ is estimated with $K_1$ linear series in the treated subsample, while $\hat{\mu}\left(x,0\right)$ is estimated separately with potentially different $K_0$ linear series in the untreated subsample. However, since we treat $D_{i}$ as a random variable, we cannot directly treat $\hat{\mu}\left(x,1\right)$ and $\hat{\mu}\left(x,0\right)$ as two completely separate nonparametric estimators with exogenously given sample sizes. Instead, we treat $\hat{\mu}$ as the least square estimator of
where $\psi^{(K_1)}(x) = (\psi_1(x), \ldots, \psi_{K_1}(x))'$ and $\psi^{(K_0)}(x) = (\psi_{K_1+1}(x), \ldots, \psi_{K_1+K_0}(x))'$ denote the B-spline basis functions used to estimate $\mu_0(x, 1)$ and $\mu_0(x, 0)$, respectively, with sieve dimensions $K_1$ and $K_0$.
Define the vector-valued function $\overline{\psi}(x)$ such that its $k$th component equals $d\psi_k(x)$ for $k \in \{1, \ldots, K_1\}$ and $(1-d)\psi_k(x)$ for $k \in \{K_1+1, \ldots, K_1+K_0\}$. Following \citet*{chen2018optimal}, specifically equations (6) and (7), the variance (or standard error) estimator in our setting takes the form
where $\hat{\Omega}$ denotes the estimated asymptotic covariance matrix of the OLS estimators in (ref) given by \[ \hat{\Omega}:=\left(\Psi^{\left(2K\right)}\Psi^{\left(2K\right)'}\right)^{-1}\left(\frac{1}{n}\sum_{i=1}^{n}u_{i}^{2}\Psi^{\left(2K\right)}\Psi^{\left(2K\right)'}\right)\left(\Psi^{\left(2K\right)}\Psi^{\left(2K\right)'}\right)^{-1} \] and the estimated directional derivative vector $D_\mu \ol{W}\left(\hat{\mu}\right)\left[\ol{\psi}\right]$ is given by \[ D_\mu \ol{W}\left(\hat{\mu}\right)\left[\ol{\psi}\right] =
, \] where the minus sign in front of the sieve terms in $\psi^{K_0}(x)$ is due to the presence of the minus sign in front of $\mu(x,0)$ in $\text{CATE}(x)=\mu(x,1)-\mu(x,0)$.
We use Sobol points to numerically compute the integral above. See the simulation section for more details.
In certain cases, such as in \citet*{kitagawa2018should}, one might be interested in the welfare functional under the original distribution of covariates $F_{0}$, which may not be known or controlled. Given the random sample of $\left(Y_{i},D_{i},X_{i}\right)_{i=1}^{n}$ used to estimate the CATE $h_0$, it is natural to use the sample mean $\frac{1}{n}\sum_{i=1}^{n}[\cdot]$ as an estimator of the population expectation $\int [\cdot]dF_0$, in which case a natural plug-in estimator of $W\left(h_{0}\right)\equiv\ol W\left(\mu_{0}\right)$ is given by \[ \hat{W}\left(\hat{h}\right)\equiv\hat{\ol W}\left(\hat{\mu}\right):=\frac{1}{n}\sum_{i=1}^{n}\left[\hat{\mu}\left(X_{i},1\right)-\hat{\mu}\left(X_{i},0\right)\right]_{+}\equiv\frac{1}{n}\sum_{i=1}^{n}\left[\hat{h}\left(X_{i}\right)\right]_{+}. \] Clearly, the approximation of the integral introduces an additional source of randomness, but the result established in the last subsection continues to be useful in this case.
The standard errors can be computed similarly, based on straightforward adaptations of the analytical formula (ref) or the sieve-based formula (ref) in Section 3.1.
We now turn to the general value functional given by
The simple plug-in estimators are defined as
Under Assumption (ref) and (ref)(b)(c) the functional derivative of $\ol V\left(\mu_{0}\right)\left[\nu\right]$ is given by
As explained in Section (ref), $V\left(h_{0}\right)\equiv \ol V\left(\mu_{0}\right)$ is generally not $\sqrt{n}$-estimable, in particular, the $D_\mu \ol V\left(\mu_{0}\right)\left[\nu\right]$ are not bounded (or continuous) linear functionals on $L^2(f)$, which means that they do not have a Riesz representer on the whole Hilbert space $L^2(f)$. Nevertheless, sieve Riesz representer is well-defined (see \citet*{chen2025semiparametric}), which will be the key ingredient for the sieve variance term for the properties of $\ol V\left(\hat{\mu}\right)-\ol V\left(\mu_{0}\right)$. In particular we need to study the asymptotic property of the plug-in estimator of the submanifold integral of form (ref), which has been studied in \citet*{chen2025semiparametric} using linear series (Bspline) first stage: in particular, Section 4.3 of \citet*{chen2025semiparametric} analyzes the integral on upper contour set of the form $V\left(h_{0}\right)$ specifically. We use the result in \citet*{chen2025semiparametric} without repeating it here, but focus on the adaptation required for the standard error computation.
The standard error estimates can be computed based on the linear sieve first stage in a way similar to that described in Section 3.1, with the following adaptions. Again, we use the formula
where the pathwise derivative $D_\mu \ol{V}\left(\hat{\mu}\right)\left[\ol{\psi}\right]$ given by \[ D_\mu \ol{V}\left(\hat{\mu}\right)\left[\ol{\psi}\right] =
, \] can be approximated via
based on the mathematical result\footnote{See Theorem 3.13.(iii) of \citet*{evans2015measure}.} that \[ \lim_{\epsilon\searrow0}\frac{1}{2\epsilon}\int_{\left\{ x\in{\cal X}:~-\epsilon<h\left(x\right)<\epsilon\right\} }\omega\left(x\right)dx=\int_{\left\{ x\in{\cal X}:~h\left(x\right)=0\right\} }\frac{\omega\left(x\right)}{\norm{\nabla_{x}h\left(x\right)}}d{\cal H}^{d-1}\left(x\right). \] Again, we use Sobol points for numerical integration. See the simulation section for details, as well as robustness checks with respect the choice of $\epsilon$ in the numerical differentiation step.
We now consider the case where $F=F_{0}$ and population expectation $\mathbb{E}[\cdot]$ is estimated by the sample average $\frac{1}{n}\sum_{i=1}^{n}[\cdot]$ , and seek to characterize the asymptotic behavior of the natural plug-in estimator of $V\left(h_{0}\right)\equiv\ol V\left(\mu_{0}\right)$ is given by
It turns out that the additional error in the approximation of $V\left(\hat{h}\right)$ by $\hat{V}\left(\hat{h}\right)$ is asymptotically negligible relative to $V\left(\hat{h}\right)-V(h_0)$, which converges at a slower-than-$\sqrt{n}$ rate.
The standard error can be computed using formula (ref), with $\hat{D}_\mu\overline{V}\left(\hat{\mu}\right)\left[\ol{\psi}\right]$ given by the following adapted estimator,
where $\hat{f}(x)$ is a nonparametric density estimator of $f_0(x)$. We then again use Sobol points to numerically evaluate the integral in (ref).
We first report the finite-sample performance of the semiparametric estimator, its associated standard error estimator, and the resulting confidence interval, based on the theoretical results in Theorem (ref), which concern the welfare functional under a known target distribution. The model specifications used in the simulations are summarized in Table (ref). For each specification, random samples of size $n$ are drawn with covariates $X_i \sim F_0$, and treatment status is assigned according to the propensity score function $p_0(X_i)$. Outcomes are then generated as $Y_i = \mu_0(X_i, D_i) + \epsilon_i$, with $\epsilon_i \sim N(0,1)$.
The semiparametric estimator of the welfare functional is constructed in two steps. In the first step, $\mu_0(x,1)$ and $\mu_0(x,0)$ are estimated separately for the treated group ($D_i = 1$) and the control group ($D_i = 0$) using B-spline regressions. In the second step, the welfare functional is approximated by numerically integrating over $M = 5{,}000$ Sobol points $\{X_j^{Sobol}\}_{j=1}^M$ drawn from the target distribution $F$: \[ \hat{\overline{W}}(\hat{\mu}) \;=\; \frac{1}{M} \sum_{j=1}^M \big[\hat{\mu}(X_j^{Sobol},1) - \hat{\mu}(X_j^{Sobol},0)\big]_+, \] where $\hat{\mu}(x,1)$ and $\hat{\mu}(x,0)$ denote the first-stage nonparametric estimators.
Let $\hat{h}(x) = \hat{\mu}(x,1) - \hat{\mu}(x,0)$ denote the estimated CATE function, and let $\hat{p}(x)$ be a B-spline sieve estimator of the propensity score. The $95\%$ confidence interval for the welfare functional takes the usual form, \[ \left[\hat{\overline{W}}(\hat{\mu}) - 1.96 \frac{\hat{\sigma}_W}{\sqrt{n}}, \;\; \hat{\overline{W}}(\hat{\mu}) + 1.96 \frac{\hat{\sigma}_W}{\sqrt{n}}\right], \] with \[ \hat{\sigma}_{W}^{2} =\frac{1}{N}\sum_{i}\frac{\mathbf{\mathbbm1}\left\{ \hat{h}\left(X_{i}\right)\geq0\right\} \lambda^{2}\left(X_{i}\right)(Y_i -\hat{h}(X_i))^2}{\hat{p}\left(X_{i}\right)\left(1-\hat{p}\left(X_{i}\right)\right)}. \]
To evaluate performance, we simulate each model 2,000 times at sample sizes $n = 1500, 3000,$ and $6000$. Table (ref) reports the true welfare functional ($W_{\text{true}}$), the average bias of the estimator $\hat{\overline{W}}(\hat{\mu})$ (Bias), its sampling standard deviation (SD), the average estimated standard error (SE), the standard deviation of SE across iterations (SD(SE)), and the empirical coverage probability of the associated 95% confidence interval (Coverage). The results\footnote{Additional simulations based on GAM are presented in the Appendix (ref).} indicate that the nominal coverage rate is attained in nearly all cases, even with relatively small samples, and improves further as sample size increases. The bias of the estimator also becomes negligible relative to its sampling variability, and the proposed SE estimator closely tracks the sampling standard deviation with high precision.
One potential drawback of the plug-in approach to estimating the asymptotic variance of the welfare functional is that it requires nonparametric estimation of the propensity score function $p_0(x)$, which can introduce additional sampling noise and numerical instability. As an alternative, one may employ the sieve variance estimator as in chen2014sieveIrregular or \citet*{chen2015sieve}, whose formula depends only on the pathwise derivative of the welfare functional and a linear regression of the outcome variable on the sieve basis. When the target distribution $F$ is known, the pathwise derivative of the welfare functional can be numerically approximated using Sobol points drawn from the target population combined with importance sampling. When the target distribution is unknown, the pathwise derivative can instead be evaluated directly at the observed data and approximated by a sample average. Simulation results and the empirical application based on the sieve variance estimator are shown in the Appendix (ref).
We now investigate the finite-sample performance of our estimation and inference procedure for the welfare functional in the case where the target distribution $F$ coincides with the population distribution $F_0$, though both remain unknown. The relevant model specifications are presented in Table (ref). In contrast to the designs considered in Section (ref), these specifications explicitly impose that $F$ and $F_0$ share identical supports.
The plug-in estimator $\hat{\overline{W}}(\hat{\mu})$ proposed here differs from that in Section (ref) only in the second step: instead of using Sobol points to approximate the integral, we take the sample average of $\big[\hat{\mu}(X_i,1) - \hat{\mu}(X_i,0)\big]_+$ over the observed data: \[ \hat{\overline{W}}(\hat{\mu}) \;=\; \frac{1}{n}\sum_{i=1}^n \big[\hat{\mu}(X_i,1) - \hat{\mu}(X_i,0)\big]_+. \]
Relative to the known-$F$ case, the asymptotic variance of $\hat{\overline{W}}(\hat{h})$ includes an additional component, $\mathrm{Var}([h_0(X_i)]_+)$. We estimate the asymptotic standard deviation $\hat{\overline{\sigma}}_W$ by substituting the B-spline sieve estimates for the nuisance functions and replacing the population mean with its sample analog. As before, we restrict attention to observations with estimated propensity scores in $[0,1]$. A 95% confidence interval is then given by \[ \left[\hat{\overline{W}}(\hat{\mu}) - 1.96 \,\frac{\hat{\overline{\sigma}}_W}{\sqrt{n}}, \;\; \hat{\overline{W}}(\hat{\mu}) + 1.96 \,\frac{\hat{\overline{\sigma}}_W}{\sqrt{n}}\right]. \]
The simulation results\footnote{Additional simulations based on GAM are presented in the Appendix (ref).} based on 2,000 iterations with sample sizes $n = 1500, 3000,$ and $6000$ are reported in Table (ref). Overall, the coverage of the proposed confidence intervals converges to the nominal 95% level as sample size increases, while the decreasing bias and standard deviation indicate a reduction in mean squared error.
To assess the finite-sample properties of our estimation inference procedure for the value functional $\overline{V}(\hat{\mu})$ under a known target distribution $F$, we analyze the model described in Table (ref):
Our parameter of interest is a scaled value functional under known $F$: \[ V(h_0) \;=\; 3^2 \int \mathbbm{1}\!\left\{ \big(1 - x_{1}^2 - x_{2}^2\big)\,\big(4 + \sin(x_{1})x_{2} + \cos(x_{2})\big) \;\geq\; 0 \right\} \, dF(x_{1},x_{2}), \] where the scaling factor $3^2$ is chosen so that the integral evaluates to $\pi$. Equivalently, this setup can be interpreted as a Monte Carlo experiment for estimating the area of the unit circle by uniformly throwing darts over a $3 \times 3$ square.
The plug-in estimator $\hat{\overline{V}}(\hat{\mu})$ for $\overline{V}(\mu_0)$ is constructed in two steps. In the first step, we estimate $\mu_0(x,1) = \mu_0(x_1,x_2,1)$ and $\mu_0(x,0) = \mu_0(x_1,x_2,0)$ separately using B-spline sieve estimators, fitted on the treated and control groups, respectively, with sieve dimensions predetermined as described in subsection (ref). In the second step, $\mathbbm{1}\{(\hat{\mu}_0(x, 1) - \hat{\mu}_0(x, 0)) \geq 0\}$ is numerically integrated using $5000$ Sobol points drawn from $F$. The resulting semiparametric two-step estimator is \[ \hat{\overline{V}}(\hat{\mu}) \;=\; \frac{1}{M} \sum_{j=1}^M \mathbbm{1}\{\hat{\mu}(X_j^{Sobol},1) - \hat{\mu}(X_j^{Sobol},0)\big\}. \]
The 95% confidence interval is then constructed as \[ \left[\hat{\overline{V}}(\hat{\mu}) - 1.96 \,\frac{\hat{\sigma}_V}{\sqrt{n}}, \;\; \hat{\overline{V}}(\hat{\mu}) + 1.96 \,\frac{\hat{\sigma}_V}{\sqrt{n}}\right], \] where the asymptotic variance estimate $\hat{\sigma}_{V}^{2}$ is given by \[ \hat{\sigma}_{V}^{2}= \hat{D}_\mu \overline{V}\left(\hat{\mu}\right)\left[\ol{\psi}\right]'\hat{\Omega} \hat{D}_\mu \overline{V} \left(\hat{\mu}\right)\left[\ol{\psi}\right], \] with $\hat{\Omega}$ being the estimated asymptotic covariance matrix for the OLS estimators in the linear regression model (ref), and \[ \hat{D}_\mu\overline{V}\left(\hat{\mu}\right)\left[\ol{\psi}\right] = 3^2
. \]
We set the tuning parameter $\epsilon = 0.005$ to mitigate bias in $\hat{\overline{V}}(\hat{\mu})$ and approximate $\hat{D}_\mu \overline{V}(\hat{\mu})[\nu]$ using the sample average over $M$ Sobol draws from $F$. Because draws from $F$ are unlikely to fall within the set $\{x \in [-1.5, 1.5]^2 : -\epsilon < \hat{h}(x) < \epsilon\}$ when $\epsilon$ is small, we use $M = 1{,}000{,}000$ Sobol points to ensure accuracy. Simulation results\footnote{Additional simulations under alternative model specifications, as well as robustness checks for different values of $\epsilon$, are presented in Appendix (ref).} based on 2,000 iterations with sample sizes $n = 1500, 3000,$ and $6000$ are reported in Table (ref). The results show that the coverage rate reaches the nominal level with relatively modest sample sizes, even though the value functional is not $\sqrt{n}$-estimable. Also, as sample size increases, both bias and standard error decrease, and the plug-in estimator for the standard error provides a close estimate of the sampling standard deviation.
We revisit the empirical application analyzed in kitagawa2018should using Job Training Parternship Act (JTPA) dataset. The JTPA study randomized whether applicants were eligible to receive a mix of training, job-search assistance, and other services provided under the program for a period of 18 months. A detailed description of the study and an assessment of average program effects for five major subgroups of the target population are provided in \citet*{BloomEtAl1997}.
We evaluate welfare using two outcome measures, following the approach in KT18. The first is total earnings over the 30 months following treatment assignment. The second adjusts this measure by subtracting \$774 for individuals assigned to treatment, thereby incorporating program costs. These outcomes are considered from an intention-to-treat perspective, meaning we focus on eligibility assignment rather than treatment effects among compliers. The available covariates include applicants' pre-program earnings, years of education, and treatment status. Our objective is to estimate and conduct inference on the first-best welfare functional and the optimal fraction of the population that should receive treatment.
As the first step in the estimation and inference procedure, we trim observations outside the common support of the treated and control groups to enforce the overlap assumption and to avoid extrapolation when applying sieve estimators. For either outcome measure, we then estimate $\mu_0(x,1)$ nonparametrically using B-spline sieves fitted on the treated sample, and $\mu_0(x,0)$ analogously on the control sample. The sieve dimensions are selected in a data-driven manner following \citet*{chen2025adaptive}. Combining these estimates on the common support yields the CATE estimate $\hat{h}(x) = \hat{\mu}(x,1) - \hat{\mu}(x,0)$. Taking sample averages of $[\hat{h}(x)]_+$ and $\mathbbm{1}\{\hat{h}(x) > 0\}$ over the trimmed dataset produces the estimators $\hat{\overline{W}}(\hat{\mu})$ and $\hat{\overline{V}}(\hat{\mu})$, respectively. Confidence intervals for $\overline{W}(\mu_0)$ and $\overline{V}(\mu_0)$ are then constructed according to their asymptotic theories: $\left(\hat{\overline{W}}(\hat{\mu}) \pm 1.96 \,\hat{\ol{\sigma}}_W/\sqrt{N}\right),$ and $\left(\hat{\overline{V}}(\hat{\mu}) \pm 1.96 \,\hat{\sigma}_V/\sqrt{N}\right)$, where $N$ denotes the sample size of the JTPA dataset.
While computing sieve estimate of $\hat{\overline{\sigma}}_W$ is straightforward, computing sieve estimate of $\hat{\sigma}_V$ involves several additional steps. Specifically, it requires a density estimate $\hat{f}(x)$ for the covariates as discussed in Remark (ref). To this end, we use a Gaussian kernel density estimator, selecting the bandwidth matrix via the smoothed cross-validation method (Hscv() in the ks package) and applying a scaling factor of $s=3$ to ensure adequate smoothness in the presence of discrete values for years of education. Furthermore, we need to specify a small hyperparameter $\epsilon$ to provide a close approximate for the pathwise derivative of the value functional in $\hat{\sigma}_V$. We set $\epsilon$ equal to a fraction $\iota = 0.01$ of the standard deviation of $\hat{h}(x)$ over the trimmed dataset. Robustness checks on the tuning parameters $s$ and $\iota$ are reported in Appendix (ref), and alternative density estimation approaches are examined in Appendix (ref).
Table (ref) presents our estimation and inference results alongside the corresponding findings from KT18. Specifically, we report our estimated welfare gain and the share of individuals to be treated, along with their confidence intervals, based on the trimmed dataset. For comparability, we also present results obtained using the untrimmed dataset-on which KT18 conducted their analysis-with the same choice of tuning parameters. The nonparametric plug-in estimates from KT18, which target the same welfare and share parameters as in our analysis, serve as an empirical benchmark. In addition, the linear rule estimates with their associated confidence intervals from KT18 are reported to assess how conservative our nonparametric inference procedure is relative to their parametric counterparts.
Across both the trimmed and untrimmed datasets, our estimates of the welfare gain and the optimal treatment share are broadly consistent with the nonparametric plug-in rule estimates reported in KT18, despite methodological difference in the first stage: we employ sieve estimators for the nuisance functions, whereas they use a Nadaraya-Watson estimator with an Epanechnikov kernel. A key contribution of our analysis is the provision of confidence intervals for both parameters of interest. By contrast, KT18 report confidence intervals only for the welfare gain under parametric rules, but not for the optimal treatment share. Although our confidence interval for the welfare gain appears wide, its length is comparable to that under the linear rule in KT18, underscoring that our inference procedure remains sharp even while relying on fully nonparametric methods.