Inference on Individual Treatment Effects in Nonseparable Triangular Models

Inference on Individual Treatment Effects in Nonseparable Triangular Models.

Introduction

Heterogeneous treatment effects have received increasing attention in the causal inference and policy evaluation literature (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., 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, vuong2017counterfactual and 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 (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., 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., Abadie2002,Abadie2003,Froelich2013,imbens1994identification among others) and the instrumental variable quantile regression (IVQR) model (see, e.g., 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 Calonico2014, we also propose bias-corrected UCBs using standard errors that incorporate additional variability from the estimated bias (Section (ref)). 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)).

The FVX estimator and our UCBs require bandwidth selection. Our results explicitly allow for random data-dependent bandwidths. Following the literature (e.g., 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 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., Mammen2012 and Ma2019. The asymptotic theory of the FVX estimator is different from the results obtained in Mammen2012 and Ma2019. E.g., among their other results, 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. Ma2019 show that in the context of first-price auctions, the estimation of latent bidders' valuations contributes to the asymptotic variance of the Guerre2000 estimator for the density of valuations. Moreover, the variance of the 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 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 Chen2020jackknife,Chernozhukov2014gaussian. Then theorems in Chen2020jackknife,Chernozhukov2014anti,Chernozhukov2016 show that the distributions of these Gaussian suprema can be approximated by bootstrapping.} This approach was recently applied, e.g., in 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) 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) shows the density estimator's uniform convergence rate and asymptotic normality. Section (ref) 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) describes the algorithm for our proposed confidence band. Section (ref) presents the results from Monte Carlo experiments. Section (ref) 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}).

Model and the FVX estimator

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

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)) 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

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” (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

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:

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:

The leave-$i$-out nonparametric estimator of $\phi_{dx}\left(y\right),d\in\left\{ 0,1\right\} $, can be constructed as

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)) with its leave-$i$-out nonparametric estimator $\widehat{\phi}_{dx}^{\left(-i\right)}\left(y\right)$:

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:

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:

where $b_{\zeta}$ is another bandwidth. Algorithm (ref) in Section (ref) provides detailed step-by-step instructions for constructing the UCB. In ((ref)) 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)) in Section (ref). The bootstrap critical value $z_{1-\alpha}^{\mathsf{jmb}}$ is defined in equation ((ref)) in Section (ref). The need for the second bandwidth $b_{\zeta}$ is discussed in Section (ref). We provide data-dependent procedures for selecting the two bandwidths $b$ and $b_{\zeta}$. See the discussions following Assumption (ref) in Section (ref) and Assumption (ref) in Section (ref) respectively. Theorem (ref) in Section (ref) 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.

Asymptotic properties of the FVX estimator

In this section, we establish two new asymptotic results for the FVX density estimator. Theorem (ref) below shows that under seemingly minimal conditions (see Assumptions (ref) and (ref) 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 Stone1982). Theorem (ref) 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) below.

Assumptions

The following assumption on the data generating process (DGP) is similar to those in VX and FVX.

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., Vytlacil2002). See 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), $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). 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). 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) 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) in Appendix (ref)). 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.

The smoothness assumption imposed by (a) with $P\geq1$ is stronger than that imposed by Assumption (ref)(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) 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) in Appendix (ref) 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 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.

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). 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)) under $P=2$. See li2007nonparametric. Theorem (ref) 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), one can show that such a choice satisfies Assumption (ref) 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$.

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)) below and serves as the basis for establishing the asymptotic properties of the FVX estimator. The latter are presented below as Theorems (ref) and (ref) (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:

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), $\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

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., Newey:1994jb) that

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)). 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:

Using the above definitions, we can write that the difference between the feasible and infeasible estimators as a U-statistic with a kernel that depends on the bandwidth:

uniformly in $\left(v,b\right)\in I_{x}\times\left[\underline{h},\overline{h}\right]$, where

By Assumption (ref)(b) and ((ref)) in Appendix (ref), $\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 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

By definition and since $\mathcal{G}_{x}^{\left[2\right]}\left(w,v;b\right)=0$, the 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 (Chen2020jackknife) if $\mathrm{Var}\left[\mathcal{G}_{x}^{\left[1\right]}\left(W,v;b\right)\right]>0$. In the proof of Lemma (ref), 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)) is dominated by the first term. See Remark (ref) below. Combining the result with ((ref)) and ((ref)), we can write the estimation error of the feasible estimator as

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)). 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{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)) 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) for more details on how the linearization is involved in the 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)) 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)) and the linearization in ((ref)) 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), ((ref)) and ((ref)) 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)) 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) below, which is a refinement of Theorem 2 in FVX. In Appendix (ref), we prove a stronger version (Theorem (ref)) of Theorem (ref). 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).

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$.

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)) and the linearization ((ref)) 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)) - ((ref))) 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)), we show that the same normality result holds if $h$ is replaced by its estimator $\widehat{h}$. The result is analogous to those in 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}$.

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) with no controls $X$. The treatment status $D$ is determined by the index model in ((ref)) 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.

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), (ref) and (ref) 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).

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). By the same arguments as those used to establish Theorems (ref) and (ref), one can show that

where

Note that $\varDelta_{j}$ can be expressed as a function of $W_{j}$ (see ((ref))). 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)) 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))), where $V\left(v,x;b\right)\coloneqq V_{1}\left(v,x;b\right)+V_{2}\left(v,x;b\right)$ and

The plug-in estimator of the $V_{1}\left(v,x;b\right)$ term is given by

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

be the reweighted kernel estimator proposed by 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 U-statistic with an estimated kernel:

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)) 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 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

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). 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).

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) 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) and (ref) 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) below provides a uniform convergence rate for $\widehat{V}\left(v\mid x;\widehat{h},\widehat{h}_{\zeta}\right)$. In Appendix (ref), Theorem (ref) 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). The stronger result of Theorem (ref) is used in the proof of Theorem (ref) below.

A pointwise $1-\alpha$ asymptotic confidence interval for $f_{\Delta\mid X}\left(v\mid x\right)$ can be constructed as

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)) over the domain are invalid in the uniform sense. In the section below, we discuss the construction of valid UCBs.

Jackknife multiplier bootstrap UCB

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:

Moreover, recall the expansion of the estimation error in ((ref)). 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)) 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) 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)), 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 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., 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:

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)). 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

We show in Appendix (ref) (the proof of Theorem (ref)) 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

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)). 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{Alternatively, a constant-width UCB (e.g., 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.

Theorem (ref) 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., Chernozhukov2014anti for discussion).

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).

Bias-corrected JMB UCB

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 Calonico2014. We assume that the third-order derivatives in Assumption (ref) with $P=2$ are Lipschitz continuous. Theorem (ref) 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)(a). The bias-corrected estimator of $f_{\Delta\mid X}\left(v\mid x\right)$ is given by

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., Xu2017) or the same bandwidth as that used for estimating the density 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 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)-(ref) with $P=2$, the bias part of the bias-corrected estimator in ((ref)) 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 (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)) 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

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) 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) (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.

Conditioning on sub-vectors of the covariates

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). 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

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

Define also

$\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). 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)) by adapting Algorithm (ref).

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

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)) 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).

Monte Carlo experiments

This section evaluates the finite-sample performance of the UCBs proposed in Section (ref) 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

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$.

Tables (ref) reports the simultaneous coverage rates of two types of UCBs: the bias-corrected JMB UCB defined in ((ref)) and the bias-corrected nonparametric bootstrap (NPB) UCB described in Appendix (ref). 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)), Lemma (ref) in the supplement and standard arguments (see the proof of 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)) 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) 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.}

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)). The additional simulation results in Section (ref) 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)) to practitioners for assessing the shape of the density of ITEs.

Empirical application: Childbearing and labor income

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 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.}

Table (ref) 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 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%.

Next, we use the FVX estimator with our bias-corrected JMB UCB to analyze the distribution of the ITE. Figure (ref) shows the unconditional PDF of the ITE together with the 95% pointwise and uniform confidence bands for the density. Following Remark (ref), 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.

Figure (ref) 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)(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)(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, 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.

\oldappendix {\appendixname \thesection\quad}

\counterwithin{thm}{section}