Treatment Effect Heterogeneity in Regression Discontinuity Designs

Online Supplement for “Treatment Effect Heterogeneity in Regression Discontinuity Designs”

\onehalfspacing

\singlespacing

Introduction

This section introduces the notation and links it with the one used in the main text, sets up the statistical framework, and enumerates the assumptions we rely on throughout this supplement.

Notation

In this Supplemental Appendix, we use $n$ to denote sample size and $h\equiv h_n$ to denote a bandwidth sequence where we omit the dependence on $n$ to ease notation. Moreover, $h$ denotes a generic bandwidth (e.g. $h=h_{-}$ or $h=h_{+}$ depending on the context).

Linear algebra. Throughout the text, $\mathbf{e}_{\nu}$ denotes a conformable vector of zeros with a 1 in its $(\nu+1)$-th element, which may take different dimensions in different places, $\mathbf{0}_k$ and $\boldsymbol{\iota}_k$ are the $k$-dimensional zero and one vectors, respectively, $\mathbf{I}_k$ and $\mathbf{0}_{k\times j}$ denote the $k\times k$ identity matrix and a $k\times j$ matrix of zeros, respectively, $\otimes$ indicates the Kronecker product, $\mathrm{tr}(\cdot)$ is the trace operator, and $\mathrm{diag}(\mathbf{x})$ yields a square diagonal matrix with the elements of $\mathbf{x}$ on its main diagonal. With a slight abuse of notation, we denote with $\mathbf{v}^k = (v_1^k,\cdots,v_n^k)^\prime$ the element-wise power for vectors $\mathbf{v}\in \operatorname{\mathbb{R}}^n$. The maximum and minimum of two real numbers $a$ and $b$ are denoted by $a \vee b$ and $a \wedge b$, respectively. We let $|\cdot|$ denote the Euclidean norm, $|\mathbf{A}|^2 = {\sum_i\sum_j |a_{ij}|^2} = {\mathrm{tr}(\mathbf{A}^\prime\mathbf{A})}$. Finally, for some $q\in\operatorname{\mathbb{N}}$, with $\mathcal{C}^q$ we denote the space of functions that are $q$-times continuously differentiable.

Asymptotic statements. For two positive sequences $\left\{a_n\right\}_n,\left\{b_n\right\}_n$, we write $a_n = O(b_n)$ if $\exists\,M\in\operatorname{\mathbb{R}}_{++}:a_n\leq M b_n$ for all large $n$, $a_n=o(b_n)$ if $\lim_{n\to\infty} a_n b_n^{-1} = 0$, and $a_n\lesssim b_n$ if there exists a constant $C$ such that $a_n\leq C b_n$ for all large $n$. For two sequences of random variables $\left\{A_n\right\}_n,\left\{B_n\right\}_n$, we write $A_n=o_{\mathbb{P}}(B_n)$ if $\forall\,\epsilon>0, \lim_{n\to\infty} \mathbb{P}[|A_n B_n^{-1}|\geq \varepsilon]= 0$ and $A_n =O_{\mathbb{P}}(B_n)$ if $\forall\,\varepsilon>0, \exists \,M,n_0\in\operatorname{\mathbb{R}}_{++} : \mathbb{P}[|A_nB_n^{-1}|>M]<\varepsilon,$ for $n>n_0$. We denote convergence in probability with $\overset{\mathbb{P}}{\longrightarrow}$ and convergence in distribution with $\rightsquigarrow$. We denote (possibly multivariate) Gaussian random variable with $\mathsf{N}(\mathbf{a}, \mathbf{B}),$ where $\mathbf{a}$ denotes the mean and $\mathbf{B}$ the variance-covariance.

Causal Model. We now describe the population causal model. The outcome variable is

with $\left(Y_i(0), Y_i(1)\right)\in\operatorname{\mathbb{R}}^2$ denoting the potential outcomes and $T_i:=\mathbbm{1}(X_i\geq c)$ denoting treatment status. We denote the vector of covariates with $\mathbf{W}_i \in \operatorname{\mathbb{R}}^d, d\in\operatorname{\mathbb{N}}$. We stress that, as the notation suggests, $\mathbf{W}_i$ is interpreted as a pretreatment vector of covariates in the sense that $\mathbf{W}_i = \mathbf{W}_i(1)=\mathbf{W}_i(0)$.

In sharp RD designs, $T_i=\mathbbm{1}\left(X_i \geq c\right)$, where $X_i\in\operatorname{\mathbb{R}}$ denotes the running variable and $c\in\operatorname{\mathbb{R}}$ is the cutoff. Throughout, $F_R(\cdot)$ denotes the cumulative distribution function (cdf) of a random variable $R$ and $f_R(\cdot)$ is the density of $R$ with respect to the Lebesgue measure. We denote with $f(\cdot)$ the density of $X_i$.

We further define

and $\boldsymbol{\sigma}^2_{W}(x) =\boldsymbol{\mu}_{WW}(x)-\boldsymbol{\mu}_{W}(x)\boldsymbol{\mu}_{W}(x)^\prime$.

Let $\mathbf{r}_q(x)=(1,x,\ldots, x^q)^\prime$ be the polynomial basis of order $q\in\operatorname{\mathbb{N}}$ and let \[\mathbf{r}_{p,s}(u, \mathbf{w}) = (\mathbf{r}_p(u)^\prime,\mathbf{w}^\prime\otimes \mathbf{r}_s(u)^\prime)^\prime.\] For a generic kernel function $k(\cdot)$ let \[K(u)=\mathbbm{1}(u<0)k(-u) + \mathbbm{1}(u\geq 0)k(u)\] and \[K_{\mathbf{h}}(u) = \mathbbm{1}(u<0)k_{h_{-}}(-u) + \mathbbm{1}(u\geq 0)k_{h_{+}}(u),\quad k_h(u) = k(u/h)/h, \quad \mathbf{h} = (h_{-},h_{+})^\prime. \] Finally, we also define the following matrices \[ {\boldsymbol{\Lambda}}_{-,p,s} := f(c)\int\limits_{-\infty}^0 K(u) \mathbf{r}_p(u)\mathbf{r}_s(u)^\prime \,\mathrm{d} u ,\qquad {\boldsymbol{\Lambda}}_{+,p,s} := f(c)\int\limits_0^\infty K(u) \mathbf{r}_p(u)\mathbf{r}_s(u)^\prime \,\mathrm{d} u, \] and \[\boldsymbol{\Xi}_{-,p,s} = f(c)\int\limits_{-\infty}^0 K^2(u)\mathbf{r}_p\left(u\right)\mathbf{r}_s\left(u\right)^\prime \,\mathrm{d} u, \qquad \boldsymbol{\Xi}_{+,p,s} =f(c)\int\limits_0^\infty K^2(u)\mathbf{r}_p\left(u\right)\mathbf{r}_s\left(u\right)^\prime \,\mathrm{d} u. \]

Setup

Formally, the RD estimators with treatment interactions can be obtained from a “long” regression, i.e.,

where

and

Note that the vector $\widehat{\boldsymbol{\lambda}}_{+,s}(h)$ (similarly, $\widehat{\boldsymbol{\lambda}}_{-,s}(h)$) stacks first all the coefficients of the polynomial basis of order $s$ for the first covariate, then for the second, and so on until the $d$-th covariate.

In what follows, we define $m:=d(1+s)$. The regression above has orthogonal regressors. Indeed, it has a set of regressors that pertains to estimation to the left of the cutoff, i.e., $(\mathbbm{1} (X_i<c)\mathbf{r}_{p}(X_i-c)^\prime, \mathbbm{1} (X_i\leq c)\mathbf{W}_i^\prime\otimes\mathbf{r}_{s}(X_i-c)^\prime)^\prime$, and a second set of regressors that refers to estimation to the right of the cutoff, i.e., $(\mathbbm{1} (X_i\geq c)\mathbf{r}_{p}(X_i-c)^\prime, \mathbbm{1} (X_i\geq c)\mathbf{W}_i^\prime\otimes\mathbf{r}_{s}(X_i-c)^\prime)^\prime$. Accordingly, by the partitioned regression theorem, the least squares coefficients in the “long” regression are equivalent to the least squares coefficients in the following “short” regressions:

which can be equivalently written in matrix form as

where the Hessian and score matrices are defined as

where \[\mathbf{R}_{p}(h)=

_{n \times (1+p)},\quad \mathbf{R}_{p,s}(h)=

_{n \times (1+p + m)},\] with the scaling and kernel matrices defined as \[\mathbf{H}_{p,s}(h) =

_{(1+p + m)\times(1+p + m)}, \qquad \mathbf{H}_q(h) = diag(h^\ell:\ell=0,\ldots,q),\]

In what follows, we mostly refer to the RD estimator to the right of the cutoff, $\widehat{\boldsymbol{\vartheta}}_{+, p, s}(h)$, but everything follows symmetrically for the other RD estimator to the left of the cutoff, $\widehat{\boldsymbol{\vartheta}}_{-, p, s}(h)$.

We focus our attention on estimators of the following form: \[\widehat{\chi}_{\nu,p,s}(\mathbf{h}):=\breve{\mathbf{e}}_{\nu}^\prime(\widehat{\boldsymbol{\vartheta}}_{+,p,s}(h_+)-\widehat{\boldsymbol{\vartheta}}_{-,p,s}(h_-)),\qquad \nu\in\{0,1,\ldots,p\wedge s\}, \] where for fixed $\mathbf{w}\in\operatorname{\mathbb{R}}^d$ \[\breve{\mathbf{e}}_\nu \equiv \breve{\mathbf{e}}_{\nu}(\mathbf{w}) :=

\] is the vector which extracts the $(\nu+1)$-th term in the polynomial basis in $(X_i-c)$ and the $(\nu+1)$-th terms in the interaction between such basis and $\mathbf{W}_i$. The identity matrices are there just to stress the different dimensions $\mathbf{e}_\nu$ takes in different places. Most of the time, we suppress the dependence of $\breve{\mathbf{e}}_\nu$ on $\mathbf{w}$ to avoid making the notation even more cumbersome and because in the case of orthogonal covariates, the only sensible evaluation point would be $\mathbf{w} = \boldsymbol{\iota}_d$.\footnote{This case corresponds to situations where each component is binary, $\mathbf{W}_i \in \{0,1\}^d$, and at most one component takes the value one, implying that $\mathbf{W}_i'\mathbf{W}_i \leq 1$.} The evaluation point $\mathbf{w}$ would only turn out to be useful when some elements in $\mathbf{W}_i$ are continuous variables, as the following example clarifies.

Moreover, let

be the fixed-$h$ best linear mean square error predictors of $\operatorname{\mathbb{E}}[Y_i(d)\mid X_i=c, \mathbf{W}_i], d\in\{0,1\}$, i.e.,

and

Assumptions

In what follows, we state all the assumptions we rely on to prove the results in this supplemental appendix. Assumption 1 in the main text contains Assumptions (ref), (ref), (ref), and (ref), whereas Assumption 2 coincides with (ref).

In words, under Assumption (ref) the sample is composed of independent draws from an underlying population of interest. Assumptions (ref), (ref), (ref), and (ref) are technical conditions on the data-generating process (DGP) we rely on when we establish convergence rates and characterize the probability limits of interest. Assumption (ref) is a standard assumption in least squares which requires the probability limit of the Gram matrix to be positive definite. Assumption (ref) also establishes that $\mathbf{W}$ can be thought of as a vector of pretreatment covariates. Assumption (ref) is necessary to Taylor expand the conditional expectation of the potential outcomes and characterize the smoothing bias. Assumptions (ref) and (ref) are standard bounded absolute higher-order moment conditions that make the Lyapunov condition satisfied and allow us to invoke a Lindeberg-Feller central limit theorem to show asymptotic normality of the RD estimator. Assumption (ref) imposes a partially linear (in $\mathbf{W}$) structure --locally at the cutoff-- on the population conditional expectation of the potential outcomes which we rely on when attaching a causal interpretation to the probability limit of the RD estimator. Assumption (ref) states standard technical conditions on the kernel used in the local polynomial regression.

Mapping between Main Text and Supplement

In the main text, we present all the results using the long regression with $p=s=1$

which aligns with the classical way of thinking about linear regressions. Indeed, this representation allows us to readily interpret the coefficients on the terms interacted with $T_i$ as specific differences between treated and control groups.

In this supplement, we set up the problem as the long regression

This is equivalent to (ref) in that the predicted values and the residuals of these two regressions are the same. However, the regression coefficients in (ref) have a different interpretation than those in (ref). For example, the coefficients on the terms interacted with $T_i$ do not describe differences between treated and control groups but rather capture only moments of the treated population. On the one hand, (ref) allows coefficients to be directly interpreted as long as differences between treated and control are of interest. On the other hand, as we already stressed in Section (ref), (ref) can be written as two separate “short" regressions with orthogonal design matrices, making statements and proofs less cumbersome from a notational standpoint. For this reason, we rely on (ref) in the main text and on (ref) in this supplement.

Table (ref) maps the notation used in this supplement with the one used in the main paper. Dependence on $h$ is omitted to simplify notation. For the MSE constants, we will use the relationship \[\widehat{\boldsymbol{\varsigma}}(\mathbf{h}) = \mathbf{M}\widehat{\boldsymbol{\vartheta}}_{1,1}(\mathbf{h}) =\mathbf{M}

, \qquad \mathbf{M} :=

.\]

We conclude this section by mapping the results in the main text with the proofs in the supplement:

Auxiliary Lemmas and Quantities

In this subsection, we introduce a series of auxiliary lemmas and quantities we rely upon in the rest of the supplemental appendix. We heavily rely on the next lemma for asymptotic statements.

The next lemma handles the Gram matrices \[\widehat{\boldsymbol{\Lambda}}_{-,p,s}(h) = \mathbf{R}_p(h)^\prime\mathbf{K}_-(h)\mathbf{R}_s(h)/n, \qquad \widehat{\boldsymbol{\Lambda}}_{+,p,s}(h) = \mathbf{R}_p(h)^\prime\mathbf{K}_+(h)\mathbf{R}_s(h)/n,\] and shows the object concentrates around in probability and the rate at which such concentration occurs.

In the next lemma, we give conditions for the asymptotic invertibility of $\widehat{\boldsymbol{\Gamma}}_{-,p,s}(h)$ and $\widehat{\boldsymbol{\Gamma}}_{+,p,s}(h)$, thereby making local polynomial estimators well-defined in large samples.

To have a more compact notation, we introduce the following quantities for integers $p,s,a\in \operatorname{\mathbb{N}}_0$:

The following two lemmas characterize the asymptotic properties of $\widehat{\boldsymbol{\zeta}}_{-,p,s,a}(h), \widehat{\boldsymbol{\zeta}}_{+,p,s,a}(h), \widehat{\boldsymbol{\varphi}}_{-,p,s,a}(h)$, and $\widehat{\boldsymbol{\varphi}}_{+,p,s,a}(h)$. These results are used below to guarantee that the “constant terms” in the MSE expansions are asymptotically well-defined.

Main Results

In this section, we mostly focus on the RD estimator to the right of the cutoff as everything holds symmetrically for the RD estimator to the left. All the lemmas and main results are reported for both estimators.

The RD estimator with interacted covariates to the right of the cutoff is defined as \[\widehat{\boldsymbol{\vartheta}}_{+,p,s}(h):= \operatorname*{arg\,min}_{\mathbf{t} \in \operatorname{\mathbb{R}}^{1+p+m}}\sum_{i=1}^n\mathbbm{1}(X_i\geq c)K\left(\frac{X_i-c}{h}\right)(Y_i - \mathbf{r}_{p,s}(X_i-c,\mathbf{W}_i)^\prime\mathbf{t})^2.\]

To analyze the statistical properties of the RD estimator with interacted covariates $\widehat{\boldsymbol{\vartheta}}_{+,p,s}(h)$, we take advantage of the following decomposition:

where $\boldsymbol{\vartheta}_{+,p,s}^\star(h)$ the fixed-$h$ best linear mean square error predictor of $\operatorname{\mathbb{E}}[Y_i(1)\mid X_i=c, \mathbf{W}_i]$ and \[u_{+,i}(h):= Y_i - \mathbf{r}_{p,s}\left(X_i-c,\mathbf{W}_i\right)^\prime \boldsymbol{\vartheta}^\star_{+,p,s}(h)\] is the fixed-$h$ $\ell_2$-projection residual. This decomposition gives us

where

Asymptotic Approximation and Asymptotic Variance

Lemma (ref) shows that $\widehat{\boldsymbol{\Gamma}}_{+,p,s}(h)$ is asymptotically invertible. This Lemma takes care of the “denominator" of the RD estimator by showing that \[\widehat{\boldsymbol{\Gamma}}_{+,p,s}(h) = {\boldsymbol{\Gamma}}_{+,p,s} + o(1) + O_{\mathbb{P}}(1/\sqrt{nh}).\]

Coming to the $\mathbf{L}_{+,p,s}(h)$ term, first note that $\operatorname{\mathbb{E}}[\mathbf{L}_{+,p,s}(h)] = 0$ because of the properties of $\ell_2$-projection residuals. Regarding the variance of this term, the next lemma shows that $\operatorname{\mathbb{V}}[\mathbf{L}_{+,p,s}(h)]$ is $O(1/\sqrt{nh})$. Then, it follows that $\mathbf{L}_{+,p,s}(h) = O_{\mathbb{P}}(1/\sqrt{nh})$ by Lemma (ref).

Typically, we are interested in estimators of the form $\breve{\mathbf{e}}_{\nu}^\prime \widehat{\boldsymbol{\vartheta}}_{+,p,s}(h_+)$. By Lemma (ref) and Lemma (ref) and Slutsky's theorem, we know that \[\sqrt{nh^{2\nu+1}}\breve{\mathbf{e}}_{\nu}^\prime(\widehat{\boldsymbol{\vartheta}}_{+,p,s}(h) - \boldsymbol{\vartheta}_{+,p,s}^\star(h)) = \sqrt{nh}\breve{\mathbf{e}}_{\nu}^\prime \widehat{\boldsymbol{\Gamma}}_{+,p,s}^{-1}(h)\mathbf{L}_{+,p,s}(h)=O_{\mathbb{P}}(1).\] The next theorem shows that the term above is not only bounded in probability but also converging in distribution to a Normal random variable.

Variance Estimation

To estimate the asymptotic variance of the RD estimator, we propose the following plug-in estimator \[\widehat{\mathcal{V}}_{+,\nu,p,s}(h) = \breve{\mathbf{e}}_{\nu}^\prime\widehat{\boldsymbol{\Gamma}}_{+,p,s}^{-1}(h) \widehat{\mathbf{V}}_{+,p,s}(h) {\widehat{\boldsymbol{\Gamma}}_{+,p,s}^{-1}(h)}^\prime\breve{\mathbf{e}}_{\nu},\] where \[\widehat{\mathbf{V}}_{+,p,s}(h):= \frac{1}{nh}\sum_{i=1}^n\sqrt{w_{+,i}(h)}\cdot\mathbbm{1}(X_i\geq c)K\left(\frac{X_i-c}{h}\right)\mathbf{r}_{p,s}\left(\frac{X_i-c}{h},\mathbf{W}_i\right)\mathbf{r}_{p,s}\left(\frac{X_i-c}{h},\mathbf{W}_i\right)^\prime {\widehat{u}_{+,i}}^2(h),\] with $\widehat{u}_{+,i}(h):= Y_i - \mathbf{r}_{p,s}(X_i-c,\mathbf{W}_i)^\prime\widehat{\boldsymbol{\vartheta}}_{+,p,s}(h)$.

The weights $w_{+,i}(h)$ allow for different HC-type estimators. Namely,

where

$\mathbf{Q}_{-, p,s}, $ and $ \mathbf{Q}_{+, p,s}$ are the “projection" matrices used to get the estimated residuals and are defined as,

and \[L_{-,i} := \mathbf{e}_i^{\prime} \mathbf{Q}_{-, p,s} \mathbf{e}_i,\qquad L_{+,i} := \mathbf{e}_i^{\prime} \mathbf{Q}_{+, p,s} \mathbf{e}_i,\] which denote the leverage of each observation.

Define the standard error of the RD estimators as \[\widehat{\sigma}_{-,\nu,p,s}(h) := \left(\frac{1}{{nh^{2\nu+1}}}\widehat{\mathcal{V}}_{-,\nu,p,s}(h)\right)^{1/2}, \qquad \widehat{\sigma}_{+,\nu,p,s}(h) := \left(\frac{1}{{nh^{2\nu+1}}}\widehat{\mathcal{V}}_{+,\nu,p,s}(h)\right)^{1/2}.\] The lemma above naturally yields the following corollary via Slutsky's theorem.

Clustered Data

In the case of clustered data, the extension of the above results is immediate. The only difference would be reflected in the form of the “meat" matrices $\mathbf{V}_{-,\nu,p,s},$ and $ \mathbf{V}_{+,\nu,p,s}$, which will ultimately depend on the particular form of clustering being used. For a review on cluster-robust inference, see Cameron2015_JHR and mackinnon2023ClusterrobustInferenceGuide.

We assume that each unit $i\in\{1,2,\cdots,n\}$ belongs to a single cluster $s(i)\in\{1,2,\cdots,G\},$ where $s:\{1,2,\cdots,n\}\to\{1,2,\cdots,G\}$. Furthermore, we consider an asymptotic regime where the number of clusters $G$ grows large, $G\to\infty$, and dominates the bandwidth $Gh\to\infty$.

To estimate the asymptotic variance of the RD estimator under clustering, we propose the following plug-in estimator \[\widehat{\mathcal{V}}^{\mathsf{CL}}_{+,\nu,p,s}(h) = \frac{(G-1)n}{(G-1)(n-p-1-d)}\breve{\mathbf{e}}_{\nu}^\prime\widehat{\boldsymbol{\Gamma}}_{+,p,s}^{-1}(h) \widehat{\mathbf{V}}^{\mathsf{CL}}_{+,p,s}(h) {\widehat{\boldsymbol{\Gamma}}_{+,p,s}^{-1}(h)}^\prime\breve{\mathbf{e}}_{\nu},\] where

This estimator, as well as many more, is implemented in our R and Stata software.

Smoothing Bias and Probability Limit

Without further assumptions on the data generating process or on $\mu_{Y-}(X_i,\mathbf{W}_i)$ and $\mu_{Y+}(X_i,\mathbf{W}_i)$, the probability limit of the RD estimator has a generic best linear mean squared error predictor interpretation. We need to leverage extra structure to attach more interpretation to this probability limit. In this spirit, we posit Assumption (ref), which we state again here below for the reader's convenience. \addtocounter{assumption}{-2}

Let

where

with $\boldsymbol{\alpha}_{-,p}\equiv \boldsymbol{\alpha}_{-,p}(c), \boldsymbol{\lambda}_{-,s}\equiv \boldsymbol{\lambda}_{-,s}(c),\boldsymbol{\alpha}_{+,p}\equiv \boldsymbol{\alpha}_{+,p}(c),$ and $ \boldsymbol{\lambda}_{+,s}\equiv \boldsymbol{\lambda}_{+,s}(c)$.

Consider the following decomposition

From Theorem (ref), we know that $\mathbf{H}_{p,s}(h)(\widehat{\boldsymbol{\vartheta}}_{+,p,s}(h) - {\boldsymbol{\vartheta}}^\star_{+,p,s}(h)) = O_{\mathbb{P}}(1/\sqrt{nh})$. In what follows, we show that the second term is of order $O_{\mathbb{P}}(h^{1+p\wedge s})$. Throughout, we consider an asymptotic regime in which this term is negligible with respect to the first term. It is precisely in this spirit that we will refer to ${\boldsymbol{\vartheta}}^\star_{+,p,s}(h) - \boldsymbol{\vartheta}_{+,p,s}$ as a “bias" term. The next lemma shows that the leading term of the bias is of order $O(h^{1+p\wedge s})$.

The following corollary follows naturally from Theorem (ref) and Lemma (ref).

Consistency

Using Lemma (ref) and Lemma (ref), we have

thus if $nh \to \infty$ and $h\to 0$, then $\widehat{\boldsymbol{\vartheta}}_{-,p,s}(h) \overset{\mathbb{P}}{\longrightarrow} \boldsymbol{\vartheta}_{-,p,s}$ and $\widehat{\boldsymbol{\vartheta}}_{+,p,s}(h) \overset{\mathbb{P}}{\longrightarrow} \boldsymbol{\vartheta}_{+,p,s}$.

Typically, we are interested in estimators of the form $\breve{\mathbf{e}}_{\nu}^\prime \widehat{\boldsymbol{\vartheta}}_{-,p,s}(h)$ or $\breve{\mathbf{e}}_{\nu}^\prime \widehat{\boldsymbol{\vartheta}}_{+,p,s}(h)$. In this case, we get

In an asymptotic regime where also $n\min\{h_-,h_+\}^{2\nu+1} \to \infty$ and $h\to 0$, then $\breve{\mathbf{e}}_{\nu}^\prime\widehat{\boldsymbol{\vartheta}}_{-,p,s}(h) \overset{\mathbb{P}}{\longrightarrow} \breve{\mathbf{e}}_{\nu}^\prime{\boldsymbol{\vartheta}}_{-,p,s}$ and $\breve{\mathbf{e}}_{\nu}^\prime\widehat{\boldsymbol{\vartheta}}_{+,p,s}(h) \overset{\mathbb{P}}{\longrightarrow} \breve{\mathbf{e}}_{\nu}^\prime{\boldsymbol{\vartheta}}_{+,p,s}$.

Define \[\widehat{\chi}_{\nu,p,s}(\mathbf{h}) := \breve{\mathbf{e}}_{\nu}^\prime\left(\widehat{\boldsymbol{\vartheta}}_{+,p,s}(h_+)-\widehat{\boldsymbol{\vartheta}}_{-,p,s}(h_-)\right), \quad {\chi}_{\nu,p,s} := \breve{\mathbf{e}}_{\nu}^\prime\left({\boldsymbol{\vartheta}}_{+,p,s}-{\boldsymbol{\vartheta}}_{-,p,s}\right).\] Then \[\widehat{\chi}_{\nu,p,s}(\mathbf{h}) - {\chi}_{\nu,p,s} = O_{\mathbb{P}}\left(\frac{1}{\sqrt{n\min\{h_-,h_+\}^{2\nu+1}}} + \max\{h_-,h_+\}^{1+p\wedge s-\nu}\right),\] and $\widehat{\chi}_{\nu,p,s}(\mathbf{h}) \overset{\mathbb{P}}{\longrightarrow} {\chi}_{\nu,p,s}.$

Causal Interpretation of the Probability Limit

Before delving into this section, we redefine the “extractor" vector as \[\breve{\mathbf{e}}_{\nu_x,\nu_w} =

\] to be the vector extracting the $(\nu_x+1)$-th term in the polynomial basis in $(X_i-c)$ and the $(\nu_w+1)$-th terms in the interaction between such basis and $\mathbf{W}_i$. Furthermore, we define $\breve{\mathbf{e}}_{\nu_x,\cdot} \equiv \mathbf{I}_{1+p+d}\mathbf{e}_{\nu_x},$ to be the vector that extracts only the $(\nu_x+1)$-th term in the polynomial basis in $(X_i-c)$. We also define accordingly \[\chi_{\nu_x,\nu_w,p,s} := \breve{\mathbf{e}}_{\nu_x,\nu_w}^\prime \left({\boldsymbol{\vartheta}}_{+,p,s}-{\boldsymbol{\vartheta}}_{-,p,s}\right).\]

The quantity $\chi_{\nu_x,\nu_w,p,s}$ has a natural causal interpretation in many cases, which depends on the elements extracted:

MSE Expansions

We now provide first-order expansions for the MSE of the RD estimator. This is crucial to obtain formulas for MSE-optimal bandwidths. We first start by giving a Nagar expansion of the unconditional MSE. Then, we provide an approximation for the bias and the variance of the RD estimator, where the approximation comes from disregarding higher-order terms.

First of all, recall that the RD estimator and the best linear mean square error predictor are \[\widehat{\boldsymbol{\vartheta}}_{+,p,s}(h) =\mathbf{H}_{p,s}^{-1}(h)\widehat{\boldsymbol{\Gamma}}^{-1}_{+,p,s}(h)\widehat{\boldsymbol{\Upsilon}}_{+,p,s}(h), \qquad \boldsymbol{\vartheta}^\star_{+,p,s}(h) =\mathbf{H}_{p,s}^{-1}(h)\widetilde{\boldsymbol{\Gamma}}^{-1}_{+,p,s}(h)\widetilde{\boldsymbol{\Upsilon}}_{+,p,s}(h).\] Consider the following decomposition for fixed $\mathbf{w}\in\operatorname{\mathbb{R}}^d$

where we refer to $\mathsf{L}_{+,\nu,p,s}(h)$ as the “linear" term, to $\mathsf{Q}_{+,\nu,p,s}(h)$ as the “quadratic" term, and to $\mathsf{B}_{+,\nu,p,s}(h)$ as the “bias" term. By Theorem (ref) and Lemma (ref), we know that $\mathsf{L}_{+,\nu,p,s}(h) = O_{\mathbb{P}}(1/\sqrt{nh^{2\nu+1}})$ and $\mathsf{Q}_{+,\nu,p,s}(h) = O_{\mathbb{P}}(1/\sqrt{n^2h^{2\nu+2}}) = o_{\mathbb{P}}(\mathsf{L}_{+,\nu,p,s}(h))$. Moreover, under Assumption (ref), we have that $\mu_+^{(\nu)}(c,\mathbf{W}_i)=\breve{\mathbf{e}}'\boldsymbol{\vartheta}_{+,p,s}(h)$, thus by Lemma (ref) we have that $\mathsf{B}_{+,\nu,p,s}(h) = O(h^{1+p\wedge s -\nu})$. In an asymptotic regime in which $nh^{2+p\wedge s}\to \infty$, we have that $\mathsf{Q}_{+,\nu,p,s}(h)=o_{\mathbb{P}}(\mathsf{B}_{+,\nu,p,s}(h))$. Therefore, under Assumptions (ref)-(ref), as $n\to\infty$ and $nh^{2+p\wedge s}\to\infty$ we have \[\breve{\mathbf{e}}_\nu'(\widehat{\boldsymbol{\vartheta}}_{+,p,s}(h) - \boldsymbol{\vartheta}_{+,p,s}(h)) = \mathsf{L}_{+,\nu,p,s}(h) + \mathsf{B}_{+,\nu,p,s}(h) + o_{\mathbb{P}}(\min\{\mathsf{L}_{+,\nu,p,s}(h),\mathsf{B}_{+,\nu,p,s}(h)\}).\] With this result at hand, we define the first-order approximation of the mean squared error as \[\mathsf{MSE}[\breve{\mathbf{e}}_{\nu}^\prime (\widehat{\boldsymbol{\vartheta}}_{+,p,s}(h) -\boldsymbol{\vartheta}_{+,p,s})]= \operatorname{\mathbb{E}}\left[(\mathsf{L}_{+,\nu,p,s}(h) + \mathsf{B}_{+,\nu,p,s}(h))^2\right] = \operatorname{\mathbb{V}}[\mathsf{L}_{+,\nu,p,s}(h)] + \mathsf{B}_{+,\nu,p,s}(h)^2,\] where the second equality follows from the fact that $ \mathsf{B}_{+,\nu,p,s}(h)$ is non-random and $\mathsf{L}_{+,\nu,p,s}(h)$ is unconditionally mean-zero. Similarly, we define

for the left side of the cutoff, whereas for the difference, we define

Under Assumptions (ref)-(ref) and relying on Lemma (ref), the bias approximations for the RD estimators are given by

where

Furthermore, we define

Under Assumptions (ref)-(ref) and relying on Corollary (ref)

The variance approximations for the RD estimators are given by \[\mathcal{V}_{-,\nu,p,s}(h)\to\mathcal{V}_{-,\nu,p,s}:=\breve{\mathbf{e}}_{\nu}^\prime{\boldsymbol{\Gamma}}_{-,p,s}^{-1}\mathbf{V}_{-,p,s} {{\boldsymbol{\Gamma}}_{-,p,s}^{-1\prime}}\breve{\mathbf{e}}_{\nu}, \qquad \mathcal{V}_{+,\nu,p,s}(h)\to\mathcal{V}_{+,\nu,p,s}:=\breve{\mathbf{e}}_{\nu}^\prime{\boldsymbol{\Gamma}}_{+,p,s}^{-1} \mathbf{V}_{+,p,s} {{\boldsymbol{\Gamma}}_{+,p,s}^{-1\prime}}\breve{\mathbf{e}}_{\nu}.\]

With these approximations, the MSE expansions of the RD estimators are given by

Accordingly, we define an MSE-optimal bandwidth as the minimizer of the Nagar expansion of the unconditional MSE.

\paragraph{One-sided Optimal Bandwidths.} The MSE-optimal bandwidths are defined as

so, under the additional assumption that $\mathcal{B}_{-,\nu,p,s}\neq0\neq \mathcal{B}_{+,\nu,p,s}$, we get

\paragraph{Two-sided Optimal Bandwidth.} In this case a single bandwidth is chosen, so $h=h_+=h_-$ and the optimal bandwidth is defined as

so, under the additional assumption that $\mathcal{B}_{+,\nu,p,s}-\mathcal{B}_{-,\nu,p,s}\neq 0$, we get \[h^\star_{\nu,p,s} =\left[\frac{1+2\nu}{2(1+(p\wedge s)-\nu)n}\frac{\mathcal{V}_{+,\nu,p,s}+\mathcal{V}_{-,\nu,p,s}}{(\mathcal{B}_{+,\nu,p,s}-\mathcal{B}_{-,\nu,p,s})^2}\right]^{\frac{1}{3+2(p\wedge s)}}.\]

Empirical Practice Investigating Covariate-Heterogeneity

To examine common empirical practices regarding heterogeneous effects, we sampled papers from AEA journals over the past decade. The papers are listed in Table (ref) below, along with a summarization of how treatment effect heterogeneity was studied in each, including the empirical setting, the type of heterogeneity, and details on how estimation and inference were conducted. For full details of each study, please refer to the original work.

\cleardoublepage

Proofs

Proof of Lemma (ref)

Proof of Lemma (ref)

Proof of Lemma (ref)

Proof of Lemma (ref)

Proof of Lemma (ref)

Proof of Theorem (ref)

Proof of Theorem (ref)

Proof of Lemma (ref)