The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
116,237 characters
Inference on Individual Treatment Effects in Nonseparable Triangular Models.
\title{\textbf{Inference on Individual Treatment Effects in Nonseparable
Triangular Models}\let\thefootnote\relax\footnotetext{\today.}\thanks{We thank our coeditor, Xiaohong Chen, the associate editor and two
anonymous referees, whose comments have greatly improved the paper.
This research is supported by the National Natural Science Foundation
of China under grant 71903190, fund for building world-class universities
(disciplines) of Renmin University of China (Ma), and Japan Society
for the Promotion of Science KAKENHI Grant Number 21K01419 (Yu).}}
\author{Jun Ma\thanks{School of Economics, Renmin University of China, P.R. China. Email:
[email removed]}\and Vadim Marmer\thanks{Vancouver School of Economics, University of British Columbia, Canada.
Email: [email removed]}\and Zhengfei Yu\thanks{Faculty of Humanities and Social Sciences, University of Tsukuba,
Japan. Email: [email removed]}}
\maketitle
\begin{abstract}
In nonseparable triangular models with a binary endogenous treatment
and a binary instrumental variable, \citet{vuong2017counterfactual}
and \citet*{feng2019estimation} respectively provide the identification
results for the individual treatment effects (ITEs) under the rank
invariance assumption and propose a uniformly consistent kernel density
estimator using estimated ITEs for the density of the ITE. This paper
establishes the asymptotic normality of the density estimator of \citet*{feng2019estimation}
and shows that the estimation error of the ITEs that vanishes at the
root $n$ rate has a non-negligible effect on the asymptotic distribution.
We propose asymptotically valid standard errors that account for estimated
ITEs, as well as a bias correction. Furthermore, we develop uniform
confidence bands for the density of the ITE using the jackknife multiplier
or nonparametric bootstrap critical values.\\
\textbf{Keywords: }Individual treatment effects, nonparametric triangular
models, two-step nonparametric estimation, bootstrap, uniform confidence
bands, labor supply and family size\\
\textbf{JEL classification: }C12, C14, C31, C36
\end{abstract}
\section{Introduction}
Heterogeneous treatment effects have received increasing attention
in the causal inference and policy evaluation literature (\citealp{angrist2004treatment,Heckman1997,Heckman2006}).
There is a vast literature studying the causal effect of ceteris paribus
change of a treatment variable using triangular models (see, e.g.,
\citealp{Chesher2003,Chesher2005,d2015identification,Imbens2009,Jun2011,Newey1999,Torgovitsky2015,Vytlacil2007}
among others). In a triangular model, the outcome variable is generated
by an outcome equation, and a selection equation determines the endogenous
treatment variable. Recently, \citet[VX, hereafter]{vuong2017counterfactual}
and \citet[FVX, hereafter]{feng2019estimation} developed nonparametric
identification and estimation methods for individual treatment effects
(ITEs) in a triangular model with a nonseparable outcome equation,
a selection equation which is a latent index model (\citealp{Vytlacil2002}),
a binary endogenous treatment variable and a binary instrument under
the rank invariance assumption. VX assumes that disturbances in both
equations are scalar-valued and that the outcome is a strictly monotone
function of the disturbance. The outcome equation in VX satisfies
the rank invariance assumption (i.e., for given covariates, the ranks
of the two potential outcomes are the same). See, e.g., \citet{chernozhukov2020instrumental}
for discussion of this assumption. The triangular model considered
in VX and FVX is also closely related to the classical local average
treatment effect (LATE) model (see, e.g., \citealp{Abadie2002,Abadie2003,Froelich2013,imbens1994identification}
among others) and the instrumental variable quantile regression (IVQR)
model (see, e.g., \citealp{Chernozhukov2005} among others).\footnote{A more detailed literature review about the triangular model with
an endogenous treatment can be found in our online supplement available
at: \href{http://ruc-econ.github.io/ITE_Supp_Rev_V13.pdf}{ruc-econ.github.io/ITE\_Supp\_Rev\_V13.pdf}.}
Since the ITE is defined and estimated for each individual, it is
natural to focus on its probability density function (PDF) when assessing
the heterogeneity of a treatment. For that purpose, FVX uses the conventional
kernel density estimator applied to estimated ITEs. They show its
uniform consistency as well as derive its rate of convergence. However,
when it comes to inference for the density of the ITE, two theoretical
problems still need to be solved. First, asymptotically valid standard
errors for the density estimator should incorporate the uncertainty
stemming from the estimation of ITEs. Second, as researchers are often
interested in the shape of the distribution of the ITE, it is important
to have asymptotically valid uniform confidence bands (UCBs) for the
density of the ITE. Our paper contributes to the literature by providing
easy-to-implement solutions to both problems.
We first provide a sharper bound for the uniform rate of convergence
of the FVX density estimator and show that it attains the optimal
rate under seemingly minimal conditions. We then show the asymptotic
normality of the FVX estimator and derive an analytical formula for
the standard error of the density estimator that incorporates the
influence of ITEs' estimation in the first step of the FVX procedure.
The asymptotic normality result in this paper is non-trivial as the
asymptotic distribution of the FVX estimator is different from that
of the infeasible estimator based on true unobserved ITEs. FVX uses
nonparametrically estimated counterfactual mappings to generate pseudo
(i.e., estimated) ITEs in the first step of their procedure. In the
second step, they apply kernel density estimation to the pseudo ITEs
to construct an estimator for the PDF of the ITE. While the estimated
counterfactual mappings converge at the root $n$ rate, we show that
the first-step estimation errors' contribution to the density estimator's
asymptotic variance is non-negligible and can substantially dominate
that of the second stage. This phenomenon is due to discontinuities
in the linearization of the first-step estimator. At the same time,
the asymptotic bias is unaffected by the first-stage estimation errors
and equal to that of the infeasible estimator.
The paper's second contribution is to propose asymptotically valid
UCBs for the density of the ITE. The proposed UCB captures the uncertainty
about the entire estimated density function and, therefore, can be
used for inference about the shape of the ITE's distribution (e.g.,
the number and locations of the modes) and for comparisons between
the distributions of the ITE in different sub-populations. UCBs can
also be used for the specification of parametric models. Our bootstrap
UCBs have the desirable property of polynomial coverage error decay
rates. Following \citet{Calonico2014}, we also propose bias-corrected
UCBs using standard errors that incorporate additional variability
from the estimated bias (Section \ref{subsec:Bias-corrected-JMB}).
This approach is common in the recent literature, as it validates
the use of conventional data-driven bandwidth selectors for inference.
In addition, we consider an extension to inference on PDF conditional
on a sub-vector of covariates (Section \ref{subsec:Conditioning-on-sub-vectors}).
The FVX estimator and our UCBs require bandwidth selection. Our results
explicitly allow for random data-dependent bandwidths. Following the
literature (e.g., \citealp{hsiao2007consistent,Li:2010dy}) and the
standard practice in applied work, we assume that the data-dependent
bandwidth consistently estimates some deterministic bandwidth. Similarly
to the existing literature, we verify that the uniform rate of convergence
and asymptotic normality results are unaffected by the bandwidth estimation.
However, we go a step further and provide an explicit estimate of
the effect of bandwidth estimation on the coverage error of the UCBs.
The result is new for this literature and adds to our understanding
of the impact of bandwidth selection on inferential procedures.
In the empirical section of the paper, we use the FVX estimator with
our UCBs to study the effect of having more than two children on their
parents' labor income using the instrument proposed in \citet{angrist1998children}.
We show that the conditional distribution of the ITE is significantly
different between households with high-school-only and college-educated
mothers. In the latter case, the ITE's distribution is more dispersed.
However, it also places significantly more weight on positive effects.
From the perspective of nonparametric inference, our paper contributes
to the literature on inference in the presence of nonparametrically
generated variables. See, e.g., \citet{Mammen2012} and \citet{Ma2019}.
The asymptotic theory of the FVX estimator is different from the results
obtained in \citet{Mammen2012} and \citet{Ma2019}. E.g., among their
other results, \citet{Mammen2012} show that the influence of variables'
estimation can be made asymptotically negligible using a proper choice
of the bandwidth. However, in the case of the FVX estimator, the influence
of ITEs' estimation has a non-negligible effect regardless of the
bandwidth choice. \citet{Ma2019} show that in the context of first-price
auctions, the estimation of latent bidders' valuations contributes
to the asymptotic variance of the \citet{Guerre2000} estimator for
the density of valuations. Moreover, the variance of the \citet{Guerre2000}
estimator has a slower decay rate than that of the infeasible estimator
constructed using the true latent valuations. Here, we show that while
the variance of the FVX estimator has the same decay rate as that
of the infeasible estimator, the estimation of ITEs cannot be ignored.
Our paper also contributes to the literature on inference for nonparametrically
estimated functions. In this paper, we take the intermediate Gaussian
approximation approach to show the asymptotic validity of our UCBs
by using tools developed by \citet{Chernozhukov2014gaussian,Chernozhukov2014anti,Chernozhukov2016,Chen2020jackknife}.\footnote{Under certain conditions, it can often be shown that the suprema of
estimation errors can be approximated by the suprema of tight Gaussian
random elements using the theorems of \citet{Chen2020jackknife,Chernozhukov2014gaussian}.
Then theorems in \citet{Chen2020jackknife,Chernozhukov2014anti,Chernozhukov2016}
show that the distributions of these Gaussian suprema can be approximated
by bootstrapping.} This approach was recently applied, e.g., in \citet{chen2018optimal,Cheng2019,Kato2019,Ma2019}
among others, to show the asymptotic validity of bootstrap UCBs for
various nonparametric curves in different contexts.
The rest of the paper is organized as follows. Section \ref{sec:Model-and-identification}
reviews the model setup and identification of the distribution of
the ITE and discusses the nonparametric estimation of the density
of the ITE. Section \ref{sec:Asymptotic-properties} shows the density
estimator's uniform convergence rate and asymptotic normality. Section
\ref{sec:Robust-inference} provides standard errors for the density
estimator that can be used for asymptotically valid inference. It
also establishes the validity of the bootstrap UCBs. Section \ref{subsec:Jackknife-multiplier-bootstrap}
describes the algorithm for our proposed confidence band. Section
\ref{sec:Monte-Carlo-experiments} presents the results from Monte
Carlo experiments. Section \ref{sec:Empirical-illustrations} applies
our inference method to study the effect of family size on labor income.
The proofs of the theorems and statements of the technical lemmas
are presented in the appendices. The proofs of the lemmas, auxiliary
results, and additional simulation evidence are collected in the online
supplement (\href{http://ruc-econ.github.io/ITE_Supp_Rev_V13.pdf}{ruc-econ.github.io/ITE\_Supp\_Rev\_V13.pdf}).
\section{Model and the FVX estimator \label{sec:Model-and-identification}}
For completeness, first, we describe the model setup of VX and FVX,
and their estimator. Let $\mathbbm{1}\left(\cdot\right)$ denote the
indicator function. The outcome and selection equations are given
respectively by
\begin{eqnarray}
Y & = & g\left(D,X,\epsilon\right)\nonumber \\
D & = & \mathbbm{1}\left(\eta\leq s\left(Z,X\right)\right),\label{eq: triangular model}
\end{eqnarray}
where $Y\in\mathbb{R}$ is a continuously distributed outcome variable,
$D\in\left\{ 0,1\right\} $ is an endogenous treatment variable, and
$X\in\mathscr{S}_{X}$ is a vector of observed explanatory variables
(covariates) with $\mathscr{S}_{V}$ denoting the support of the distribution
of a random vector $V$ (i.e., the smallest closed set $C$ satisfying
$\mathrm{Pr}\left[V\in C\right]=1$). $Z\in\left\{ 0,1\right\} $
is a binary instrumental variable that has no direct effect on $Y$
and, therefore, are excluded from the outcome equation. $\left(\epsilon,\eta\right)$
are unobserved scalar-valued disturbances conditionally independent
of $Z$ given $X$. $g$ and $s$ are unknown functions.
The functions $g\left(d,x,\cdot\right)$ and $s\left(\cdot,x\right)$
are assumed to be strictly increasing. The selection equation in (\ref{eq: triangular model})
has the form of a latent index selection model: treatment is assigned
if some latent index or utility $s\left(Z,X\right)$ crosses the threshold
$\eta$. The ITE is defined as
\begin{equation}
\varDelta\coloneqq g\left(1,X,\epsilon\right)-g\left(0,X,\epsilon\right),\label{eq:Delta through Y}
\end{equation}
where ``$a\coloneqq b$'' is understood as ``$a$ is defined by
$b$''. Note that $\varDelta$ is random conditionally on $X$ due
to the unobserved $\epsilon$, i.e., the treatment effect varies among
individuals with the same observed characteristics. The unobserved
disturbances $\epsilon$ and $\eta$ are allowed to be correlated
conditionally on $X$.\footnote{The model allows the ITEs to be ``essentially heterogeneous'' (\citealp{Heckman2006})
since whether or not individuals who have the same observed characteristics
select into treatment can be correlated with the gain from treatment.} Denote $d'\coloneqq1-d$, and let $g^{-1}\left(d',x,\cdot\right)$
be the inverse function of $g\left(d',x,\cdot\right)$. Let $\mathscr{S}_{V\mid W=w}$
denote the support of the conditional distribution of $V$ given $W=w$.
For $y\in\mathscr{S}_{g\left(d',x,\epsilon\right)\mid X=x}$, define
the corresponding counterfactual mapping $\phi_{dx}\left(y\right)\coloneqq g\left(d,x,g^{-1}\left(d',x,y\right)\right)$,
i.e., $\phi_{dx}\left(y\right)$ is the outcome one would observe
instead of $y$ if the treatment status $d'$ were switched to $d$.
Using the counterfactual mappings $\left(\phi_{0x},\phi_{1x}\right)$,
we can write the ITE as
\begin{equation}
\varDelta=D\left(Y-\phi_{0X}\left(Y\right)\right)+\left(1-D\right)\left(\phi_{1X}\left(Y\right)-Y\right).\label{eq:Delta}
\end{equation}
VX shows constructive nonparametric identification of the counterfactual
mappings. This result establishes the identification of the entire
distribution of $\varDelta$. FVX shows that $\phi_{dx}\left(y\right)$
is the unique minimizer of the strictly convex function $Q_{dx}\left(\cdot;y\right)$
defined below:
\begin{multline}
Q_{dx}\left(t;y\right)\coloneqq\\
\left(\mathrm{E}\left[\mathbbm{1}\left(D=d\right)\left|Y-t\right|\mid Z=d,X=x\right]-\mathrm{E}\left[\mathbbm{1}\left(D=d'\right)\mathrm{sgn}\left(Y-y\right)\mid Z=d,X=x\right]\times t\right)\\
-\left(\mathrm{E}\left[\mathbbm{1}\left(D=d\right)\left|Y-t\right|\mid Z=d',X=x\right]-\mathrm{E}\left[\mathbbm{1}\left(D=d'\right)\mathrm{sgn}\left(Y-y\right)\mid Z=d',X=x\right]\times t\right),\label{eq:Q_dx definition}
\end{multline}
where $\mathrm{sgn}\left(u\right)\coloneqq2\times\mathbbm{1}\left(u>0\right)-1$
denotes the left continuous sign function.
The econometrician observes $\{\left(Y_{i},D_{i},X_{i},Z_{i}\right):i=1,\ldots,n\}$,
a sample of observations on $\left(Y,D,X^{\top},Z\right)^{\top}$
generated by the model. Let $\widehat{Q}_{dx}^{\left(-i\right)}\left(t;y\right)$
denote the leave-$i$-out sample analogue of $Q_{dx}\left(t;y\right)$
constructed under the FVX assumption that $X$ is discretely distributed:
\begin{multline}
\widehat{Q}_{dx}^{\left(-i\right)}\left(t;y\right)\coloneqq\\
\frac{\sum_{j\neq i}\left\{ \mathbbm{1}\left(D_{j}=d,Z_{j}=d,X_{j}=x\right)\left|Y_{j}-t\right|-\mathbbm{1}\left(D_{j}=d',Z_{j}=d,X_{j}=x\right)\mathrm{sgn}\left(Y_{j}-y\right)t\right\} }{\sum_{j\neq i}\mathbbm{1}\left(Z_{j}=d,X_{j}=x\right)}\\
-\frac{\sum_{j\neq i}\left\{ \mathbbm{1}\left(D_{j}=d,Z_{j}=d',X_{j}=x\right)\left|Y_{j}-t\right|-\mathbbm{1}\left(D_{j}=d',Z_{j}=d',X_{j}=x\right)\mathrm{sgn}\left(Y_{j}-y\right)t\right\} }{\sum_{j\neq i}\mathbbm{1}\left(Z_{j}=d',X_{j}=x\right)}.\label{eq:Q_hat definition}
\end{multline}
The leave-$i$-out nonparametric estimator of $\phi_{dx}\left(y\right),d\in\left\{ 0,1\right\} $,
can be constructed as
\begin{equation}
\widehat{\phi}_{dx}^{\left(-i\right)}\left(y\right)\coloneqq\underset{t\in\left[\underline{y}_{dx},\overline{y}_{dx}\right]}{\mathrm{argmin}}\widehat{Q}_{dx}^{\left(-i\right)}\left(t;y\right),\label{eq:phi_hat definition}
\end{equation}
where we write $\mathscr{S}_{g\left(d,x,\epsilon\right)\mid X=x}=\left[\underline{y}_{dx},\overline{y}_{dx}\right]$.\footnote{As FVX, we assume that $\underline{y}_{dx}$ and $\overline{y}_{dx}$
are known. Lemma 1 of VX shows that the supports of the potential
outcomes, $\mathscr{S}_{g\left(d,x,\epsilon\right)\mid X=x}=\left[\underline{y}_{dx},\overline{y}_{dx}\right]$,
are identified by $\mathscr{S}_{g\left(d,x,\epsilon\right)\mid X=x}=\mathscr{S}_{Y\mid D=d,X=x}$.
In practical implementation, $\underline{y}_{dx}$ and $\overline{y}_{dx}$
can be estimated. See Section 3 of FVX for discussion.} One can now estimate the ITEs by replacing $\phi_{dx}(y)$ in (\ref{eq:Delta})
with its leave-$i$-out nonparametric estimator $\widehat{\phi}_{dx}^{\left(-i\right)}\left(y\right)$:
\begin{equation}
\widehat{\varDelta}_{i}=D_{i}\left(Y_{i}-\widehat{\phi}_{0X_{i}}^{\left(-i\right)}\left(Y_{i}\right)\right)+\left(1-D_{i}\right)\left(\widehat{\phi}_{1X_{i}}^{\left(-i\right)}\left(Y_{i}\right)-Y_{i}\right).\label{eq:pseudo ITE definition}
\end{equation}
The FVX estimator of $f_{\varDelta\mid X}\left(v\mid x\right)$, the
conditional density of $\varDelta$ given $X=x$, is the kernel density
estimator that uses $\widehat{\varDelta}_{i}$ in place of the true
unobserved ITEs:
\begin{equation}
\widehat{f}_{\varDelta\mid X}\left(v\mid x;b\right)\coloneqq\frac{\sum_{i=1}^{n}\frac{1}{b}K\left(\frac{\widehat{\varDelta}_{i}-v}{b}\right)\mathbbm{1}\left(X_{i}=x\right)}{\sum_{i=1}^{n}\mathbbm{1}\left(X_{i}=x\right)},\label{eq:f_hat definition}
\end{equation}
where $K\left(\cdot\right)$ and $b>0$ denote the kernel function
and bandwidth, respectively.
Let $\left[a\pm b\right]$ denote the interval $\left[a-b,a+b\right]$.
In this paper, we propose a bootstrap UCB defined by the following
collection of random intervals:
\begin{equation}
\mathit{CB}_{\mathsf{jmb}}\left(v\mid x;b,b_{\zeta}\right)\coloneqq\left[\widehat{f}_{\varDelta\mid X}\left(v\mid x;b\right)\pm z_{1-\alpha}^{\mathsf{jmb}}\sqrt{\frac{\widehat{V}\left(v\mid x;b,b_{\zeta}\right)}{nb}}\right],\label{eq:confidence band definition}
\end{equation}
where $b_{\zeta}$ is another bandwidth. Algorithm \ref{alg:multiplier}
in Section \ref{subsec:Jackknife-multiplier-bootstrap} provides detailed
step-by-step instructions for constructing the UCB. In (\ref{eq:confidence band definition})
above, the standard error $\sqrt{\widehat{V}\left(v\mid x;b,b_{\zeta}\right)/\left(nb\right)}$
uses the variance estimator $\widehat{V}\left(v\mid x;b,b_{\zeta}\right)$
defined in equation (\ref{eq:asy var}) in Section \ref{subsec:Standard-errors}.
The bootstrap critical value $z_{1-\alpha}^{\mathsf{jmb}}$ is defined
in equation (\ref{eq:z_pound definition}) in Section \ref{subsec:Jackknife-multiplier-bootstrap}.
The need for the second bandwidth $b_{\zeta}$ is discussed in Section
\ref{subsec:Standard-errors}. We provide data-dependent procedures
for selecting the two bandwidths $b$ and $b_{\zeta}$. See the discussions
following Assumption \ref{assu:h_hat} in Section \ref{subsec:Assumptions}
and Assumption \ref{assu:btw zeta} in Section \ref{subsec:Standard-errors}
respectively. Theorem \ref{thm:confidence band} in Section \ref{subsec:Jackknife-multiplier-bootstrap}
is our main result. It establishes that the proposed UCB covers $f_{\varDelta\mid X}\left(v\mid x\right)$
simultaneously over a range of $v$ values with a pre-specified confidence
level in large samples.
\section{Asymptotic properties of the FVX estimator\label{sec:Asymptotic-properties}}
In this section, we establish two new asymptotic results for the FVX
density estimator. Theorem \ref{thm:uniform rate} below shows that
under seemingly minimal conditions (see Assumptions \ref{assu: DGP1}
and \ref{assu: DGP4} ahead), the FVX estimator has the same uniform
rate of convergence as that of the infeasible kernel density estimator
that uses true ITEs, and attains the optimal uniform rate of convergence
(see \citealp{Stone1982}). Theorem \ref{thm:asymptotic normality}
shows that the FVX estimator is asymptotically normal. However, its
asymptotic variance is larger than that of the infeasible estimator.
We show that these results hold under either a deterministic bandwidth
or a data-dependent bandwidth that satisfies Assumption \ref{assu:h_hat}
below.
\subsection{Assumptions\label{subsec:Assumptions}}
The following assumption on the data generating process (DGP) is similar
to those in VX and FVX.
\begin{assumption}[DGP]
\label{assu: DGP1} (a) For all $\left(d,x\right)\in\mathscr{S}_{\left(D,X\right)}$,
$g\left(d,x,\cdot\right)$ is continuously differentiable and strictly
increasing. (b) $Z\in\left\{ 0,1\right\} $ is independent of $\left(\epsilon,\eta\right)$
conditionally on $X$. (c) For all $x\in\mathscr{S}_{X}$, $\mathrm{Pr}\left[D=1\mid Z=1,X=x\right]\ne\mathrm{Pr}\left[D=1\mid Z=0,X=x\right]$.
(d) The conditional distribution of $(\epsilon,\eta)$ given $X$
is absolutely continuous with respect to the Lebesgue measure, has
a compact support, and its PDF is continuous and bounded. (e) The
supports $\mathscr{S}_{\left(D,X\right)}$ and $\mathscr{S}_{\left(Z,X\right)}$
are equal to $\{0,1\}\times\mathscr{S}_{X}$. (f) For all $x\in\mathscr{S}_{X}$,
$s\left(0,x\right)<s\left(1,x\right)$. (g) For $D_{zx}\coloneqq\mathbbm{1}\left(\eta\leq s\left(z,x\right)\right)$,
the complier group is given by $D_{0x}<D_{1x}$. We assume that for
all $\left(d,x\right)\in\mathscr{S}_{\left(D,X\right)}$, $\mathscr{S}_{g\left(d,x,\epsilon\right)\mid X=x,D_{0x}<D_{1x}}=\mathscr{S}_{g\left(d,x,\epsilon\right)\mid X=x}$.
(h) For all $d\in\{0,1\}$, the conditional distribution of $g\left(d,x,\epsilon\right)$
given $X=x$ and $D_{0x}<D_{1x}$ has a bounded away from zero density
$f_{dx\mid C_{x}}$. (i) For all $x\in\mathscr{S}_{X}$ and $d\in\left\{ 0,1\right\} $,
the conditional distributions of $g\left(d,x,\epsilon\right)$ has
the support $\mathscr{S}_{g\left(d,x,\epsilon\right)\mid X=x}=\left[\underline{y}_{dx},\overline{y}_{dx}\right]$
with known boundaries $-\infty<\underline{y}_{dx}<\overline{y}_{dx}<+\infty$.
(j) The data $\left\{ W_{i}\coloneqq\left(Y_{i},D_{i},X_{i}^{\top},Z_{i}\right)^{\top}:i=1,\ldots,n\right\} $
are i.i.d. observations on $W=\left(Y,D,X^{\top},Z\right)^{\top}$.
(k) $X$ is discretely distributed and $\mathscr{S}_{X}$ is finite.
\end{assumption}
In the above assumption, the continuity and monotonicity conditions
in (a), the standard instrument exogeneity assumption in (b), the
instrument relevance condition in (c), the absolute continuity condition
in (d), (e,f), as well as the equality of the supports condition in
(g) are imposed for identification. The assumption that $\epsilon$
is scalar-valued and the condition in (a) impose rank invariance on
the potential outcomes. See Section 2.1 of VX. Parts (d,e) of the
assumption are mild regularity conditions. Parts (c,f) and the latent
index assumption on the selection equation imply that $\mathrm{Pr}\left[D=1\mid Z=1,X=x\right]>\mathrm{Pr}\left[D=1\mid Z=0,X=x\right]$.
Under (c,f) and the latent index assumption, we have $D_{0x}\leq D_{1x}$,
for all $x\in\mathscr{S}_{X}$. Clearly, the model satisfies the LATE
independence and monotonicity assumptions (see, e.g., \citealp[Section 4]{Vytlacil2002}).
See \citet{Kitagawa2015} for testable implications.
In part (g), it is assumed that conditionally on $X=x$, the support
of the conditional distribution of $g\left(d,x,\epsilon\right)$ in
the complier group $D_{0x}<D_{1x}$ is the same as that of the conditional
distribution of $g\left(d,x,\epsilon\right)$ given $X=x$. VX argues
that (g) is satisfied if the conditional distribution of $\left(\epsilon,\eta\right)$
given $X=x$ has a rectangular support, for all $x\in\mathscr{S}_{X}$.
Note that the identification result $\mathscr{S}_{g\left(d,x,\epsilon\right)\mid X=x}=\mathscr{S}_{Y\mid D=d,X=x}$
in VX, together with (a), implies that $\mathscr{S}_{\epsilon\mid X=x}=\mathscr{S}_{\epsilon\mid D=d,X=x}$,
for all $d\in\left\{ 0,1\right\} $. The rest of the conditions are
imposed for estimation. As in FVX, part (k) restricts the estimation
framework to discretely distributed covariates $X$. Under these assumptions,
the conditional distribution of $Y$ or $\epsilon$ given $\left(D,X,Z\right)$
is absolutely continuous with respect to the Lebesgue measure and
admits a continuous and bounded Lebesgue density. Under Assumption
\ref{assu: DGP1}, $f_{dx\mid C_{x}}$ is also continuous and bounded.\footnote{Let $F_{dx\mid C_{x}}$ denote the conditional CDF of $g\left(d,x,\epsilon\right)$
given $X=x$ and $D_{0x}<D_{1x}$. Then it is clear that $F_{dx\mid C_{x}}\left(y\right)=\mathrm{Pr}\left[\epsilon\leq g^{-1}\left(d,x,y\right)\mid s\left(0,x\right)<\eta\leq s\left(1,x\right),X=x\right]$.}
Theorem 1 of FVX shows the asymptotic properties of the estimated
counterfactual mappings $\widehat{\phi}_{dx}$ under Assumption \ref{assu: DGP1}.
Let $\varDelta_{x}\left(\cdot\right)\coloneqq g\left(1,x,\cdot\right)-g\left(0,x,\cdot\right)$,
which is continuously differentiable under Assumption \ref{assu: DGP1}.
FVX assumes in their Assumption 5(i) that the conditional density
of $\varDelta=\varDelta_{x}\left(\epsilon\right)$ given $X=x$ exists
and is $P$-times continuously differentiable. Without imposing further
restrictions, Assumptions \ref{assu: DGP1} alone does not guarantee
that the distribution of $\varDelta=\varDelta_{x}\left(\epsilon\right)$
is absolutely continuous with respect to the Lebesgue measure.\footnote{E.g., the distribution has a mass point if $\varDelta_{x}\left(\cdot\right)$
is constant on some sub-interval of $\mathscr{S}_{\epsilon\mid X=x}$.} The assumption below provides mild sufficient conditions for the
existence and differentiability of the Lebesgue density of the ITE
$\varDelta$ given $X=x$ (see Lemma \ref{lem:differentiability of the PDF}
in Appendix \ref{sec:Appendix A}). Let $F_{V\mid W}\left(\cdot\mid w\right)$
and $f_{V\mid W}\left(\cdot\mid w\right)$ denote the conditional
cumulative distribution function (CDF) and PDF of $V$ given $W=w$,
respectively.
\begin{assumption}
[Existence and differentiability of the conditional PDF of the ITE]\label{assu: DGP4}(a)
For all $x\in\mathscr{S}_{X}$, the conditional CDF $F_{\epsilon\eta\mid X}\left(\cdot\mid x\right)$
of $\left(\epsilon,\eta\right)$ given $X=x$ and $g\left(d,x,\cdot\right)$
are both $\left(P+1\right)$-times continuously differentiable. (b)
There is a partition of $\mathscr{S}_{\epsilon\mid X=x}=\left[\underline{\epsilon}_{x},\overline{\epsilon}_{x}\right]$,
$\underline{\epsilon}_{x}=\epsilon_{x,0}<\epsilon_{x,1}<\cdots<\epsilon_{x,m}=\overline{\epsilon}_{x}$
with $\mathscr{S}_{\epsilon\mid X=x}=\bigcup_{j=1}^{m}\left[\epsilon_{x,j-1},\epsilon_{x,j}\right]$,
such that $\varDelta_{x}\left(\cdot\right)$ is piecewise monotone:
for all $j=1,...,m$, the restriction of $\varDelta_{x}\left(\cdot\right)$
on $\left[\epsilon_{x,j-1},\epsilon_{x,j}\right]$, $\left.\varDelta_{x,j}\left(\cdot\right)\coloneqq\varDelta_{x}\right|_{\left[\epsilon_{x,j-1},\epsilon_{x,j}\right]}$
is strictly monotone. (c) Let $\varDelta_{x}\left(\left(\epsilon_{x,j-1},\epsilon_{x,j}\right)\right)\coloneqq\left\{ \varDelta_{x}\left(e\right):e\in\left(\epsilon_{x,j-1},\epsilon_{x,j}\right)\right\} $
denote the image of $\varDelta_{x,j}(\cdot)$. We assume that $\varDelta_{x}\left(\left(\epsilon_{x,0},\epsilon_{x,1}\right)\right)=\cdots=\varDelta_{x}\left(\left(\epsilon_{x,m-1},\epsilon_{x,m}\right)\right)$.
\end{assumption}
The smoothness assumption imposed by (a) with $P\geq1$ is stronger
than that imposed by Assumption \ref{assu: DGP1}(a,d). Under (a),
$f_{dx\mid C_{x}}$ is $P$-times continuously differentiable. The
piecewise monotonicity condition in (b) is easily satisfied if $\varDelta_{x}$
has finitely many local extrema on $\left[\underline{\epsilon}_{x},\overline{\epsilon}_{x}\right]$.\footnote{Since it was assumed in Assumption \ref{assu: DGP1} that $\varDelta_{x}$
is continuously differentiable, this condition is satisfied if the
set of zeros of the continuous derivative function $\varDelta_{x}'$,
$\left\{ e\in\left[\underline{\epsilon}_{x},\overline{\epsilon}_{x}\right]:\varDelta_{x}'\left(e\right)=0\right\} $,
contains only isolated points.} Parts (a,b) of the assumption guarantee the existence of the Lebesgue
density $f_{\varDelta\mid X}\left(\cdot\mid x\right)$. Note that
the knowledge of the partition in (b) is not required for estimation
or inference. Part (c) rules out discontinuities in the interior of
$\mathscr{S}_{\varDelta\mid X=x}$. See the proof of Lemma \ref{lem:differentiability of the PDF}
in Appendix \ref{sec:Appendix A} for more details. We are unaware
of any weaker conditions that could be imposed on $\varDelta_{x}$
to guarantee the existence and differentiability of the conditional
PDF of $\varDelta=\varDelta_{x}\left(\epsilon\right)$ given $X=x$.
Application of kernel-based nonparametric techniques is complicated
by the bandwidth selection issue. A common practice in applied work
is using a data-dependent bandwidth approximating some underlying
deterministic bandwidth. We allow the bandwidth $\widehat{h}=\widehat{h}_{n}$
used in the implementation to be data-dependent and, following \citet{Li:2010dy},
assume that $\widehat{h}$ is a consistent estimator of some deterministic
bandwidth sequence $h=h_{n}\downarrow0$ in the sense that $\widehat{h}_{n}/h_{n}\rightarrow_{p}1$.
To simplify the notation, we suppress the dependence of the bandwidths
on $n$. Formally, we make the following assumption.
\begin{assumption}
\label{assu:h_hat}$\mathrm{Pr}\left[\left|\widehat{h}/h-1\right|>\varepsilon_{n}\right]\leq\delta_{n}$
for some deterministic bandwidth $h$ and positive sequences $\varepsilon_{n},\delta_{n}\downarrow0$.\footnote{Assumption \ref{assu:h_hat} is equivalent to requiring $\widehat{h}/h\rightarrow_{p}1$.
It is clear that Assumption \ref{assu:h_hat} implies $\widehat{h}/h\rightarrow_{p}1$.
On the other hand, if $\widehat{h}/h\rightarrow_{p}1$, $\varepsilon_{n}=\delta_{n}\downarrow0$
can be taken to be the Ky Fan metric between $\widehat{h}/h$ and
$1$, which converges to 0 as $n\uparrow\infty$. See \citet[Theorem 9.2.2]{dudley2002real}.}
\end{assumption}
The deterministic bandwidth assumption ($\varepsilon_{n}=\delta_{n}=0$)
is nested as a special case. Clearly, $\widehat{h}/h-1=O_{p}\left(\varepsilon_{n}\right)$,
under Assumption \ref{assu:h_hat}. E.g., as in FVX, one can consider
a feasible version of the Silverman rule-of-thumb (ROT) bandwidth
by setting $\widehat{h}=\widehat{h}^{\mathsf{rot}}\coloneqq C_{K}\cdot\widehat{\sigma}_{\varDelta\mid X=x}\cdot n_{x}^{-1/5}$,
where $\widehat{\sigma}_{\varDelta\mid X=x}$ is the sample analogue
of the standard deviation $\sigma_{\varDelta\mid X=x}\coloneqq\sqrt{\mathrm{Var}\left[\varDelta\mid X=x\right]}$
computed using the pseudo (estimated) ITEs, $C_{K}$ is a known constant
that depends only on the kernel function $K$, and $n_{x}\coloneqq\sum_{i=1}^{n}\mathbbm{1}\left(X_{i}=x\right)$.\footnote{The ROT bandwidth is a parametric estimator of the asymptotic mean
integrated squared error (AMISE) optimal bandwidth for $\widetilde{f}_{\varDelta\mid X}\left(v\mid x;h\right)$
defined by (\ref{eq:f_til definition}) under $P=2$. See \citet{li2007nonparametric}.
Theorem \ref{thm:asymptotic normality} shows that the asymptotic
mean squared error (AMSE) of $\widehat{f}_{\varDelta\mid X}\left(v\mid x;h\right)$
is given by $\mathrm{AMSE}\left(v\mid x\right)\coloneqq f_{\varDelta\mid X}^{\left(P\right)}\left(v\mid x\right)^{2}\mu_{K,P}^{2}h^{2P}+\mathscr{V}\left(v\mid x\right)/\left(nh\right)$.
Hence, the bandwidth that minimizes the AMISE $\int_{I_{x}}\mathrm{AMSE}\left(v\mid x\right)\mathrm{d}v$
is also a multiple of $n^{-1/\left(2P+1\right)}$.} Let $a_{n}\wasypropto b_{n}$ denote $a_{n}=C\cdot b_{n}$ for some
constant $C>0$. Let $p_{x}\coloneqq\mathrm{Pr}\left[X=x\right]$.
By Lemma \ref{lem:Lemma 1 phi linear representation}, one can show
that such a choice satisfies Assumption \ref{assu:h_hat} with $h=C_{K}\cdot\left(\sigma_{\varDelta\mid X=x}p_{x}^{-1/5}\right)\cdot n^{-1/5}$,
$\varepsilon_{n}\wasypropto\sqrt{\mathrm{log}\left(n\right)/n}$ and
$\delta_{n}\wasypropto n^{-1}$. We also assume that the kernel function
$K$ is of order $P\geq2$.
\begin{assumption}
\label{assu:kernel}(a) $K$ is symmetric, compactly supported on
$\left[-1,1\right]$ and twice continuously differentiable on $\mathbb{R}$
with Lipschitz derivatives. (b) $\int K\left(u\right)\mathrm{d}u=1$.
(c) $\int u^{k}K\left(u\right)\mathrm{d}u=0$, for all $k=1,2,...,P-1$
($P\geq2$), and $\int u^{P}K\left(u\right)\mathrm{d}u\ne0$.
\end{assumption}
\subsection{Rate of convergence and asymptotic distribution}
In this section, we first derive a linearization for the FVX estimator.
The result is given by equation (\ref{eq:linearization f_hat - f three terms})
below and serves as the basis for establishing the asymptotic properties
of the FVX estimator. The latter are presented below as Theorems \ref{thm:uniform rate}
and \ref{thm:asymptotic normality} (the rate of convergence and asymptotic
distribution, respectively).
Let $\widetilde{f}_{\varDelta\mid X}\left(v\mid x;b\right)$ denote
the infeasible estimator of the density $f_{\varDelta\mid X}\left(v\mid x\right)$
that uses the true latent ITEs:
\begin{equation}
\widetilde{f}_{\varDelta\mid X}\left(v\mid x;b\right)\coloneqq\frac{\sum_{i=1}^{n}\frac{1}{b}K\left(\frac{\varDelta_{i}-v}{b}\right)\mathbbm{1}\left(X_{i}=x\right)}{\sum_{i=1}^{n}\mathbbm{1}\left(X_{i}=x\right)}.\label{eq:f_til definition}
\end{equation}
Let $I_{x}$ denote an inner closed sub-interval of $\mathscr{S}_{\Delta\mid X=x}$,
$\underline{h}\coloneqq\left(1-\varepsilon_{n}\right)h$ and $\overline{h}\coloneqq\left(1+\varepsilon_{n}\right)h$.
Under Assumption \ref{assu:h_hat}, $\widehat{h}\in\left[\underline{h},\overline{h}\right]$
with probability $1-\delta_{n}$. For $\left(v,b\right)\in I_{x}\times\left[\underline{h},\overline{h}\right]$,
we decompose the estimation error $\widehat{f}_{\varDelta\mid X}\left(v\mid x;b\right)-f_{\varDelta\mid X}\left(v\mid x\right)$
into that of the infeasible estimator $\widetilde{f}_{\varDelta\mid X}\left(v\mid x;b\right)-f_{\varDelta\mid X}\left(v\mid x\right)$
and the difference between the feasible and infeasible estimators
$\widehat{f}_{\varDelta\mid X}\left(v\mid x;b\right)-\widetilde{f}_{\varDelta\mid X}\left(v\mid x;b\right)$.
We show that the former satisfies
\begin{multline}
\widetilde{f}_{\varDelta\mid X}\left(v\mid x;b\right)-f_{\varDelta\mid X}\left(v\mid x\right)=p_{x}^{-1}\left(\widetilde{f}_{\varDelta X}\left(v,x;b\right)-m_{\varDelta X}\left(v,x;b\right)\right)\\
+p_{x}^{-1}\left(m_{\varDelta X}\left(v,x;b\right)-f_{\varDelta X}\left(v,x\right)\right)+O_{p}\left(\sqrt{\frac{\mathrm{log}\left(n\right)}{n}}\right),\label{eq:f_til - f linearization}
\end{multline}
where $\widetilde{f}_{\varDelta X}\left(v,x;b\right)\coloneqq\widetilde{f}_{\varDelta\mid X}\left(v\mid x;b\right)\widehat{p}_{x}$
with $\widehat{p}_{x}\coloneqq n^{-1}\sum_{i=1}^{n}\mathbbm{1}\left(X_{i}=x\right)$
is the infeasible estimator of the joint density function $f_{\varDelta X}\left(v,x\right)\coloneqq f_{\varDelta\mid X}\left(v\mid x\right)p_{x}$,
$m_{\varDelta X}\left(v,x;b\right)\coloneqq\mathrm{E}\left[\widetilde{f}_{\varDelta X}\left(v,x;b\right)\right]$,
and the remainder term is uniform in $\left(v,b\right)\in I_{x}\times\left[\underline{h},\overline{h}\right]$.
Note that $m_{\varDelta X}\left(v,x;b\right)-f_{\varDelta X}\left(v,x\right)$
is the bias term that depends on the bandwidth $b$. Let $\mu_{K,P}\coloneqq\left(\int u^{P}K\left(u\right)\mathrm{d}u\right)/P!$,
and let $f_{\varDelta\mid X}^{\left(P\right)}\left(v\mid x\right)\coloneqq\partial^{P}f_{\varDelta\mid X}\left(v\mid x\right)/\partial v^{P}$
denote the derivatives of the conditional PDF. It follows from standard
arguments for kernel density estimators (see, e.g., \citealp{Newey:1994jb})
that
\begin{equation}
m_{\varDelta X}\left(v,x;b\right)-f_{\varDelta X}\left(v,x\right)=f_{\varDelta\mid X}^{\left(P\right)}\left(v\mid x\right)p_{x}\mu_{K,P}b^{P}+o\left(h^{P}\right),\label{eq:infeasible decompose}
\end{equation}
uniformly in $\left(v,b\right)\in I_{x}\times\left[\underline{h},\overline{h}\right]$.
For a deterministic bandwidth sequence $h$ such that $nh\uparrow\infty$,
it follows from standard arguments that $\sqrt{nh}\left(\widetilde{f}_{\varDelta X}\left(v,x;h\right)-m_{\varDelta X}\left(v,x;h\right)\right)$
is asymptotically normal.
Let $K'$ denote the derivative of the kernel function $K$. Denote
$\widehat{f}_{\varDelta X}\left(v,x;b\right)\coloneqq\widehat{f}_{\varDelta X}\left(v\mid x;b\right)\widehat{p}_{x}$.
We approximate $\widehat{f}_{\varDelta X}\left(v,x;b\right)-\widetilde{f}_{\varDelta X}\left(v,x;b\right)$
by $\left(nb^{2}\right)^{-1}\sum_{i=1}^{n}K'\left(\left(\varDelta_{i}-v\right)/b\right)\left(\widehat{\varDelta}_{i}-\varDelta_{i}\right)\mathbbm{1}\left(X_{i}=x\right)$.
The first-stage estimation errors $\widehat{\phi}_{dX_{i}}^{\left(-i\right)}\left(Y_{i}\right)-\phi_{dX_{i}}\left(Y_{i}\right)$
in $\widehat{\varDelta}_{i}-\varDelta_{i}$ can be approximated using
its linear representation (see Theorem 1 of FVX and Lemma \ref{lem:Lemma 1 phi linear representation}).
After recalling that $f_{dx\mid C_{x}}$ is the conditional PDF of
$g\left(d,x,\epsilon\right)$ given $X=x$ in the complier group,
we define:
\begin{eqnarray}
\zeta_{dx}\left(y\right) & \coloneqq & f_{dx\mid C_{x}}\left(y\right)\left(\mathrm{Pr}\left[D=d\mid Z=1,X=x\right]-\mathrm{Pr}\left[D=d\mid Z=0,X=x\right]\right),\nonumber \\
R_{d'x}\left(y\right) & \coloneqq & \mathrm{Pr}\left[Y\leq\phi_{dx}\left(y\right),D=d\mid X=x\right]+\mathrm{Pr}\left[Y\leq y,D=d'\mid X=x\right],\nonumber \\
q_{dx}\left(W_{i},W_{j}\right) & \coloneqq & \frac{\mathbbm{1}\left(D_{i}=d',X_{i}=x\right)}{\zeta_{dx}\left(\phi_{dx}\left(Y_{i}\right)\right)}\nonumber \\
& & \times\left\{ \mathbbm{1}\left(Y_{j}\leq\phi_{dx}\left(Y_{i}\right),D_{j}=d\right)+\mathbbm{1}\left(Y_{j}\leq Y_{i},D_{j}=d'\right)-R_{d'x}\left(Y_{i}\right)\right\} ,\label{eq:q_dx definition}\\
q_{x}\left(W_{i},W_{j}\right) & \coloneqq & q_{1x}\left(W_{i},W_{j}\right)-q_{0x}\left(W_{i},W_{j}\right),\nonumber \\
p_{zx} & \coloneqq & \mathrm{Pr}\left[Z=z,X=x\right],\nonumber \\
\pi_{x}\left(Z_{i},X_{i}\right) & \coloneqq & \frac{\mathbbm{1}\left(Z_{i}=0,X_{i}=x\right)}{p_{0x}}-\frac{\mathbbm{1}\left(Z_{i}=1,X_{i}=x\right)}{p_{1x}}.\label{eq:pi_x definition}
\end{eqnarray}
Using the above definitions, we can write that the difference between
the feasible and infeasible estimators as a \textit{U}-statistic with
a kernel that depends on the bandwidth:
\begin{equation}
\widehat{f}_{\varDelta\mid X}\left(v\mid x;b\right)-\widetilde{f}_{\varDelta\mid X}\left(v\mid x;b\right)=p_{x}^{-1}\frac{1}{n\left(n-1\right)}\sum_{i=1}^{n}\sum_{j\neq i}\mathcal{G}_{x}\left(W_{i},W_{j},v;b\right)+O_{p}\left(\frac{\mathrm{log}\left(n\right)}{nh^{2}}+\frac{\mathrm{log}\left(n\right)^{3/4}}{n^{3/4}h}\right),\label{eq:f_hat - f_til U stat expansion}
\end{equation}
uniformly in $\left(v,b\right)\in I_{x}\times\left[\underline{h},\overline{h}\right]$,
where
\begin{equation}
\mathcal{G}_{x}\left(W_{i},W_{j},v;b\right)\coloneqq\frac{1}{b^{2}}K'\left(\frac{\varDelta_{i}-v}{b}\right)q_{x}\left(W_{i},W_{j}\right)\pi_{x}\left(Z_{j},X_{j}\right).\label{eq:G kernel definition}
\end{equation}
By Assumption \ref{assu: DGP1}(b) and (\ref{eq:q_dx simplify}) in
Appendix \ref{sec:Appendix A}, $\mathcal{G}_{x}^{\left[2\right]}\left(w,v;b\right)\coloneqq\mathrm{E}\left[\mathcal{G}_{x}\left(w,W,v;b\right)\right]=0$,
for all $w$ and $\mathrm{E}\left[\mathcal{G}_{x}\left(W_{i},W_{j},v;b\right)\right]=0$,
for all $\ensuremath{i\neq j}$. The leading term (or the H$\acute{\mathrm{a}}$jek
projection) in the Hoeffding decomposition of the \textit{U}-statistic
is given by $\mathcal{G}_{x}^{\left[1\right]}\left(w,v;b\right)\coloneqq\mathrm{E}\left[\mathcal{G}_{x}\left(W,w,v;b\right)\right]$.
Therefore, the Hoeffding decomposition is given by
\begin{multline}
\frac{1}{n\left(n-1\right)}\sum_{i=1}^{n}\sum_{j\neq i}\mathcal{G}_{x}\left(W_{i},W_{j},v;b\right)=\frac{1}{n}\sum_{i=1}^{n}\mathcal{G}_{x}^{\left[1\right]}\left(W_{i},v;b\right)\\
+\frac{1}{n\left(n-1\right)}\sum_{i=1}^{n}\sum_{j\neq i}\left\{ \mathcal{G}_{x}\left(W_{i},W_{j},v;b\right)-\mathcal{G}_{x}^{\left[1\right]}\left(W_{j},v;b\right)\right\} .\label{eq:G Hoeffding}
\end{multline}
By definition and since $\mathcal{G}_{x}^{\left[2\right]}\left(w,v;b\right)=0$,
the \textit{U}-statistic $\left(n\left(n-1\right)\right)^{-1}\sum_{i=1}^{n}\sum_{j\neq i}\mathcal{G}_{x}\left(W_{i},W_{j},v;b\right)$
is non-degenerate (\citealp{Chen2020jackknife}) if $\mathrm{Var}\left[\mathcal{G}_{x}^{\left[1\right]}\left(W,v;b\right)\right]>0$.
In the proof of Lemma \ref{lem:lemma 3}, we show that the condition
holds for all $\left(v,b\right)\in I_{x}\times\left[\underline{h},\overline{h}\right]$
when $n$ is sufficiently large. We also show that the second term
on the right-hand side of (\ref{eq:G Hoeffding}) is dominated by
the first term. See Remark \ref{rmk: proof of normality} below. Combining
the result with (\ref{eq:f_til - f linearization}) and (\ref{eq:f_hat - f_til U stat expansion}),
we can write the estimation error of the feasible estimator as
\begin{multline}
\widehat{f}_{\varDelta\mid X}\left(v\mid x;b\right)-f_{\varDelta\mid X}\left(v\mid x\right)=p_{x}^{-1}\left(\widetilde{f}_{\varDelta X}\left(v,x;b\right)-m_{\varDelta X}\left(v,x;b\right)\right)+p_{x}^{-1}\frac{1}{n}\sum_{i=1}^{n}\mathcal{G}_{x}^{\left[1\right]}\left(W_{i},v;b\right)\\
+p_{x}^{-1}\left(m_{\varDelta X}\left(v,x;b\right)-f_{\varDelta X}\left(v,x\right)\right)+O_{p}\left(\sqrt{\frac{\mathrm{log}\left(n\right)}{n}}+\frac{\mathrm{log}\left(n\right)}{nh^{2}}+\frac{\mathrm{log}\left(n\right)^{3/4}}{n^{3/4}h}\right),\label{eq:linearization f_hat - f three terms}
\end{multline}
where the remainder term is uniform in $\left(v,b\right)\in I_{x}\times\left[\underline{h},\overline{h}\right]$.
\sloppy The asymptotic variance of the FVX estimator $\widehat{f}_{\varDelta\mid X}\left(v\mid x;h\right)$
under the deterministic bandwidth sequence is determined by the sum
of $\left(\widetilde{f}_{\varDelta X}\left(v,x;h\right)-m_{\varDelta X}\left(v,x;h\right)\right)/p_{x}$
and $n^{-1}\sum_{i=1}^{n}\mathcal{G}_{x}^{\left[1\right]}\left(W_{i},v;h\right)/p_{x}$.
The first term comes from the infeasible estimator $\widetilde{f}_{\varDelta X}\left(v,x;h\right)$
as in equation (\ref{eq:f_til - f linearization}). The second term
reflects the influence of the estimation of ITEs in the first stage.
We show that these terms are both of order $O_{p}\left(\sqrt{\mathrm{log}\left(n\right)/\left(nh\right)}\right)$
uniformly in $v\in I_{x}$ and independent of each other. Consequently,
the first-stage estimation errors unambiguously add to the asymptotic
variance of $\widehat{f}_{\varDelta\mid X}\left(v\mid x;h\right)$
and their contribution must be taken into account for valid inference.\footnote{\label{fn:phi_hat - phi}Let $\widehat{\phi}_{dx}\left(y\right)$
be the leave-in version of $\widehat{\phi}_{dx}^{\left(-i\right)}\left(y\right)$
(i.e., $\widehat{\phi}_{dx}\left(y\right)$ minimizes the sample analogue
of $Q_{dx}\left(\cdot;y\right)$). The crucial observation is that
the linearization of $\widehat{\phi}_{dx}\left(y\right)-\phi_{dx}\left(y\right)$
derived in FVX (also see Lemma \ref{lem:Lemma 1 phi linear representation})
is discontinuous in both $W_{i}$ and $y$. As a result, the difference
between the feasible and infeasible estimators $\widehat{f}_{\varDelta\mid X}\left(v\mid x;h\right)-\widetilde{f}_{\varDelta\mid X}\left(v\mid x;h\right)$
converges at a rate slower than $n^{-1/2}$. One can show that the
difference would be of order $O_{p}(n^{-1/2})$ if the linearization
were smooth. See the proof of Lemma \ref{lem:lemma 2} for more details
on how the linearization is involved in the \textit{U}-statistic representation
given by $\left(n\left(n-1\right)\right)^{-1}\sum_{i=1}^{n}\sum_{j\neq i}\mathcal{G}_{x}\left(W_{i},W_{j},v;b\right)$.} Equation (\ref{eq:infeasible decompose}) implies that $m_{\varDelta X}\left(v,x;h\right)-f_{\varDelta X}\left(v,x\right)=O\left(h^{P}\right)$
uniformly in $v\in I_{x}$. These results provide the uniform rate
of convergence of $\widehat{f}_{\varDelta\mid X}\left(\cdot\mid x;h\right)$.
The bias expansion in (\ref{eq:infeasible decompose}) and the linearization
in (\ref{eq:linearization f_hat - f three terms}) are also valid
for a continuum $\left[\underline{h},\overline{h}\right]$ of bandwidths.
Since $\mathrm{Pr}\left[\widehat{h}\in\left[\underline{h},\overline{h}\right]\right]>1-\delta_{n}$
under Assumption \ref{assu:h_hat}, (\ref{eq:infeasible decompose})
and (\ref{eq:linearization f_hat - f three terms}) with $b$ replaced
by $\widehat{h}$ still hold. We show that the first two terms on
the right-hand side of the decomposition in (\ref{eq:linearization f_hat - f three terms})
are of the same $O_{p}\left(\sqrt{\mathrm{log}\left(n\right)/\left(nh\right)}\right)$
order uniformly in $\left(v,b\right)\in I_{x}\times\left[\underline{h},\overline{h}\right]$.
These results show that the uniform convergence rate remains the same
if $h$ is replaced by its estimator $\widehat{h}$.
We summarize the above results in Theorem \ref{thm:uniform rate}
below, which is a refinement of Theorem 2 in FVX. In Appendix \ref{sec:Appendix A},
we prove a stronger version (Theorem \ref{thm: uniform rate appendix})
of Theorem \ref{thm:uniform rate}. The latter establishes the non-asymptotic
deviation bounds for the uniform estimation errors of $\widehat{f}_{\varDelta\mid X}\left(\cdot\mid x;h\right)$
and $\widehat{f}_{\varDelta\mid X}\left(\cdot\mid x;\widehat{h}\right)$
and is used in the proof of Theorem \ref{thm:confidence band}.
For a function $f:A\rightarrow\mathbb{R}$, let $\left\Vert f\right\Vert _{A}\coloneqq\mathrm{sup}_{x\in A}\left|f\left(x\right)\right|$
denote the sup-norm of $f$ on $A$. For a subset $A$ in the Euclidean
space, let $\mathrm{Int}\left(A\right)$ denote the interior of $A$.
\begin{thm}
\label{thm:uniform rate}Suppose that Assumptions \ref{assu: DGP1}-\ref{assu:kernel}
hold, and the deterministic bandwidth $h$ is such that $\mathrm{log}\left(n\right)/\left(nh^{3}\right)\downarrow0$.
Then, for any $x\in\mathscr{S}_{X}$ and compact $I_{x}\subseteq\mathrm{Int}\left(\mathscr{S}_{\Delta\mid X=x}\right)$,
\[
\left\Vert \widehat{f}_{\varDelta\mid X}\left(\cdot\mid x;h\right)-f_{\varDelta\mid X}\left(\cdot\mid x\right)\right\Vert _{I_{x}}=O_{p}\left(\sqrt{\frac{\mathrm{log}\left(n\right)}{nh}}+h^{P}\right)
\]
and
\[
\left\Vert \widehat{f}_{\varDelta\mid X}\left(\cdot\mid x;\widehat{h}\right)-f_{\varDelta\mid X}\left(\cdot\mid x\right)\right\Vert _{I_{x}}=O_{p}\left(\sqrt{\frac{\mathrm{log}\left(n\right)}{nh}}+h^{P}\right).
\]
\end{thm}
\begin{rembold}In comparison, Theorem 2 of FVX has a slower $O_{p}\left(\sqrt{\mathrm{log}\left(n\right)/\left(nh^{2}\right)}+h^{P}\right)$
convergence rate. Theorem \ref{thm:uniform rate} implies that the
FVX and infeasible estimators of $f_{\varDelta\mid X}\left(\cdot\mid x\right)$
have the same uniform convergence rate. Moreover, the convergence
rate is unaffected by the estimation of the bandwidth. The optimal
bandwidth rate that leads to the fastest possible convergence rate
is of order $\left(\mathrm{log}\left(n\right)/n\right)^{1/\left(2P+1\right)}$.
Hence, both the FVX and infeasible estimators attain the optimal uniform
convergence rate $\left(\mathrm{log}\left(n\right)/n\right)^{P/\left(2P+1\right)}$.
Note that under our smoothness conditions, any uniformly consistent
estimator cannot converge uniformly at a rate faster than $\left(\mathrm{log}\left(n\right)/n\right)^{P/\left(2P+1\right)}$
(see \citealp{Stone1982}).\end{rembold}
The next theorem establishes the asymptotic normality of the FVX estimator
and quantifies the contribution of the first-stage estimation errors
to the asymptotic variance. By using (\ref{eq:infeasible decompose})
and the linearization (\ref{eq:linearization f_hat - f three terms})
for a single bandwidth $h$, we show that for any fixed $v\in I_{x}$,
asymptotic normality holds for $\sqrt{nh}\left(\widehat{f}_{\varDelta\mid X}\left(v\mid x;h\right)-f_{\varDelta\mid X}\left(v\mid x\right)-f_{\varDelta\mid X}^{\left(P\right)}\left(v\mid x\right)\mu_{K,P}h^{P}\right)$.
By using the uniform-in-bandwidth approximation ((\ref{eq:f_til - f linearization})
- (\ref{eq:f_hat - f_til U stat expansion})) of $\widehat{f}_{\varDelta\mid X}\left(v\mid x;b\right)-f_{\varDelta\mid X}\left(v\mid x\right)$
and an asymptotic equivalence result (Lemma \ref{lem:lemma random btw}),
we show that the same normality result holds if $h$ is replaced by
its estimator $\widehat{h}$. The result is analogous to those in
\citet{Li:2010dy}.
Let $f_{\epsilon DX}\left(e,d,x\right)\coloneqq f_{\epsilon\mid DX}\left(e\mid d,x\right)\mathrm{Pr}\left[D=d,X=x\right]$
denote the joint density of $(\epsilon,D,X^{\top})^{\top}$.
\begin{thm}
\label{thm:asymptotic normality}Suppose that Assumptions \ref{assu: DGP1}-\ref{assu:kernel}
hold, and $h\wasypropto n^{-\lambda}$ with $1/\left(2P+1\right)\leq\lambda<1/3$.
Then, for any $x\in\mathscr{S}_{X}$ and $v$ in a compact sub-interval
$I_{x}$ of $\mathscr{S}_{\Delta\mid X=x}$,
\[
\sqrt{nh}\left(\widehat{f}_{\varDelta\mid X}\left(v\mid x;h\right)-f_{\varDelta\mid X}\left(v\mid x\right)-f_{\varDelta\mid X}^{\left(P\right)}\left(v\mid x\right)\mu_{K,P}h^{P}\right)\to_{d}\mathrm{N}\left(0,\mathscr{V}\left(v\mid x\right)\right)
\]
and if $\varepsilon_{n}=o\left(\mathrm{log}\left(n\right)^{-1/2}\right)$
in Assumption \ref{assu:h_hat},
\[
\sqrt{n\widehat{h}}\left(\widehat{f}_{\varDelta\mid X}\left(v\mid x;\widehat{h}\right)-f_{\varDelta\mid X}\left(v\mid x\right)-f_{\varDelta\mid X}^{\left(P\right)}\left(v\mid x\right)\mu_{K,P}\widehat{h}^{P}\right)\to_{d}\mathrm{N}\left(0,\mathscr{V}\left(v\mid x\right)\right),
\]
where $\mathscr{V}\left(v\mid x\right)\coloneqq p_{x}^{-2}\left(\mathscr{V}_{1}\left(v,x\right)+\mathscr{V}_{2}\left(v,x\right)\right)$,
\begin{eqnarray*}
\mathscr{V}_{1}\left(v,x\right) & \coloneqq & f_{\varDelta X}\left(v,x\right)\int K\left(u\right)^{2}\mathrm{d}u,\\
\mathscr{V}_{2}\left(v,x\right) & \coloneqq & \left\{ \sum_{j=1}^{m}\left(\frac{f_{\epsilon DX}\left(\varDelta_{x,j}^{-1}\left(v\right),0,x\right)}{\zeta_{1x}\left(g\left(1,x,\varDelta_{x,j}^{-1}\left(v\right)\right)\right)}-\frac{f_{\epsilon DX}\left(\varDelta_{x,j}^{-1}\left(v\right),1,x\right)}{\zeta_{0x}\left(g\left(0,x,\varDelta_{x,j}^{-1}\left(v\right)\right)\right)}\right)^{2}\frac{f_{\epsilon\mid X}\left(\varDelta_{x,j}^{-1}\left(v\right)\mid x\right)}{\left|\varDelta_{x,j}'\left(\varDelta_{x,j}^{-1}(v)\right)\right|^{3}}\right\} \\
& & \times\left(p_{1x}^{-1}+p_{0x}^{-1}\right)\int K\left(u\right)^{2}\mathrm{d}u,
\end{eqnarray*}
and $\varDelta_{x,j}'$ is the derivative of $\varDelta_{x,j}$ defined
in Assumption \ref{assu: DGP4}.
\end{thm}
\begin{rembold}Under the same assumptions, the infeasible kernel
estimator that uses the true ITEs satisfies
\[
\sqrt{nh}\left(\widetilde{f}_{\varDelta\mid X}\left(v\mid x;h\right)-f_{\varDelta\mid X}\left(v\mid x\right)-f_{\varDelta\mid X}^{\left(P\right)}\left(v\mid x\right)\mu_{K,P}h^{P}\right)\rightarrow_{d}\mathrm{N}\left(0,p_{x}^{-2}\mathscr{V}_{1}\left(v,x\right)\right).
\]
Note that the estimation of ITEs does not affect the leading bias
term.\end{rembold}
\begin{rembold}\label{rmk: proof of normality}In the proof of Lemma
\ref{lem:lemma 3}, we show that $\mathrm{Var}\left[\mathcal{G}_{x}^{\left[1\right]}\left(W,v;b\right)\right]=b^{-1}\left(\mathscr{V}_{2}\left(v,x\right)+o\left(1\right)\right)$
uniformly in $\left(v,b\right)\in I_{x}\times\left[\underline{h},\overline{h}\right]$.
The proof of this result and derivation of the form of $\mathscr{V}_{2}\left(v,x\right)$
crucially rely on Assumption \ref{assu: DGP4}(b,c). It is clear from
the definition of $\zeta_{dx}$, the fact that $\mathscr{S}_{\epsilon\mid X=x}=\mathscr{S}_{\epsilon\mid D=d,X=x}$,
and Assumption \ref{assu: DGP1}(g) that for all $j$,
\begin{equation}
\left(\frac{f_{\epsilon DX}\left(\varDelta_{x,j}^{-1}\left(v\right),0,x\right)}{\zeta_{1x}\left(g\left(1,x,\varDelta_{x,j}^{-1}\left(v\right)\right)\right)}-\frac{f_{\epsilon DX}\left(\varDelta_{x,j}^{-1}\left(v\right),1,x\right)}{\zeta_{0x}\left(g\left(0,x,\varDelta_{x,j}^{-1}\left(v\right)\right)\right)}\right)^{2}\frac{f_{\epsilon\mid X}\left(\varDelta_{x,j}^{-1}\left(v\right)\mid x\right)}{\left|\varDelta_{x,j}'\left(\varDelta_{x,j}^{-1}(v)\right)\right|^{3}}>0\label{eq:V_2 > 0}
\end{equation}
and therefore, $\mathscr{V}_{2}\left(v,x\right)>0$. Moreover, in
the proof of Lemma \ref{lem:lemma random btw}, we show that $\sqrt{nb}\left(n^{-1}\sum_{i=1}^{n}\mathcal{G}_{x}^{\left[1\right]}\left(W_{i},v;b\right)\right)-\sqrt{nh}\left(n^{-1}\sum_{i=1}^{n}\mathcal{G}_{x}^{\left[1\right]}\left(W_{i},v;h\right)\right)=o_{p}\left(1\right)$,
uniformly in $b\in\left[\underline{h},\overline{h}\right]$. Further,
in the proof of Theorem \ref{thm:asymptotic normality}, we show that
$\sqrt{nh}\left(n^{-1}\sum_{i=1}^{n}\mathcal{G}_{x}^{\left[1\right]}\left(W_{i},v;h\right)\right)\rightarrow_{d}\mathrm{N}\left(0,\mathscr{V}_{2}\left(v,x\right)\right)$
and the second term on the right-hand side of (\ref{eq:G Hoeffding})
is $O_{p}\left(\left(nh^{2}\right)^{-1}\right)$. Therefore, in (\ref{eq:G Hoeffding}),
$n^{-1}\sum_{i=1}^{n}\mathcal{G}_{x}^{\left[1\right]}\left(W_{i},v;b\right)$
dominates the second term under our assumption on the rate of $h$.\end{rembold}
The $p_{x}^{-2}\mathscr{V}_{1}\left(v,x\right)$ term in the asymptotic
variance of the FVX estimator is the asymptotic variance of the infeasible
estimator. The $p_{x}^{-2}\mathscr{V}_{2}\left(v,x\right)$ term is
due to the estimation of the ITEs. Thus, the estimation of the ITEs
increases the variance (but not the bias). To illustrate the effect
of estimation of the ITEs numerically, consider the DGP used for the
Monte Carlo simulations in Section \ref{sec:Monte-Carlo-experiments}
with no controls $X$. The treatment status $D$ is determined by
the index model in (\ref{eq:model simulation}) with coefficients
$(\gamma_{0},\gamma_{1})=(-0.5,0.5)$. The kernel function $K$ is
taken to be the triweight kernel. In this case for $v=2$, $\mathscr{V}_{1}\left(v\right)=0.16$
and $\mathscr{V}_{2}\left(v\right)=2.30$. Hence, the contribution
of the ITE estimation errors to the asymptotic variance of the FVX
estimator can be substantial and even exceed the asymptotic variance
of the infeasible estimator.
\section{Inference\label{sec:Robust-inference}}
In this section, we discuss the construction of asymptotically valid
standard errors as well as construction of asymptotically valid UCBs
for $\left\{ f_{\varDelta\mid X}\left(v\mid x\right):v\in I_{x}\right\} $.
We maintain Assumptions \ref{assu: DGP1}, \ref{assu: DGP4} and \ref{assu:kernel}
with $P=2$. We also maintain the assumption that, as in the practical
implementation of many nonparametric econometric methods, the bandwidth
is data-driven and satisfies Assumption \ref{assu:h_hat}.
\subsection{Standard errors\label{subsec:Standard-errors}}
Inference for $f_{\varDelta\mid X}\left(v\mid x\right)$ requires
a consistent estimator of the asymptotic variance term $\mathscr{V}\left(v\mid x\right)$
defined in Theorem \ref{thm:asymptotic normality}. By the same arguments
as those used to establish Theorems \ref{thm:uniform rate} and \ref{thm:asymptotic normality},
one can show that
\begin{multline}
\sqrt{n\widehat{h}}\left(\widehat{f}_{\varDelta\mid X}\left(v\mid x;\widehat{h}\right)-f_{\varDelta\mid X}\left(v\mid x\right)\right)=\frac{1}{\sqrt{n}}\sum_{i=1}^{n}p_{x}^{-1}\left(\mathcal{U}_{x}^{\left[1\right]}\left(W_{i},v;h\right)-\mu_{\mathcal{U}_{x}}\left(v;h\right)\right)\\
+O_{p}\left(\varepsilon_{n}\sqrt{\mathrm{log}\left(n\right)}+\frac{\mathrm{log}\left(n\right)}{\sqrt{nh^{3}}}+\left(\frac{\mathrm{log}\left(n\right)^{3}}{nh^{2}}\right)^{1/4}+\sqrt{\mathrm{log}\left(n\right)h}+\sqrt{nh^{5}}\right),\label{eq:bias corrected linearization}
\end{multline}
where
\begin{eqnarray*}
\mathcal{U}_{x}^{\left[1\right]}\left(w,v;b\right) & \coloneqq & \mathrm{E}\left[\mathcal{U}_{x}\left(W,w,v;b\right)\right],\\
\mathcal{U}_{x}\left(W_{i},W_{j},v;b\right) & \coloneqq & \frac{1}{\sqrt{b}}K\left(\frac{\varDelta_{j}-v}{b}\right)\mathbbm{1}\left(X_{i}=x\right)+\sqrt{b}\cdot\mathcal{G}_{x}\left(W_{i},W_{j},v;b\right),\\
\mu_{\mathcal{U}_{x}}\left(v;b\right) & \coloneqq & \mathrm{E}\left[\mathcal{U}_{x}\left(W_{1},W_{2},v;b\right)\right].
\end{eqnarray*}
Note that $\varDelta_{j}$ can be expressed as a function of $W_{j}$
(see (\ref{eq:Delta})). Also note that the second H$\acute{\mathrm{a}}$jek
projection term $\mathcal{U}_{x}^{\left[2\right]}\left(w,v;b\right)\coloneqq\mathrm{E}\left[\mathcal{U}_{x}\left(w,W,v;b\right)\right]$
is constant and equal to $\mu_{\mathcal{U}_{x}}\left(v;b\right)=\sqrt{b}\cdot m_{\varDelta X}\left(v,x;b\right)$.
Since $\epsilon$ is conditionally independent of $Z$ given $X$,
one can show that the finite-sample variance of the right-hand side
term in (\ref{eq:bias corrected linearization}) is given by $V\left(v\mid x;h\right)\coloneqq p_{x}^{-2}\mathrm{Var}\left[\mathcal{U}_{x}^{\left[1\right]}\left(W,v;h\right)\right]=p_{x}^{-2}V\left(v,x;h\right)$
(see (\ref{eq:cross term zero})), where $V\left(v,x;b\right)\coloneqq V_{1}\left(v,x;b\right)+V_{2}\left(v,x;b\right)$
and
\begin{eqnarray*}
V_{1}\left(v,x;b\right) & \coloneqq & \mathrm{E}\left[\frac{1}{b}K\left(\frac{\varDelta-v}{b}\right)^{2}\mathbbm{1}\left(X=x\right)\right]-b\cdot m_{\varDelta X}\left(v,x;b\right)^{2},\\
V_{2}\left(v,x;b\right) & \coloneqq & \mathrm{E}\left[\frac{1}{b^{3}}K'\left(\frac{\varDelta_{3}-v}{b}\right)q_{x}\left(W_{3},W_{1}\right)K'\left(\frac{\varDelta_{2}-v}{b}\right)q_{x}\left(W_{2},W_{1}\right)\mathbbm{1}\left(X_{1}=x\right)\right]\\
& & \times p_{x}^{-1}\left(p_{1x}^{-1}+p_{0x}^{-1}\right).
\end{eqnarray*}
The plug-in estimator of the $V_{1}\left(v,x;b\right)$ term is given
by
\begin{equation}
\widehat{V}_{1}\left(v,x;b\right)\coloneqq\frac{1}{n}\sum_{i=1}^{n}\frac{1}{b}K\left(\frac{\widehat{\varDelta}_{i}-v}{b}\right)^{2}\mathbbm{1}\left(X_{i}=x\right)-b\cdot\widehat{f}_{\varDelta X}\left(v,x;b\right)^{2}.\label{eq:V_hat_dagger definition}
\end{equation}
Denote $\widehat{p}_{zx}\coloneqq n^{-1}\sum_{i=1}^{n}\mathbbm{1}\left(Z_{i}=z,X_{i}=x\right)$
and $\widehat{p}_{z\mid x}\coloneqq\widehat{p}_{zx}/\widehat{p}_{x}$.
Let
\begin{eqnarray}
\widehat{\zeta}_{1x}\left(y;b_{\zeta}\right) & \coloneqq & \frac{\sum_{i=1}^{n}\frac{1}{b_{\zeta}}K_{\zeta}\left(\frac{Y_{i}-y}{b_{\zeta}}\right)D_{i}\left(\frac{Z_{i}-\widehat{p}_{1\mid x}}{\widehat{p}_{1\mid x}\widehat{p}_{0\mid x}}\right)\mathbbm{1}\left(X_{i}=x\right)}{\sum_{i=1}^{n}\mathbbm{1}\left(X_{i}=x\right)}\nonumber \\
\widehat{\zeta}_{0x}\left(y;b_{\zeta}\right) & \coloneqq & \frac{\sum_{i=1}^{n}\frac{1}{b_{\zeta}}K_{\zeta}\left(\frac{Y_{i}-y}{b_{\zeta}}\right)\left(1-D_{i}\right)\left(\frac{\widehat{p}_{0\mid x}-\left(1-Z_{i}\right)}{\widehat{p}_{1\mid x}\widehat{p}_{0\mid x}}\right)\mathbbm{1}\left(X_{i}=x\right)}{\sum_{i=1}^{n}\mathbbm{1}\left(X_{i}=x\right)}\label{eq:reweighted kernel estimators}
\end{eqnarray}
be the reweighted kernel estimator proposed by \citet{Abadie2002},
where $b_{\zeta}>0$ is the bandwidth and $K_{\zeta}\left(\cdot\right)$
is a second-order kernel. Let
\[
\widehat{R}_{d'x}\left(y\right)\coloneqq\frac{\sum_{i=1}^{n}\left\{ \mathbbm{1}\left(Y_{i}\leq\widehat{\phi}_{dx}\left(y\right),D_{i}=d,X_{i}=x\right)+\mathbbm{1}\left(Y_{i}\leq y,D_{i}=d',X_{i}=x\right)\right\} }{\sum_{i=1}^{n}\mathbbm{1}\left(X_{i}=x\right)}
\]
be the plug-in nonparametric estimator of $R_{d'x}$. The $V_{2}\left(v,x;b\right)$
term can be estimated by a \textit{U}-statistic with an estimated
kernel:
\begin{multline}
\widehat{V}_{2}\left(v,x;b,b_{\zeta}\right)\coloneqq\frac{1}{n\left(n-1\right)\left(n-2\right)}\sum_{i=1}^{n}\sum_{j\neq i}\sum_{k\neq i,\,k\neq j}\frac{1}{b^{3}}K'\left(\frac{\widehat{\varDelta}_{j}-v}{b}\right)\widehat{q}_{x}\left(W_{j},W_{i};b_{\zeta}\right)\\
\times K'\left(\frac{\widehat{\varDelta}_{k}-v}{b}\right)\widehat{q}_{x}\left(W_{k},W_{i};b_{\zeta}\right)\mathbbm{1}\left(X_{i}=x\right)\widehat{p}_{x}^{-1}\left(\widehat{p}_{1x}^{-1}+\widehat{p}_{0x}^{-1}\right),\label{eq:V_hat_ddager definition}
\end{multline}
where $\widehat{q}_{x}\left(W_{i},W_{j};b_{\zeta}\right)\coloneqq\widehat{q}_{1x}\left(W_{i},W_{j};b_{\zeta}\right)-\widehat{q}_{0x}\left(W_{i},W_{j};b_{\zeta}\right)$,
$\widehat{q}_{dx}\left(W_{i},W_{j};b_{\zeta}\right)$ is the plug-in
nonparametric estimator of $q_{dx}\left(W_{i},W_{j}\right)$ defined
in (\ref{eq:q_dx definition}) constructed by replacing $\zeta_{dx}$,
$\phi_{dx}$, and $R_{d'x}$ with their nonparametric estimators $\widehat{\zeta}_{dx}\left(\cdot;b_{\zeta}\right)$,
$\widehat{\phi}_{dx},$ and $\widehat{R}_{d'x}$ respectively.\footnote{It is known that the kernel estimator $\widehat{\zeta}_{dx}\left(y;b_{\zeta}\right)$
is asymptotically biased if $y$ is near the boundaries of the support
$\left[\underline{y}_{dx},\overline{y}_{dx}\right]$. As \citet{Guerre2000},
we can trim off the estimated counterfactual outcomes $\widehat{\phi}_{dx}\left(Y_{i}\right)$
that lie in the boundary region $\left[\underline{y}_{dx},\underline{y}_{dx}+b_{\zeta}\right)\cup\left(\overline{y}_{dx}-b_{\zeta},\overline{y}_{dx}\right]$
by multiplying $\widehat{\zeta}_{dx}\left(\widehat{\phi}_{dx}\left(Y_{i}\right);b_{\zeta}\right)^{-1}$
in $\widehat{q}_{dx}\left(W_{j},W_{i};b_{\zeta}\right)$ by a trimming
factor $\mathbbm{1}\left(\underline{y}_{dx}+b_{\zeta}\leq\widehat{\phi}_{dx}\left(Y_{i}\right)\leq\overline{y}_{dx}-b_{\zeta}\right)$.
It can be shown that the effect of the trimming factor is asymptotically
negligible. All of our asymptotic results remain true, and the finite-sample
performances may improve when trimming is used.}
Let
\begin{equation}
\widehat{V}\left(v\mid x;b,b_{\zeta}\right)\coloneqq\widehat{p}_{x}^{-2}\left(\widehat{V}_{1}\left(v,x;b\right)+\widehat{V}_{2}\left(v,x;b,b_{\zeta}\right)\right).\label{eq:asy var}
\end{equation}
For estimating $V\left(v\mid x;h\right)$, we set $b=\widehat{h}$
in $\widehat{V}\left(v\mid x;b,b_{\zeta}\right)$, where $\widehat{h}$
satisfies Assumption \ref{assu:h_hat}. Similarly, we set the second
bandwidth $b_{\zeta}=\widehat{h}_{\zeta}$, where $\widehat{h}_{\zeta}$
is a random bandwidth that satisfies the following assumption similar
to Assumption \ref{assu:h_hat}.
\begin{assumption}
\label{assu:btw zeta}$\mathrm{Pr}\left[\left|\widehat{h}_{\zeta}/h_{\zeta}-1\right|>\varepsilon_{n}^{\zeta}\right]\leq\delta_{n}^{\zeta}$
for some deterministic bandwidth $h_{\zeta}$ and positive sequences
$\varepsilon_{n}^{\zeta},\delta_{n}^{\zeta}\downarrow0$.
\end{assumption}
Suppose $\widehat{h}_{\zeta}$ is the Silverman ROT bandwidth of the
form $\widehat{h}_{\zeta}^{\mathsf{rot}}=C_{K_{\zeta}}\cdot\widehat{\sigma}_{Y\mid X=x}\cdot n_{x}^{-1/5}$,
where $\widehat{\sigma}_{Y\mid X=x}$ is the sample analogue of $\sigma_{Y\mid X=x}\coloneqq\sqrt{\mathrm{Var}\left[Y\mid X=x\right]}$
and $C_{K_{\zeta}}$ is a constant that depends on $K_{\zeta}$. In
this case, Assumption \ref{assu:btw zeta} is satisfied with $\varepsilon_{n}^{\zeta}\wasypropto\sqrt{\mathrm{log}\left(n\right)/n}$
and $\delta_{n}^{\zeta}\wasypropto n^{-1}$.\footnote{One may use estimators $\left(\widehat{\zeta}_{0x}\left(y;b_{\zeta,0}\right),\widehat{\zeta}_{1x}\left(y;b_{\zeta,1}\right)\right)$with
different bandwidths $\left(b_{\zeta,0},b_{\zeta,1}\right)$. By easily
modifying the proofs, we get results similar to Theorems \ref{thm:variance estimator}
and \ref{thm:confidence band} under two data-dependent bandwidths
$\left(\widehat{h}_{\zeta,0},\widehat{h}_{\zeta,1}\right)$ that satisfy
the same assumption for $\widehat{h}_{\zeta}$. The ROT bandwidths
can be set as $\widehat{h}_{\zeta,d}^{\mathsf{rot}}=C_{K_{\zeta}}\cdot\widehat{\sigma}_{Y\mid X=x,D=d}\cdot n_{dx}^{-1/5}$,
where $n_{dx}\coloneqq\sum_{i=1}^{n}\mathbbm{1}\left(D_{i}=d,X_{i}=x\right)$,
for $d=0,1$, and $\widehat{\sigma}_{Y\mid X=x,D=d}$ denotes the
sample analogue of $\sigma_{Y\mid X=x,D=d}\coloneqq\sqrt{\mathrm{Var}\left[Y\mid X=x,D=d\right]}$.} Theorem \ref{thm:variance estimator} below provides a uniform convergence
rate for $\widehat{V}\left(v\mid x;\widehat{h},\widehat{h}_{\zeta}\right)$.
In Appendix \ref{sec:Proofs B}, Theorem \ref{thm:variance theorem appendix}
presents a non-asymptotic deviation bound for the uniform estimation
error of $\widehat{V}\left(v\mid x;\widehat{h},\widehat{h}_{\zeta}\right)$,
which implies the result of Theorem \ref{thm:variance estimator}.
The stronger result of Theorem \ref{thm:variance theorem appendix}
is used in the proof of Theorem \ref{thm:confidence band} below.
\begin{thm}
\label{thm:variance estimator}Suppose that Assumptions \ref{assu: DGP1}-\ref{assu:btw zeta}
hold with $P=2$, the third-order derivative functions in Assumption
\ref{assu: DGP4}(a) are Lipschitz continuous, $h\wasypropto n^{-\lambda}$
with $0<\lambda<1/4$, and $h_{\zeta}\wasypropto n^{-\lambda_{\zeta}}$
with $0<\lambda_{\zeta}<1$. Then, for any $x\in\mathscr{S}_{X}$
and compact $I_{x}\subseteq\mathrm{Int}\left(\mathscr{S}_{\Delta\mid X=x}\right)$,
\[
\left\Vert \widehat{V}\left(\cdot\mid x;\widehat{h},\widehat{h}_{\zeta}\right)-V\left(\cdot\mid x;h\right)\right\Vert _{I_{x}}=O_{p}\left(\sqrt{\frac{\mathrm{log}\left(n\right)}{nh_{\zeta}}}+h_{\zeta}^{2}+\frac{\mathrm{log}\left(n\right)}{nh^{4}}+\sqrt{\frac{\mathrm{log}\left(n\right)}{nh^{2}}}+\varepsilon_{n}h\right).
\]
\end{thm}
\begin{rembold}While the estimator $\widehat{V}_{2}\left(v,x;b,b_{\zeta}\right)$
may be negative in finite samples, its modification $
\savestack{\tmpbox}{\stretchto{
\scaleto{
\scalerel*[\widthof{\ensuremath{V}}]{\kern-.6pt\bigwedge\kern-.6pt}
{\rule[-\textheight/2]{1ex}{\textheight}}
}{\textheight}
}{0.5ex}}
\stackon[1pt]{V}{\scalebox{-1}{\tmpbox}}
_{2}\left(v,x;b,b_{\zeta}\right)$
defined below is always non-negative:
\[
\savestack{\tmpbox}{\stretchto{
\scaleto{
\scalerel*[\widthof{\ensuremath{V}}]{\kern-.6pt\bigwedge\kern-.6pt}
{\rule[-\textheight/2]{1ex}{\textheight}}
}{\textheight}
}{0.5ex}}
\stackon[1pt]{V}{\scalebox{-1}{\tmpbox}}
_{2}\left(v,x;b,b_{\zeta}\right)\coloneqq\frac{1}{n}\sum_{i=1}^{n}\frac{1}{b}\left\{ \frac{1}{n}\sum_{j=1}^{n}\frac{1}{b}K'\left(\frac{\widehat{\varDelta}_{j}-v}{b}\right)\widehat{q}_{x}\left(W_{j},W_{i};b_{\zeta}\right)\right\} ^{2}\mathbbm{1}\left(X_{i}=x\right)\widehat{p}_{x}^{-1}\left(\widehat{p}_{1x}^{-1}+\widehat{p}_{0x}^{-1}\right).
\]
One can show that the difference $\left\Vert
\savestack{\tmpbox}{\stretchto{
\scaleto{
\scalerel*[\widthof{\ensuremath{V}}]{\kern-.6pt\bigwedge\kern-.6pt}
{\rule[-\textheight/2]{1ex}{\textheight}}
}{\textheight}
}{0.5ex}}
\stackon[1pt]{V}{\scalebox{-1}{\tmpbox}}
_{2}\left(\cdot,x;\widehat{h},\widehat{h}_{\zeta}\right)-\widehat{V}_{2}\left(\cdot,x;\widehat{h},\widehat{h}_{\zeta}\right)\right\Vert _{I_{x}}$
is of a smaller order than $\left\Vert \widehat{V}_{2}\left(\cdot,x;\widehat{h},\widehat{h}_{\zeta}\right)-V_{2}\left(\cdot,x;h\right)\right\Vert _{I_{x}}$.\end{rembold}
A pointwise $1-\alpha$ asymptotic confidence interval for $f_{\Delta\mid X}\left(v\mid x\right)$
can be constructed as
\begin{equation}
\left[\widehat{f}_{\varDelta\mid X}\left(v\mid x;\widehat{h}\right)\pm z_{1-\alpha/2}\sqrt{\frac{\widehat{V}\left(v\mid x;\widehat{h},\widehat{h}_{\zeta}\right)}{n\widehat{h}}}\right],\label{eq:pointwise CI}
\end{equation}
where $z_{1-\alpha/2}$ is the $1-\alpha/2$ quantile of the standard
normal distribution, and $h$ satisfies $nh^{5}\downarrow0$. However,
if one is interested in constructing valid confidence bands for the
density function, the $z_{1-\alpha/2}$ critical value must be replaced
with a bigger one determined by the distribution of the supremum of
the estimation errors along the domain, as interpolations of the pointwise
confidence intervals (\ref{eq:pointwise CI}) over the domain are
invalid in the uniform sense. In the section below, we discuss the
construction of valid UCBs.
\subsection{Jackknife multiplier bootstrap UCB\label{subsec:Jackknife-multiplier-bootstrap}}
Let $S\left(v\mid x;b\right)$ and $Z\left(v\mid x;b,b_{\zeta}\right)$
denote the non-studentized and studentized estimation errors, respectively:
\begin{equation}
S\left(v\mid x;b\right)\coloneqq\sqrt{nb}\left(\widehat{f}_{\varDelta\mid X}\left(v\mid x;b\right)-f_{\varDelta\mid X}\left(v\mid x\right)\right)\textrm{ and }Z\left(v\mid x;b,b_{\zeta}\right)\coloneqq\frac{S\left(v\mid x;b\right)}{\sqrt{\widehat{V}\left(v\mid x;b,b_{\zeta}\right)}}.\label{eq:S Z definition}
\end{equation}
Moreover, recall the expansion of the estimation error in (\ref{eq:bias corrected linearization}).
An asymptotically valid $1-\alpha$ UCB simultaneously covers $\left\{ f_{\varDelta\mid X}\left(v\mid x\right):v\in I_{x}\right\} $
with a pre-specified asymptotic coverage $1-\alpha$. To construct
a valid UCB, one has to replace the standard normal quantile $z_{1-\alpha/2}$
in the pointwise confidence interval (\ref{eq:pointwise CI}) with
a critical value approximating the $1-\alpha$ quantile of the distribution
of $\left\Vert Z\left(\cdot\mid x;\widehat{h},\widehat{h}_{\zeta}\right)\right\Vert _{I_{x}}$.
In this section, we discuss the validity of the computationally fast
jackknife multiplier bootstrap (JMB). Appendix \ref{sec:Nonparametric-bootstrap}
provides the algorithm and theoretical results for the nonparametric
bootstrap.
We consider the problem of estimating the distribution of $\left\Vert S\left(\cdot\mid x;\widehat{h}\right)\right\Vert _{I_{x}}$
or $\left\Vert Z\left(\cdot\mid x;\widehat{h},\widehat{h}_{\zeta}\right)\right\Vert _{I_{x}}$
using the linearization (\ref{eq:bias corrected linearization}),
under the ``undersmoothing'' assumption $nh^{5}\downarrow0$ to
ensure that the bias is asymptotically negligible in comparison to
the standard deviation. The conventional bandwidth selectors that
estimate the AMISE-optimal bandwidth such as the Silverman ROT method
or cross-validation (using the pseudo ITEs in our case) violate the
undersmoothing assumption. The undersmoothing assumption requires
that the selected bandwidth should vanish at a faster rate. In practical
implementation of undersmoothing for many nonparametric econometric
techniques, a commonly used strategy is to shrink a conventional approximately
AMISE-optimal data-driven bandwidth by an ad hoc amount.
The JMB approach of \citet[ CK hereafter,]{Chen2020jackknife} approximates
the distribution of the supremum (with respect to $v$) of $n^{-1/2}\sum_{i=1}^{n}\left(\mathcal{U}_{x}^{\left[1\right]}\left(W_{i},v;h\right)-\mu_{\mathcal{U}_{x}}\left(v;h\right)\right)$
with that of the Gaussian multiplier process (e.g., \citealp{Chernozhukov2014anti})
that uses the jackknife estimator of $\mathcal{U}_{x}^{\left[1\right]}\left(W_{i},v;h\right)$.
However unlike in CK, in our case the kernel $\mathcal{U}_{x}$ involves
the unknown nonparametric objects ($\phi_{dx},\zeta_{dx},R_{d'x}$),
unknown probabilities ($p_{0x}$ and $p_{1x}$), and latent ITEs.
Therefore, we use the estimated version of $\mathcal{U}_{x}$ that
replaces the unknown objects with their nonparametric estimators:
\begin{equation}
\widehat{\mathcal{U}}_{x}\left(W_{j},W_{i},v;b,b_{\zeta}\right)\coloneqq\frac{1}{\sqrt{b}}K\left(\frac{\widehat{\varDelta}_{i}-v}{b}\right)\mathbbm{1}\left(X_{i}=x\right)+\frac{1}{b^{3/2}}K'\left(\frac{\widehat{\varDelta}_{j}-v}{b}\right)\widehat{q}_{x}\left(W_{j},W_{i};b_{\zeta}\right)\widehat{\pi}_{x}\left(Z_{i},X_{i}\right),\label{eq:estimated kernel}
\end{equation}
where $\widehat{\pi}_{x}\left(Z_{i},X_{i}\right)$ is constructed
by replacing $\left(p_{0x},p_{1x}\right)$ with $\left(\widehat{p}_{0x},\widehat{p}_{1x}\right)$
in the definition of $\pi_{x}\left(Z_{i},X_{i}\right)$ in (\ref{eq:pi_x definition}).
Let $\left(\nu_{1},...,\nu_{n}\right)$ denote i.i.d. standard normal
random variables that are drawn independently from the data. Let $\left\{ \widehat{S}_{\mathsf{jmb}}\left(\cdot\mid x;b,b_{\zeta}\right):v\in I_{x}\right\} $
be the feasible JMB process, where
\begin{eqnarray}
\widehat{\mathcal{U}}_{x}^{\left[1\right]}\left(W_{i},v;b,b_{\zeta}\right) & \coloneqq & \frac{1}{n-1}\sum_{j\neq i}\widehat{\mathcal{U}}_{x}\left(W_{j},W_{i},v;b,b_{\zeta}\right),\nonumber \\
\widehat{S}_{\mathsf{jmb}}\left(v\mid x;b,b_{\zeta}\right) & \coloneqq & \frac{1}{\sqrt{n}}\sum_{i=1}^{n}\nu_{i}\widehat{p}_{x}^{-1}\left\{ \widehat{\mathcal{U}}_{x}^{\left[1\right]}\left(W_{i},v;b,b_{\zeta}\right)-\sqrt{b}\cdot\widehat{f}_{\varDelta X}\left(v,x;b\right)\right\} .\label{eq:Z bootstrap analogue}
\end{eqnarray}
We show in Appendix \ref{sec:Proofs B} (the proof of Theorem \ref{thm:confidence band appendix})
that the distribution of $\left\Vert Z\left(\cdot\mid x;\widehat{h},\widehat{h}_{\zeta}\right)\right\Vert _{I_{x}}$
can be approximated by the conditional distribution of $\left\Vert \widehat{Z}_{\mathsf{jmb}}\left(\cdot\mid x;\widehat{h},\widehat{h}_{\zeta}\right)\right\Vert _{I_{x}}$
given the original sample $W_{1}^{n}\coloneqq\left\{ W_{1},...,W_{n}\right\} $,
where
\[
\widehat{Z}_{\mathsf{jmb}}\left(v\mid x;b,b_{\zeta}\right)\coloneqq\frac{\widehat{S}_{\mathsf{jmb}}\left(v\mid x;b,b_{\zeta}\right)}{\sqrt{\widehat{V}\left(v\mid x;b,b_{\zeta}\right)}}.
\]
Let $\mathrm{Pr}_{\mid W_{1}^{n}}\left[\cdot\right]$ and $\mathrm{E}_{\mid W_{1}^{n}}\left[\cdot\right]$
denote the conditional probability and expectation respectively given
$W_{1}^{n}$, and
\begin{equation}
z_{1-\alpha}^{\mathsf{jmb}}\coloneqq\mathrm{inf}\left\{ t\in\mathbb{R}:\mathrm{Pr}_{\mid W_{1}^{n}}\left[\left\Vert \widehat{Z}_{\mathsf{jmb}}\left(\cdot\mid x;\widehat{h},\widehat{h}_{\zeta}\right)\right\Vert _{I_{x}}\leq t\right]\geq1-\alpha\right\} \label{eq:z_pound definition}
\end{equation}
be the $1-\alpha$ quantile of the conditional distribution of $\left\Vert \widehat{Z}_{\mathsf{jmb}}\left(\cdot\mid x;\widehat{h},\widehat{h}_{\zeta}\right)\right\Vert _{I_{x}}$
given $W_{1}^{n}$. Recall that $\mathit{CB}_{\mathsf{jmb}}\left(v\mid x;b,b_{\zeta}\right)$
is defined by (\ref{eq:confidence band definition}). The JMB confidence
band is given by the family of random intervals $\left\{ \mathit{CB}_{\mathsf{jmb}}\left(v\mid x;\widehat{h},\widehat{h}_{\zeta}\right):v\in I_{x}\right\} $.
Note that one can approximate $z_{1-\alpha}^{\mathsf{jmb}}$ to any
degree of accuracy by Monte Carlo simulations, and that the width
of $\mathit{CB}_{\mathsf{jmb}}\left(v\mid x;\widehat{h},\widehat{h}_{\zeta}\right)$
varies with $\widehat{V}\left(v\mid x;\widehat{h},\widehat{h}_{\zeta}\right)$.\footnote{\label{fn:constant width}Alternatively, a constant-width UCB (e.g.,
\citealp{Cheng2019}) $\left[\widehat{f}_{\varDelta\mid X}\left(v\mid x;\widehat{h}\right)\pm s_{1-\alpha}^{\mathsf{jmb}}/\sqrt{n\widehat{h}}\right]$
is based on the critical value $s_{1-\alpha}^{\mathsf{jmb}}\coloneqq\mathrm{inf}\left\{ t\in\mathbb{R}:\mathrm{Pr}_{\mid W_{1}^{n}}\left[\left\Vert \widehat{S}_{\mathsf{jmb}}\left(\cdot\mid x;\widehat{h},\widehat{h}_{\zeta}\right)\right\Vert _{I_{x}}\leq t\right]\geq1-\alpha\right\} $
that approximates the $1-\alpha$ quantile of the distribution of
$\left\Vert S\left(\cdot\mid x;\widehat{h}\right)\right\Vert _{I_{x}}$.
However, such a UCB cannot exploit the fact that the tails of a density
function approach zero, and therefore the variable-width UCB is preferred
in this context.} The following algorithm summarizes the construction of the JMB confidence
band for the density of the ITE.
\begin{lyxalgorithm}[JMB confidence band]
\label{alg:multiplier}\textbf{Step 1}: Compute the pseudo ITEs using
(\ref{eq:phi_hat definition})-(\ref{eq:pseudo ITE definition}).
\textbf{Step 2}: Select the covariates' value $x$, the number of
grid points $G$, and a grid $I_{x}^{G}\coloneqq\left\{ v_{1}<v_{2}<\ldots<v_{G}\right\} $
over which the density is estimated. \textbf{Step 3}: Select a kernel
$K$, use the ROT bandwidth $\widehat{h}^{\mathsf{rot}}$ with undersmoothing
(e.g., $\widehat{h}=\widehat{h}^{\mathsf{rot}}n^{-1/40}$), and for
all $v\in I_{x}^{G}$ compute the kernel density estimator $\widehat{f}_{\varDelta\mid X}\left(v\mid x;\widehat{h}\right)$
using (\ref{eq:f_hat definition}) and the pseudo ITEs from Step 1.
\textbf{Step 4}: Select a kernel $K_{\zeta}$, use the ROT bandwidth
$\widehat{h}_{\zeta}^{\mathsf{rot}}$, and for all $v\in I_{x}^{G}$
compute the variance estimator $\widehat{V}\left(v\mid x;\widehat{h},\widehat{h}_{\zeta}\right)$
using (\ref{eq:V_hat_dagger definition}), (\ref{eq:V_hat_ddager definition}),
and (\ref{eq:asy var}). \textbf{Step 5}: Select the number of bootstrap
repetitions $B$, for $r=1,...,B$ generate i.i.d. standard normal
random variables $\left(\nu_{1}^{\left(r\right)},...,\nu_{n}^{\left(r\right)}\right)$,
and for all $v\in I_{x}^{G}$ compute $\widehat{Z}_{\mathsf{jmb}}^{\left(r\right)}\left(v\mid x;\widehat{h},\widehat{h}_{\zeta}\right)=\widehat{S}_{\mathsf{jmb}}^{\left(r\right)}\left(v\mid x;\widehat{h},\widehat{h}_{\zeta}\right)/\sqrt{\widehat{V}\left(v\mid x;\widehat{h},\widehat{h}_{\zeta}\right)}$
using (\ref{eq:Z bootstrap analogue}) with $\left(\nu_{1},...,\nu_{n}\right)$
replaced by $\left(\nu_{1}^{\left(r\right)},...,\nu_{n}^{\left(r\right)}\right)$.
\textbf{Step 6}: Select the coverage level $1-\alpha$ and compute
the critical value $z_{1-\alpha}^{\mathsf{jmb}}$:
\begin{equation}
z_{1-\alpha}^{\mathsf{jmb}}=\mathrm{inf}\left\{ t\in\mathbb{R}:\frac{1}{B}\sum_{r=1}^{B}\mathbbm{1}\left(\underset{v\in I_{x}^{G}}{\mathrm{max}}\left|\widehat{Z}_{\mathsf{jmb}}^{\left(r\right)}\left(v\mid x;\widehat{h},\widehat{h}_{\zeta}\right)\right|\leq t\right)\geq1-\alpha\right\} .\label{eq:JMB critical value}
\end{equation}
\textbf{Step 7}: Compute the JMB confidence band $\mathit{CB}_{\mathsf{jmb}}$
using (\ref{eq:confidence band definition}) over $v\in I_{x}^{G}$.
\end{lyxalgorithm}
Theorem \ref{thm:confidence band} below shows that the proposed JMB
confidence band is asymptotically valid and its coverage error decays
at a polynomial rate. The result rules out coverage probability errors
with logarithmic decay rates (see, e.g., \citealp{Chernozhukov2014anti}
for discussion).
\begin{thm}
\label{thm:confidence band}Suppose that Assumptions \ref{assu: DGP1}-\ref{assu:btw zeta}
hold with $P=2$, the third-order derivative functions in Assumption
\ref{assu: DGP4}(a) are Lipschitz continuous, $h\wasypropto n^{-\lambda}$
with $1/5<\lambda<1/4$, $h_{\zeta}\wasypropto n^{-\lambda_{\zeta}}$
with $1/8<\lambda_{\zeta}<1/2$, $\varepsilon_{n}=O\left(n^{-\varrho_{\varepsilon}}\right)$,
$\delta_{n}=O\left(n^{-\varrho_{\delta}}\right)$, and $\delta_{n}^{\zeta}=O\left(n^{-\varrho_{\delta}^{\zeta}}\right)$
for some $\varrho_{\varepsilon},\varrho_{\delta},\varrho_{\delta}^{\zeta}>0$.
Then, for all $x\in\mathscr{S}_{X}$ and any compact $I_{x}\subseteq\mathrm{Int}\left(\mathscr{S}_{\Delta\mid X=x}\right)$,
\[
\mathrm{Pr}\left[f_{\varDelta\mid X}\left(v\mid x\right)\in\mathit{CB}_{\mathsf{jmb}}\left(v\mid x;\widehat{h},\widehat{h}_{\zeta}\right),\,\textrm{for all }v\in I_{x}\right]=\left(1-\alpha\right)+O\left(n^{-\varrho}\right)
\]
for some $\varrho>0$.
\end{thm}
\begin{rembold}In the proof of Theorem \ref{thm:confidence band},
we explicitly derive an estimate of the coverage probability error,
which is presented in Theorem \ref{thm:confidence band appendix}
in Appendix \ref{sec:Proofs B}. We show that replacing the deterministic
bandwidths $\left(h,h_{\zeta}\right)$ with estimators satisfying
Assumptions \ref{assu:h_hat} and \ref{assu:btw zeta} incurs an additional
$O\left(\mathrm{log}\left(n\right)\varepsilon_{n}+\delta_{n}+\delta_{n}^{\zeta}\right)$
coverage probability error term. This result shows how the noise in
the estimated bandwidth translates into an error in the coverage probability.
The JMB confidence band achieves a potentially higher coverage accuracy
if $\left(\widehat{h},\widehat{h}_{\zeta}\right)$ converge to their
population counterparts at a fast rate, which is the case for the
simple and feasible Silverman ROT approach we recommended.\end{rembold}
\begin{rembold}\label{Rmk: truncation} Since density functions are
non-negative, the lower bound of the UCB can be truncated to zero
to avoid negative values. Thus, the lower bound of $\mathit{CB}_{\mathsf{jmb}}\left(v\mid x;\widehat{h},\widehat{h}_{\zeta}\right)$
can be replaced with $\mathrm{max}\left\{ 0,\widehat{f}_{\varDelta\mid X}\left(v\mid x;\widehat{h}\right)-z_{1-\alpha}^{\mathsf{jmb}}\sqrt{\widehat{V}\left(v\mid x;\widehat{h},\widehat{h}_{\zeta}\right)/\left(n\widehat{h}\right)}\right\} $
without affecting the coverage properties of the UCB. Hence, the result
of Theorem \ref{thm:confidence band} continues to hold with the modified
lower bound.\footnote{We thank the associate editor for suggesting the modification to us.}\end{rembold}
Our main focus is on the multiplier bootstrap approach, as it is computationally
fast even with large sample sizes. The more commonly used nonparametric
bootstrap would require re-calculation of the bootstrap versions of
the estimated ITEs at every bootstrap repetition, which can be computationally
burdensome. However, the constant-width version of the nonparametric
bootstrap confidence band has the advantage of fewer tuning parameters
as it does not require estimation of $\zeta_{dx}$ and therefore does
not need the second bandwidth $b_{\zeta}$. The validity of the nonparametric
bootstrap approach is discussed in Appendix \ref{sec:Nonparametric-bootstrap}.
\subsection{Bias-corrected JMB UCB\label{subsec:Bias-corrected-JMB}}
In this section, we discuss the bias correction approach to inference
that can accommodate conventional bandwidth selectors such as ROT
bandwidths that decay at the $n^{-1/5}$ rate as in \citet{Calonico2014}.
We assume that the third-order derivatives in Assumption \ref{assu: DGP4}
with $P=2$ are Lipschitz continuous. Theorem \ref{thm:asymptotic normality}
implies that in large samples and when $P=2$, $\widehat{f}_{\varDelta\mid X}\left(v\mid x;h\right)-f_{\varDelta\mid X}^{(2)}\left(v\mid x\right)\mu_{K,2}h^{2}$
is approximately distributed as $\mathrm{N}\left(f_{\varDelta\mid X}\left(v\mid x\right),\mathscr{V}\left(v\mid x\right)/\left(nh\right)\right)$.
We use an estimator of the density derivative $f_{\varDelta\mid X}^{(2)}\left(v\mid x\right)$
to remove the bias. Let
\[
\widehat{f}_{\varDelta\mid X}^{(2)}\left(v\mid x;b_{\mathsf{b}}\right)\coloneqq\frac{\sum_{i=1}^{n}\frac{1}{b_{\mathsf{b}}^{3}}K_{\mathsf{b}}^{(2)}\left(\frac{\widehat{\varDelta}_{i}-v}{b_{\mathsf{b}}}\right)\mathbbm{1}\left(X_{i}=x\right)}{\sum_{i=1}^{n}\mathbbm{1}\left(X_{i}=x\right)}
\]
denote the kernel estimator of $f_{\varDelta\mid X}^{(2)}\left(v\mid x\right)$
using a bandwidth $b_{\mathsf{b}}$ and a second-order kernel $K_{\mathsf{b}}$
($K_{\mathsf{b}}^{(2)}$ denotes the second derivative of $K_{\mathsf{b}}$).
We assume that $K_{\mathsf{b}}^{(2)}$ satisfies Assumption \ref{assu:kernel}(a).
The bias-corrected estimator of $f_{\Delta\mid X}\left(v\mid x\right)$
is given by
\begin{equation}
\widehat{f}_{\varDelta\mid X}^{\mathsf{bc}}\left(v\mid x;b,b_{\mathsf{b}}\right)\coloneqq\widehat{f}_{\varDelta\mid X}\left(v\mid x;b\right)-\widehat{f}_{\varDelta\mid X}^{(2)}\left(v\mid x;b_{\mathsf{b}}\right)\mu_{K,2}b^{2}=\frac{\sum_{i=1}^{n}\frac{1}{b}M\left(\frac{\widehat{\varDelta}_{i}-v}{b};b,b_{\mathsf{b}}\right)\mathbbm{1}\left(X_{i}=x\right)}{\sum_{i=1}^{n}\mathbbm{1}\left(X_{i}=x\right)},\label{eq:debiased estimator definition}
\end{equation}
where $M\left(u;b,b_{\mathsf{b}}\right)\coloneqq K\left(u\right)-\left(b/b_{\mathsf{b}}\right)^{3}\mu_{K,2}K_{\mathsf{b}}^{(2)}\left(\left(b/b_{\mathsf{b}}\right)u\right)$.
Typical choices for the bandwidth $b_{\mathsf{b}}$ used for bias
correction include the estimation-optimal bandwidth for the second-order
density derivative (e.g., \citealp{Xu2017}) or the same bandwidth
as that used for estimating the density \citep[e.g.,][]{Cheng2019}.
We assume that $b_{\mathsf{b}}$ is some random bandwidth $\widehat{h}_{\mathsf{b}}$
that consistently estimates some target bandwidth $h_{\mathsf{b}}$
in the sense that $\widehat{h}_{\mathsf{b}}/h_{\mathsf{b}}\rightarrow_{p}1$.
We assume that $h/h_{\mathsf{b}}\rightarrow\varsigma\in\left[0,\infty\right)$
as in \citet{Calonico2014}. In practice, we can take $\widehat{h}_{\mathsf{b}}$
to be the feasible Silverman ROT bandwidth for the second-order density
derivative that uses the pseudo ITEs and aims to approximate the AMISE-optimal
bandwidth (in this case, $h_{\mathsf{b}}\wasypropto n^{-1/9}$). Alternatively,
one can set $\widehat{h}_{\mathsf{b}}=\widehat{h}$.
One can show that the bias of $\widehat{f}_{\varDelta\mid X}^{\mathsf{bc}}\left(v\mid x;h,h_{\mathsf{b}}\right)$
is of the order $O\left(h^{2}h_{\mathsf{b}}\right)$, and its standard
deviation is of the order $\left(nh\right)^{-1/2}$.\footnote{Under Assumptions \ref{assu: DGP1}-\ref{assu:kernel} with $P=2$,
the bias part of the bias-corrected estimator in (\ref{eq:debiased estimator definition})
is $o\left(h^{2}\right)$. Lipschitz continuity ensures that the bias
part is $O\left(h^{2}h_{\mathsf{b}}\right)$, if $h=O\left(h_{\mathsf{b}}\right)$.} Hence, when the bandwidth $h$ is chosen AMISE-optimally so that
$h\wasypropto n^{-1/5}$, the bias of the bias-corrected estimator
$\widehat{f}_{\varDelta\mid X}^{\mathsf{bc}}\left(v\mid x;h,h_{\mathsf{b}}\right)$
is of smaller order than its standard deviation. Robust standard errors
(\citealp{Calonico2014}) have to take into account the additional
stochastic variability coming from the bias correction, i.e., estimation
of $f_{\varDelta\mid X}^{(2)}\left(v\mid x\right)$.
\sloppy Let $\widehat{V}^{\mathsf{bc}}\left(v\mid x;b,b_{\zeta},b_{\mathsf{b}}\right)$
and $\widehat{S}_{\mathsf{jmb}}^{\mathsf{bc}}\left(v\mid x;b,b_{\zeta},b_{\mathsf{b}}\right)$
be defined similarly to $\widehat{V}\left(v\mid x;b,b_{\zeta}\right)$
and $\widehat{S}_{\mathsf{jmb}}\left(v\mid x;b,b_{\zeta}\right)$
respectively with $K\left(\cdot\right)$ replaced by the bias-correcting
kernel $M\left(\cdot;b,b_{\mathsf{b}}\right)$. Define $\widehat{Z}_{\mathsf{jmb}}^{\mathsf{bc}}\left(v\mid x;b,b_{\zeta},b_{\mathsf{b}}\right)\coloneqq\widehat{S}_{\mathsf{jmb}}^{\mathsf{bc}}\left(v\mid x;b,b_{\zeta},b_{\mathsf{b}}\right)/\sqrt{\widehat{V}^{\mathsf{bc}}\left(v\mid x;b,b_{\zeta},b_{\mathsf{b}}\right)}$,
and let $z_{1-\alpha}^{\mathsf{jmb},\mathsf{bc}}$ be as in (\ref{eq:z_pound definition})
with $\left\Vert \widehat{Z}_{\mathsf{jmb}}\left(\cdot\mid x;\widehat{h},\widehat{h}_{\zeta}\right)\right\Vert _{I_{x}}$
replaced by $\left\Vert \widehat{Z}_{\mathsf{jmb}}^{\mathsf{bc}}\left(\cdot\mid x;\widehat{h},\widehat{h}_{\zeta},\widehat{h}_{\mathsf{b}}\right)\right\Vert _{I_{x}}$.
The bias-corrected JMB UCB is given by
\begin{equation}
\mathit{CB}_{\mathsf{jmb}}^{\mathsf{bc}}\left(v\mid x;\widehat{h},\widehat{h}_{\zeta},\widehat{h}_{\mathsf{b}}\right)\coloneqq\left[\widehat{f}_{\varDelta\mid X}^{\mathsf{bc}}\left(v\mid x;\widehat{h},\widehat{h}_{\mathsf{b}}\right)\pm z_{1-\alpha}^{\mathsf{jmb},\mathsf{bc}}\sqrt{\frac{\widehat{V}^{\mathsf{bc}}\left(v\mid x;\widehat{h},\widehat{h}_{\zeta},\widehat{h}_{\mathsf{b}}\right)}{n\widehat{h}}}\right].\label{eq:bias corrected CB}
\end{equation}
The bias-corrected UCB can be computed by replacing $\widehat{f}_{\varDelta\mid X}\left(v\mid x;\widehat{h}\right)$,
$\widehat{V}\left(v\mid x;\widehat{h},\widehat{h}_{\zeta}\right)$,
and $\widehat{Z}_{\mathsf{jmb}}\left(v\mid x;\widehat{h},\widehat{h}_{\zeta}\right)$
in Algorithm \ref{alg:multiplier} with their respective bias-corrected
versions $\widehat{f}_{\varDelta\mid X}^{\mathsf{bc}}\left(v\mid x;\widehat{h},\widehat{h}_{\mathsf{b}}\right)$,
$\widehat{V}^{\mathsf{bc}}\left(v\mid x;\widehat{h},\widehat{h}_{\zeta},\widehat{h}_{\mathsf{b}}\right)$,
and $\widehat{Z}_{\mathsf{jmb}}^{\mathsf{bc}}\left(v\mid x;\widehat{h},\widehat{h}_{\zeta},\widehat{h}_{\mathsf{b}}\right)$.
We show that the conclusion of Theorem \ref{thm:confidence band}
(asymptotic validity with a polynomial rate) holds for $\mathit{CB}_{\mathsf{jmb}}^{\mathsf{bc}}$
under the assumption that $h_{\mathsf{b}}\wasypropto n^{-\lambda_{\mathsf{b}}}$,
$\lambda_{\mathsf{b}}\leq\lambda$, $\left(1-2\lambda_{\mathsf{b}}\right)/5<\lambda<1/4$,
and $\mathrm{Pr}\left[\left|\widehat{h}_{\mathsf{b}}/h_{\mathsf{b}}-1\right|>\varepsilon_{n}^{\mathsf{b}}\right]\leq\delta_{n}^{\mathsf{b}}$,
where $\varepsilon_{n}^{\mathsf{b}}=O\left(n^{-\varrho_{\varepsilon}^{\mathsf{b}}}\right)$
and $\delta_{n}^{\mathsf{b}}=O\left(n^{-\varrho_{\delta}^{\mathsf{b}}}\right)$
for some $\varrho_{\varepsilon}^{\mathsf{b}},\varrho_{\delta}^{\mathsf{b}}>0$.
Note that now the ROT bandwidth without undersmoothing satisfies the
rate requirement and thus can be used in Step 3. In the online supplement,
we sketch how to show this result.
\subsection{Conditioning on sub-vectors of the covariates\label{subsec:Conditioning-on-sub-vectors}}
In applications, researchers are often interested in the unconditional
PDF of the ITE, or the conditional PDF of the ITE after conditioning
only on some of the covariates. See, e.g., the application in Section
\ref{sec:Empirical-illustrations}. This section discusses how our
results can be applied in such cases.
\sloppy Partition the vector of covariates as $X=\left(X_{1}^{\top},X_{2}^{\top}\right)^{\top}$
(similarly, $X_{i}=\left(X_{1,i}^{\top},X_{2,i}^{\top}\right)^{\top}$).
Let $f_{\varDelta\mid X_{1}}\left(\cdot\mid x_{1}\right)$ denote
the conditional density of the ITE for some fixed $x_{1}\in\mathscr{S}_{X_{1}}$.
The estimator of the conditional density is given by
\begin{equation}
\widehat{f}_{\varDelta\mid X_{1}}\left(v\mid x_{1};b\right)\coloneqq\frac{\sum_{i=1}^{n}\frac{1}{b}K\left(\frac{\widehat{\varDelta}_{i}-v}{b}\right)\mathbbm{1}\left(X_{1,i}=x_{1}\right)}{\sum_{i=1}^{n}\mathbbm{1}\left(X_{1,i}=x_{1}\right)}.\label{eq:f_hat x_dag definition}
\end{equation}
Let $\widehat{q}_{X_{i}}\left(W_{i},W_{j};b_{\zeta}\right)$ be defined
similarly to $\widehat{q}_{x}\left(W_{i},W_{j};b_{\zeta}\right)$
by the same formula with $x$ replaced by $X_{i}$. Note that in the
definition of $\widehat{q}_{X_{i}}\left(W_{i},W_{j};b_{\zeta}\right)$,
$b_{\zeta}$ may depend on $X_{i}$. Next, let $\widehat{q}_{x_{1}}\left(W_{i},W_{j};b_{\zeta}\right)\coloneqq\mathbbm{1}\left(X_{1,i}=x_{1}\right)\widehat{q}_{X_{i}}\left(W_{i},W_{j};b_{\zeta}\right)$,
$\widehat{p}_{x_{1}}\coloneqq n^{-1}\sum_{i=1}^{n}\mathbbm{1}\left(X_{1,i}=x_{1}\right)$,
$\widehat{f}_{\varDelta X_{1}}\left(v,x_{1};b\right)\coloneqq\widehat{f}_{\varDelta\mid X_{1}}\left(v\mid x_{1};b\right)\widehat{p}_{x_{1}}$,
and
\begin{eqnarray*}
\widehat{\mathcal{U}}_{x_{1}}\left(W_{j},W_{i},v;b,b_{\zeta}\right) & \coloneqq & \frac{1}{\sqrt{b}}K\left(\frac{\widehat{\varDelta}_{i}-v}{b}\right)\mathbbm{1}\left(X_{1,i}=x_{1}\right)\\
& & +\frac{1}{b^{3/2}}K'\left(\frac{\widehat{\varDelta}_{j}-v}{b}\right)\widehat{q}_{x_{1}}\left(W_{j},W_{i};b_{\zeta}\right)\widehat{\pi}_{X_{j}}\left(Z_{i},X_{i}\right),\\
\widehat{\mathcal{U}}_{x_{1}}^{\left[1\right]}\left(W_{i},v;b,b_{\zeta}\right) & \coloneqq & \frac{1}{n-1}\sum_{j\neq i}\widehat{\mathcal{U}}_{x_{1}}\left(W_{j},W_{i},v;b,b_{\zeta}\right),\\
\widehat{S}_{\mathsf{jmb}}\left(v\mid x_{1};b,b_{\zeta}\right) & \coloneqq & \frac{1}{\sqrt{n}}\sum_{i=1}^{n}\nu_{i}\widehat{p}_{x_{1}}^{-1}\left\{ \widehat{\mathcal{U}}_{x_{1}}^{\left[1\right]}\left(W_{i},v;b,b_{\zeta}\right)-\sqrt{b}\cdot\widehat{f}_{\varDelta X_{1}}\left(v,x_{1};b\right)\right\} .
\end{eqnarray*}
Define also
\begin{eqnarray}
\widehat{V}_{1}\left(v,x_{1};b\right) & \coloneqq & \frac{1}{n}\sum_{i=1}^{n}\frac{1}{b}K\left(\frac{\widehat{\varDelta}_{i}-v}{b}\right)^{2}\mathbbm{1}\left(X_{1,i}=x_{1}\right)-b\cdot\widehat{f}_{\varDelta X_{1}}\left(v,x_{1};b\right)^{2},\nonumber \\
\widehat{V}_{2}\left(v,x_{1};b,b_{\zeta}\right) & \coloneqq & \frac{1}{n\left(n-1\right)\left(n-2\right)}\sum_{i=1}^{n}\sum_{j\neq i}\sum_{k\neq i,\,k\neq j}\frac{1}{b^{3}}K'\left(\frac{\widehat{\varDelta}_{j}-v}{b}\right)\widehat{q}_{x_{1}}\left(W_{j},W_{i};b_{\zeta}\right)\nonumber \\
& & \times K'\left(\frac{\widehat{\varDelta}_{k}-v}{b}\right)\widehat{q}_{x_{1}}\left(W_{k},W_{i};b_{\zeta}\right)\left(\frac{1}{\widehat{p}_{0X_{i}}}+\frac{1}{\widehat{p}_{1X_{i}}}\right)\frac{\mathbbm{1}\left(X_{i}=X_{j}=X_{k}\right)}{\widehat{p}_{X_{i}}},\label{eq:V_hat x_dag definition}
\end{eqnarray}
$\widehat{V}\left(v\mid x_{1};b,b_{\zeta}\right)\coloneqq\widehat{p}_{x_{1}}^{-2}\left(\widehat{V}_{1}\left(v,x_{1};b\right)+\widehat{V}_{2}\left(v,x_{1};b,b_{\zeta}\right)\right)$.
For bandwidth selection, we can set $b$ to be the pseudo-ITE-based
version of the Silverman ROT bandwidth computed using the subsample
corresponding to $X_{1,i}=x_{1}$, and set $b_{\zeta}$ to be the
same ROT bandwidth as in Algorithm \ref{alg:multiplier}. A JMB UCB
for $\left\{ f_{\varDelta\mid X_{1}}\left(v\mid x_{1}\right):v\in I_{x_{1}}\right\} $,
where $I_{x_{1}}$ denotes an inner closed sub-interval of $\mathscr{S}_{\varDelta\mid X_{1}=x_{1}}$,
can be constructed similarly to $\mathit{CB}_{\mathsf{jmb}}$ in (\ref{eq:confidence band definition})
by adapting Algorithm \ref{alg:multiplier}.
For the unconditional PDF $f_{\varDelta}$ of the ITE, the estimator
is given by $\widehat{f}_{\varDelta}\left(v;b\right)\coloneqq\left(nb\right)^{-1}\sum_{i=1}^{n}K\left(\left(\widehat{\varDelta}_{i}-v\right)/b\right)$.
Let
\begin{eqnarray*}
\widehat{\mathcal{U}}\left(W_{j},W_{i},v;b,b_{\zeta}\right) & \coloneqq & \frac{1}{\sqrt{b}}K\left(\frac{\widehat{\varDelta}_{i}-v}{b}\right)+\frac{1}{b^{3/2}}K'\left(\frac{\widehat{\varDelta}_{j}-v}{b}\right)\widehat{q}_{X_{j}}\left(W_{j},W_{i};b_{\zeta}\right)\widehat{\pi}_{X_{j}}\left(Z_{i},X_{i}\right),\\
\widehat{\mathcal{U}}^{\left[1\right]}\left(W_{i},v;b,b_{\zeta}\right) & \coloneqq & \frac{1}{n-1}\sum_{j\neq i}\widehat{\mathcal{U}}\left(W_{j},W_{i},v;b,b_{\zeta}\right),\\
\widehat{S}_{\mathsf{jmb}}\left(v;b,b_{\zeta}\right) & \coloneqq & \frac{1}{\sqrt{n}}\sum_{i=1}^{n}\nu_{i}\left\{ \widehat{\mathcal{U}}^{\left[1\right]}\left(W_{i},v;b,b_{\zeta}\right)-\sqrt{b}\cdot\widehat{f}_{\varDelta}\left(v;b\right)\right\} .
\end{eqnarray*}
Let $I$ denote an inner closed sub-interval of $\mathscr{S}_{\varDelta}$,
$\widehat{V}_{1}\left(v;b\right)\coloneqq\left(nb\right)^{-1}\sum_{i=1}^{n}K\left(\left(\widehat{\varDelta}_{i}-v\right)/b\right)^{2}-b\cdot\widehat{f}_{\varDelta}\left(v;b\right)^{2}$,
and define $\widehat{V}_{2}\left(v;b,b_{\zeta}\right)$ as in (\ref{eq:V_hat x_dag definition})
but with $\widehat{q}_{x_{1}}\left(W_{j},W_{i};b_{\zeta}\right)$
and $\widehat{q}_{x_{1}}\left(W_{k},W_{i};b_{\zeta}\right)$ replaced
by $\widehat{q}_{X_{j}}\left(W_{j},W_{i};b_{\zeta}\right)$ and $\widehat{q}_{X_{k}}\left(W_{k},W_{i};b_{\zeta}\right)$
respectively. Further, define $\widehat{V}\left(v;b,b_{\zeta}\right)\coloneqq\widehat{V}_{1}\left(v;b\right)+\widehat{V}_{2}\left(v;b,b_{\zeta}\right)$.
A JMB UCB can be constructed similarly to the above. We set $b$ to
be the pseudo-ITE-based version of the Silverman ROT bandwidth computed
using the entire sample, and set $b_{\zeta}$ to be the same ROT bandwidth
as in Algorithm \ref{alg:multiplier}.
\section{Monte Carlo experiments\label{sec:Monte-Carlo-experiments}}
This section evaluates the finite-sample performance of the UCBs proposed
in Section \ref{sec:Robust-inference} for the density $f_{\varDelta}\left(v\right)$
of the ITE. We consider the following experiment design based on FVX.
The outcome and treatment status variables are generated according
to
\begin{eqnarray}
Y & = & \left(\epsilon+1\right)^{2+D}\nonumber \\
D & = & \mathbbm{1}\left(\gamma_{0}+\gamma_{1}\cdot Z+\eta\geq0\right),\label{eq:model simulation}
\end{eqnarray}
where $\left(\epsilon,\eta\right)=\left(\varPhi\left(U\right),\varPhi\left(V\right)\right)$,
$\left(U,V\right)$ has a mean-zero bivariate normal distribution
with $\mathrm{Var}\left[U\right]=\mathrm{Var}\left[V\right]=1$ and
$\mathrm{Cov}\left[U,V\right]=0.3$, and $\varPhi$ is the standard
normal CDF. The instrument is generated according to $Z=\mathbbm{1}\left(N>0\right)$,
where $N$ is a standard normal random variable independent of $\left(\epsilon,\nu\right)$.
In this design, the ITE satisfies $\varDelta=\epsilon\left(\epsilon+1\right)^{2}$,
where $\epsilon$ is uniformly distributed on $\left[0,1\right]$
and $\varDelta$ is supported on $\left[0,4\right]$. We consider
two sets of values for $(\gamma_{0},\gamma_{1}):$ $(-0.5,0.5)$ and
$(-0.4,0.6).$ We use the triweight kernel for $K$, $K_{\zeta}$
and $K_{\mathsf{b}}$, and the Silverman ROT bandwidths for $\widehat{h}$,
$\widehat{h}_{\mathsf{b}}$ and $\widehat{h}_{\zeta}.$\footnote{To be specific, we take $\widehat{h}=3.15\cdot\widehat{\sigma}_{\varDelta}\cdot n^{-1/5}$,
$\widehat{h}_{\mathsf{b}}=2.7\cdot\widehat{\sigma}_{\varDelta}\cdot n^{-1/9}$.
When it comes to $\widehat{h}_{\zeta}$, we distinguish between the
treated and control subsamples, i.e., $\widehat{h}_{\zeta,d}=3.15\cdot\widehat{\sigma}_{Y\mid D=d}\cdot n_{d}^{-1/5}$
with $n_{d}\coloneqq\sum_{i=1}^{n}\mathbbm{1}\left(D_{i}=d\right)$.} The number of Monte Carlo replications is set to $1,000$.
\begin{table}[t]
\caption{Simultaneous coverage rates for UCBs and the interpolated bootstrap
percentile CIs}
\label{table: uniform_CP}\resizebox{\textwidth}{!}{
\begin{tabular}{cccccccccccccc}
\toprule
\multicolumn{1}{c}{} & \multicolumn{1}{c}{} & \multicolumn{6}{c}{$\gamma_{0}=-0.5,\gamma_{1}=0.5$} & \multicolumn{6}{c}{$\gamma_{0}=-0.4,\gamma_{1}=0.6$}\tabularnewline
\cmidrule{3-14} \cmidrule{4-14} \cmidrule{5-14} \cmidrule{6-14} \cmidrule{7-14} \cmidrule{8-14} \cmidrule{9-14} \cmidrule{10-14} \cmidrule{11-14} \cmidrule{12-14} \cmidrule{13-14} \cmidrule{14-14}
\multicolumn{1}{c}{} & \multicolumn{1}{c}{} & \multicolumn{3}{c}{$v\in\left[0.5,3.5\right]$} & \multicolumn{3}{c}{$v\in\left[0.8,3.2\right]$} & \multicolumn{3}{c}{$v\in\left[0.5,3.5\right]$} & \multicolumn{3}{c}{$v\in\left[0.8,3.2\right]$}\tabularnewline
\cmidrule{3-14} \cmidrule{4-14} \cmidrule{5-14} \cmidrule{6-14} \cmidrule{7-14} \cmidrule{8-14} \cmidrule{9-14} \cmidrule{10-14} \cmidrule{11-14} \cmidrule{12-14} \cmidrule{13-14} \cmidrule{14-14}
\multicolumn{1}{c}{$n$} & \multicolumn{1}{c}{Methods} & \multicolumn{1}{c}{0.90} & \multicolumn{1}{c}{0.95} & \multicolumn{1}{c}{0.99} & \multicolumn{1}{c}{0.90} & \multicolumn{1}{c}{0.95} & \multicolumn{1}{c}{0.99} & \multicolumn{1}{c}{0.90} & \multicolumn{1}{c}{0.95} & \multicolumn{1}{c}{0.99} & \multicolumn{1}{c}{0.90} & \multicolumn{1}{c}{0.95} & \multicolumn{1}{c}{0.99}\tabularnewline
\midrule
2000 & Bias-corrected JMB & 0.823 & 0.891 & 0.953 & 0.836 & 0.893 & 0.957 & 0.817 & 0.887 & 0.943 & 0.820 & 0.883 & 0.946\tabularnewline
& Bias-corrected NPB & 0.852 & 0.921 & 0.980 & 0.869 & 0.930 & 0.980 & 0.851 & 0.930 & 0.986 & 0.861 & 0.926 & 0.987\tabularnewline
& Interpolated CIs & 0.507 & 0.715 & 0.927 & 0.605 & 0.795 & 0.958 & 0.579 & 0.769 & 0.949 & 0.663 & 0.829 & 0.963\tabularnewline
\midrule
4000 & Bias-corrected JMB & 0.831 & 0.891 & 0.956 & 0.847 & 0.905 & 0.953 & 0.819 & 0.885 & 0.946 & 0.835 & 0.891 & 0.946\tabularnewline
& Bias-corrected NPB & 0.868 & 0.932 & 0.985 & 0.876 & 0.936 & 0.987 & 0.875 & 0.936 & 0.989 & 0.870 & 0.932 & 0.989\tabularnewline
& Interpolated CIs & 0.405 & 0.628 & 0.914 & 0.533 & 0.729 & 0.954 & 0.476 & 0.703 & 0.934 & 0.581 & 0.782 & 0.957\tabularnewline
\midrule
6000 & Bias-corrected JMB & 0.837 & 0.902 & 0.966 & 0.864 & 0.909 & 0.966 & 0.844 & 0.900 & 0.960 & 0.849 & 0.904 & 0.958\tabularnewline
& Bias-corrected NPB & 0.860 & 0.944 & 0.994 & 0.894 & 0.954 & 0.996 & 0.876 & 0.946 & 0.991 & 0.884 & 0.944 & 0.990\tabularnewline
& Interpolated CIs & 0.364 & 0.580 & 0.885 & 0.548 & 0.744 & 0.956 & 0.431 & 0.655 & 0.930 & 0.586 & 0.773 & 0.967\tabularnewline
\bottomrule
\end{tabular}}
\end{table}
Tables \ref{table: uniform_CP} reports the simultaneous coverage
rates of two types of UCBs: the bias-corrected JMB UCB defined in
(\ref{eq:bias corrected CB}) and the bias-corrected nonparametric
bootstrap (NPB) UCB described in Appendix \ref{sec:Nonparametric-bootstrap}.
FVX interpolates pointwise (nonparametric) bootstrap percentile confidence
intervals (CIs) to construct a confidence band for the density of
the ITE. It follows from our results that such intervals are valid
in the pointwise sense.\footnote{Validity of pointwise bootstrap percentile confidence intervals follows
from (\ref{eq:bias corrected linearization}), Lemma \ref{lem:S_star linearization}
in the supplement and standard arguments (see the proof of \citealp[Theorem 4.2]{Ma2019}).} We also report the coverage probability of the confidence band constructed
by interpolating the bootstrap percentile pointwise confidence intervals.
The nominal coverage rates are $0.90$, $0.95$, and $0.99$. We consider
two ranges of $v$: a longer interval $I=[0.5,3.5]$ and a shorter
one $I=[0.8,3.2]$. We use grid search to solve the one-dimensional
optimization problems in (\ref{eq:phi_hat definition}) when estimating
the pseudo ITEs and calculating the supremum of the bootstrap process.
The number of bootstrap replications is set to $5,000$. Table \ref{table: uniform_length}
reports the average widths of the bias-corrected JMB and NPB UCBs
relative to the interpolated pointwise CIs.\footnote{The average width of bias-corrected JMB (or NPB) UCB is computed by
first averaging the widths of the confidence band over all grid points
in the given range $I$ and then averaging over all simulation replications.
The reported number is the ratio of the average width of the UCBs
to that of the interpolated pointwise CIs.}
\begin{table}[t]
\caption{Average widths of the $95\%$ UCBs relative to the interpolated bootstrap
percentile CIs}
\label{table: uniform_length} \centering{}
\begin{tabular}{cccccc}
\toprule
\multicolumn{1}{c}{} & \multicolumn{1}{c}{} & \multicolumn{2}{c}{$\gamma_{0}=-0.5,\gamma_{1}=0.5$} & \multicolumn{2}{c}{$\gamma_{0}=-0.4,\gamma_{1}=0.6$}\tabularnewline
\cmidrule{3-6} \cmidrule{4-6} \cmidrule{5-6} \cmidrule{6-6}
\multicolumn{1}{c}{} & \multicolumn{1}{c}{Methods} & \multicolumn{1}{c}{$v\in\left[0.5,3.5\right]$} & \multicolumn{1}{c}{$v\in\left[0.8,3.2\right]$} & \multicolumn{1}{c}{$v\in\left[0.5,3.5\right]$} & \multicolumn{1}{c}{$v\in\left[0.8,3.2\right]$}\tabularnewline
\midrule
2000 & Bias-corrected JMB & 1.424 & 1.527 & 1.516 & 1.634\tabularnewline
& Bias-corrected NPB & 1.454 & 1.595 & 1.576 & 1.754\tabularnewline
\midrule
4000 & Bias-corrected JMB & 1.384 & 1.480 & 1.430 & 1.527\tabularnewline
& Bias-corrected NPB & 1.416 & 1.543 & 1.492 & 1.642\tabularnewline
\midrule
6000 & Bias-corrected JMB & 1.526 & 1.631 & 1.413 & 1.509\tabularnewline
& Bias-corrected NPB & 1.550 & 1.681 & 1.458 & 1.595\tabularnewline
\bottomrule
\end{tabular}
\end{table}
We make the following observations regarding the simulation results.
First, as expected, interpolation of pointwise CIs exhibits substantial
under-coverage, especially for the nominal coverage probabilities
$0.90$ and $0.95.$ Therefore, appealing as it looks to practitioners,
interpolation of pointwise CIs fails to cover the true density curve
with the desired coverage probability even in large samples. Second,
the bias-corrected JMB and NPB UCBs yield reasonably good coverage
rates across different setups and sample sizes. The JMB UCBs are slightly
narrower than the NPB UCBs; however they are also less accurate. Third,
we focus on the studentized UCB because it has variable width and
thus is narrower than the non-studentized counterpart (Footnote \ref{fn:constant width}).
The additional simulation results in Section \ref{sec:Additional-Monte-Carlo}
of the online supplement confirm that the non-studentized UCBs are
on average wider than the studentized ones. The computation of the
JMB is faster than the NPB, as the former avoids the estimation of
ITEs for each bootstrap sample. Therefore, we recommend the bias-corrected
JMB UCB defined in (\ref{eq:bias corrected CB}) to practitioners
for assessing the shape of the density of ITEs.
\section{Empirical application: Childbearing and labor income\label{sec:Empirical-illustrations}}
In this section, we apply the FVX estimator for the density of the
ITE and our bias-corrected JMB UCB to investigate the effect of family
size on labor income. Understanding the relationship between the two
variables is important for policymakers; however, estimation of the
effect can be complicated due to the simultaneity between the labor
supply and fertility decisions \citep[AE hereafter]{angrist1998children}.
We revisit the 1980 Census Public Use Micro Samples (PUMS) previously
used by AE and other authors. Following AE, we focus on married women
aged 21-35 with at least two children. The focus on households with
at least two children is due to the identification strategy developed
in AE, as explained below. Our outcome variable $Y$ is the sum of
the mother's and father's 1979 labor incomes (in thousands of dollars).
The binary treatment variable $D$ takes the value one if the mother
has more than two children. The instrument proposed in AE is the ``same-sex''
dummy variable that takes the value one when the first two children
are of the same sex. This identification strategy relies on parental
preferences for a mixed sibling-sex composition: parents whose first
two children are of the same sex are more likely to have an additional
child. AE shows that in 1980, the estimated probabilities of having
a third child for women with same-sex and mixed-sex children were
0.432 and 0.372, respectively. Moreover, the difference between the
two groups is highly significant.
Our covariates $X$ include the mother's education level (less than
high school, high school, some college, college), the quartile in
the age distribution, race (white, black, Hispanic, others), and the
sex of the first birth.\footnote{The three quartiles of the distribution of mother's age in our sample
are 28, 31, and 33.} Following FVX, we drop observations in $X$-defined cells that contain
less than $1\%$ of the sample. Our remaining sample has $224,962$
observations.
Using 2SLS with a linear IV regression model, the estimated effect
of having more than two children on parents' labor income is $-3.54$
with a standard error of $1.51$. The estimates suggest that in 1979,
having more than two children reduced parents' labor income by \$3,540
(or \$12,624 in 2020 dollars).\footnote{We used the CPI series from FRED, Federal Reserve Bank of St. Louis,
for the conversion.} The effect is substantial and corresponds to 7.7\% of the average
household labor income in our sample.\footnote{The average household labor income in our sample is \$45,829 in 1979
dollars.}
\begin{table}[t]
\caption{Summary statistics for ITE estimates of the effect of having more
than two children on parents' labor income (in thousands of 1979 dollars).}
\label{table: ITE_stats}
\centering{}\resizebox{\textwidth}{!}{
\begin{tabular}{ccccccccc}
\toprule
& Mean & Std.dev & 1st decile & 1st quartile & Median & 3rd quartile & 9th decile & $\Pr[\varDelta>0]$\tabularnewline
\midrule
Full sample & -2.67 & 36.72 & -16.99 & -8.71 & -4.14 & 0.25 & 10.37 & 0.264\tabularnewline
($n=224,692$) & & & & & & & & \tabularnewline
\midrule
High School, Age> 31 & -5.01 & 30.75 & -10.93 & -7.59 & -5.45 & -3.00 & 1.03 & 0.122\tabularnewline
($n=48,129$) & & & & & & & & \tabularnewline
College, Age> 31 & 10.88 & 64.20 & -39.99 & -14.00 & -1.42 & 21.69 & 53.22 & 0.462\tabularnewline
($n=19,103$) & & & & & & & & \tabularnewline
\midrule
High School, Age $\leq$ 31 & -7.93 & 23.56 & -12.89 & -8.96 & -4.70 & -2.06 & 0 & 0.097\tabularnewline
($n=65,447$) & & & & & & & & \tabularnewline
College, Age $\leq$ 31 & 25.93 & 68.24 & -16.08 & -5.58 & 9.95 & 27.07 & 178.43 & 0.667\tabularnewline
($n=8,547$) & & & & & & & & \tabularnewline
\bottomrule
\end{tabular}}
\end{table}
Table \ref{table: ITE_stats} reports summary statistics for the ITE
estimates. According to the results, estimated ITEs display substantial
heterogeneity. For example, the median ITE in our sample is -4.14
with an interquartile range of 9.03.\footnote{The results are consistent with the findings in \citet{Froelich2013}
who also report substantial heterogeneity using the 2000 PUMS data
and quantile treatment effects.} Conditioning on the above median age and college-level education
produces an even wide range of estimated ITEs: -1.42 for the median
effect with the interquartile range of 35.69. While for the below
median age, college-educated mothers, the median effect is positive
(9.95), the corresponding interquartile range is similarly wide (32.65).
The table also shows that in the case of mothers with only high-school-level
education, ITEs tend to be more negative. For example, conditional
on the above median age and only high-school-level education, the
median ITE is -5.45, with an interquartile range of 7.93. Only 12.2\%
of the households in this group have positive estimated ITEs, compared
to 46.2\% of the households with college-educated mothers from the
same age group. The group with the largest fraction of households
with positive estimated ITEs is the below-median-age mothers with
a college education: 66.7\%.
\begin{figure}[t]
\caption{The unconditional PDF of the ITE with the 95\% pointwise and uniform
confidence bands\label{fig:unconditional-distribution}}
\centering{}\includegraphics[scale=0.8]{unconditional}
\end{figure}
Next, we use the FVX estimator with our bias-corrected JMB UCB to
analyze the distribution of the ITE. Figure \ref{fig:unconditional-distribution}
shows the unconditional PDF of the ITE together with the 95\% pointwise
and uniform confidence bands for the density. Following Remark \ref{Rmk: truncation},
the lower bounds of the confidence bands are truncated to zero. One
can see that while the UCB developed in this paper is somewhat wider
than the pointwise, it is still informative. The estimated mode of
the unconditional distribution is -4.08, and according to the UCB
the mode is located between -5.60 and -3.52.
\begin{figure}[t]
\caption{The conditional PDFs of the ITE (solid lines) and their 95\% UCBs
(shaded areas) conditional on the mother's age, high school (red),
or college (black) education: the left panel is (a) Above the median
age, and the right panel is (b) Below the median age}
\noindent \centering{}\label{fig:conditional PDFs }\includegraphics[scale=0.42]{HS_vs_Col_Age_gr_2}\includegraphics[scale=0.42]{HS_vs_Col_Age_leq_2}
\end{figure}
Figure \ref{fig:conditional PDFs } shows the conditional PDFs of
the ITE conditional on the mother's age (above or below the median
age in our sample) and education (high school only or college levels)
with their 95\% UCBs. Figure \ref{fig:conditional PDFs }(a) displays
the results conditional on the above median age for high school only
and college education levels. As there are regions where the two UCBs
do not intersect, we can conclude with at least 90.25\% confidence
that the two densities conditional on high school and college are
different.\footnote{The two 95\% UCBs are independent since they are computed on non-overlapping
samples. Hence, the overall confidence level when comparing the two
distributions is $0.95^{2}=0.9025$.} In particular, while the distribution of the ITE conditional on college
is much more dispersed, it also has more probability mass in the positive
range. E.g., on the interval corresponding to ITEs between \$20,000
and \$30,000, the confidence band for the density conditional on high
school is very narrow and close to zero. At the same time, the density
conditional on college is significantly different from zero. Hence,
a non-negligible fraction of households with college-educated mothers
above the median age received a substantial positive effect of a magnitude
between 43.6\% and 65.5\% of the average household labor income in
our sample. There is no evidence that households with only high-school-educated
mothers experienced ITEs of this magnitude.
In the case of high school only, the UCB does not rule out a bimodal
density with the two modes at -8.80 and -2.56. According to these
results, households with high-school-educated mothers above the median
age are likely to experience either a strong negative effect around
19.2\% of the average household labor income or a more moderate negative
effect around 5.6\% of the average household labor income.
Figure \ref{fig:conditional PDFs }(b) shows similar results conditional
on the below median age. We can again conclude with a 90.25\% confidence
that the conditional distributions by education level (high school
or college) are different. The ITE distribution is more dispersed
for households with college-educated mothers than for households with
high-school-educated mothers. Similarly to the previous case, the
distribution conditional on college has more mass in the positive
range than the distribution conditional on high school. The results
conditional on high school again cannot rule out bimodality; however,
this time the first mode at -6.56 corresponds to a more moderate effect,
and the second mode at 0.13 occurs in the positive range.
We conclude that there are significant differences in the distributions
of the ITE across the education levels. For households with high-school-only-educated
mothers, the distribution of the ITE is heavily concentrated in the
negative range and potentially bimodal. Households with college-educated
mothers have a wider range of ITEs. However, such households can also
experience positive effects of a large magnitude. Predicting the effect
of having more than two children on labor income for such households
is difficult as the distribution is thinly spread from large negative
to large positive values.
Recently, \citet{abrevaya2021estimation} studied the distributional
effect of having a third child on female labor supply by applying
a weakly nonseparable model (equipped with the mean-variance-effect
structure) to the PUMS dataset in 2000. Similar to our findings, they
also documented a large amount of heterogeneity in the ITE distributions
(See their Figures 2 to 5). In terms of how the ITE distribution varies
across mothers' education levels, they found less variation with the
2000 data than we document in this paper with the 1980 data.
\bibliographystyle{chicago}
\bibliography{ITE}
\oldappendix
\newcommand{\appendixname~\thesection\quad}{\appendixname~\thesection\quad}
\counterwithin{thm}{section}