EconBase
← Back to paper

Extrapolation in Regression Discontinuity Design Using Comonotonicity

Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.

77,129 characters · 11 sections · 38 citation commands

Rendered from LaTeX for readability, not typeset faithfully. Citation keys are highlighted; maths is left as source; figures, tables and equation environments are summarised rather than reproduced; unrecognised commands are greyed out so nothing is silently dropped. Email addresses are removed.

Extrapolation in Regression Discontinuity Design Using Comonotonicity

abstractWe present a novel approach for extrapolating causal effects away from the margin between treatment and non-treatment in sharp regression discontinuity designs with multiple covariates. Our methods apply both to settings in which treatment is a function of multiple observables and settings in which treatment is determined based on a single running variable. Our key identifying assumption is that conditional average treated and untreated potential outcomes are comonotonic: covariate values associated with higher average untreated potential outcomes are also associated with higher average treated potential outcomes. We provide an estimation method based on local linear regression. Our estimands are weighted average causal effects, even if comonotonicity fails. We apply our methods to evaluate counterfactual mandatory summer school policies.

Regression discontinuity design (RDD) provides a means for identifying causal effects without the need to isolate a source of exogenous variation in treatment. Instead, causal effects are identified under an a priori continuity assumption which is plausible in many empirical settings. Unfortunately, while causal inference in RDD design may be credible, its scope is limited. Standard RDD analysis can identify average causal effects only within the sub-population whose covariates place them on the margin between treatment and non-treatment. As such, identification of the causal impact of counterfactual policies is limited to those under which there is at most a marginal change in the treatment rule.

To be more precise, consider the sharp RDD framework which is the focus of the present work. In sharp RDD, the binary treatment $D$ is deterministic given some observable individual characteristics $X$. For example, treatment may be assigned if and only if a particular observable covariate known as the `running variable' exceeds a certain threshold. In other cases treatment may be assigned by a more complicated rule that depends on multiple observables. Sharp RDD settings arise when receipt of an intervention is decided by some central authority purely on the basis of individual characteristics that are observed by both the authority and researcher. For example, eligibility for a particular welfare program may be decided according to a strict rule based on an individual's age, income, household size, and marital status. Let $Y$ be an outcome of interest and $Y(d)$ the potential outcome from treatment level $d$. In RDD, treatment is entirely explained by the covariates, that is, the value of $D$ depends entirely on whether $X$ lies in a `treatment region' or a `no-treatment region'. Under a continuity assumption, the conditional average treatment effect (CATE), given by $E[Y(1)-Y(0)|X=x]$, is identified for values of $x$ that lie on the `frontier' between the treatment and no-treatment regions. Without further assumptions, the CATE is not identified away from the frontier.

In this work we provide a novel means for identifying causal effects away from the frontier in sharp RDD settings in which $X$ is multivariate: that is, more than one covariate is available. We show that our approach is valid when a particular comonotonicity (in the sense of e.g., Schmeidler1989) condition holds. In its simplest form, our comonotonicity condition states that the mean untreated potential outcome is weakly greater conditional on a covariate value $x$ that conditional on a value $x'$, if and only if the conditional average treated potential outcome is also greater for value $x$ than $x'$. For example, if wealthier and older individuals have higher potential outcomes on average when they are not treated, they must also have higher average potential outcomes when treated. This assumption has testable implications. Under the standard RDD assumptions, both treated and untreated conditional average potential outcomes are identified along the frontier, and so one can assess whether the comonotonicity condition holds for covariate values on the frontier.

The comonotonicity assumption suggests a straightforward means of extrapolating conditional average treatment effects away from the frontier. Given a vector of covariate values $x$ that lie in say, the interior of the treatment region, we estimate the conditional average treated potential outcome for these values $E[Y(1)|X=x]$, which is equal to the conditional mean of the factual outcome $E[Y|X=x]$. We then find a point $x^*$ on the frontier with the same associated conditional mean treated potential outcome. The conditional means of both treated and untreated potential outcomes are identified on the frontier under standard RDD assumptions. If such a point exists, then the conditional average treatment effect at $x$ is identical to that at $x^*$.

We develop a novel empirical approach based on our identification results. In particular, we develop a non-parametric method based on local linear regression. Notably, the methods do not require applied researchers to specify the frontier between the treatment and non-treatment region and thus our techniques are applicable even when treatment status is observed but the assignment rule is not a priori known. This facilitates implementation of our methods, an applied researcher need only input the outcomes, covariates, and treatments into a provided software package <>, with no need to encode the formula for the treatment rule.

Our estimator targets a particular weighted average of causal effects for individuals along the frontier. Thus our estimand is a weighted average of causal effects, regardless of whether comonotonicity holds. Comonotonicity then ensures that this weighted average is equal to the CATE at a chosen covariate value away from the frontier.

We apply our methods to analyse the impact of mandatory summer school attendance. This empirical setting was previously analysed by Matsudaira2008. We use our methods to examine heterogeneity in treatment effects.

Related Literature

This paper contributes to the extensive literature on RDD (see Cattaneo2022 for a recent review). A number of works use multivariate RDD to test a kind of selection-on-observables assumption. These include Battistin2008 and Mealli2012. A closely related approach is that of Angrist2015. Angrist2015 extrapolate treatment effects away from the cut-off by assuming that, conditional on some covariates, a scalar running variable is unrelated to potential outcomes. Extrapolation from the cut-off using higher order derivatives is considered in DiNardo2011. A full statistical analysis using first derivatives and extension to fuzzy RDD is given in Dong2015 who identify the impact of a marginal change in the threshold. Cattaneo2021 uses the existence of different cut-offs for certain sub-populations to extrapolate causal effects in univariate RDD under a type of parallel trends assumption. Another strand of the literature identifies causal effects in RDD away from the frontier by using auxiliary information, for example, Grembi2016 and Wing2013.

A number of other works consider RDD settings with one or more covariates in addition to a running variable. Calonico2019 and Noack2024 develop inference methods that can be applied to RDD setting with covariates, where the primary role of covariates is efficiency gain. In our work, we exploit the availability of additional covariates in order to identify richer causal objects rather than for statistical precision. Our paper also relates to the literature that considers multivariate running variables. For example, papers such as Papay2011 in the canonical RD design or Keele2015 in the geographic RDD setting. The methods by armstrong2018optimal and Imbens2019 that put a bound on the smoothness and work out the optimal bias-variance trade-off naturally apply to the multivariate setting as well.

Comonotonicity is closely related to, but distinct from, rank invariance and rank similarity of potential outcomes. An extensive literature employs rank invariance or similarity assumptions in order to achieve causal inference using instrumental variables. For example, Chernozhukov2005, Chernozhukov2007, Horowitz2007, Imbens2009, Torgovitsky2015, and DHaultfoeuille2015. In contrast to the rank invariance/similarity assumptions employed in these works, comonotonicity restricts conditional mean treated and untreated potential outcomes rather than the joint distribution of treated and untreated potential outcomes.

We provide further comparison with Angrist2015 and IV methods using rank-invariance in Appendix B.4.

A Simple Demonstration

Before we present formal results, let us begin with a worked example. For concreteness, we frame this in a simplified version of the empirical setting in Section 5, however we have not explicitly calibrated the example to that data.

We suppose the vector of covariates $X$ is two-dimensional, that is $X=(X_{1},X_{2})'$ where $X_{1}$ and $X_{2}$ are both scalars. In line with the application in Section 5, let us suppose $X_1$ and $X_2$ are, respectively, students scores on a math and reading test. We assume that an individual receives treatment, which we indicate by $D=1$, if and only if a linear combination of the individual's test scores falls below a threshold. In our empirical application, treatment is mandatory attendance of summer school. Thus treatment is deterministic given the covariates: the individual receives treatment if the vector of test scores $X$ falls in some region $\mathcal{X}_{1}$. Conversely, the individual does not receive treatment ($D=0$) if $X$ is in the region $\mathcal{X}_{0}$. Together $\mathcal{X}_{1}$ and $\mathcal{X}_{0}$ make up the support of $X$.

In the Figure 1 we have normalized $X_{1}$ and $X_{2}$ so that their joint support is the unit square. An individual is treated when $X$ lies in the turquoise region labelled $\mathcal{X}_{1}$ and is not treated if $X$ falls into the magenta region marked $\mathcal{X}_{0}$. The intersection of the boundaries of these two regions is marked with a blue line labelled $\mathcal{F}$ which we refer to as the `frontier'. To be precise, the figure indicates that an individual is treated if $0.4 X_{1}+X_{2}\leq0.7$.

figure[figure omitted — 128 chars of source]

Consider an outcome, $Y$ that is an individual's normalized score on a math test taken a year after the initial tests in $X$ and thus after any treatment has been administered. The dashed curves in the figure are contour curves (level sets) of the conditional average outcome $E[Y|X=x]$ within the treated or untreated regions. That is, if two points $x_1,x_2\in\mathcal{X}_{1}$ lie on the same contour, then $E[Y|X=x_1]=E[Y|X=x_2]$ and likewise if $x_1,x_2\in\mathcal{X}_{0}$, but not if say, $x_1\in\mathcal{X}_{1}$ and $x_2\in\mathcal{X}_{0}$.

In Figure 1, $(x_{1}^{*},x_{2}^{*})$ lies on the frontier $\mathcal{F}$. As such, both the conditional average treated and untreated potential outcomes $E[Y(0)|X=x^{*}]$ and $E[Y(1)|X=x^{*}]$, are identified at this point under standard RDD assumptions. The point $(x_{1}^{\circ},x_{2}^{\circ})$ is in the interior of the treatment region $\mathcal{X}_{1}$ and as such, $E[Y(1)|X=x^{\circ}]$ is identified, indeed this quantity simply equals $E[Y|X=x^{\circ}]$. However, without further assumptions we cannot identify $E[Y(0)|X=x^{\circ}]$ nor the conditional average treatment effect $E[\tau|X=x^{\circ}]$.

To progress, we make the following `comonotonicity' assumption. Consider two vectors of initial math and reading test scores $x=(x_1,x_2)$ and $x'=(x_1',x_2')$. We suppose $E[Y(0)|X=x]$ is greater than $E[Y(0)|X=x']$ if and only if $E[Y(1)|X=x]$ is greater than $E[Y(1)|X=x']$. That is, the functions $E[Y(0)|X=\cdot]$ and $E[Y(1)|X=\cdot]$ are `comonotonic' in the sense of e.g., Schmeidler1989. In the present context, this assumption states that initial test scores associated with higher average math outcomes under the counterfactual of no-treatment, are also associated with higher average math scores under the counterfactual of treatment. We provide further discussion of the assumption in the context of our empirical application in Section 5. In addition, Appendix B.1 specifies a model of skill accumulation motivated by the test score application, and we show that the assumption holds under this model.

Note that $x^{\circ}$ lies on the same contour as $x^{*}$ which is indicated by the solid black curve. For all values of $x\in\mathcal{X}_{1}$ along this contour, $E[Y|X=x]$ and thus $E[Y(1)|X=x]$, is constant. It follows from the comonotonicity condition that $E[Y(0)|X=x^{*}]$ is identical to $E[Y(0)|X=x^{\circ}]$. Therefore, the conditional average treatment effect is the same at $x^{\circ}$ as at $x^{*}$, and so $E[\tau|X=x^{\circ}]$ is identified.

In short, we extrapolate conditional average treatment effects away from the frontier and into the treatment region by matching values of $x^{\circ}$ with a value $x^{*}$ on the frontier with the same conditional average treated potential outcome. Similarly, we can extrapolate into the interior of the no-treatment region by finding covariate values on the frontier with the same conditional average untreated potential outcomes.

Note that this approach does not always allow us to extrapolate to all points in the support of $X$. For example, for some $x\in\mathcal{X}_{1}$ there may be no $x'\in\mathcal{F}$ such that $E[Y(1)|X=x']=E[Y(1)|X=x]$. Indeed, this is true for the point $(x_{1}^{+},x_{2}^{+})$ in Figure 1. Consider the extreme case in which the contours of $E[Y|X=x]$ are parallel to the frontier. Then we cannot extrapolate anywhere. Conversely, if contours are perpendicular to the frontier and $X$ has sufficiently large support, then we can extrapolate everywhere. Which case is closer to the truth depends on the context. Under comonotonicity the contours are contours of treatment effects. If a social planner wishes to treat individuals with high causal effects, it is optimal to set the frontier equal to a contour. In practise, the social planner may be constrained to employ simple treatment rules, may not be able to accurately assess conditional average causal effects, or may care about more than one outcome or outcomes that differ from those of interest to the researcher. Such cases may be intermediate to these two extreme examples allowing us to extrapolate causal effects to some but not all points in the support of $X$. Indeed, this appears to be the case in our empirical application.

Recall that along the frontier, both the treated and untreated conditional average potential outcomes can be identified under standard RDD assumptions. On the left panel in Figure 1.2, we plot conditional average potential outcomes at each point along the frontier as a function of the $x_1$-coordinate of that point. The conditional average treated potential outcome is in turquoise and the conditional average untreated potential outcome in magenta. The vertical dashed blue line is at the location $x^*_1$ and intersects each of the curves at the corresponding dashed horizontal line. The elevation of the dashed horizontal turquoise line is equal to $E[Y(1)|X=x^*]$ which is also equal to $E[Y(1)|X=x^\circ]$. The elevation of the dashed magenta line is $E[Y(0)|X=x^*]$.

On the right panel in Figure 1.2, we plot the pairs $(E[Y(1)|X=x],E[Y(0)|X=x])$ for all values of $x$ along the frontier $\mathcal{F}$. This curve is shown in blue. The curve is strictly increasing which suggests that along the frontier values of $x$ associated with a greater conditional average treated potential outcomes also have greater conditional average untreated potential outcomes. We can also see this on the left panel, an $x_1$ with a larger associated value of $E[Y(0)|X=x]$ also has a larger value of $E[Y(1)|X=x]$. Note that the comonotonicity condition implies that the blue curve in the right panel must be strictly increasing, and this can be used as the basis for a test of this assumption.

figure[figure omitted — 208 chars of source]

The vertical dashed line (shown in turquois) in the right panel of Figure 1.2 is located at the value $E[Y(1)|X=x^{\circ}]$. As we discuss above, this quantity is trivially identified. This is the same value as $E[Y(1)|X=x^{*}]$, which is located along the frontier. The corresponding value $E[Y(0)|X=x^{*}]$ is thus equal to the elevation of the horizontal dashed line, shown in magenta. Under the comonotonicity condition, this is equal to the value of $E[Y(1)|X=x^{\circ}]$.

More generally, for a given point $x$ in the treated region, we can find the value of $E[Y(0)|X=x]$ from the right panel by the $y$-coordinate of the point on the blue curve with $x$-coordinate equal to $E[Y|X=x]$. Conversely, if $x$ is in the no-treatment region, then $E[Y(1)|X=x]$ is the $x$-coordinate of the point on the blue curve whose $y$-coordinate is $E[Y|X=x]$. For some values of $x$ there may be no corresponding point on the blue curve in which case we cannot extrapolate conditional average treatment effects to that value of $x$.

We can understand the blue curve in Figure 1.2 as a function that returns the value of $E[Y(0)|X=x]$ given the value of $E[Y(1)|X=x]$. We denote this function by $q_0$ and conversely, the function that returns $E[Y(1)|X=x]$ in given $E[Y(0)|X=x]$ by $q_1$.

Formal Identification Results

We now state our formal results, we begin by re-stating the standard sharp RDD identification result that conditional average treatment effects are identified for covariate values on the frontier between the treatment and no-treatment regions. Again we let $Y$ be a scalar outcome of interest, $D$ a binary treatment indicator. We denote by $Y(d)$ as the potential outcome under treatment level $d$, $\tau=Y(1)-Y(0)$ as the individual level treatment effect, and $X$ a vector of covariates with support $\mathcal{X} \subset \mathbb{R}^{d_{X}}$. The support $\mathcal{X}$ is partitioned into a treatment region $\mathcal{X}_{1}$ and no-treatment region $\mathcal{X}_{0}$.

The frontier consists of those covariate values that lie on the margin between the treated and untreated regions. To define this formally, let $int(\mathcal{X}_1)$ denote the interior of the treated region and $int(\mathcal{X}_0)$ the interior of the untreated region. Then the frontier $\mathcal{F}$ is the intersection of the boundaries between these two sets. Equivalently, letting $cl(\cdot)$ return the closure of a set, \[ \mathcal{F}=cl\big(int(\mathcal{X}_{1})\big)\cap cl\big(int(\mathcal{X}_{0})\big). \]

If $\mathcal{X}_{1}$ and $\mathcal{X}_{0}$ are sufficiently regular then the above is simply the intersection of their boundaries.\footnote{In particular, if the closures of $\mathcal{X}_{1}$ and $\mathcal{X}_{0}$ are `regular closed sets'.} The use of the interiors in the definition allows us to avoid cluttering the analysis with regularity conditions on $\mathcal{X}_{1}$ and $\mathcal{X}_{0}$.

Standard RDD analysis identifies conditional average treatment effects at the frontier under the conditions below.

\theoremstyle{definition} \newtheorem*{A01}{Assumption 1.1} \begin{A01}[Deterministic Treatment] $D=1\{X\in\mathcal{X}_{1}\}$, where $\mathcal{X}_{1}$ and $\mathcal{X}_{0}$ do not overlap and their union is the support of $X$. \end{A01}

\theoremstyle{definition} \newtheorem*{A02}{Assumption 1.2} \begin{A02}[Continuity] The functions $x\mapsto E[Y(0)|X=x]$ and $x\mapsto E[Y(1)|X=x]$ are continuous. \end{A02}

Assumption 1.1 simply states that treatment is deterministic given $X$ and defines the treatment and no-treatment regions $\mathcal{X}_{1}$ and $\mathcal{X}_{0}$. Assumption 1.2 imposes that the conditional average potential outcomes are continuous. Strictly speaking, in order to identify conditional average treatment effects at the frontier, it would suffice to replace Assumption 1.2 with the weaker condition that the conditional expectations are continuous at each $x$ on the frontier. However, the stronger continuity condition allows us to sidestep some technicalities.\footnote{In particular, when $X$ is continuously distributed, conditional expectations like $E[Y(0)|X=x]$ and $E[Y(1)|X=x]$ are not uniquely defined. From a measure-theoretic perspective, $x\mapsto E[Y(d)|X=x]$ is a measurable function $g$ on the support of $X$ so that for any measurable set $\mathcal{S}$ with $P(X\in\mathcal{S})>0$, $E[g(X)|X\in\mathcal{S}]=E[Y(d)|X\in\mathcal{S}]$, but such a function $g$ is not unique. However, if there exists a continuous $g$ with this property, then it is unique. Thus we understand $x\mapsto E[Y(d)|X=x]$ for $d=0,1$ to refer to the unique continuous conditional expectation functions. We employ this convention throughout the text.}

Theorem 1 below is standard, it establishes identification of the conditional average potential outcomes, and thus conditional average treatment effect $E[\tau|X=x]$ for each $x$ on the frontier $\mathcal{F}$.

\theoremstyle{plain} \newtheorem*{T1}{Theorem 1} \begin{T1} Suppose Assumptions 1.1 and 1.2 hold. Then $E[Y(1)|X=x]$ and $E[Y(0)|X=x]$ are identified for every $x\in\mathcal{F}$.

More precisely, for each $d=0,1$ there is a function $g_d$ on $\mathcal{X}_d\cup \mathcal{F}$ so that for all $x\in \mathcal{X}_d$, $g_d(x)=E[Y|X=x]$, and $g_d$ is continuous at each point in $\mathcal{F}$. For each $x\in\mathcal{F}$, conditional average potential outcomes are identified as follows:

equation[equation omitted — 54 chars of source]

\end{T1}

Extrapolation Under Comonotonicity

Theorem 1 identifies conditional average potential outcomes at the frontier. In order to identify effects away from the frontier we require an additional identifying assumption. \theoremstyle{definition} \newtheorem*{A2}{Assumption 2} \begin{A2}[Comonotonicity] for any $x_{2},x_{1}\in\mathcal{X}$:

equation*[equation* omitted — 94 chars of source]

\end{A2}

Assumption 2 states that values of the covariates associated with higher average untreated potential outcomes also correspond to higher treated potential outcomes.Before we discuss Assumption 2 in greater depth and provide sufficient conditions, let us first state a key consequence of Assumption 2 which constitutes our main identification result. \theoremstyle{plain} \newtheorem*{T2}{Theorem 2} \begin{T2} Suppose Assumptions 1 and 2 hold and define $g_0$ and $g_1$ as in Theorem 1. Suppose that for some $x\in\mathcal{X}_{d}$, there exists an $x^{*}\in\mathcal{F}$ so that $E[Y|X=x]=g_{d}(x^{*})$. Then $E[Y(1)|X=x]=g_{1}(x^{*})$ and $E[Y(0)|X=x]=g_{0}(x^{*})$. Moreover, for any $x_1,x_2\in\mathcal{F}$, $g_0(x_1)\geq g_0(x_2)\iff g_1(x_1)\geq g_1(x_2)$. \end{T2}

The first statement of Theorem 2 shows how comonotonicity allows us to impute conditional average potential outcomes away from the frontier for certain values of $x$. The second statement notes that under our assumptions, $g_1(x)$ is increasing in the corresponding value of $g_0(x)$ along the frontier and vice versa. The latter condition suggests that comonotonicity is testable because $g_0(x)$ and $g_1(x)$ are identified on the frontier.\footnote{An implication (in fact, a sharp testable implication) of comonotonicity is $$\inf_{x,x' \in \mathcal{F}} (g_{d}(x) - g_{d}(x'))(g_{1-d}(x) - g_{1-d}(x')) \geq 0.$$ Hence, one can directly apply the methods of chernozhukov2013intersection to falsify comonotonicity. However, since $g_{d}$ and $g_{1-d}$ can be estimated independently, we suspect that more powerful tests tailored to our setting can be developed. We view this as an interesting direction for future research.}

commentWhen we consider estimators motivated by Theorem 2 it is helpful to re-frame the result as follows. Define the function $q_{1-d}$ so that for any $x$ on the frontier, \[q_{1-d}\big(E[Y(d)|X=x]\big)=E[Y(1-d)|X=x].\] The domain of this function consists of all the values of $E[Y(d)|X=x]$ for some $x$. Under regularity conditions on the frontier this set is an interval with lower and upper end points $\underline{y}_{1-d}$ and and $\bar{y}_{1-d}$ (which can be positive or negative infinity) defined below: \begin{equation} y_{1-d}=\inf_{x\in\mathcal{F}}E[Y(d)|X=x],\,\,\,\,\,\bar{y}_{1-d} =\sup_{x\in\mathcal{F}}E[Y(d)|X=x] \end{equation} Under the conditions of Theorem 2, the function $q_{1-d}$ is well-defined and weakly increasing. Note that the curve in Figure 1.2 is precisely the function $q_0$. Now, consider some $x\in\mathcal{X}_d$ not necessarily on the frontier. Under the conditions of Theorem 2, if $E[Y|X=x]$ is in the domain of $q_{1-d}$, then the conditional mean potential outcomes at $x$ are identified as follows: \begin{align} E[Y(d)|X=x]&=E[Y|X=x]\nonumber\\ E[Y(1-d)|X=x]&=q_{1-d}\big(E[Y|X=x]\big) \end{align} That is, the conditional average treated and untreated potential outcomes can be written in terms of the conditional mean outcome and the function $q_{1-d}$ which is identified by standard RDD arguments.

When we consider estimators motivated by Theorem 2 it is helpful to re-frame the result as follows. Under the conditions of Theorem 2, there is a unique, increasing function $q_{1-d}$ so that for every $x\in\mathcal{F}$,

equation[equation omitted — 71 chars of source]

The domain of $q_{1-d}$ consists of all the values of $E[Y(d)|X=x]$ for some $x\in\mathcal{F}$. Under regularity conditions on the frontier this set is an interval with lower and upper end points $\underline{y}_{1-d}$ and and $\bar{y}_{1-d}$ (which can be positive or negative infinity) with formulae

equation[equation omitted — 149 chars of source]

Under the conditions of Theorem 2, $q_{1-d}$ is weakly increasing. The curve in Figure 1.2 is precisely the function $q_0$.

Now, consider some $x\in\mathcal{X}_d$ not necessarily on the frontier. Under the conditions of Theorem 2, if $E[Y|X=x]$ is in the domain of $q_{1-d}$, then the conditional mean potential outcomes at $x$ are identified as follows.

align[align omitted — 104 chars of source]

That is, the conditional average treated and untreated potential outcomes can be written in terms of the conditional mean outcome and the function $q_{1-d}$, which are both identified. Conditional average treatment effects may thus be identified at a point $x\in\mathcal{X}_d$ by

equation[equation omitted — 90 chars of source]

If comonotonicity does not hold, then there may be no function that satisfies ((ref)) and so it is helpful to use the following definition of $q_{1-d}$, which coincides with ((ref)) under the conditions of Theorem 2.

equation[equation omitted — 90 chars of source]

This function remains well-defined even if comonotoncity fails. Moreover, using this definition, the right hand side of ((ref)) is a weighted average of conditional average treatment effects, regardless of whether comonotonicity holds.

Proposition 1 states two sufficient, but by no means necessary, conditions for Assumption 2. The first sufficient condition states that if average untreated potential outcomes are higher within some covariate stratum, then average treatment effects are also greater in this stratum. The second condition states that, in a certain sense, there is less heterogeneity in treatment effects between strata of $X$ than heterogeneity in average baseline potential outcomes.

\theoremstyle{plain} \newtheorem*{P1}{Proposition 1} \begin{P1} Either of the following two conditions implies Assumption 2.

a. For all $x_{2},x_{1}\in\mathcal{X}$ \[ E[Y(0)|X=x_{1}]\geq E[Y(0)|X=x_{2}] \implies E[\tau|X=x_{1}]\geq E[\tau|X=x_{2}]. \]

b. For all $x_{2},x_{1}\in\mathcal{X}$

equation[equation omitted — 90 chars of source]

and the inequality is strict when the quantity on the RHS is non-zero. \end{P1}

An immediate corollary of Proposition 1 is that, if average treatment effects do not differ between strata of $X$, then comonotonicity holds. Of course, homogeneity in treatment effects makes extrapolation away from the frontier trivial. Proposition 1 demonstrates two ways in which our assumptions weaken treatment effect homogeneity. In particular, it suffices for Assumption 2 that differences in average treatment effects between covariate strata have the same sign as differences in average untreated potential outcomes, or that the difference in treatment effects is smaller than the difference in baseline potential outcomes.

In a randomized controlled trial (RCT), both conditional average treated and untreated potential outcomes are directly identified. Thus one can formally test or informally assess whether the assumption of comonotonicity holds in the setting of a given RCT. In Appendix B.5 we provide an informal empirical assessment of whether comonotonicity is credible in the context of the RCT in Alan. The empirical setting in Alan is apt in that it shares key features with the setting for our empirical application, in particular both the outcomes and covariates consist of students' reading and math test scores as in our application.

In particular settings, comonotonicity may be justified by primitive modelling conditions. In Appendix B.1 we show that the assumption holds under a model of skill formation and noisy test scores that is motivated by our empirical application.

Comonotoncity is related to, but distinct from, rank invariance of potential outcomes. This rank invariance condition states that if we draw two individuals $i$ and $j$ iid from the population, so that $(Y_i(0),Y_i(1))$ and $(Y_j(0),Y_j(1))$ are independent copies of $(Y(0),Y(1))$, then $Y_{i}(0)\geq Y_{j}(0)\iff Y_{i}(1)\geq Y_{j}(1)$. Note that rank-invariance has to hold for every pair of individuals in the population. By contrast, Assumption 2 concerns only averages of potential outcomes between subgroups. Moreover, rank invariance of potential outcomes is not verifiable, even with an ideal randomized controlled trial, because it involves a comparison between an individual's outcomes in two distinct counterfactual worlds, whereas all of the objects in Assumption 2 can be identified in an ideal randomized experiment.

In a formal sense, rank-invariance of potential outcomes is neither stronger nor weaker than Assumption 2 (that is, neither one implies the other). However, one can show that rank-invariance implies a conditional quantile version of Assumption 2:

equation[equation omitted — 111 chars of source]

where $Q_{Y(d)|X}(t|x)$ is the $q$-th conditional quantile of $Y(d)$ given $X=x$. Thus if the conditional distributions of potential outcomes are symmetric, then rank-invariance implies Assumption 2 (because the conditional means and medians are identical). Note that if we were to assume rank invariance in place of Assumption 2, we could use ((ref)) to identify conditional quantiles of potential outcomes and in some cases the entire distribution of individual treatment effects within some covariate strata. We discuss this further in Appendix B along with formal results.

Conditional Comonotoncity

In cases with strictly more than two covariates, the comonotonicity condition is generally stronger than necessary for identification. Suppose $X$ has three or more components. Let $X^{(1)}$ be the subvector containing precisely two of the characteristics in $X$ and let $X^{(2)}$ be the subvector that contains the remaining characteristics. Assumption 3 below weakens Assumption 2 in that it only requires comonotonicity to hold within each stratum of the characteristics in $X^{(2)}$ and not between strata.

\theoremstyle{definition} \newtheorem*{A3}{Assumption 3} \begin{A3}[Conditional Comonotonicity] For any $x_{1},x_{2}\in\mathcal{X}$ with $x_{1}^{(2)}=x_{2}^{(2)}$, we have:

equation*[equation* omitted — 95 chars of source]

\end{A3}

Theorem 3 is similar to Theorem 2. Again we extrapolate conditional average potential outcomes from a point on the frontier $x^*$ to a point $x$ away from the frontier. However, Theorem 3 strengthens the requirement on $x^*$ so that its subvector $x^{*(2)}$ is identical to the corresponding subvector $x^{(2)}$ of $x$.

\theoremstyle{plain} \newtheorem*{T3}{Theorem 3} \begin{T3} Suppose Assumptions 1 and 3 hold and define $g_0$ and $g_1$ as in Theorem 1. Suppose that for some $x\in\mathcal{X}_{d}$, there exists an $x^{*}\in\mathcal{F}$ so that $x^{*(2)}=x^{(2)}$ and $E[Y|X=x]=g_{d}(x^{*})$. Then $E[Y(1)|X=x]=g_{1}(x^{*})$ and $E[Y(0)|X=x]=g_{0}(x^{*})$. Moreover, for any $x_1,x_2\in\mathcal{F}$,with $x_1^{(2)}=x_2^{(2)}$, $g_0(x_1)\geq g_0(x_2)\iff g_1(x_1)\geq g_1(x_2)$. \end{T3}

The second statement in Theorem 3 states the weaker condition that Assumption 3 places on $g_0$ and $g_1$ along the frontier. Whereas under full comonotonicity they must move in the same direction, under the conditional comonotonicity assumption, this co-movement need only apply within each fixed stratum of $X^{(2)}$.

In the appendix we provide further alternatives that weaken Assumption 2 and which may still be sufficient for identification. In particular, we consider local comonotonicity. In this case, we require that the vector of first derivatives of the conditional average treated and untreated potential outcomes are identical. Moreover, we provide an even weaker, conditional version of local comonotonicity, akin to the conditioning in Assumption 3.

Estimation

Our analysis motivates estimates of conditional average potential outcomes and treatment effects. From these, one may obtain estimates of the causal effect of an alternative treatment allocation policy. For example, the average effect of treating all individuals. We first consider estimation under full comonotonicity based on Theorem 2. We then consider estimation under the conditional comonotonicity assumption based on the identification result in Theorem 3.

Conditional Average Potential Outcome Estimates

We obtain conditional average potential outcome estimates by the plug-in principle. That is, we replace the functions $E[Y|X=\cdot]$, $q_{0}$, and $q_{1}$ with estimates in ((ref)). Let $\hat{q}_{d}$ be an estimate of $q_{d}$ for $d=0,1$. We estimate $E[Y|X=x]$ separately for $x$ in the treated and untreated regions, using only the data for individuals within the corresponding region. Let $\hat{g}_d(x)$ denote the estimate of $E[Y|X=x]$ for $x\in\mathcal{X}_d$. Consider an $x$ in the treated region $\mathcal{X}_1$. We estimate the conditional average potential outcomes as follows.

align[align omitted — 133 chars of source]

Taking the difference between these two estimates yields an estimate of the conditional average treatment effect at $x$: \[ \hat{\tau}(x)=\hat{g}_1(x)-\hat{q}_{0}\big(\hat{g}_1(x)\big) \] Conversely, for $x$ in the untreated region $\mathcal{X}_0$, we simply switch $0$ with $1$ in the formula above and estimate conditional average causal effects accordingly.

align*[align* omitted — 199 chars of source]

We obtain $\hat{g}_1$ by non-parametric regression of $Y$ on $X$, using only data on individuals in the treated region. Similarly, for $x$ in the untreated region, we obtain $\hat{g}_0(x)$ by regression using only data in the untreated region. We suggest local linear regression for this purpose but one could use alternative non-parametric regression methods. The local linear regression estimates take the form below for $d=0,1$, where $K$ is a kernel (e.g., uniform, Gaussian, or triangular) and $h$ is a bandwidth:

align*[align* omitted — 200 chars of source]

where we define $\mathcal{I}_{d} := \{i:X_{i}\in\mathcal{X}_{d}\}$ set of indices for the observations with treatment status $d$.

Estimation of $q_{0}$ and $q_1$ is a little more involved. For ease of exposition, we focus below on estimation of conditional average untreated potential outcomes in the treated region and thus on the estimation of $q_0$. Estimation of conditional average treated potential outcomes in the untreated region proceeds symmetrically.

Recall the definition of $q_{1-d}$ given in equation ((ref)). We see $q_0(y)$ is the mean untreated potential outcome among individuals on the frontier whose conditional mean treated potential outcome is $y$. This characterization motivates a nonparametric regression-based approach.\footnote{Note that under comonotonicity, $q_{1-d}$ is also uniquely defined by ((ref)). However, if we replace $E[Y(1)|X=x]$ and $E[Y(0)|X=x]$ with noisy regression estimates $\hat{g}_1(x)$ and $\hat{g}_0(x)$, there is no guarantee that there is a functional relationship between them on the frontier. That is, there may not exist any function $\tilde{q}_0$ so that $\tilde{q}_0\big(\hat{g}_1(x)\big)=\hat{g}_0(x)$ for all $x\in\mathcal{F}$. This motivates the use of the characterization in ((ref)) which is remains well-defined even if we replace conditional average potential outcomes on the frontier with noisy estimates, and even if comonotoncity fails.} In brief, we take the untreated individuals who lie within a small shrinking neighbourhood of the frontier and we regress their outcomes on estimates of $E[Y(1)|X]$. To assess whether an untreated individual $i$ is within a small neighbourhood of the frontier, we check whether $i$'s covariates $X_i$ are sufficiently close to those of $i$'s nearest treated neighbour. Let $NN(x)$ be the index of this nearest neighbour for $x \in \mathcal{X}_{0}$ (in the event of ties, one can choose among them at random), formally \[ NN(x):=\arg\min_{j\in\mathcal{I}_1}\|X_j-x\|, \] and $NN_{i} := NN(X_{i})$. If $\|X_i-X_{NN_{i}}\|\leq d$, and treated and untreated regions satisfy some regularity conditions, it follows that $X_i$ is within distance $d$ of a point on the frontier. This approach allows researchers to apply our method in settings in which treatment statuses are observed but the location of the frontier itself is unknown. It also obviates the need for applied practitioners to calculate distances of points to the frontier which could complicate practical implementation.

For untreated individuals within a shrinking distance from the frontier, the conditional mean treated potential outcome can be consistently estimated. A crude estimate is $\hat{g}_1(X_{NN_{i}})$, which is an estimate of the conditional mean outcome of the nearest treated neighbour. However, if the conditional mean potential outcome is twice differentiable, one can reduce bias by estimating the derivative of the conditional mean outcome at $X_{NN_{i}}$. Let $\triangledown\hat{g}_1(x)$ be an estimate of the vector of partial derivatives $\frac{\partial}{\partial x} E[Y|X=x]$, then for $x \in \mathcal{X}_{0}$, we estimate $E[Y(1)|X=x]$ by \[\tilde{g}_1(x)=\hat{g}_1(X_{NN(x)})+\triangledown\hat{g}_1(X_{NN(x)})'(X_i-X_{NN(x)}).\] In the case of local linear regression, the local linear regression coefficients $\hat{\beta}_{X_{NN(x)}}$ (absent the intercept) are a natural choice for $\triangledown\hat{g}_1(X_{NN(x)})$. The resulting estimator can be written succinctly as follows:

align*[align* omitted — 234 chars of source]

As described above, we estimate $q_{0}(y)$ by nonparametric regression of $Y_i$ on $\tilde{g}_1(X_i)$ using only treated individuals $i$ for whom $\|X_i-X_{NN_{i}}\|$ is sufficiently small. Again we suggest local linear regression. Let $b$ be a bandwidth and $W_i=1\{\|X_i-X_{NN_{i}}\|\leq \varepsilon\}$. The estimate is given below.

align[align omitted — 256 chars of source]

The bandwidth $b$ may be chosen by cross-validation. One could replace $Y_i$ in the above with $\hat{g}_0(X_i)$, which may increase precision but at the expense of greater bias. The weighting by the binary indicator $W_i$ ensures that the local linear regression only includes data on individuals sufficiently close to the frontier. As a rule of thumb, we suggest setting $d=\omega h$, where $\omega$ is twice the $0.75$ quantile of the kernel $K$ and $h$ is the bandwidth used in the estimate $\hat{g}_1$. This choice is designed to ensure that the bias from using points not exactly at the frontier is of a similar magnitude to the bias in the estimation of $\hat{g}_1$.

Recall that under regularity conditions, the domain of $q_0$ is an interval with end points $\underline{y}_0$ and $\bar{y}_0$ defined in ((ref)). We estimate these end points by $\hat{\underline{y}}_0$ and $\hat{\bar{y}}_0$ defined below. \[ \hat{\underline{y}}_0=\min_{i:i\in\mathcal{I}_1,W_i=1}\tilde{g}_{1}(X_{i}),\,\,\,\,\,\hat{\bar{y}}_0=\max_{i: i \in \mathcal{I}_{1},W_i=1}\tilde{g}_{1}(X_{i}) \]

Figure 3.1 provides a visual demonstration of the method in which the non-parametric regressions are carried out with a uniform kernel. Figure 3.1.a contains simulated data from the expository model in Section 1. The crosses represent the covariate values of different data points. Let $i$ be the index of the untreated individual with covariate values indicated by the large cross in Figure 3.1.a. The covariate vales of $i$'s nearest treated neighbour $X_{NN_{i}}$, are indicated by the large plus sign. The large black circle is centered at $X_{NN_{i}}$ and has radius $h$ (which was here chosen arbitrarily to equal $0.2$). Those data points in the treated region within distance $h$ of $X_{NN_{i}}$ are indicated by the small circles. The estimate $\tilde{g}_0(X_i)$ is the fitted value at $X_i$ from linear regression of the outcomes on the covariates using only those data points whose covariates are indicated by small circles. The distance between $X_{i}$ and $X_{NN_{i}}$ is the length of the black dashed line between these two points. Because we use a uniform kernel $\omega=1$, and thus the suggested value of $d$ described above is simply $h$. Therefore, $W_i=1$ because $\|X_{i}-X_{NN_{i}}\|\leq h$ as can be seen from the fact $X_i$ lies within the black circle.

figure[figure omitted — 197 chars of source]

Using the same simulated data, Figure 3.1.b plots values of $Y_i$ for untreated individuals with $W_i=1$ against corresponding values of $\tilde{g}_1(X_i)$ which is evaluated as in the previous paragraph. Regressing these outcomes on the corresponding values of $\tilde{g}_1(X_i)$ by local linear regression with a uniform kernel yields the solid black curve which is the estimate of $q_0$ (here we simply use the same ad hoc bandwidth of $0.2$). The true value of $q_0$ is indicated by the dashed blue curve. The values of $\hat{\underline{y}}_0$ and $\hat{\bar{y}}_0$ are indicated by the vertical dashed lines.

Effects of Counterfactual Policies

Given estimates of conditional average potential outcomes, one can estimate the causal effects of counterfactual treatment regimes. Consider a policy in which an individual with $X=x$ is treated with probability $p(x)$. Note that if $p(x)\in\{0,1\}$ for all $x$, then treatment remains deterministic under this regime. The mean causal impact of this counterfactual policy on the outcomes of individuals with values of $X$ in a set $\mathcal{S}$ denoted $\theta_0$, is given below.

align[align omitted — 191 chars of source]

The restriction $X\in\mathcal{S}$ allows us to restrict attention to those values of $x$ at which both conditional average potential outcomes are identified. Recall that under regularity conditions on the frontier, for $d=0,1$ we identify $E[Y(0)|X=x]$ and $E[Y(1)|X=x]$ for each $x$ in $\mathcal{X}_{d}$ such that that $\underline{y}_{1-d}<E[Y|X=x]<\bar{y}_{1-d}$. As such, we focus on estimation for $\mathcal{S}$ the set that contains these values of $x$. In order to estimate ((ref)) given this value of $\mathcal{S}$ we use the plug-in principle. The resulting estimate $\hat{\theta}$ is given by

align[align omitted — 314 chars of source]

where $S_i$ is a binary indicator defined below. \[S_i=(1-D_i)1\{\hat{\underline{y}}_{1}\leq\hat{g}_{0}(X_{i})\leq\hat{\bar{y}}_{1}\}+D_i 1\{\hat{\underline{y}}_{0}\leq\hat{g}_{1}(X_{i})\leq\hat{\bar{y}}_{0}\}\]

Figure 3.2 plots the same simulated data as in Figure 3.1. In Figure 3.2.a, the covariate values of individuals with $S_i=1$ are indicated by crosses, whereas those for whom $S_i=0$ (i.e., those for whom conditional average causal effects are not identified) are indicated by circles. The dashed blue line in the figure indicates a counterfactual treatment rule, individuals are treated if and only if their covariate values are to the right of this dashed line (That is, if their math score exceeds $0.5$). Note that all individuals whose treatment under this rule differs from their factual treatment (those with covariate values in the orange and deep blue triangles) have $S_i=1$, and thus the average treatment effect of this change is fully identified.

Using the simulated data, Figure 3.2.b plots estimated mean causal effects of a counterfactual treatment rule in which individuals are treated if and only if their math score is greater than some cut-off. The x-axis gives the value of this cut-off and the y-value of the solid blue curve is the corresponding causal effect estimate calculated from the formula ((ref)).

figure[figure omitted — 219 chars of source]

It may be helpful to compare $\hat{\theta}$ to the expected change in treatment status under the counterfactual policy among those individuals in the sample with $S_i=1$. That is, to compare it with $\frac{1}{\sum_{i=1}^{n}S_{i}} \sum_{i=1}^{n}S_{i}[(1-D_{i})p(X_{i})-D_{i}(1-p(X_{i}))]$. If the costs of a policy are proportional to the number of individuals treated, then this quantity is proportional to the expected cost of the counterfactual treatment policy on the subsample with $S_i=1$ less the cost of the factual treatment assignments. This quantity is plotted by the dashed red curve in Figure 3.2.b with units given by the y-axis on the right of the figure.

Conditional Comonotonicity

In settings with more than two covariates in $X_i$ we suggest researchers estimate effects using the weaker conditional comonotonicity condition in Section 2.2. In these settings $X_i$ can be decomposed into two subvectors $X_i^{(1)}$ and $X_i^{(2)}$. If $X_i^{(2)}$ is discrete (with a finite number of support points), then one can apply the estimation method described above separately within the different strata of $X_i^{(2)}$. More generally, if $X_i^{(2)}$ is continuously distributed, one needs only slightly adjust the method in the previous section.

In particular, $\tilde{g}_1$ is estimated exactly as in the previous subsection. However, in place of $q_{0}(y)$ we estimate a value of this quantity within each stratum of $X_i^{(2)}$. That is, for each $x^{(2)}$ in the support of $X_i^{(2)}$ we estimate a stratum-specific value of $q_{0}(y)$ by $\hat{q}_{0}(y,x^{(2)}):=(1,y,{x^{(2)}}')\hat{\gamma}_{y,x^{(2)}}$, where $\hat{\gamma}_{y,x^{(2)}}$ is given by

align*[align* omitted — 211 chars of source]

That is, $\hat{q}_{0}(y,x^{(2)})$ is a fitted value from regressing outcomes $Y_i$ for treated individuals near the frontier on both an estimate of $E[Y(1)|X=x]$ and also $X^{(2)}$. The corresponding estimated conditional average treated and untreated potential outcome estimates at $x=({x^{(1)}}',{x^{(2)}}')'$ are then

align*[align* omitted — 128 chars of source]

Asymptotic Analysis

We now provide asymptotic properties of the estimation procedure introduced in Section (ref). We refer to the estimation of $g_{d}$ as the first stage and estimation of $q_{1-d}$ as the second stage. We first show a (strong) uniform convergence result for the first stage, which is rather standard in the literature of local polynomial estimation (e.g., stone1982optimal, masry1996multivariate, hansen2008uniform). A minor deviation is that we consider a local extrapolation to estimate $g_{d}$ for observations are “slightly outside” of $\mathcal{X}_{d}$.

Then, we establish asymptotic normality of $\hat{q}_{1-d}(y)$. As evident from its definition given in (ref), the estimator $\hat{q}_{1-d}(y)$ is obtained by a local linear regression on regressors that are also generated by a local linear regression. Hence, the estimation method falls within the general framework of mammen2012nonparametric, whose arguments we closely follow, with two key differences. First, we incorporate weights $W_i$ to ensure that only observations near the boundary are used in estimating $g_d$ (for units in $\mathcal{X}_{1-d}$). These weights do not affect the convergence rate of the first stage but do enter the variance of the second stage, introducing a trade-off absent in MRS. Second, the first stage estimation error is independent of the data used for the second stage since the estimation of $\tilde{g}_{d}(x)$ only uses observations in $\mathcal{X}_{d}$, whereas $\hat{q}_{1-d}(y)$ only uses observations in $\mathcal{X}_{1-d}$.

We start with showing strong uniform consistency of $\tilde{g}_{d}(x)$ over $$\mathcal{X}_{d, \varepsilon} := \mathcal{X}_{1-d} \cap \{x: \min_{X_{i} \in \mathcal{X}_{d}}\lVert x - X_{i} \rVert \leq \varepsilon \},$$ which is the set of $x \in \mathcal{X}_{1-d}$ such that there exists at least one observation in $ \mathcal{X}_{d}$ that lies within an $\varepsilon$-ball centered at $x$. By the triangle inequality and mean-value theorem, we have

align*[align* omitted — 202 chars of source]

where $x_{t} = tx + (1-t)NN(x)$ for some $t \in [0, 1]$. The convergence of the second term in the last line is standard. By Cauchy-Schwarz and another application of the mean-value theorem, the first term can be bounded uniformly over $\mathcal{X}_{d,\varepsilon}$ by $$\norm{\hat{\beta}_{NN(x)} - \nabla g_{d}(NN(x))}\varepsilon + M\varepsilon^{2},$$ which shows that the convergence rate is essentially driven by the convergence rate by $\norm{\hat{\beta}_{NN(x)} - \nabla g_{d}(NN(x))},$ which is again standard in the local polynomial literature.

We assume the following to establish the said uniform convergence result. \setcounter{assumption}{3}

assumption[Support and Treatment Region] The support $\mathcal{X}$ of $X$ is compact. Assume that for $d \in \{0,1 \}$, $\pi_{d} :=P(X \in \mathcal{X}_{d}) \in [\delta_{\pi}, 1-\delta_{\pi}]$ for some $\delta_{\pi} > 0.$

The compactness of $\mathcal{X}$ guarantees that the closures of both $\mathcal{X}_{1}$ and $\mathcal{X}_{0}$ are compact. Such compactness conditions are used to derive uniform convergence results for $\tilde{g}_{d}(\cdot)$. These can be relaxed by imposing tail conditions on $f_{X}(\cdot)$ as in hansen2008uniform.

We make the following smoothness and boundedness conditions on the kernel, density of $X$, and conditional expectation functions $g_{d}(\cdot)$, which are typical in the nonparametric literature.

assumption[Smoothness and boundedness] (i) The kernel $K$ is a symmetric density function of order 2, Lipschitz continuous and has bounded support, (ii) $E[Y^{2+\delta_{y}} | X] < \infty$ for some $\delta_{y} > 0$, (iii) the density $f_{X}(\cdot)$ is twice differentiable and satisfies $0 < \inf_{x \in \mathcal{X}} f_{X}(x) $ and (iv) $g_{d}$ is twice differentiable with Lipschitz continuous second derivatives.

Define $a_{n} = O\left( \left(\frac{\log n}{n h^{k}} \right)^{1/2} + h^{2}\right)$, which is the usual uniform convergence rate of local linear estimators. The following shows that the convergence rate of the local linear estimator is unaffected under our local extrapolation. The result is intuitive given the boundary properties of the local linear estimator and that points that lie within $h$ distance from the support are essentially the same with being on the boundary.

\setcounter{thm}{3}

thm[Strong Uniform Consistency of $\tilde{g}_{d}(\cdot)$] Suppose Assumptions (ref) and (ref) hold, and set $\varepsilon \asymp h$. Then, \begin{equation} \sup_{x \in \mathcal{X}_{d,\varepsilon}} \abs{\tilde{g}_{d}(x) - g_{d}(x)} = O(a_{n}), a.s. \end{equation}
proofThe result follows immediately by the decomposition and bounds provided above, and then applying Theorem 6 of masry1996multivariate to our i.i.d. setting.

We now move on to establishing asymptotic results for $\hat{q}_{1-d}(y)$ for $y \in \mathcal{Y}:= \{g_{d}(x): x \in \mathcal{F}\}$. We first derive an asymptotic normality result for $q^{\ast}_{1-d}$, which is the infeasible estimator that uses $g_{1-d}$ instead of the estimated $\tilde{g}_{d}$.\footnote{Formally, ${q}^{\ast}_{0}(y)=(1,y){\gamma}^{\ast}_{y},$ where ${\gamma}^{\ast}_{y}=\arg\min_{\gamma}\sum_{i \in \mathcal{I}_{_{1-d}}}W_{i}K\left(\frac{|{g}_{d}(X_{i})-y|}{b}\right)[Y_{i}-\big(1,{g}_{d}(X_{i})\big)\gamma]^{2}$. } Then, we follow MRS and first derive a stochastic expansion of $\hat{q}_{1-d}(y)$ around $q^{\ast}_{1-d}(y)$. Combining these results we provide conditions under which $\hat{q}_{1-d}(y)$ is asymptotically normal.

Due to the fact that $\tilde{g}_{d}$ is independent of $(X_{i})_{i\in \mathcal{I}_{1-d}}$, the more involved assumptions of MRS (such as Assumptions 3 and 4) are unnecessary. Since we rely on local linear estimation, we need a smoothness assumption on $q_{1-d}$. We assume the following.

assumption[Smoothness of $q_{1-d}$] The function $q_{1-d}(\cdot)$ is twice differentiable with uniformly bounded second derivatives.

Due to the “local to the boundary” nature of the estimation procedure, we impose some regularity condition on the boundary. These assumptions are not restrictive and satisfied in most empirical contexts. Let $v_{\varepsilon}$ denote the volume of an $\varepsilon$-ball in $\mathbb{R}^{k}$.

assumption[Regularity of $\mathcal{F}$] (i) $\mathcal{F}$ is piecewise $C^{2}$ and continuous. (ii) For each $x \in \mathcal{F}$ and $d \in \{0, 1\}$, $\mathrm{vol}\, (\mathcal{X}_{d} \cap B_{\varepsilon}(x))/v_{\varepsilon} \to \kappa_{d}(x) $ as $\varepsilon \to 0$ where $\kappa_{d}(x) \in [\delta, 1-\delta]$ for some $\delta > 0$.
assumption[Locally nonvanishing $\nabla g_{d}(\cdot)$] The conditional expectation function $g_{d}$ satisfies $\inf_{x \in B_{\delta_{g}}(y)} \norm{\nabla g_{d}(x)} > 0 $ for some $\delta_{g} > 0$.

Assumption (ref)(i) and (ref) are used to perform a change of variables to calculate the density of $g_{d}(X)$ local to the boundary. They can be relaxed at the cost of a longer proof. Assumption (ref)(ii) ensures that there is enough volume in either side of the boundary, which ensures that the probabilty that $X$ lies on either side of the boundary scales like $\varepsilon^{d}$ near the boundary. These assumptions ensure that $P(W_{i}=1|X_{i}=1)$ scales linearly in $\varepsilon$.

Define $f_{g}(u) = v_{1} \int_{g^{-1}_{d}(u) \cap \mathcal{F}} \kappa_{d}(z) \frac{f_{X}(z)}{\lVert\nabla g_d(z)\rVert} d\mathcal H^{k-2}(z), $\footnote{$\mathcal{H}^{k-2}$ denotes the $k-2$ dimensional Hausdorff measure.} which is effectively the density of $g_{d}(X_{i})$ that accounts for the weight $W_{i}$ at a boundary point, the conditional variance $\sigma^{2}(y) := E[(Y_{i}-E[Y_{i}|g_{x}(X_{i})])^{2}| g_{d}(X_{i}) = y]$ and second moment of the kernel $\mu_{2} = \int u^{2}K(u)du$. The following theorem establishes the asymptotic normality of $q^{\ast}_{1-d}(y)$.

thm[Asymptotic normality of $q^*_{1-d}(y)$] Suppose Assumptions (ref)-(ref) hold, $a_{n}/b \to 0$, $n\varepsilon^{d} \to \infty$, and $n \varepsilon b \to \infty$. Then, we have \begin{equation*} (n \varepsilon b)^{1/2} \left({q}^{\ast}_{1-d}(y) - q_{1-d}(y) - \frac{b^{2}}{2} q”(y)\mu_{2}\right) \to N\left(0,\frac{ \sigma^2(y) R(K)}{f_g(y)}\right). \end{equation*}

The asymptotic distribution is essentially the same as what one would expect from a local regression of $Y_{i}$ on $g_{d}(X_{i})$. The difference is that the weights affect the variance through the normalizing factor $(n\varepsilon b)^{1/2}$ and density $f_{g}(y)$. However, the weights do not affect the bias term.

We now analyze how close $q^{\ast}_{1-d}(y)$ is to $\hat{q}_{1-d}(y)$, and derive conditions under which the difference is negligible. Similar to MRS, define

equation*[equation* omitted — 151 chars of source]

To derive uniform results, we need tail conditions on the regression error $\eta_{i} := Y_{i}-g_{d}(X_{i})$. We impose a sub-exponential moment condition on the error term as follows, which is the same with Assumption 1(iv) of MRS.

assumption[Tail behavior of $\eta_i$] The regression error $\eta_{i}$ satisfies $E[\exp(\ell |\eta_{i}|) \mid g_{d}(X_{i})] \leq C$ for some $\ell, C > 0$.

The following expansion, which is a minor modification of that by mammen2012nonparametric, characterizes the difference between the oracle and feasible estimators.

thm[Expansion of $\hat{q}_{1-d}(y) - {q}^{*}_{1-d}(y)$] Suppose Assumptions (ref)-(ref) hold. Then, we have \begin{align*} \sup_{y} \big| \hat{q}_{1-d}(y) - {q}^{*}_{1-d}(y) + q_{1-d}'(y) \hat{\Delta}(y) \big| &= O_P(c_{n}), \end{align*} where $c_{n} =(\log n\varepsilon)^{1/2}((n\varepsilon)^{-1/2}a_{n}b^{-3/2} \vee a_{n}b \vee a_{n}^{2}/b))$.

Compared with MRS, the main differences are the presence of the weights $W_{i}$ and the fact that the estimation of the first stage is independent of the second stage data. The former requires minor modifications to the proof and latter in fact makes the proof easier. This difference is also reflected in the convergence rate; in MRS, the first component of the rate $((n\varepsilon)^{-1/2}a_{n}b^{-3/2}$ includes an additional correction term to account for the possible correlation between the two stages, which leads to slower convergence rates.

Also, under possible violations of comonotonicity, the expansion includes an additional bias term $\hat{\Gamma}(y) := E[ (W_{i}/\varepsilon)\,K_b'({g}_{d}(X_{i}) - y) (\tilde{g}_{d}(X_i) - g_{d}(X_i) )\xi_{i}]$. This is because $\xi_{i}:=E[Y_{i}- q_{1-d}(g_{d}(X_{i}))|X_{i}] \neq 0 $ in general when comonotonicity is violated. However, we note that dropping the assumption that $\xi_{i} = 0$, does not affect the convergence rate above or the following corollary.

If $(n \varepsilon b)^{1/2}c_{n} = o(1)$, the first stage estimation is negligible so that the asymptotic distribution derived in Theorem (ref) holds with $\hat{q}_{1-d}(y)$ in place of $q^{\ast}_{1-d}(y)$ as well. The following corollary shows when this is possible under $k=2$.

corSuppose Assumptions (ref)-(ref) hold, and let $k=2$, $\varepsilon \asymp h \asymp n^{-r_{h}}$ and $b \asymp n^{-r_{b}}$. Then, if $r_h \in (1/8, 4/15)$ and $r_{b}\in (\{\max\{ (3/4)r_h,\, {(1 - 5r_h)}/{3} \}, \min\{ 1/2 - r_h,\, 2r_h,\, 1 - 3r_h,\, 9r_h - 1 \})$, then \begin{equation*} (n \varepsilon b)^{1/2} \left(\hat{q}_{1-d}(y) - q_{1-d}(y) - \frac{b^{2}}{2} q”(y)\mu_{2}\right) \to N\left(0,\frac{ \sigma^2(y) R(K)}{f_g(y)}\right). \end{equation*} Moreover, setting $(1 - r_{h})/5 < r_{b}$ corresponding to the undersmoothing regime where the $O(b^{2})$ bias term can be dropped.

For example, when the mean squared error optimal bandwidth of $r_{h} = 1/6$ is chosen for the first stage, the bandwidth choice of $r_{b} \in (1/6, 1/3)$ is feasible for the second stage. That is, one can choose $r_{b} = 1/5$ yet “undersmooth.” This seems counterintuitive but in fact it is due to the fact that the first stage bandwidth enters the second stage variance that reduces the effect of the bias in the second stage.

We leave the development of a full inferential theory under minimal conditions for future work, including the problem of testing comonotonicity (see the discussion following Theorem 2). However, the results of this section show that for appropriate bandwidth choices, the estimate $\hat{q}_{1-d}$ is asymptotically equivalent to an oracle estimator which has the form of a standard local linear regression estimate. It is shown in, for example, chernozhukov2014gaussian, that for such estimates multiplier bootstrap confidence bands are valid under appropriate assumptions, which we implement in our empirical application. Moreover, extending the robust bias correction method by calonico2014robust to our setting with generated regressors seems to be a promising direction for future research.

Application: The Impact of Mandatory Summer School

We apply the methods to the empirical setting of Matsudaira2008 who investigate the impact of mandatory summer schooling on future test scores in reading and mathematics. The authors use data on 5th-grade students in a large school district in the northeastern US in 2001. Students in this district scored below a threshold on year-end reading or math tests were required to either repeat a grade or attend mandatory summer schooling. The primary outcomes considered are reading and math test scores one year after the initial tests. It is important to note that some students attended summer school regardless, leading to an intent-to-treat interpretation of the results. The test scores that determine mandatory summer school are observed by the researcher, and so conditional on these observed covariates, treatment is deterministic.

Following Imbens2019, we include only students with scores 40 points above or below the cutoffs, leaving a sample size of $n = 30,741$. The distribution of test scores in this restricted sample is presented in Figure (ref). The blue curve is the frontier: students whose scores were below either cutoff were faced with the option of summer school or grade repetition.

figure[figure omitted — 164 chars of source]

Below we provide the results from regressing the outcome on covariates using only observations in either the treated or untreated.

figure[figure omitted — 523 chars of source]

Figure (ref) plots quadratic regression results for the math scores. The contours in the magenta untreated region correspond to regression using only untreated individuals while the contours in the turquoise treated region from regression only on treated individuals. Figure (ref) presents results from the same exercise using reading scores. In both cases, the contours in the two regions appear approximately aligned, consistent with the hypothesis of comonotonicity. Figures (ref) and (ref) plot results from the same exercise but with local linear regression instead of quadratic and gives similar results to quadratic.

In this setting, we may justify the comonotonicity condition by considering the underlying skills in reading and math measured by the test scores. In Appendix B.1 we specify a formal model of skill formation in which test scores are noisy measures of underlying skills and show that the model implies comonotonicity. Here we provide a verbal summary. Suppose students begin with some levels of skill in reading and math and the initial tests are noisy measures of these skills. If reading and math skills are positively correlated in the population, then both reading and math test scores are informative about the reading skill, even if only the reading test directly measures this quantity. Likewise, both reading and math scores are informative about math skills.

The outcomes, reading and math scores a year after initial tests, may be understood to measure reading and math skills at this later date. The reading outcome would thus reflect reading skill at this date which may be depend on reading skill at the time of the initial test, the (possibly heterogeneous) impact of treatment, and exogenous day-of-test noise. One may then expect that the reading outcome is positively correlated with both initial reading and math scores due to the mutual positive association with initial reading skill, and similarly for the math outcome. Indeed, in Figures (ref) and (ref) we see that this is the case. A higher initial math score is associated with a higher reading outcome, but this association is weaker than that between initial reading scores and the reading outcome, which presumably reflects that the initial reading tests directly measure initial reading skill and thus have a higher association with this quantity than the math scores. A similar pattern is evident for the math outcome.

Now suppose some initial test scores $x$ in the treated region are associated with higher average reading outcomes than some scores $x'$. Suppose we accept the premises above, then this observation suggests individuals with scores $x$ tend to have higher initial reading skills than those with scores $x'$. Both the treated and untreated potential outcomes are increasing in initial reading skills. As such, we may expect individuals with scores $x$ to also have higher average reading outcomes than those with $x'$ under a counterfactual in which none are treated. Indeed we show formally that comonotonicity applies under the skill formation model in Appendix B.1.

figure[figure omitted — 754 chars of source]

Figures (ref) and (ref) present estimates of $\hat{q}_0$ and $\hat{q}_1$ respectively for math scores using the methods in Section 3. $\hat{q}_0$ and $\hat{q}_1$ are respectively below and above the 45 degree line for all but the largest values of the conditional average potential outcomes. Thus the estimates provide evidence of a positive treatment effect for individuals with low and moderate conditional average baseline potential outcomes. More generally the slope of $\hat{q}_0$ ($\hat{q}_1$) is greater (lower) than that of the 45 degree line, which suggests the treatment effect is lower on average for individuals whose covariate values are associated with greater untreated potential outcomes. That is, that individuals whose scores would be worse without treatment, tend to benefit more from treatment. In all figures, green curves are upper and lower $90\%$ pointwise confidence bands evaluated using the multiplier bootstrap.\footnote{To calculate the multiplier boostrap intervals, for each bootstrap simulation $s$ we draw $n$ independent standard exponential random variables. In each of the kernel sums in the first and second stages of the estimator $\hat{q}_d$ we multiply the $i$-th term by the corresponding exponential random variable to obtain a bootstrap estimate $\hat{q}_{d,s}$. The band at $y$ is then $\hat{q}_{d}(y)$ plus and minus the $90$-th percentile of $|\hat{q}_{d,s}(y)-\hat{q}_{d}(y)|$ over the bootstrap draws. We use $100$ bootstrap draws.}

Figure (ref) compares these estimates. Under comonotonicity the population versions of these functions should align, and indeed in (ref) these appear close.

Figures (ref) and (ref) present estimates of $\hat{q}_0$ and $\hat{q}_1$ for reading scores, respectively. In this case, conditional average treatment effect estimates are positive for all individuals and increasing in conditional average baseline potential outcomes. Figure (ref) compares these estimates which are closely aligned.

figure[figure omitted — 277 chars of source]

Figures (ref) and (ref) respectively display the covariate values at which we are able to extrapolate conditional average treatment effects. Crosses represent values at which identification is achieved under comonotonicity and circles indicate points to which we cannot extrapolate. In terms of the notation introduced in Section 3, points marked with crosses represent values of $X_i$ for which we estimate $S_i=1$ and those with circles, values for which $S_i=0$.

Figures (ref) and (ref) indicate which counterfactual objects one can identify using comonotonicity in this setting. For example, Figure (ref) suggests that we can identify the average effect on math outcomes of counterfactual treatment policies in which the math score threshold for mandatory summer school is moderately increased. However, the average effect of a moderate decrease in this cut-off is not identified. Similarly, Figure (ref) suggests that the average counterfactual effect on reading outcomes of a moderate increase in the reading score is identified.

figure[figure omitted — 303 chars of source]

We estimate average causal effects on math scores from counterfactual increases in the math cut-off and the effects on reading outcomes of increases in the reading threshold. The methods in Section 3 allow us to estimate the average effects among those individuals for whom the CATE is identified, i.e., for whom $S_i=1$. From Figures (ref) and (ref) we see that for all but very large increases in the threshold, every individual impacted by the increase in the cut-off has $S_i=1$ and so for all but the very largest increases in the cut-off these conditional (on $S_i=1$) effects are exactly the unconditional average effects. The solid blue curves in (ref) and (ref) plot these causal effects for math and reading respectively. The dashed red curves show the number of individuals in the sample with $S_i=1$ whose treatment status is impacted by the counterfactual policy.

Figure (ref) shows that as the math threshold increases, the causal effect increases initially but the rate at a rate that slows with the level of the threshold. When the threshold is raised to $0.75$ increases in the threshold lead to a slight decrease in the effect which then levels out. The number of individuals with $S_i=1$ impacted by the policy increases roughly linearly with the threshold until the threshold reaches around $0.9$ before it levels off, which may reflect the lack of identification of causal effects for individuals with very high initial math scores seen in (ref). The slow down and eventual reversal of the increase in the causal effect up to threshold level $0.9$ can be understood from Figures (ref) and (ref). When the threshold is high, any marginal increase impacts only individuals with high initial math scores and from (ref) we see these individuals have high average untreated potential outcomes. (ref) suggests that this is associated with a smaller CATE, and thus a smaller increase in the causal effect of the policy.

By contrast, (ref) shows a steadily increasing relationship between the CATE and conditional average untreated potential outcome and the causal effect shown in (ref) increases at a similar rate to the number of individuals impacted by the policy. Note that in both cases, the multiplier bootstrap confidence bands exclude zero average impact of some of the included counterfactual policy rules.