The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
113,860 characters
Online Supplement for ``Treatment Effect Heterogeneity in Regression Discontinuity Designs''
\title{Online Supplement for ``Treatment Effect Heterogeneity in Regression Discontinuity Designs''}
\author{Sebastian Calonico\thanks{Graduate School of Management, UC Davis.} \and
Matias D. Cattaneo\thanks{Department of Operations Research and Financial Engineering, Princeton University.} \and
Max H. Farrell\thanks{Department of Economics, UC Santa Barbara.} \and
Filippo Palomba\thanks{Department of Economics, Princeton University.} \and
Roc{\'i}o Titiunik\thanks{Department of Politics, Princeton University.}}
\maketitle
\onehalfspacing
\begin{abstract}
\noindent This supplement contains all proofs, additional results, and other technical details about estimation, inference, and identification. Section \ref{sec: introduction} describes setup and notation, states the assumptions we rely on, and introduces some auxiliary results. Section \ref{sec: main results} illustrates the main technical results. Section \ref{sec: empirical app} discusses in detail the common strategies used in empirical analyses. Section \ref{sec: proofs} contains all the proofs.
\end{abstract}
\vspace{-.1in}
\clearpage
\tableofcontents
\singlespacing
\clearpage
\section{Introduction}
\label{sec: 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.
\subsection{Notation}
\label{sec: 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).
\textbf{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.
\textbf{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.
\textbf{Causal Model.} We now describe the population causal model. The outcome variable is
\begin{align*}
Y_i&=T_i \cdot Y_i(1)+\left(1-T_i\right) \cdot Y_i(0),
\end{align*}
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
\begin{align*}
\mathbf{Y}&=\left[Y_{1}, \cdots, Y_{n}\right]^{\prime},&
\mathbf{W} &= [\mathbf{W}_1,\mathbf{W}_2,\cdots,\mathbf{W}_n]^\prime,\\
\boldsymbol{\mu}_{W}(x)&=\operatorname{\mathbb{E}}[\mathbf{W}_i|X_i=x], & \boldsymbol{\mu}_{WW}(x)&=\operatorname{\mathbb{E}}[\mathbf{W}_i\mathbf{W}_i^\prime|X_i=x], \\
\mu_{W_\ell}(x)&=\operatorname{\mathbb{E}}[W_{i\ell}|X_i=x], & \mu_{W_\ell W_j}(x)&=\operatorname{\mathbb{E}}[W_{i\ell}W_{ij}|X_i=x], \quad \ell,j \in\{ 1,\cdots, d\}, \\
\mu_{Y-}(x,\mathbf{w}) &= \operatorname{\mathbb{E}}[Y_i(0)\mid X_i=x, \mathbf{W}_i=\mathbf{w}], & \mu_{Y+}(x,\mathbf{w}) &= \operatorname{\mathbb{E}}[Y_i(1)\mid X_i=x, \mathbf{W}_i=\mathbf{w}], \\
\mu^{(\nu)}_{Y-}(x,\mathbf{w}) &= \frac{\partial^\nu}{\partial x^\nu}\operatorname{\mathbb{E}}[Y_i(0)\mid X_i=x, \mathbf{W}_i=\mathbf{w}], & \mu_{Y+}(x,\mathbf{w}) &= \frac{\partial^\nu}{\partial x^\nu}\operatorname{\mathbb{E}}[Y_i(1)\mid X_i=x, \mathbf{W}_i=\mathbf{w}], \\
\sigma^2_{Y-}(x,\mathbf{w}) &= \operatorname{\mathbb{V}}[Y_i(0)\mid X_i=x, \mathbf{W}_i=\mathbf{w}], & \sigma^2_{Y+}(x,\mathbf{w}) &= \operatorname{\mathbb{V}}[Y_i(1)\mid X_i=x, \mathbf{W}_i=\mathbf{w}],
\end{align*}
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. \]
\subsection{Setup}
\label{subsec: setup}
Formally, the RD estimators with treatment interactions can be obtained from a ``long'' regression, i.e.,
\begin{align*}
\widehat{\boldsymbol{\vartheta}}_{p, s}(\mathbf{h})
&:=\operatorname*{arg\,min}_{\substack{\mathbf{a}_-,\mathbf{a}_+ \in \mathbb{R}^{1+p}, \\ \boldsymbol{\ell}_-,\boldsymbol{\ell}_+ \in \mathbb{R}^{d(1+s)}}} \sum_{i=1}^n \left(Y_i-
\begin{bmatrix}\mathbf{r}_{p}\left(X_i-c\right)\\\mathbf{r}_{p}\left(X_i-c\right)\end{bmatrix}^{\prime}
\begin{bmatrix}\mathbf{a}_-\\\mathbf{a}_+\end{bmatrix} -
\begin{bmatrix}(\mathbf{W}_i\otimes\mathbf{r}_{s}(X_i-c)) \\ (\mathbf{W}_i\otimes\mathbf{r}_{s}(X_i-c))\end{bmatrix}^{\prime}
\begin{bmatrix}\boldsymbol{\ell}_- \\ \boldsymbol{\ell}_+ \end{bmatrix}\right)^2 K_{\mathbf{h}}\left(X_i-c\right) \\
\widehat{\boldsymbol{\vartheta}}_{ p, s}(\mathbf{h})
&=
\begin{bmatrix}
\widehat{\boldsymbol{\vartheta}}_{-, p, s}(h_-)\\
\widehat{\boldsymbol{\vartheta}}_{+, p, s}(h_+)
\end{bmatrix}
=
\begin{bmatrix}
\widehat{\boldsymbol{\alpha}}_{-, p}(h_-) \\
\widehat{\boldsymbol{\lambda}}_{-, s}(h_-) \\
\widehat{\boldsymbol{\alpha}}_{+, p}(h_+) \\
\widehat{\boldsymbol{\lambda}}_{+, s}(h_+)
\end{bmatrix}, \qquad p,s\in\operatorname{\mathbb{N}},
\end{align*}
where
\begin{align*}
\widehat{\boldsymbol{\alpha}}_{-,p}(h) =(\widehat{\alpha}_{-,0}, \cdots, \widehat{\alpha}_{-,p})',\qquad
\widehat{\boldsymbol{\alpha}}_{+,p}(h) =(\widehat{\alpha}_{+,0}, \cdots, \widehat{\alpha}_{+,p})',
\end{align*}
and
\begin{align*}
\widehat{\boldsymbol{\lambda}}_{-,s}(h) &= (\widehat{\lambda}_{-,10}(h),\cdots, \widehat{\lambda}_{-,1s}(h), \cdots, \widehat{\lambda}_{-,d0}(h),\cdots, \widehat{\lambda}_{-,ds}(h))', \\
\widehat{\boldsymbol{\lambda}}_{+,s}(h) &= (\widehat{\lambda}_{+,10}(h),\cdots, \widehat{\lambda}_{+,1s}(h), \cdots, \widehat{\lambda}_{+,d0}(h),\cdots, \widehat{\lambda}_{+,ds}(h))'.
\end{align*}
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:
\begin{align*}
\widehat{\boldsymbol{\vartheta}}_{-, p,s}(h)=\left[\begin{array}{c}
\widehat{\boldsymbol{\alpha}}_{-, p}(h) \\
\widehat{\boldsymbol{\lambda}}_{-, s}(h)
\end{array}\right]=\operatorname*{arg\,min}_{\substack{\mathbf{a} \in \mathbb{R}^{1+p},\\ \boldsymbol{\ell} \in \mathbb{R}^m}} \sum_{i=1}^n \mathbbm{1}\left(X_i<c\right)\left(Y_i-\mathbf{r}_p\left(X_i-c\right)^{\prime} \mathbf{a}-\left(\mathbf{W}_i\otimes\mathbf{r}_s(X_i-c)\right)^{\prime} \boldsymbol{\ell}\right)^2 K_h\left(X_i-c\right), \\
\widehat{\boldsymbol{\vartheta}}_{+, p, s}(h)=\left[\begin{array}{c}
\widehat{\boldsymbol{\alpha}}_{+, p}(h) \\
\widehat{\boldsymbol{\lambda}}_{+, s}(h)
\end{array}\right]=\operatorname*{arg\,min}_{\substack{\mathbf{a} \in \mathbb{R}^{1+p},\\ \boldsymbol{\ell} \in \mathbb{R}^m}} \sum_{i=1}^n \mathbbm{1}\left(X_i \geq c\right)\left(Y_i-\mathbf{r}_p\left(X_i-c\right)^{\prime} \mathbf{a}-(\mathbf{W}_i\otimes\mathbf{r}_s(X_i-c))^{\prime} \boldsymbol{\ell}\right)^2 K_h\left(X_i-c\right),
\end{align*}
which can be equivalently written in matrix form as
\begin{align*}
\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
\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),
\end{align*}
where the Hessian and score matrices are defined as
\begin{align*}
\widehat{\boldsymbol{\Gamma}}_{-,p,s}(h)=\mathbf{R}_{p,s}(h)^\prime\mathbf{K}_-(h)\mathbf{R}_{p,s}(h)/n, \qquad &\widehat{\boldsymbol{\Upsilon}}_{-,p,s}(h) =\mathbf{R}_{p,s}(h)^\prime\mathbf{K}_-(h)\mathbf{Y}/n, \\
\widehat{\boldsymbol{\Gamma}}_{+,p,s}(h)=\mathbf{R}_{p,s}(h)^\prime\mathbf{K}_+(h)\mathbf{R}_{p,s}(h)/n, \qquad &\widehat{\boldsymbol{\Upsilon}}_{+,p,s}(h) =\mathbf{R}_{p,s}(h)^\prime\mathbf{K}_+(h)\mathbf{Y}/n, \\
\end{align*}
where
\[\mathbf{R}_{p}(h)=\begin{bmatrix}
\mathbf{r}_p\left(\frac{X_1-c}{h}\right)^\prime \\
\mathbf{r}_p\left(\frac{X_2-c}{h}\right)^\prime \\
\vdots \\
\mathbf{r}_p\left(\frac{X_n-c}{h}\right)^\prime \\
\end{bmatrix}_{n \times (1+p)},\quad \mathbf{R}_{p,s}(h)=\begin{bmatrix}
\mathbf{r}_p\left(\frac{X_1-c}{h}\right)^\prime & \mathbf{W}_1'\otimes\mathbf{r}_s\left(\frac{X_1-c}{h}\right)^\prime \\
\mathbf{r}_p\left(\frac{X_2-c}{h}\right)^\prime & \mathbf{W}_2'\otimes\mathbf{r}_s\left(\frac{X_2-c}{h}\right)^\prime \\
\vdots & \vdots \\
\mathbf{r}_p\left(\frac{X_n-c}{h}\right)^\prime & \mathbf{W}_n'\otimes\mathbf{r}_s\left(\frac{X_n-c}{h}\right)^\prime \\
\end{bmatrix}_{n \times (1+p + m)},\]
with the scaling and kernel matrices defined as
\[\mathbf{H}_{p,s}(h) =
\begin{bmatrix}
\mathbf{H}_p(h) & \mathbf{0}_{(1+p)\times m} \\
\mathbf{0}_{m\times (1+p)} & \mathbf{I}_d\otimes \mathbf{H}_s(h) \\
\end{bmatrix}_{(1+p + m)\times(1+p + m)}, \qquad \mathbf{H}_q(h) = \text{diag}(h^\ell:\ell=0,\ldots,q),\]
\begin{align*}
\mathbf{K}_-(h) &= \text{diag}\left(\mathbbm{1} ( X_i< c)k_h(X_i-c): i = 1, \ldots, n\right), \quad
\mathbf{K}_+(h) = \text{diag}\left(\mathbbm{1} ( X_i\geq c)k_h(X_i-c): i = 1, \ldots, n\right).
\end{align*}
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}) := \begin{bmatrix}
\nu\mathbf{I}_{p+1}\mathbf{e}_{\nu} \\ \nu!\mathbf{w} \otimes \mathbf{I}_{s+1}\mathbf{e}_{{\nu}}
\end{bmatrix}\]
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.
\begin{example}
Suppose $d=1, W_i \in\{0,1\}$ and $\nu=0$, then
\[\breve{\mathbf{e}}_{0}^\prime\widehat{\boldsymbol{\vartheta}}_{+,p,s}(h) = \mathbf{e}_{0}^\prime\widehat{\boldsymbol{\alpha}}_{+,p}(h) + \mathbf{e}_{0}^\prime\widehat{\boldsymbol{\lambda}}_{+,s}(h)= \underbrace{\widehat{\alpha}_{+,0}(h)}_{\text{coefficient on } \mathbbm{1}(X_i\geq c)} + \underbrace{\widehat{\lambda}_{+,10}(h)}_{\text{coefficient on } W_i\cdot \mathbbm{1}(X_i\geq c) },\]
and
\[\breve{\mathbf{e}}_{0}^\prime\widehat{\boldsymbol{\vartheta}}_{-,p,s}(h) = \mathbf{e}_{0}^\prime\widehat{\boldsymbol{\alpha}}_{-,p}(h) + \mathbf{e}_{0}^\prime\widehat{\boldsymbol{\lambda}}_{-,s}(h)= \underbrace{\widehat{\alpha}_{-,0}(h)}_{\text{coefficient on } \mathbbm{1}(X_i< c)} + \underbrace{\widehat{\lambda}_{-,10}(h)}_{\text{coefficient on } W_i\cdot \mathbbm{1}(X_i< c) }.\]
If instead $\nu=1$, then
\[\breve{\mathbf{e}}_{1}^\prime\widehat{\boldsymbol{\vartheta}}_{+,p,s}(h) = \mathbf{e}_{1}^\prime\widehat{\boldsymbol{\alpha}}_{+,p}(h) + \mathbf{e}_{1}^\prime\widehat{\boldsymbol{\lambda}}_{+,s}(h)= \underbrace{\widehat{\alpha}_{+,1}(h)}_{\text{coefficient on } \mathbbm{1}(X_i\geq c)\cdot(X_i-c)} + \underbrace{\widehat{\lambda}_{+,11}(h)}_{\text{coefficient on } W_i\cdot \mathbbm{1}(X_i\geq c)\cdot(X_i-c) },\]
and
\[\breve{\mathbf{e}}_{1}^\prime\widehat{\boldsymbol{\vartheta}}_{-,p,s}(h) = \mathbf{e}_{1}^\prime\widehat{\boldsymbol{\alpha}}_{-,p}(h) + \mathbf{e}_{1}^\prime\widehat{\boldsymbol{\lambda}}_{-,s}(h)= \underbrace{\widehat{\alpha}_{-,1}(h)}_{\text{coefficient on } \mathbbm{1}(X_i< c)\cdot(X_i-c)} + \underbrace{\widehat{\lambda}_{-,11}(h)}_{\text{coefficient on } W_i\cdot \mathbbm{1}(X_i< c)\cdot(X_i-c) }.\]
Now, suppose that $W_i$ is a continuous covariate and $\nu = 0$. In this case, the evaluation point $w=w_0$ does matter, and so
\[\breve{\mathbf{e}}_{0}^\prime\widehat{\boldsymbol{\vartheta}}_{+,p,s}(h) = \mathbf{e}_{0}^\prime\widehat{\boldsymbol{\alpha}}_{+,p}(h) + w_0\mathbf{e}_{0}^\prime\widehat{\boldsymbol{\lambda}}_{+,s}(h)= \underbrace{\widehat{\alpha}_{+,0}(h)}_{\text{coefficient on } \mathbbm{1}(X_i\geq c)} + \underbrace{\widehat{\lambda}_{+,10}(h)}_{\substack{\text{coefficient on } \\ W_i\cdot \mathbbm{1}(X_i\geq c)} }\cdot\, w_0,\]
and
\[\breve{\mathbf{e}}_{0}^\prime\widehat{\boldsymbol{\vartheta}}_{-,p,s}(h) = \mathbf{e}_{0}^\prime\widehat{\boldsymbol{\alpha}}_{-,p}(h) + \mathbf{e}_{0}^\prime\widehat{\boldsymbol{\lambda}}_{-,s}(h)= \underbrace{\widehat{\alpha}_{-,0}(h)}_{\text{coefficient on } \mathbbm{1}(X_i< c)} + \underbrace{\widehat{\lambda}_{-,10}(h)}_{\substack{\text{coefficient on } \\ W_i\cdot \mathbbm{1}(X_i< c)} }\cdot\, w_0.\]
\end{example}
Moreover, let
\begin{align*}
\boldsymbol{\vartheta}_{-,p,s}^\star(h)&:= \operatorname*{arg\,min}_{\mathbf{t}\in\operatorname{\mathbb{R}}^{1+p+m}}\operatorname{\mathbb{E}}\left[\mathbbm{1}(X_i< 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\right],\\
\boldsymbol{\vartheta}_{+,p,s}^\star(h)&:= \operatorname*{arg\,min}_{\mathbf{t}\in\operatorname{\mathbb{R}}^{1+p+m}}\operatorname{\mathbb{E}}\left[\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\right],
\end{align*}
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.,
\begin{align*}
\boldsymbol{\vartheta}^\star_{-,p,s}(h) &= \mathbf{H}_{p,s}^{-1}(h)\operatorname{\mathbb{E}}\left[\mathbbm{1}(X_i< 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\right]^{-1} \times \\
&\qquad\qquad\qquad\qquad\qquad\qquad \operatorname{\mathbb{E}}\left[\mathbbm{1}(X_i< c) K\left(\frac{X_i-c}{h}\right)\mathbf{r}_{p,s}\left(\frac{X_i-c}{h},\mathbf{W}_i\right)\mu_{Y-}(X_i,\mathbf{W}_i)\right] \\
&=\mathbf{H}_{p,s}^{-1}(h)\operatorname{\mathbb{E}}[\widehat{\boldsymbol{\Gamma}}_{-,p,s}(h)]^{-1}\operatorname{\mathbb{E}}[\widehat{\boldsymbol{\Upsilon}}_{-,p,s}(h)],
\end{align*}
and
\begin{align*}
\boldsymbol{\vartheta}^\star_{+,p,s}(h) &= \mathbf{H}_{p,s}^{-1}(h)\operatorname{\mathbb{E}}\left[\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\right]^{-1} \times \\
&\qquad\qquad\qquad\qquad\qquad\qquad \operatorname{\mathbb{E}}\left[\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)\mu_{Y+}(X_i,\mathbf{W}_i)\right] \\
&=\mathbf{H}_{p,s}^{-1}(h)\operatorname{\mathbb{E}}[\widehat{\boldsymbol{\Gamma}}_{+,p,s}(h)]^{-1}\operatorname{\mathbb{E}}[\widehat{\boldsymbol{\Upsilon}}_{+,p,s}(h)].
\end{align*}
\subsection{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{ass: sampling}, \ref{ass: density running variable}, \ref{ass: DGP}, and \ref{ass: kernel}, whereas Assumption 2 coincides with \ref{ass: conditional expectation}.
\begin{assumption}[Sampling]\label{ass: sampling}
$\{(Y_i,X_i,\mathbf{W}_i^\prime)^\prime\}_{i=1}^n$ is a random sample.
\end{assumption}
\begin{assumption}[Density of Running Variable] \label{ass: density running variable}
The Lebesgue density of $X_i$ is continuous, bounded, and bounded away from zero.
\end{assumption}
\begin{assumption}[DGP]\label{ass: DGP}
For $ t\in \{0,1\}$ and $\mathsf{k}\in \operatorname{\mathbb{R}}_{++}, $ with $-\mathsf{k}<c<\mathsf{k}$ and for all $x\in [-\mathsf{k},\mathsf{k}], \mathbf{w}\in\operatorname{\mathbb{R}}^d:$
\begin{enumerate}[label=(\alph*), ref={SA3\alph*}]
\item $\operatorname{\mathbb{E}}[\mathbf{W}_i\mid X_i=x]$ is continuous and $\operatorname{\mathbb{E}}[\mathbf{W}_i\mathbf{W}_i'\mid X_i=x]$ is continuous and invertible. \label{ass: SRD continuous W}
\item $\operatorname{\mathbb{E}}[Y_i(t)\mid X_i=x, \mathbf{W}_i=\mathbf{w}]$is $q$-times continuously differentiable in $x$ and continuous in $\mathbf{w}$ for $q\geq 1$. \label{ass: SRD continuous Y}
\item $\operatorname{\mathbb{V}}[Y_i(t)|X_i=x]$ and $\operatorname{\mathbb{V}}[Y_i(t)|X_i=x, \mathbf{W}_i=\mathbf{w}]$ are continuous and bounded away from zero. \label{ass: SRD continuous variance}
\item $\operatorname{\mathbb{E}}[|Y_i(t)|^4|X_i=x, \mathbf{W}_i=\mathbf{w}]$ is continuous in both arguments. \label{ass: continuous fourth moment of Y}
\item $\operatorname{\mathbb{E}}[|\mathbf{W}_i|^4|X_i=x]$ is continuous. \label{ass: SRD continuous fourth moment of W}
\item $\operatorname{\mathbb{E}}[\mathbf{W}_i \operatorname{\mathbb{V}}[Y_i(t)|X_i,\mathbf{W}_i]\mid X_i=x]$ and $\operatorname{\mathbb{E}}[\mathbf{W}_i\mathbf{W}_i^\prime\operatorname{\mathbb{V}}[Y_i(t)|X_i,\mathbf{W}_i]\mid X_i=x]$ are continuous. \label{ass: SRD continuous covariances}
\item ${\boldsymbol{\Gamma}}_{-,p,s}$ and ${\boldsymbol{\Gamma}}_{+,p,s}$ in Equation \eqref{eq: plim hessian} are positive definite matrices. \label{ass: pd gram matrix}
\end{enumerate}
\end{assumption}
\begin{assumption}[Identification]
\label{ass: conditional expectation} For $\mathsf{k}\in \operatorname{\mathbb{R}}, $ with $-\mathsf{k}<c<\mathsf{k}$ and for all $x\in [-\mathsf{k},\mathsf{k}], \mathbf{w}\in\operatorname{\mathbb{R}}^d$ the expectation of $Y_i(t), t\in\{0,1\}$ conditional on $(X_i,\mathbf{W}_i)$ is
\begin{align*}
\mu_-(x, \mathbf{w}) = \alpha_{-}(x) + \boldsymbol{\lambda}_{-}(x)^\prime \mathbf{w}, \qquad
\mu_+(x, \mathbf{w}) = \alpha_{+}(x) + \boldsymbol{\lambda}_{+}(x)^\prime \mathbf{w},
\end{align*}
where $\boldsymbol{\lambda}_{-}(x):=(\lambda_{-,1}(x), \cdots, \lambda_{-,d}(x))^\prime, \boldsymbol{\lambda}_{+}(x):=(\lambda_{+,1}(x), \cdots, \lambda_{+,d}(x))^\prime,$ and $\alpha_{-}(x),\alpha_{+}(x)\in\mathcal{C}^{q},$ and $\lambda_{-,\ell}(x),\lambda_{+,\ell}(x),\in\mathcal{C}^{q}$ for all $\ell\in\{1,\cdots,d\}$ for some $q\in\operatorname{\mathbb{N}}$.
\end{assumption}
\begin{assumption}[Kernel] \label{ass: kernel}
The kernel function $k(\cdot):[0,1]\to \mathbb{R}$ is continuous and nonnegative.
\end{assumption}
In words, under Assumption \ref{ass: sampling} the sample is composed of independent draws from an underlying population of interest. Assumptions \ref{ass: density running variable}, \ref{ass: SRD continuous W}, \ref{ass: SRD continuous variance}, and \ref{ass: SRD continuous covariances} 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{ass: pd gram matrix} is a standard assumption in least squares which requires the probability limit of the Gram matrix to be positive definite. Assumption \ref{ass: SRD continuous W} also establishes that $\mathbf{W}$ can be thought of as a vector of pretreatment covariates. Assumption \ref{ass: SRD continuous Y} is necessary to Taylor expand the conditional expectation of the potential outcomes and characterize the smoothing bias. Assumptions \ref{ass: continuous fourth moment of Y} and \ref{ass: SRD continuous fourth moment of W} 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{ass: conditional expectation} 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{ass: kernel} states standard technical conditions on the kernel used in the local polynomial regression.
\subsection{Mapping between Main Text and Supplement}
In the main text, we present all the results using the long regression with $p=s=1$
\begin{align}\label{eq: long regression interacted}
Y_i\quad \text{onto}\quad (\mathbf{r}_p(X_i-c)', T_i\mathbf{r}_p(X_i-c)', \mathbf{W}_i'\otimes\mathbf{r}_s(X_i-c)', T_i\mathbf{W}_i'\otimes\mathbf{r}_s(X_i-c)')',
\end{align}
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
\begin{align}\label{eq: long regression interacted 2}
Y_i\quad \text{onto}\quad((1-T_i)\mathbf{r}_p(X_i-c)', T_i\mathbf{r}_p(X_i-c)', (1-T_i)\mathbf{W}_i'\otimes\mathbf{r}_s(X_i-c)', T_i\mathbf{W}_i'\otimes\mathbf{r}_s(X_i-c)')'.
\end{align}
This is equivalent to \eqref{eq: long regression interacted} in that the predicted values and the residuals of these two regressions are the same. However, the regression coefficients in \eqref{eq: long regression interacted 2} have a different interpretation than those in \eqref{eq: long regression interacted}. 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, \eqref{eq: long regression interacted} 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{subsec: setup}, \eqref{eq: long regression interacted 2} 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 \eqref{eq: long regression interacted} in the main text and on \eqref{eq: long regression interacted 2} in this supplement.
Table \ref{tab: mapping notation} 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}
\begin{bmatrix}
\widehat{\boldsymbol{\vartheta}}_{-, p, s}(h_-)\\
\widehat{\boldsymbol{\vartheta}}_{+, p, s}(h_+)
\end{bmatrix}, \qquad
\mathbf{M} :=
\begin{bmatrix}
-\mathbf{e}_0'\mathbf{I}_{p+1} & \mathbf{0}'_{m} \\
\mathbf{0}_{d\times (p+1)} & -(\boldsymbol{\iota}_d\otimes\mathbf{I}_{s+1}\mathbf{e}_0)' \\
\mathbf{e}_0'\mathbf{I}_{p+1} & \mathbf{0}'_{m} \\
\mathbf{0}_{d\times (p+1)} & (\boldsymbol{\iota}_d\otimes\mathbf{I}_{s+1}\mathbf{e}_0)'
\end{bmatrix}.\]
\begin{table}[!ht]
\centering
\caption{Mapping between main text and supplement notation.}
\begin{tabular}{ccc}
\toprule\toprule
& Main & Supplement \\
\midrule
\textit{Estimands} & & \\
& $\alpha(\cdot)$ & $\alpha_-(\cdot)$ \\
& $\theta(\cdot)$ & $\alpha_+(\cdot)-\alpha_-(\cdot)$ \\
& $\boldsymbol{\lambda}(\cdot)$ & $\boldsymbol{\lambda}_-(\cdot)$ \\
& $\boldsymbol{\xi}(\cdot)$ & $\boldsymbol{\lambda}_+(\cdot)-\boldsymbol{\lambda}_-(\cdot)$ \\
& $\boldsymbol{\varsigma}$ & $(\alpha_+(0)-\alpha_-(0),\boldsymbol{\lambda}_+(0)'-\boldsymbol{\lambda}_-(0)')'$\\
\textit{Estimators} & & \\
& $\widehat{\alpha}$ & $\mathbf{e}_0'\widehat{\boldsymbol{\alpha}}_-$ \\
& $\widehat{\theta}$ & $\mathbf{e}_0'(\widehat{\boldsymbol{\alpha}}_+-\widehat{\boldsymbol{\alpha}}_-)$ \\
& $\widehat{\boldsymbol{\lambda}}$ & $(\boldsymbol{\iota}_d\otimes\mathbf{I}_{s+1}\mathbf{e}_0)'\widehat{\boldsymbol{\lambda}}_-$\\
& $\widehat{\boldsymbol{\xi}}$ & $(\boldsymbol{\iota}_d\otimes\mathbf{I}_{s+1}\mathbf{e}_0)'(\widehat{\boldsymbol{\lambda}}_+-\widehat{\boldsymbol{\lambda}}_-)$\\
& $\widehat{\boldsymbol{\varsigma}}$ & $\breve{\mathbf{e}}_0'(\widehat{\boldsymbol{\vartheta}}_{+,1,1}-\widehat{\boldsymbol{\vartheta}}_{-,1,1})$\\
& $\widehat{\omega}_1$ & $\mathbf{e}_1'\widehat{\boldsymbol{\alpha}}_-$ \\
& $\widehat{\omega}_2$ & $\mathbf{e}_1'(\widehat{\boldsymbol{\alpha}}_+-\widehat{\boldsymbol{\alpha}}_-)$ \\
& $\widehat{\boldsymbol{\omega}}_3$ & $(\boldsymbol{\iota}_d\otimes\mathbf{I}_{s+1}\mathbf{e}_1)'\widehat{\boldsymbol{\lambda}}_-$ \\
& $\widehat{\boldsymbol{\omega}}_4$ & $(\boldsymbol{\iota}_d\otimes\mathbf{I}_{s+1}\mathbf{e}_1)'(\widehat{\boldsymbol{\lambda}}_+-\widehat{\boldsymbol{\lambda}}_-)$ \\
\textit{Constants} & & \\
& $\mathsf{B}_{\boldsymbol{\mathsf{s}}}$ & $\mathsf{s}'\mathbf{M}(\mathbf{B}_{+,1,1}^{[0]} + \mathbf{B}_{+,1,1}^{[1]}-\mathbf{B}_{-,1,1}^{[0]} + \mathbf{B}_{-,1,1}^{[1]})$ \\
& $\mathsf{V}_{\boldsymbol{\mathsf{s}}}$ & $\mathsf{s}'\mathbf{M}({\boldsymbol{\Gamma}}_{+,1,1}^{-1} \mathbf{V}_{+,1,1}
{{\boldsymbol{\Gamma}}_{+,1,1}^{-1\prime}} + {\boldsymbol{\Gamma}}_{-,1,1}^{-1} \mathbf{V}_{-,1,1}
{{\boldsymbol{\Gamma}}_{-,1,1}^{-1\prime}})\mathbf{M}'\mathsf{s}$ \\
\bottomrule\bottomrule
\end{tabular}
\label{tab: mapping notation}
\end{table}
We conclude this section by mapping the results in the main text with the proofs in the supplement:
\begin{itemize}
\item Theorem 1 is a particular case of the results proved in Section \ref{subsec: consistency} with $p=s=1$ and $\nu=0$;
\item Theorem 2 follows from the results in Section \ref{subsec: mse expansions} with $p=s=1, \nu=0,$ and $\breve{\mathbf{e}}_\nu=\mathbf{s}$;
\item Theorem 3 is a particular case of Corollary \ref{cor: asympt distr under sa5} with $p=s=1$ and $\nu=0$.
\end{itemize}
\subsection{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.
\begin{lemma}\label{lemma: markov}
Let $\mathbf{A}_n$ be a sequence of random matrices with finite first two moments. Then
\[\mathbf{A}_n = \operatorname{\mathbb{E}}[\mathbf{A}_n] + O_{\mathbb{P}}(|\operatorname{\mathbb{V}}[\mathbf{A}_n]|^{1/2}).\]
\end{lemma}
\begin{flushright}
[\hyperref[proof: markov]{Proof}]
\end{flushright}
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.
\begin{lemma}\label{lemma: hessian}
Let Assumptions \ref{ass: sampling}, \ref{ass: density running variable}, and \ref{ass: kernel} hold with $\mathsf{k} \geq h$.
If $nh\to\infty$ and $h\to 0,$ then
$$
\widehat{\boldsymbol{\Lambda}}_{-, p, s}(h)=\widetilde{\boldsymbol{\Lambda}}_{-, p, s}(h)+O_{\mathbb{P}}(1/\sqrt{nh}), \qquad \widehat{\boldsymbol{\Lambda}}_{+, p, s}(h)=\widetilde{\boldsymbol{\Lambda}}_{+, p, s}(h)+O_{\mathbb{P}}(1/\sqrt{nh}),
$$
with
\begin{align*}
\widetilde{\boldsymbol{\Lambda}}_{-, p, s}(h)&:= \mathbb{E}[\widehat{\boldsymbol{\Lambda}}_{-,p, s}(h)] = \int\limits_{-\infty}^0K(u)\mathbf{r}_{p}(u)\mathbf{r}_s(u)^\prime f(uh+c)\,\mathrm{d} u, \\
\widetilde{\boldsymbol{\Lambda}}_{+, p, s}(h) &:= \mathbb{E}[\widehat{\boldsymbol{\Lambda}}_{+, p, s}(h)] =\int\limits_0^\infty K(u)\mathbf{r}_{p}(u)\mathbf{r}_s(u)^\prime f(uh+c)\,\mathrm{d} u.
\end{align*}
\end{lemma}
\begin{flushright}
[\hyperref[proof: hessian]{Proof}]
\end{flushright}
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.
\begin{lemma}
\label{lemma: grand hessian}
Let Assumptions \ref{ass: sampling}, \ref{ass: density running variable}, \ref{ass: DGP}, and \ref{ass: kernel} hold with $\mathsf{k} \geq h$.
If $nh\to\infty$ and $h\to0$, then
\begin{align*}
\widehat{\boldsymbol{\Gamma}}_{-,p,s}(h) &= \widetilde{\boldsymbol{\Gamma}}_{-,p,s}(h) + O_{\mathbb{P}}(1/\sqrt{nh}),\qquad
\widehat{\boldsymbol{\Gamma}}_{+,p,s}(h) = \widetilde{\boldsymbol{\Gamma}}_{+,p,s}(h) + O_{\mathbb{P}}(1/\sqrt{nh}),
\end{align*}
where
\begin{align*}
\widetilde{\boldsymbol{\Gamma}}_{-,p,s}(h):=
\begin{bmatrix}
\widetilde{\boldsymbol{\Lambda}}_{-,p,p}(h) & \widetilde{\mathbf{G}}_{-,2}(h) \\
\widetilde{\mathbf{G}}_{-,2}(h)^\prime & \widetilde{\mathbf{G}}_{-, 3}(h)
\end{bmatrix}, \qquad
\widetilde{\boldsymbol{\Gamma}}_{+,p,s}(h):=\begin{bmatrix}
\widetilde{\boldsymbol{\Lambda}}_{+,p,p}(h) & \widetilde{\mathbf{G}}_{+,2}(h) \\
\widetilde{\mathbf{G}}_{+,2}(h)^\prime & \widetilde{\mathbf{G}}_{+, 3}(h)
\end{bmatrix},
\end{align*}
where
\begin{align*}
\widetilde{\mathbf{G}}_{-,2}(h)&=\int\limits_{-\infty}^0 K\left(u\right)\left[ \boldsymbol{\mu}_{W}(uh+c)^\prime\otimes \mathbf{r}_p\left(u\right)\mathbf{r}_s\left(u\right)^\prime\right]f(uh+c)\,\mathrm{d} u,\\
\widetilde{\mathbf{G}}_{+,2}(h)&=\int\limits_{0}^\infty K(u)\left[ \boldsymbol{\mu}_{W}(uh+c)^\prime\otimes \mathbf{r}_p\left(u\right)\mathbf{r}_s\left(u\right)^\prime\right]f(uh+c)\,\mathrm{d} u,\\
\widetilde{\mathbf{G}}_{-,3}(h)&=\int\limits_{-\infty}^0 K\left(u\right)\left[ \boldsymbol{\mu}_{WW}(uh+c)\otimes \mathbf{r}_s\left(u\right)\mathbf{r}_s\left(u\right)^\prime\right]f(uh+c)\,\mathrm{d} u,\\
\widetilde{\mathbf{G}}_{+,3}(h)&=\int\limits_{0}^\infty K\left(u\right) \left[\boldsymbol{\mu}_{WW}(uh+c)\otimes \mathbf{r}_s\left(u\right)\mathbf{r}_s\left(u\right)^\prime\right]f(uh+c)\,\mathrm{d} u.
\end{align*}
Last,
\[\widetilde{\boldsymbol{\Gamma}}_{-,p,s}(h) = \boldsymbol{\Gamma}_{-,p,s}\{1+o(1)\}, \qquad \widetilde{\boldsymbol{\Gamma}}_{+,p,s}(h) = \boldsymbol{\Gamma}_{+,p,s}\{1+o(1)\}, \]
where
\begin{align}
\label{eq: plim hessian}
\boldsymbol{\Gamma}_{-,p,s} = \begin{bmatrix}
\boldsymbol{\Lambda}_{-,p,p} & \mathbf{G}_{-,2} \\
\mathbf{G}_{-,2}' & \mathbf{G}_{-,3}
\end{bmatrix},\qquad \boldsymbol{\Gamma}_{+,p,s} = \begin{bmatrix}
\boldsymbol{\Lambda}_{+,p,p} & \mathbf{G}_{+,2} \\
\mathbf{G}_{+,2}' & \mathbf{G}_{+,3}
\end{bmatrix},
\end{align}
with
\begin{alignat*}{3}
{\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 \\
{\mathbf{G}}_{-,2} &= \boldsymbol{\mu}_{W}^\prime \otimes {\boldsymbol{\Lambda}}_{-,p,s},\qquad &{\mathbf{G}}_{+,2} &= \boldsymbol{\mu}_{W}^\prime \otimes {\boldsymbol{\Lambda}}_{+,p,s}, \\
{\mathbf{G}}_{-,3} &= \boldsymbol{\mu}_{WW} \otimes {\boldsymbol{\Lambda}}_{-,s,s},\qquad &{\mathbf{G}}_{+,3} &= \boldsymbol{\mu}_{WW} \otimes {\boldsymbol{\Lambda}}_{+,s,s}.
\end{alignat*}
\end{lemma}
\begin{flushright}
[\hyperref[proof: grand hessian]{Proof}]
\end{flushright}
To have a more compact notation, we introduce the following quantities for integers $p,s,a\in \operatorname{\mathbb{N}}_0$:
\begin{align*}
\widehat{\boldsymbol{\zeta}}_{-,p,s,a}(h) &= \frac{1}{nh}\sum_{i=1}^n\mathbbm{1}(X_i< c)K\left(\frac{X_i-c}{h}\right)\mathbf{r}_{p,s}\left(\frac{X_i-c}{h},\mathbf{W}_i\right)\left(\frac{X_i-c}{h}\right)^{a+1}, \\
\widehat{\boldsymbol{\zeta}}_{+,p,s,a}(h) &= \frac{1}{nh}\sum_{i=1}^n\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)\left(\frac{X_i-c}{h}\right)^{a+1}, \\
\widehat{\boldsymbol{\varphi}}_{-,p,s,a}(h) &= \frac{1}{nh}\sum_{i=1}^n\mathbbm{1}(X_i< c)K\left(\frac{X_i-c}{h}\right)\mathbf{r}_{p,s}\left(\frac{X_i-c}{h},\mathbf{W}_i\right)\mathbf{W}_i^\prime\left(\frac{X_i-c}{h}\right)^{a+1},\\
\widehat{\boldsymbol{\varphi}}_{+,p,s,a}(h) &= \frac{1}{nh}\sum_{i=1}^n\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{W}_i^\prime\left(\frac{X_i-c}{h}\right)^{a+1}.
\end{align*}
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.
\begin{lemma}\label{lemma: vartheta hat order}
Let Assumptions \ref{ass: sampling}, \ref{ass: density running variable}, \ref{ass: DGP}, and \ref{ass: kernel} hold with $\mathsf{k} \geq h$. If $nh\to\infty$ and $h\to 0 $, then
\begin{alignat*}{3}
\widehat{\boldsymbol{\zeta}}_{-,p,s,a}(h) &= \widetilde{\boldsymbol{\zeta}}_{-,p,s,a}(h) + O_{\mathbb{P}}(1/\sqrt{nh}), \qquad &\widehat{\boldsymbol{\zeta}}_{+,p,s,a}(h) &= \widetilde{\boldsymbol{\zeta}}_{+,p,s,a}(h) + O_{\mathbb{P}}(1/\sqrt{nh}),\\
\widehat{\boldsymbol{\varphi}}_{-,p,s,a}(h) &= \widetilde{\boldsymbol{\varphi}}_{-,p,s,a}(h) + O_{\mathbb{P}}(1/\sqrt{nh}), \qquad &\widehat{\boldsymbol{\varphi}}_{+,p,s,a}(h) &= \widetilde{\boldsymbol{\varphi}}_{+,p,s,a}(h) + O_{\mathbb{P}}(1/\sqrt{nh}),
\end{alignat*}
with
\begin{align*}
\widetilde{\boldsymbol{\zeta}}_{-,p,s,a}(h)&:=\operatorname{\mathbb{E}}[\widehat{\boldsymbol{\zeta}}_{-,p,s,a}(h)] =\int\limits_{-\infty}^0 K(u)\mathbf{r}_{p,s}(u,\boldsymbol{\mu}_{W}(uh+c))u^{a+1}f(uh+c)\,\mathrm{d} u, \\
\widetilde{\boldsymbol{\zeta}}_{+,p,s,a}(h)&:=\operatorname{\mathbb{E}}[\widehat{\boldsymbol{\zeta}}_{+,p,s,a}(h)]=\int\limits_0^\infty K(u)\mathbf{r}_{p,s}(u,\boldsymbol{\mu}_{W}(uh+c))u^{a+1} f(uh+c)\,\mathrm{d} u,\\
\widetilde{\boldsymbol{\varphi}}_{-,p,s,a}(h)&:= E[\widehat{\boldsymbol{\varphi}}_{-,p,s,a}(h)]=\int\limits_{-\infty}^0 K(u) \begin{bmatrix}
\mathbf{r}_p(u)\boldsymbol{\mu}_{W}(uh+c)^\prime \\
\boldsymbol{\mu}_{WW}(uh+c)\otimes \mathbf{r}_s(u)
\end{bmatrix} u^{a+1} f(uh+c)\,\mathrm{d} u, \\
\widetilde{\boldsymbol{\varphi}}_{+,p,s,a}(h)&:= \operatorname{\mathbb{E}}[\widehat{\boldsymbol{\varphi}}_{+,p,s,a}(h)]=\int\limits_0^\infty K(u) \begin{bmatrix}
\mathbf{r}_p(u)\boldsymbol{\mu}_{W}(uh+c)^\prime \\
\boldsymbol{\mu}_{WW}(uh+c)\otimes \mathbf{r}_s(u)
\end{bmatrix} u^{a+1} f(uh+c)\,\mathrm{d} u.
\end{align*}
Last
\begin{alignat*}{3}
\widetilde{\boldsymbol{\zeta}}_{-,p,s,a}(h)&={\boldsymbol{\zeta}}_{-,p,s,a}\{1+o(1)\}, \qquad &\widetilde{\boldsymbol{\zeta}}_{+,p,s,a}(h)&={\boldsymbol{\zeta}}_{+,p,s,a}\{1+o(1)\}, \\
\widetilde{\boldsymbol{\varphi}}_{-,p,s,a}(h)&={\boldsymbol{\varphi}}_{-,p,s,a}\{1+o(1)\}, \qquad &\widetilde{\boldsymbol{\varphi}}_{+,p,s,a}(h)&={\boldsymbol{\varphi}}_{+,p,s,a}\{1+o(1)\},
\end{alignat*}
where
\begin{alignat*}{3}
{\boldsymbol{\zeta}}_{-,p,s,a} &= f(c)\int\limits_{-\infty}^0 K(u)\mathbf{r}_{p,s}(u,\boldsymbol{\mu}_{W})u^{a+1} \,\mathrm{d} u, &{\boldsymbol{\zeta}}_{+,p,s,a}&=f(c)\int\limits_0^\infty K(u)\mathbf{r}_{p,s}(u,\boldsymbol{\mu}_{W})u^{a+1} \,\mathrm{d} u, \\
{\boldsymbol{\varphi}}_{-,p,s,a} &= f(c)\int\limits_{-\infty}^0 K(u) \begin{bmatrix}
\mathbf{r}_p(u)\boldsymbol{\mu}_{W}^\prime \\
\boldsymbol{\mu}_{WW}\otimes \mathbf{r}_s(u)
\end{bmatrix} u^{a+1} \,\mathrm{d} u, &\quad
{\boldsymbol{\varphi}}_{+,p,s,a}&=f(c)\int\limits_0^\infty K(u)\begin{bmatrix}
\mathbf{r}_p(u)\boldsymbol{\mu}_{W}^\prime \\
\boldsymbol{\mu}_{WW}\otimes \mathbf{r}_s(u)
\end{bmatrix}u^{a+1} \,\mathrm{d} u.
\end{alignat*}
\end{lemma}
\begin{flushright}
[\hyperref[proof: vartheta hat order]{Proof}]
\end{flushright}
\section{Main Results}\label{sec: 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:
\begin{align*}
\widehat{\boldsymbol{\vartheta}}_{+,p,s}(h)&= \mathbf{H}_{p,s}^{-1}(h)\widehat{\boldsymbol{\Gamma}}_{+,p,s}^{-1}(h) \frac{1}{nh}\sum_{i=1}^n\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)Y_i, \\
&= \boldsymbol{\vartheta}_{+,p,s}^\star(h) + \mathbf{H}_{p,s}^{-1}(h)\widehat{\boldsymbol{\Gamma}}_{+,p,s}^{-1}(h) \frac{1}{nh}\sum_{i=1}^n\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)u_{+,i}(h),
\end{align*}
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
\begin{align*}
\widehat{\boldsymbol{\vartheta}}_{+,p,s}(h)- \boldsymbol{\vartheta}^\star_{+,p,s}(h) = \mathbf{H}_{p,s}^{-1}(h)\widehat{\boldsymbol{\Gamma}}_{+,p,s}^{-1}(h)\mathbf{L}_{+,p,s}(h),
\end{align*}
where
\begin{align*}
\mathbf{L}_{+,p,s}(h)&=\frac{1}{nh}\sum_{i=1}^n\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)u_{+,i}(h).
\end{align*}
\subsection{Asymptotic Approximation and Asymptotic Variance}
Lemma \ref{lemma: grand hessian} 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{lemma: markov}.
\begin{lemma}
\label{lemma: first order variance}
Let Assumptions \ref{ass: sampling}, \ref{ass: density running variable}, \ref{ass: DGP}, and \ref{ass: kernel} hold with $\mathsf{k} \geq h$. Then, for fixed $h>0$
\begin{align*}
\operatorname{\mathbb{V}}[\mathbf{L}_{-,p,s}(h)] = \frac{1}{nh}\mathbf{V}_{-,p,s}(h), \qquad \operatorname{\mathbb{V}}[\mathbf{L}_{+,p,s}(h)] = \frac{1}{nh}\mathbf{V}_{+,p,s}(h).
\end{align*}
Furthermore, as $nh\to\infty$ and $h\to 0$
\begin{alignat*}{3}
\mathbf{V}_{-,p,s}(h) &=\mathbf{V}_{-,p,s}\{1+o(1)\}, \quad &\mathbf{V}_{+,p,s}(h) &=\mathbf{V}_{+,p,s}\{1+o(1)\},
\end{alignat*}
where $|\mathbf{V}_{-,p,s}|<\infty$ and $|\mathbf{V}_{+,p,s}|<\infty$.
\end{lemma}
\begin{flushright}
[\hyperref[proof: first order variance]{Proof}]
\end{flushright}
Typically, we are interested in estimators of the form $\breve{\mathbf{e}}_{\nu}^\prime \widehat{\boldsymbol{\vartheta}}_{+,p,s}(h_+)$. By Lemma \ref{lemma: grand hessian} and Lemma \ref{lemma: first order variance} 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.
\begin{theorem}\label{thm: asymptotic distribution}
Let Assumptions \ref{ass: sampling}, \ref{ass: density running variable}, \ref{ass: DGP}, and \ref{ass: kernel} hold with $\mathsf{k} \geq h$. If $nh\to \infty $ and $ h\to 0$, then
\begin{align*}
\sqrt{nh}\mathbf{H}_{p,s}(h)(\widehat{\boldsymbol{\vartheta}}_{-,p,s}(h) - \boldsymbol{\vartheta}^\star_{-,p,s}(h)) &\rightsquigarrow \mathsf{N}(\mathbf{0}_{1+p+d},\boldsymbol{\Omega}_{-,p,s}), \\
\sqrt{nh}\mathbf{H}_{p,s}(h)(\widehat{\boldsymbol{\vartheta}}_{+,p,s}(h) - \boldsymbol{\vartheta}^\star_{+,p,s}(h)) &\rightsquigarrow \mathsf{N}(\mathbf{0}_{1+p+d},\boldsymbol{\Omega}_{+,p,s}),
\end{align*}
with
\begin{align*}
\boldsymbol{\Omega}_{-,p,s}&={\boldsymbol{\Gamma}}_{-,p,s}^{-1}\mathbf{V}_{-,p,s}
{{\boldsymbol{\Gamma}}_{-,p,s}^{-1\prime}}, \qquad
\boldsymbol{\Omega}_{+,p,s}={\boldsymbol{\Gamma}}_{+,p,s}^{-1}\mathbf{V}_{+,p,s}
{{\boldsymbol{\Gamma}}_{+,p,s}^{-1\prime}},
\end{align*}
where $\boldsymbol{\Omega}_{-,p,s}$ and $\boldsymbol{\Omega}_{+,p,s}$ are positive definite matrices.
\end{theorem}
\begin{flushright}
[\hyperref[proof: asymptotic distribution]{Proof}]
\end{flushright}
\begin{corollary}
Let the conditions in Theorem \ref{thm: asymptotic distribution} hold, then for $\nu\in\{0,\cdots,p\wedge s\}$:
\begin{align*}
\sqrt{nh^{2\nu+1}}\breve{\mathbf{e}}_\nu^\prime(\widehat{\boldsymbol{\vartheta}}_{-,p,s}(h) - \boldsymbol{\vartheta}^\star_{-,p,s}(h)) &\rightsquigarrow \mathsf{N}({0},\mathcal{V}_{-,\nu,p,s}), \\
\sqrt{nh^{2\nu+1}}\breve{\mathbf{e}}_\nu^\prime(\widehat{\boldsymbol{\vartheta}}_{+,p,s}(h) - \boldsymbol{\vartheta}^\star_{+,p,s}(h)) &\rightsquigarrow \mathsf{N}({0},\mathcal{V}_{+,\nu,p,s}),
\end{align*}
where
\[\mathcal{V}_{-,\nu,p,s}:= \breve{\mathbf{e}}_\nu^\prime \boldsymbol{\Omega}_{-,p,s}\breve{\mathbf{e}}_\nu,\qquad\mathcal{V}_{+,\nu,p,s}:= \breve{\mathbf{e}}_\nu^\prime \boldsymbol{\Omega}_{+,p,s}\breve{\mathbf{e}}_\nu.\]
\end{corollary}
\subsection{Variance Estimation}
To estimate the asymptotic variance of the RD estimator, we propose the following \textit{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,
\begin{table}[!ht]\centering
\begin{tabular}{ccccc}
\toprule & $\mathrm{HC}_0$ & $\mathrm{HC}_1$ & $\mathrm{HC}_2$ & $\mathrm{HC}_3$ \\
\midrule$w_{-, i}(h)$ & 1 & $\frac{N_{-}}{N_{-}-2 \operatorname{tr}\left(\mathbf{Q}_{-, p,s}\right)+\operatorname{tr}\left(\mathbf{Q}_{-, p,s} \mathbf{Q}_{-, p,s}\right)}$ & $\frac{1}{1-L_{+,i}}$ & $\frac{1}{\left(1-L_{+,i}\right)^2}$ \\
$w_{+, i}(h)$ & 1 & $\frac{N_{+}}{N_{+}-2 \operatorname{tr}\left(\mathbf{Q}_{+, p,s}\right)+\operatorname{tr}\left(\mathbf{Q}_{+, p,s} \mathbf{Q}_{+, p,s}\right)}$ & $\frac{1}{1-L_{+,i}}$ & $\frac{1}{\left(1-L_{+,i}\right)^2}$ \\
\bottomrule
\end{tabular}
\end{table}
where
\begin{align*}
N_{-}:=\sum_{i=1}^n \mathbbm{1}\left(X_i<c\right) \quad \text { and } \quad N_{+}:=\sum_{i=1}^n \mathbbm{1}\left(X_i \geq c\right),
\end{align*}
$\mathbf{Q}_{-, p,s}, $ and $ \mathbf{Q}_{+, p,s}$ are the ``projection" matrices used to get the estimated residuals and are defined as,
\begin{align*}
\mathbf{Q}_{-, p,s}:=\mathbf{R}_{p,s}(h) \widehat{\boldsymbol{\Gamma}}_{-, p,s}^{-1} \mathbf{R}_{p,s}(h)^{\prime} \mathbf{K}_{-}(h) / n, \quad \mathbf{Q}_{+, p,s}=\mathbf{R}_{p,s}(h) \widehat{\boldsymbol{\Gamma}}_{+, p,s}^{-1} \mathbf{R}_{p,s}(h)^{\prime} \mathbf{K}_{+}(h) / n,
\end{align*}
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.
\begin{theorem}
\label{theorem: variance estimation plug-in}
Let the assumptions of Theorem \ref{thm: asymptotic distribution} hold, then
\begin{align*}
\widehat{\mathcal{V}}_{-,\nu,p,s}(h) \overset{\mathbb{P}}{\longrightarrow} {\mathcal{V}}_{-,\nu,p,s}, \qquad \widehat{\mathcal{V}}_{+,\nu,p,s}(h) \overset{\mathbb{P}}{\longrightarrow}{\mathcal{V}}_{+,\nu,p,s}.
\end{align*}
\end{theorem}
\begin{flushright}
[\hyperref[proof: variance estimation plug-in]{Proof}]
\end{flushright}
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.
\begin{corollary}
Let Assumptions \ref{ass: sampling}, \ref{ass: density running variable}, \ref{ass: DGP}, and \ref{ass: kernel} hold with $\mathsf{k} \geq h$. If $nh^{2\nu+1}\to \infty$ and $h\to 0$ then for $\nu\in\{0,1,\ldots,p\wedge s\}$
\begin{align*}
\widehat{\sigma}_{-,\nu,p,s}^{-1}(h)\breve{\mathbf{e}}_{\nu}^\prime(\widehat{\boldsymbol{\vartheta}}_{-,p,s}(h) - \boldsymbol{\vartheta}^\star_{-,p,s}(h)) &\rightsquigarrow \mathsf{N}(0,1), \\
\widehat{\sigma}_{+,\nu,p,s}^{-1}(h)\breve{\mathbf{e}}_{\nu}^\prime(\widehat{\boldsymbol{\vartheta}}_{+,p,s}(h) - \boldsymbol{\vartheta}^\star_{+,p,s}(h)) &\rightsquigarrow \mathsf{N}(0,1).
\end{align*}
\end{corollary}
\subsubsection*{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 \cite{Cameron2015_JHR} and \cite{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 \textit{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
\begin{align*}
\widehat{\mathbf{V}}^{\mathsf{CL}}_{+,p,s}(h)&:= \frac{1}{Gh}\sum_{g=1}^G\sum_{i,j=1}^n\mathbbm{1}(s(i)=s(j))\cdot\mathbbm{1}(X_i\geq c, X_j\geq c)\times \\
&\qquad\quad K\left(\frac{X_i-c}{h}\right)K\left(\frac{X_j-c}{h}\right)\mathbf{r}_{p,s}\left(\frac{X_i-c}{h},\mathbf{W}_i\right)\mathbf{r}_{p,s}\left(\frac{X_j-c}{h},\mathbf{W}_j\right)^\prime {\widehat{u}_{+,i}}(h){\widehat{u}_{+,j}}(h).
\end{align*}
This estimator, as well as many more, is implemented in our \texttt{R} and \texttt{Stata} software.
\subsection{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{ass: conditional expectation}, which we state again here below for the reader's convenience.
\addtocounter{assumption}{-2}
\begin{assumption}
The expectation of $Y_i(t), t\in\{0,1\}$ conditional on $(X_i,\mathbf{W}_i)$ is
\begin{align*}
\mu_-(X_i, \mathbf{W}_i) = \alpha_{-}(X_i) + \boldsymbol{\lambda}_{-}(X_i)^\prime \mathbf{W}_i, \qquad
\mu_+(X_i, \mathbf{W}_i) = \alpha_{+}(X_i) + \boldsymbol{\lambda}_{+}(X_i)^\prime \mathbf{W}_i,
\end{align*}
where $\boldsymbol{\lambda}_{-}(x):=(\lambda_{-,1}(x), \cdots, \lambda_{-,d}(x))^\prime, \boldsymbol{\lambda}_{+}(x):=(\lambda_{+,1}(x), \cdots, \lambda_{+,d}(x))^\prime,$ and $\alpha_{-}(x),\alpha_{+}(x)\in\mathcal{C}^{q+2},$ and $\lambda_{-,\ell}(x),\lambda_{+,\ell}(x),\in\mathcal{C}^{q+2}$ for all $\ell\in\{1,\cdots,d\}$ for some $q\in\operatorname{\mathbb{N}}$.
\end{assumption}
Let
\begin{align}\label{eq: xi star under partial linear model}
\boldsymbol{\vartheta}_{-,p,s} := \begin{bmatrix}
\boldsymbol{\alpha}_{-,p}(c) \\
\boldsymbol{\lambda}_{-,s}(c)
\end{bmatrix}, \qquad
\boldsymbol{\vartheta}_{+,p,s} := \begin{bmatrix}
\boldsymbol{\alpha}_{+,p}(c) \\
\boldsymbol{\lambda}_{+,s}(c)
\end{bmatrix},
\end{align}
where
\begin{align*}
\boldsymbol{\alpha}_{-,p}(x) &= \begin{bmatrix}
\alpha_{-}(x) \\ \alpha_{-}^{(1)}(x) \\ \vdots \\
\frac{\alpha_{-}^{(p)}(x)}{p!}
\end{bmatrix}, \quad &\boldsymbol{\lambda}_{-,s}(x) &=
\begin{bmatrix}
\lambda_{-,1}(x) \\ \vdots \\
\frac{\lambda_{-,1}^{(s)}(x)}{s!} \\ \vdots \\
\lambda_{-,d}(x) \\ \vdots \\
\frac{\lambda_{-,d}^{(s)}(x)}{s!}
\end{bmatrix}, \quad
\boldsymbol{\alpha}_{+,p}(x) &= \begin{bmatrix}
\alpha_{+}(x) \\ \alpha_{+}^{(1)}(x) \\ \vdots \\
\frac{\alpha_{+}^{(p)}(x)}{p!}
\end{bmatrix}, \quad &\boldsymbol{\lambda}_{+,s}(x) &=
\begin{bmatrix}
\lambda_{+,1}(x) \\ \vdots \\
\frac{\lambda_{+,1}^{(s)}(x)}{s!} \\ \vdots \\
\lambda_{+,d}(x) \\ \vdots \\
\frac{\lambda_{+,d}^{(s)}(x)}{s!}
\end{bmatrix},
\end{align*}
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
\begin{align*}
\mathbf{H}_{p,s}(h)(\widehat{\boldsymbol{\vartheta}}_{+,p,s}(h) - \boldsymbol{\vartheta}_{+,p,s}) = \mathbf{H}_{p,s}(h)(\widehat{\boldsymbol{\vartheta}}_{+,p,s}(h) - {\boldsymbol{\vartheta}}^\star_{+,p,s}(h)) + \mathbf{H}_{p,s}(h)({\boldsymbol{\vartheta}}^\star_{+,p,s}(h) - \boldsymbol{\vartheta}_{+,p,s}).
\end{align*}
From Theorem \ref{thm: asymptotic distribution}, 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})$.
\begin{lemma}[Smoothing Bias]\label{lemma: smoothing bias}
Let Assumptions \ref{ass: sampling}-\ref{ass: kernel} hold for some $q\geq 2 + p\vee s$ and with $\mathsf{k} \geq h$. If $nh\to 0$ and $h\to 0$, then
\begin{align*}
{\boldsymbol{\vartheta}}_{-,p,s}^\star(h) &- \boldsymbol{\vartheta}_{-,p,s} =\\
&\mathbf{H}_{p,s}^{-1}(h)\left(h^{1+p}\widetilde{\mathbf{B}}^{[0]}_{-,p}(h) + h^{2+p}\widetilde{\mathbf{B}}^{[0]}_{-,p+1}(h) + h^{1+s}\widetilde{\mathbf{B}}^{[1]}_{-,s}(h) + h^{2+s}\widetilde{\mathbf{B}}^{[1]}_{-,s+1}(h) + o(h^{2+p\wedge s})\right),\\
{\boldsymbol{\vartheta}}_{+,p,s}^\star(h) &- \boldsymbol{\vartheta}_{+,p,s} =\\
&\mathbf{H}_{p,s}^{-1}(h)\left(h^{1+p}\widetilde{\mathbf{B}}^{[0]}_{+,p}(h) + h^{2+p}\widetilde{\mathbf{B}}^{[0]}_{+,p+1}(h) + h^{1+s}\widetilde{\mathbf{B}}^{[1]}_{+,s}(h) + h^{2+s}\widetilde{\mathbf{B}}^{[1]}_{+,s+1}(h) + o(h^{2+p\wedge s})\right),\\
\end{align*}
where
\begin{alignat*}{3}
\widetilde{\mathbf{B}}^{[0]}_{-,a}(h) =\widetilde{\boldsymbol{\Gamma}}_{-,p,s}^{-1}(h)\widetilde{\boldsymbol{\zeta}}_{-,p,s,a}(h) \frac{\alpha^{(a+1)}_{-}(c)}{(a+1)!}\to {\boldsymbol{\Gamma}}_{-,p,s}^{-1}{\boldsymbol{\zeta}}_{-,p,s,a} \frac{\alpha^{(a+1)}_{+}(c)}{(a+1)!}=:\mathbf{B}^{[0]}_{-,a}, \\
\widetilde{\mathbf{B}}^{[0]}_{+,a}(h) =\widetilde{\boldsymbol{\Gamma}}_{+,p,s}^{-1}(h)\widetilde{\boldsymbol{\zeta}}_{+,p,s,a}(h) \frac{\alpha^{(a+1)}_{+}(c)}{(a+1)!}\to {\boldsymbol{\Gamma}}_{+,p,s}^{-1}{\boldsymbol{\zeta}}_{+,p,s,a} \frac{\alpha^{(a+1)}_{+}(c)}{(a+1)!}=:\mathbf{B}^{[0]}_{+,a}, \\
\widetilde{\mathbf{B}}^{[1]}_{-,a}(h) = \widetilde{\boldsymbol{\Gamma}}_{-,p,s}^{-1}(h)\widetilde{\boldsymbol{\varphi}}_{-,p,s,a}(h) \frac{\boldsymbol{\lambda}^{(a+1)}_{-}(c)}{(a+1)!}\to {\boldsymbol{\Gamma}}_{-,p,s}^{-1}{\boldsymbol{\varphi}}_{-,p,s,a}\frac{\boldsymbol{\lambda}^{(a+1)}_{-}(c)}{(a+1)!}=:\mathbf{B}^{[1]}_{-,a}, \\
\widetilde{\mathbf{B}}^{[1]}_{+,a}(h) = \widetilde{\boldsymbol{\Gamma}}_{+,p,s}^{-1}(h)\widetilde{\boldsymbol{\varphi}}_{+,p,s,a}(h) \frac{\boldsymbol{\lambda}^{(a+1)}_{+}(c)}{(a+1)!}\to {\boldsymbol{\Gamma}}_{+,p,s}^{-1}{\boldsymbol{\varphi}}_{+,p,s,a}\frac{\boldsymbol{\lambda}^{(a+1)}_{+}(c)}{(a+1)!}=:\mathbf{B}^{[1]}_{+,a}.
\end{alignat*}
\end{lemma}
\begin{flushright}
[\hyperref[proof: smoothing bias]{Proof}]
\end{flushright}
The following corollary follows naturally from Theorem \ref{thm: asymptotic distribution} and Lemma \ref{lemma: smoothing bias}.
\begin{corollary}\label{cor: asympt distr under sa5}
Let the assumptions of Lemma \ref{lemma: smoothing bias} hold. If $nh^{2\nu+1}\to \infty, h\to 0,$ and $nh^{2(p\wedge s) + 3}\to 0$, then
\begin{align*}
\sqrt{nh^{2\nu+1}}\breve{\mathbf{e}}_{\nu}^\prime(\widehat{\boldsymbol{\vartheta}}_{-,p,s}(h) - \boldsymbol{\vartheta}_{-,p,s})&\rightsquigarrow \mathsf{N}(0,\mathcal{V}_{-,\nu,p,s}), \\
\sqrt{nh^{2\nu+1}}\breve{\mathbf{e}}_{\nu}^\prime(\widehat{\boldsymbol{\vartheta}}_{+,p,s}(h) - \boldsymbol{\vartheta}_{+,p,s})&\rightsquigarrow \mathsf{N}(0,\mathcal{V}_{+,\nu,p,s}),
\end{align*}
where
\[\mathcal{V}_{-,\nu,p,s}:= \breve{\mathbf{e}}_\nu^\prime \boldsymbol{\Omega}_{-,p,s}\breve{\mathbf{e}}_\nu,\qquad\mathcal{V}_{+,\nu,p,s}:= \breve{\mathbf{e}}_\nu^\prime \boldsymbol{\Omega}_{+,p,s}\breve{\mathbf{e}}_\nu,\]
with
\begin{align*}
\boldsymbol{\Omega}_{-,p,s}&={\boldsymbol{\Gamma}}_{-,p,s}^{-1}\mathbf{V}_{-,p,s}
{{\boldsymbol{\Gamma}}_{-,p,s}^{-1\prime}}, \qquad
\boldsymbol{\Omega}_{+,p,s}={\boldsymbol{\Gamma}}_{+,p,s}^{-1}\mathbf{V}_{+,p,s}
{{\boldsymbol{\Gamma}}_{+,p,s}^{-1\prime}},
\end{align*}
and where
\small
\begin{align*}
\mathbf{V}_{-,p,s} &= f(c) \int\limits_{-\infty}^0 K^2(u)\begin{bmatrix}
\mathbf{r}_p\left(u\right)\mathbf{r}_p\left(u\right)^\prime \operatorname{\mathbb{E}}[\sigma^2_{Y-}(X_i,\mathbf{W}_i)\mid X_i=c] & \operatorname{\mathbb{E}}[\mathbf{W}_i^\prime\sigma^2_{Y-}(X_i,\mathbf{W}_i)\mid X_i=c] \otimes \mathbf{r}_p\left(u\right)\mathbf{r}_s\left(u\right)^\prime \\
\operatorname{\mathbb{E}}[\mathbf{W}_i \sigma^2_{Y-}(X_i,\mathbf{W}_i)\mid X_i=c] \otimes \mathbf{r}_s\left(u\right)\mathbf{r}_p\left(u\right)^\prime & \operatorname{\mathbb{E}}[\mathbf{W}_i\mathbf{W}_i^\prime\sigma^2_{Y-}(X_i,\mathbf{W}_i)\mid X_i=c] \otimes \mathbf{r}_s\left(u\right)\mathbf{r}_s\left(u\right)^\prime
\end{bmatrix} \,\mathrm{d} u,\\
\mathbf{V}_{+,p,s} &= f(c) \int\limits_0^\infty K^2(u)\begin{bmatrix}
\mathbf{r}_p\left(u\right)\mathbf{r}_p\left(u\right)^\prime \operatorname{\mathbb{E}}[\sigma^2_{Y+}(X_i,\mathbf{W}_i)\mid X_i=c] & \operatorname{\mathbb{E}}[\mathbf{W}_i^\prime\sigma^2_{Y+}(X_i,\mathbf{W}_i)\mid X_i=c] \otimes \mathbf{r}_p\left(u\right)\mathbf{r}_s\left(u\right)^\prime \\
\operatorname{\mathbb{E}}[\mathbf{W}_i \sigma^2_{Y+}(X_i,\mathbf{W}_i)\mid X_i=c] \otimes \mathbf{r}_s\left(u\right)\mathbf{r}_p\left(u\right)^\prime & \operatorname{\mathbb{E}}[\mathbf{W}_i\mathbf{W}_i^\prime\sigma^2_{Y+}(X_i,\mathbf{W}_i)\mid X_i=c] \otimes \mathbf{r}_s\left(u\right)\mathbf{r}_s\left(u\right)^\prime
\end{bmatrix} \,\mathrm{d} u
\end{align*}
are positive-definite matrices.
\end{corollary}
\subsection{Consistency} \label{subsec: consistency}
Using Lemma \ref{lemma: first order variance} and Lemma \ref{lemma: smoothing bias}, we have
\begin{align*}
\mathbf{H}_{p,s}(h)\left(\widehat{\boldsymbol{\vartheta}}_{-,p,s}(h) - \boldsymbol{\vartheta}_{-,p,s}\right) &= O_{\mathbb{P}}\left(\frac{1}{\sqrt{nh}} + h^{1+p\wedge s}\right),\\
\mathbf{H}_{p,s}(h)\left(\widehat{\boldsymbol{\vartheta}}_{+,p,s}(h) - \boldsymbol{\vartheta}_{+,p,s}\right)
&= O_{\mathbb{P}}\left(\frac{1}{\sqrt{nh}} + h^{1+p\wedge s}\right),
\end{align*}
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
\begin{align*}
\breve{\mathbf{e}}_{\nu}^\prime\left(\widehat{\boldsymbol{\vartheta}}_{-,p,s}(h) - \boldsymbol{\vartheta}_{-,p,s}\right) = O_{\mathbb{P}}\left(\frac{1}{\sqrt{nh^{2\nu+1}}} + h^{1+p\wedge s-\nu}\right),\\
\breve{\mathbf{e}}_{\nu}^\prime\left(\widehat{\boldsymbol{\vartheta}}_{+,p,s}(h) - \boldsymbol{\vartheta}_{+,p,s}\right) = O_{\mathbb{P}}\left(\frac{1}{\sqrt{nh^{2\nu+1}}} + h^{1+p\wedge s-\nu}\right).
\end{align*}
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}.$
\subsection{Causal Interpretation of the Probability Limit}
Before delving into this section, we redefine the ``extractor" vector as
\[\breve{\mathbf{e}}_{\nu_x,\nu_w} = \begin{bmatrix}
\nu_x\mathbf{I}_{p+1}\mathbf{e}_{\nu_x} \\ \nu_w!\boldsymbol{\iota}_{d} \otimes \mathbf{I}_{s+1}\mathbf{e}_{{\nu_w}}
\end{bmatrix}\]
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:
\begin{enumerate}
\item When $d=1$ and $W_i$ is binary, we have that
\[\chi_{\nu_x,\cdot,p,s} = \frac{\partial^{\nu_x} }{\partial x^{\nu_x}}\operatorname{\mathbb{E}}[Y_i(1)-Y_i(0)\mid X_i=x, W_i=0] \Big|_{x=c},\]
and
\[\chi_{\nu,\nu,p,s} = \frac{\partial^{\nu} }{\partial x^{\nu}}\left(\operatorname{\mathbb{E}}[Y_i(1)-Y_i(0)\mid X_i=x, W_i=1]-\operatorname{\mathbb{E}}[Y_i(1)-Y_i(0)\mid X_i=x, W_i=0]\right) \Big|_{x=c}.\]
In other words, if a researcher is interested in the conditional average treatment effect (CATE) for the two sub-populations defined by $W_i\in\{0,1\},$ then the first element of $\widehat{\boldsymbol{\vartheta}}_{+,p,s}(h_+)-\widehat{\boldsymbol{\vartheta}}_{-,p,s}(h_-)$ identifies the CATE for the baseline ($W_i=0$) group, whereas the $(2+p)$-th element identifies the difference in CATEs between in the sub-population with $W_i=1$ and the baseline one.
\item When there is more than one binary covariate, say $U_i$ and $V_i$, then $\mathbf{W}_i=(U_i,V_i,U_i\cdot V_i)^\prime$ to correctly identify the CATEs in all the sub-populations.
\item If $d=1$ and the covariate of interest $W_i$ is categorical (e.g., race) or multi-valued discrete (e.g., age) --i.e., takes on $J$ distinct values $(1,\ldots,J)$-- then the same interpretation offered above holds as long as the covariate is dummied out. In other words, it means that
$$\mathbf{W}_i = (I_2,I_3,\ldots,I_J)^\prime, \qquad I_j:=\mathbbm{1}(W_i=j), j=2,\ldots, J .$$
\item If $d=1$ and $W_i$ is a continuous covariate (e.g., parental income), then we still have
\[\chi_{\nu_x,\cdot,p,s} = \frac{\partial^{\nu_x} }{\partial x^{\nu_x}}\operatorname{\mathbb{E}}[Y_i(1)-Y_i(0)\mid X_i=x, W_i=0] \Big|_{x=c},\]
but the interpretation of the coefficients on the interaction terms between $(X_i-c)$ and $W_i$ changes as follows
\[\chi_{\nu,p,s} = \frac{\partial^{\nu_x} }{\partial x^{\nu_x}}\frac{\partial^{\nu_w} }{\partial w^{\nu_w}}\operatorname{\mathbb{E}}[Y_i(1)-Y_i(0)\mid X_i=x, W_i=w] \Big|_{x=c},\]
where the evaluation point of the derivative in the $W$-dimension needs not to be specified as this derivative is constant in the $W$-dimension by Assumption \ref{ass: conditional expectation}. The coefficient on the interaction terms simply captures the change in the CATEs due to a marginal change in $W_i$. This change is assumed to be linear in virtue of Assumption \ref{ass: conditional expectation}.
\end{enumerate}
\subsection{MSE Expansions} \label{subsec: 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$
\begin{align*}
\breve{\mathbf{e}}_\nu'&\widehat{\boldsymbol{\vartheta}}_{+,p,s}(h) - \mu_+^{(\nu)}(c,\mathbf{w})
\\
&=h^{-\nu}\breve{\mathbf{e}}_\nu'\widehat{\boldsymbol{\Gamma}}^{-1}_{+,p,s}(h)\left(\widehat{\boldsymbol{\Upsilon}}_{+,p,s}(h) - \widehat{\boldsymbol{\Gamma}}_{+,p,s}(h)\mathbf{H}_{p,s}(h)\boldsymbol{\vartheta}_{+,p,s}^\star(h)\right) + \breve{\mathbf{e}}_\nu'\boldsymbol{\vartheta}_{+,p,s}^\star(h) - \mu_+^{(\nu)}(c,\mathbf{w}) \\
&= \underbrace{h^{-\nu}\breve{\mathbf{e}}_\nu'\widetilde{\boldsymbol{\Gamma}}^{-1}_{+,p,s}(h)\mathbf{L}_{+,p,s}(h)}_{\mathsf{L}_{+,\nu,p,s}(h)} + \underbrace{h^{-\nu}\breve{\mathbf{e}}_\nu'\left(\widehat{\boldsymbol{\Gamma}}^{-1}_{+,p,s}(h)-\widetilde{\boldsymbol{\Gamma}}^{-1}_{+,p,s}(h)\right)\mathbf{L}_{+,p,s}(h)}_{\mathsf{Q}_{+,\nu,p,s}(h)} + \underbrace{\breve{\mathbf{e}}_\nu'\boldsymbol{\vartheta}_{+,p,s}^\star(h) - \mu_+^{(\nu)}(c,\mathbf{w})}_{\mathsf{B}_{+,\nu,p,s}(h)},
\end{align*}
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{thm: asymptotic distribution} and Lemma \ref{lemma: grand hessian}, 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{ass: conditional expectation}, we have that $\mu_+^{(\nu)}(c,\mathbf{W}_i)=\breve{\mathbf{e}}'\boldsymbol{\vartheta}_{+,p,s}(h)$, thus by Lemma \ref{lemma: smoothing bias} 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{ass: sampling}-\ref{ass: conditional expectation}, 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
\begin{align*}
\mathsf{MSE}[\breve{\mathbf{e}}_{\nu}^\prime (\widehat{\boldsymbol{\vartheta}}_{-,p,s}(h) -\boldsymbol{\vartheta}_{-,p,s})]&= \operatorname{\mathbb{V}}[\mathsf{L}_{-,\nu,p,s}(h)] + \mathsf{B}_{-,\nu,p,s}(h)^2
\end{align*}
for the left side of the cutoff, whereas for the difference, we define
\begin{align*}
\mathsf{MSE}[\breve{\mathbf{e}}_{\nu}^\prime (\widehat{\boldsymbol{\vartheta}}_{+,p,s}(h)-\widehat{\boldsymbol{\vartheta}}_{-,p,s}(h) &-({\boldsymbol{\vartheta}}_{+,p,s}(h) -\boldsymbol{\vartheta}_{-,p,s}))] \\
&= \operatorname{\mathbb{V}}[\mathsf{L}_{+,\nu,p,s}(h)] +\operatorname{\mathbb{V}}[\mathsf{L}_{-,\nu,p,s}(h)] + (\mathsf{B}_{+,\nu,p,s}(h)-\mathsf{B}_{-,\nu,p,s}(h))^2.
\end{align*}
Under Assumptions \ref{ass: sampling}-\ref{ass: kernel} and relying on Lemma \ref{lemma: smoothing bias}, the bias approximations for the RD estimators are given by
\begin{align*}
\mathsf{B}_{-,\nu,p,s}(h) = h^{1+p-\nu} \mathcal{B}^{[0]}_{-,\nu,p,s}(h) + h^{1+s-\nu} \mathcal{B}^{[1]}_{-,\nu,p,s}(h) + o(h^{1+(p\wedge s)-\nu}), \\
\mathsf{B}_{+,\nu,p,s}(h) = h^{1+p-\nu} \mathcal{B}^{[0]}_{+,\nu,p,s}(h) + h^{1+s-\nu} \mathcal{B}^{[1]}_{+,\nu,p,s}(h) + o(h^{1+(p\wedge s)-\nu}),
\end{align*}
where
\begin{alignat*}{3}
\mathcal{B}^{[0]}_{-,\nu,p,s}(h) =\breve{\mathbf{e}}_{\nu}^\prime\mathbf{B}_{-,p,s}^{[0]}(h)&\to \breve{\mathbf{e}}_{\nu}^\prime\mathbf{B}_{-,p,s}^{[0]}=:\mathcal{B}^{[0]}_{-,\nu,p,s}, \quad \mathcal{B}^{[1]}_{-,\nu,p,s}(h) = \breve{\mathbf{e}}_{\nu}^\prime\mathbf{B}_{-,p,s}^{[1]}(h)&\to \breve{\mathbf{e}}_{\nu}^\prime\mathbf{B}_{-,p,s}^{[1]}(h)=:\mathcal{B}^{[1]}_{-,\nu,p,s}, \\
\mathcal{B}^{[0]}_{+,\nu,p,s}(h) =\breve{\mathbf{e}}_{\nu}^\prime\mathbf{B}_{+,p,s}^{[0]}(h) &\to \breve{\mathbf{e}}_{\nu}^\prime\mathbf{B}_{-,p,s}^{[0]} =:\mathcal{B}^{[0]}_{+,\nu,p,s}, \quad
\mathcal{B}^{[1]}_{+,\nu,p,s}(h) = \breve{\mathbf{e}}_{\nu}^\prime\mathbf{B}_{+,p,s}^{[1]}(h)&\to \breve{\mathbf{e}}_{\nu}^\prime\mathbf{B}_{+,p,s}^{[1]}(h)=:\mathcal{B}^{[1]}_{+,\nu,p,s}.
\end{alignat*}
Furthermore, we define
\begin{alignat*}{3}
\mathcal{B}_{-,\nu,p,s}(h) &:=\mathbbm{1}(p\leq s)\mathcal{B}^{[0]}_{-,\nu,p,s}(h)+\mathbbm{1}(p\geq s)\mathcal{B}^{[1]}_{-,\nu,p,s}(h),\qquad
&\mathcal{B}_{-,\nu,p,s} &:=\mathbbm{1}(p\leq s)\mathcal{B}^{[0]}_{-,\nu,p,s}+\mathbbm{1}(p\geq s)\mathcal{B}^{[1]}_{-,\nu,p,s},\\
\mathcal{B}_{+,\nu,p,s}(h) &:=\mathbbm{1}(p\leq s)\mathcal{B}^{[0]}_{+,\nu,p,s}(h)+\mathbbm{1}(p\geq s)\mathcal{B}^{[1]}_{+,\nu,p,s}(h), \qquad
&\mathcal{B}_{+,\nu,p,s} &:=\mathbbm{1}(p\leq s)\mathcal{B}^{[0]}_{+,\nu,p,s}+\mathbbm{1}(p\geq s)\mathcal{B}^{[1]}_{+,\nu,p,s}.
\end{alignat*}
Under Assumptions \ref{ass: sampling}-\ref{ass: kernel} and relying on Corollary \ref{cor: asympt distr under sa5}
\begin{align*}
\operatorname{\mathbb{V}}[\mathsf{L}_{-,\nu,p,s}(h)] &= \frac{1}{nh^{2\nu+1}}\mathcal{V}_{-,\nu,p,s}(h), \qquad \mathcal{V}_{-,\nu,p,s}(h):=\breve{\mathbf{e}}_{\nu}^\prime\widetilde{\boldsymbol{\Gamma}}_{-,p,s}^{-1}(h)\mathbf{V}_{-,p,s}(h)
{\widetilde{\boldsymbol{\Gamma}}_{-,p,s}^{-1\prime}(h)}\breve{\mathbf{e}}_{\nu},\\
\operatorname{\mathbb{V}}[\mathsf{L}_{+,\nu,p,s}(h)] &= \frac{1}{nh^{2\nu+1}}\mathcal{V}_{+,\nu,p,s}(h),\qquad \mathcal{V}_{+,\nu,p,s}(h):=\breve{\mathbf{e}}_{\nu}^\prime\widetilde{\boldsymbol{\Gamma}}_{+,p,s}^{-1}(h)\mathbf{V}_{+,p,s}(h)
{\widetilde{\boldsymbol{\Gamma}}_{+,p,s}^{-1\prime}(h)}\breve{\mathbf{e}}_{\nu}.
\end{align*}
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
\begin{align*}
\mathsf{MSE}[\breve{\mathbf{e}}_{\nu}^\prime (\widehat{\boldsymbol{\vartheta}}_{-,p,s}(h) -\boldsymbol{\vartheta}_{-,p,s})] &= \frac{1}{nh^{2\nu+1}}\mathcal{V}_{-,\nu,p,s} + h^{2(1+p\wedge s-\nu)}\mathcal{B}^2_{-,\nu,p,s}, \\
\mathsf{MSE}[\breve{\mathbf{e}}_{\nu}^\prime (\widehat{\boldsymbol{\vartheta}}_{+,p,s}(h) -\boldsymbol{\vartheta}_{+,p,s})] &= \frac{1}{nh^{2\nu+1}}\mathcal{V}_{+,\nu,p,s} + h^{2(1+p\wedge s-\nu)}\mathcal{B}^2_{+,\nu,p,s}.
\end{align*}
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
\begin{align*}
h^\star_{-,\nu,p,s} &:= \operatorname*{arg\,min}_{h>0} \left[\frac{1}{nh^{2\nu+1}}\mathcal{V}_{-,\nu,p,s} + h^{2((p\wedge s)+1-\nu)}\mathcal{B}_{-,\nu,p,s}^2\right], \\
h^\star_{+,\nu,p,s} &:= \operatorname*{arg\,min}_{h>0} \left[\frac{1}{nh^{2\nu+1}}\mathcal{V}_{+,\nu,p,s} + h^{2((p\wedge s)+1-\nu)}\mathcal{B}_{+,\nu,p,s}^2\right],
\end{align*}
so, under the additional assumption that $\mathcal{B}_{-,\nu,p,s}\neq0\neq \mathcal{B}_{+,\nu,p,s}$, we get
\begin{align*}
h^\star_{-,\nu,p,s} =\left[\frac{1+2\nu}{2(1+(p\wedge s)-\nu)n}\frac{\mathcal{V}_{-,\nu,p,s}}{\mathcal{B}_{-,\nu,p,s}^2}\right]^{\frac{1}{3+2(p\wedge s)}}, \qquad h^\star_{+,\nu,p,s} =\left[\frac{1+2\nu}{2(1+(p\wedge s)-\nu)n}\frac{\mathcal{V}_{+,\nu,p,s}}{\mathcal{B}_{+,\nu,p,s}^2}\right]^{\frac{1}{3+2(p\wedge s)}}.
\end{align*}
\paragraph{Two-sided Optimal Bandwidth.} In this case a single bandwidth is chosen, so $h=h_+=h_-$ and the optimal bandwidth is defined as
\begin{align*}
h^\star_{\nu,p,s} := \operatorname*{arg\,min}_{h>0} \left[\frac{1}{nh^{2\nu+1}}(\mathcal{V}_{+,\nu,p,s}+\mathcal{V}_{-,\nu,p,s}) + h^{2((p\wedge s)+1-\nu)}(\mathcal{B}_{+,\nu,p,s}-\mathcal{B}_{-,\nu,p,s})^2\right],
\end{align*}
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)}}.\]
\section{Empirical Practice Investigating Covariate-Heterogeneity}
\label{sec: empirical app}
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{table:papers} 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.
\begin{sidewaystable}[!ht]
\centering
\caption{Empirical Practice for RD-HTE}
\resizebox{\textwidth}{!}{\small
\begin{tabular}{p{5cm} p{4.5cm} p{3.5cm} p{2cm} | p{0.6cm} p{0.6cm} p{0.6cm} | m{1.6cm} m{0.6cm} m{0.6cm} m{1.2cm} | m{0.7cm} m{0.6cm} }
\toprule\toprule
& & & & \multicolumn{3}{c}{\textit{Heterogeneity}} & \multicolumn{4}{c}{\textit{Estimation}} & \multicolumn{2}{c}{\textit{Std. Errs.}} \\
\textit{Paper} & \textit{Outcome} & \textit{Treatment} & \textit{Running Var.} & \textit{Disc.} & \textit{Time} & \textit{Cont.} & \textit{Bandwidth} & \textit{Local} & \textit{Joint} & \textit{Covs/FE} & \textit{Robust} & \textit{Cluster} \\ \midrule
\cite{Adams} & Healthcare utilization & Financial Assistance & Poverty level & & \checkmark & & None & & & & & \\
\cite{Akhtari-Moreira-Trucco2022_AER} & Municipal Bureaucracy & Political turnover & Vote Margin & \checkmark & & & MSE-HTE & \checkmark & & \checkmark & & \checkmark \\
\cite{Asher} & Economic development & Road construction & Population & \checkmark & & & MSE & \checkmark & & \checkmark & & \checkmark \\
\cite{Brollo} & Political corruption, quality & Government revenues & Population & \checkmark & & \checkmark & None & \checkmark & \checkmark & & \checkmark \\
\cite{Dell} & Drug-related violence & Drug enforcement & Vote Margin & \checkmark & \checkmark & & Manual & \checkmark & \checkmark & \checkmark & & \checkmark \\
\cite{Miralles} & Savings, Retirement behavior & Public pension benefits & Age & \checkmark & \checkmark & & MSE & \checkmark & & \checkmark & \checkmark & \\
\cite{Han} & Healthcare Utilization & Patient Cost-Sharing & Age & \checkmark & & & Manual & \checkmark & & & & \checkmark \\
\cite{Huh} & Mortality, risky behaviors & Teenage Driving & Age & \checkmark & & & MSE & \checkmark & & & \checkmark & \\
\cite{Jones} & Irrigation Adoption, Profits & Access to water & Spatial & \checkmark & & & Manual & \checkmark & & \checkmark & & \checkmark \\
\cite{Lindo} & Academic Performance & Academic probation & Test Scores & \checkmark & & & Manual & \checkmark & & \checkmark & & \checkmark \\
\cite{McEwan} & Economics Major Choice & Higher letter grade & Test Scores & \checkmark & & & MSE & \checkmark & \checkmark & \checkmark & \checkmark & \\
\cite{Miglino} & Health outcomes & Financial Assistance & Age & \checkmark & \checkmark & & Manual & \checkmark & & \checkmark & & \checkmark \\
\cite{Pop-Eleches-Urquiola2013_AER} & Academic Performance & Access to better schools & Test Scores & \checkmark & \checkmark & & Manual & \checkmark & & \checkmark & & \checkmark \\
\cite{Shigeoka} & Utilization, health & Patient Cost Sharing & Age & \checkmark & \checkmark & & Manual & \checkmark & & \checkmark & & \checkmark \\
\cite{Silliman} & Labor market returns & Vocational Education & Test Scores & \checkmark & \checkmark & & Manual & \checkmark & & \checkmark & & \checkmark \\
\cite{Zimmerman} & Better jobs, Income & Access to better schools & Test Scores & \checkmark & & & Manual & \checkmark & \checkmark & \checkmark & & \checkmark \\ \bottomrule\bottomrule
\end{tabular}
}
\footnotesize{
\begin{flushleft}Notes:\end{flushleft}\vspace{-0.2in}
\begin{enumerate}[label=\normalfont(\roman*),noitemsep,leftmargin=*]
\item The left four columns give the citation and the empirical setting in terms of the outcome variable, treatment studied, and running variable.
\item The type of heterogeneity is studied is listed in columns 5--7, showing, respectively, if the covariates used for heterogeneity are discrete (``Disc.''), time-related (``Time''), or continuous (``Cont.'').
\item For estimation, the ``Bandwidth'' column reports the method of bandwidth selection used, where ``Manual'' indicates a hand-selected bandwidth, ``MSE'' for a mean squared error optimal bandwidth, and ``MSE-HTE'' if the bandwidth was chosen particularly for the heterogeneity analysis. The ``Local'' indicates that the analysis is conducted local to the cutoff; ``Joint'' if the heterogeneity was conducted jointly across levels of the variable, and ``Covs/FE'' if additional covariates or fixed effects are included for efficiency.
\item Standard errors are described in the final two columns, indicating if robust standard errors are used or if the standard errors are clustered.
\end{enumerate}
}
\label{table:papers}
\end{sidewaystable}
\cleardoublepage
\section{Proofs}
\label{sec: proofs}
\subsection{Proof of Lemma \ref{lemma: markov}}
\begin{proof}\label{proof: markov}
Here, for simplicity, we just prove the lemma for the scalar case. By Markov inequality, we have
\[\forall\, M>0,\quad \mathbb{P}(|A_n - \operatorname{\mathbb{E}}[A_n]| \geq M) \leq \frac{\operatorname{\mathbb{E}}[|A_n-\operatorname{\mathbb{E}}[A_n]|]}{M} \leq \frac{\operatorname{\mathbb{V}}[A_n]^{1/2}}{M},\]
which implies that $\operatorname{\mathbb{V}}[A_n]^{-1/2}(A_n - \operatorname{\mathbb{E}}[A_n]) = O_{\mathbb{P}}(1)$ which was to be shown. The matrix case follows using $\mathbb{P}(\mathbf{A}_n\preceq \mathbf{I} ) \leq \mathrm{tr}(\operatorname{\mathbb{E}}[\mathbf{A}_n])$, where $\mathbf{A}_n\preceq \mathbf{I}$ means that $\mathbf{A}_n- \mathbf{I}$ is negative semi-definite.
\end{proof}
\subsection{Proof of Lemma \ref{lemma: hessian}}
\begin{proof}\label{proof: hessian}
The proof covers the case to the right of the cutoff. Everything follows symmetrically for the other case. First, for any $p,s\in\operatorname{\mathbb{N}}$
\begin{align*}
\operatorname{\mathbb{E}}[\widehat{\boldsymbol{\Lambda}}_{+,p,s}(h)] &= \operatorname{\mathbb{E}}\left[\frac{1}{n}\sum_{i=1}^n\mathbbm{1} ( X_i\geq c)\mathbf{r}_p\left(\frac{X_i-c}{h}\right)\mathbf{r}_s\left(\frac{X_i-c}{h}\right)^\prime K\left(\frac{X_i-c}{h}\right)\frac{1}{h}\right] \tag{definition}\\
&=\frac{1}{h}\operatorname{\mathbb{E}}\left[\mathbbm{1} ( X_i\geq c)\mathbf{r}_p\left(\frac{X_i-c}{h}\right)\mathbf{r}_s\left(\frac{X_i-c}{h}\right)^\prime K\left(\frac{X_i-c}{h}\right)\right] \tag{Assumption \ref{ass: sampling}} \\
&= \frac{1}{h}\int\limits_{c}^\infty \mathbf{r}_p\left(\frac{x-c}{h}\right)\mathbf{r}_s\left(\frac{x-c}{h}\right)^\prime K\left(\frac{x-c}{h}\right)f(x)\,\mathrm{d} x \\
&= \int\limits_0^\infty \mathbf{r}_p(u)\mathbf{r}_s(u)^\prime K(u)f(hu+c) \,\mathrm{d} u \tag{$u = (x-c)/h$}
=: \widetilde{\boldsymbol{\Gamma}}_{+,p,s}(h).
\end{align*}
By Lemma \ref{lemma: markov}
\[\widehat{\boldsymbol{\Lambda}}_{+,p,s}(h) = \operatorname{\mathbb{E}}[\widehat{\boldsymbol{\Lambda}}_{+,p,s}(h)] + O_{\mathbb{P}}(|\operatorname{\mathbb{V}}[\widehat{\boldsymbol{\Lambda}}_{+,p,s}(h)]|).\]
The generic element of $\widehat{\boldsymbol{\Lambda}}_{+,p,s}(h)$ is
\[\frac{1}{nh}\sum_{i=1}^n\mathbbm{1}(X_i\geq c)\left(\frac{X_i-c}{h}\right)^{j} K\left(\frac{X_i-c}{h}\right), \quad j\in\{0,\cdots,p+s\}.\]
Fix $j\in\{0,\cdots,p+s\}$, then
\begin{align*}
&\operatorname{\mathbb{V}}\left[\frac{1}{nh}\sum_{i=1}^n\mathbbm{1} ( X_i\geq c)\left(\frac{X_i-c}{h}\right)^{j} K\left(\frac{X_i-c}{h}\right)\right] \\
&\qquad\qquad= \frac{1}{nh^2}\operatorname{\mathbb{V}}\left[\mathbbm{1} ( X_i\geq c)\left(\frac{X_i-c}{h}\right)^{j} K\left(\frac{X_i-c}{h}\right)\right] \tag{Assumption \ref{ass: sampling}} \\
&\qquad \qquad \leq \frac{1}{nh^2}\operatorname{\mathbb{E}}\left[\mathbbm{1} ( X_i\geq c)\left(\frac{X_i-c}{h}\right)^{2j} K\left(\frac{X_i-c}{h}\right)^2\right]\\
&\qquad \qquad = \frac{1}{nh}\int\limits_0^\infty u^{2j} K(u)^2 f(uh+c)\,\mathrm{d} u = O(n^{-1}h^{-1}), \tag{$u = (x-c)/h$}
\end{align*}
where the last equality follows because $u^{2j} K(u)^2 f(uh+c)$ is integrable due to Assumptions \ref{ass: density running variable}, \ref{ass: kernel}, and $\mathsf{k}\geq h$. Hence $|\operatorname{\mathbb{V}}[\widehat{\boldsymbol{\Lambda}}_{+,p,s}(h)]| = O(n^{-1}h^{-1})$. Finally, we have $\widehat{\boldsymbol{\Lambda}}_{+,p,s}(h) = \operatorname{\mathbb{E}}[\widehat{\boldsymbol{\Lambda}}_{+,p,s}(h)] + O_{\mathbb{P}}(1/\sqrt{nh})$, which was to be shown.
\end{proof}
\subsection{Proof of Lemma \ref{lemma: grand hessian}}
\begin{proof}\label{proof: grand hessian}
The proof covers the case to the right of the cutoff. Throughout the proof, we maintain that $n$ is large enough so that $\mathsf{k}\leq c+h$ in Assumption \ref{ass: DGP}. Everything follows symmetrically for the other case. First, note that
\[ \widehat{\boldsymbol{\Gamma}}_{+,p,s}(h) = \begin{bmatrix}
\widehat{\boldsymbol{\Lambda}}_{+,p,p}(h) & \widehat{\mathbf{G}}_{+,2}(h) \\ \widehat{\mathbf{G}}_{+,2}(h)^\prime & \widehat{\mathbf{G}}_{+,3}(h)
\end{bmatrix} ,\]
where
\begin{align*}
\widehat{\boldsymbol{\Lambda}}_{+,p,p}(h) &= \frac{1}{nh}\sum_{i=1}^n\mathbbm{1} ( X_i\geq c) \mathbf{r}_p\left(\frac{X_i-c}{h}\right)\mathbf{r}_p\left(\frac{X_i-c}{h}\right)^\prime K\left(\frac{X_i-c}{h}\right), \\
\widehat{\mathbf{G}}_{+,2}(h)&=\frac{1}{nh}\sum_{i=1}^n\mathbbm{1} ( X_i\geq c) \left[\mathbf{W}_i^\prime\otimes K\left(\frac{X_i-c}{h}\right)\mathbf{r}_p\left(\frac{X_i-c}{h}\right)\mathbf{r}_s\left(\frac{X_i-c}{h}\right)^\prime\right], \\
\widehat{\mathbf{G}}_{+,3}(h)&=\frac{1}{nh}\sum_{i=1}^n\mathbbm{1} ( X_i\geq c) \left[\mathbf{W}_i\mathbf{W}_i^\prime\otimes K\left(\frac{X_i-c}{h}\right)\mathbf{r}_s\left(\frac{X_i-c}{h}\right)\mathbf{r}_s\left(\frac{X_i-c}{h}\right)^\prime\right].
\end{align*}
Then, from Lemma \ref{lemma: hessian} we immediately have
\[\widehat{\boldsymbol{\Lambda}}_{+,p,p}(h) = \widetilde{\boldsymbol{\Lambda}}_{+,p,p}(h) + O_{\mathbb{P}}(1/\sqrt{nh}).\]
Consider the expectation of $\widehat{\mathbf{G}}_{+,2}(h)$. Using Assumption \ref{ass: sampling} and the change of variable $u = (x-c)/h$:
\begin{align*}
\operatorname{\mathbb{E}}[\widehat{\mathbf{G}}_{+,2}(h)] &= \operatorname{\mathbb{E}}\left[\frac{1}{h}\mathbbm{1} ( X_i\geq c)\left\{\boldsymbol{\mu}_{W}(X_i)^\prime\otimes K\left(\frac{X_i-c}{h}\right)\mathbf{r}_p\left(\frac{X_i-c}{h}\right)\mathbf{r}_s\left(\frac{X_i-c}{h}\right)^\prime\right\}\right] \\
&=\int\limits_{0}^\infty \boldsymbol{\mu}_{W}(uh+c)^\prime\otimes \mathbf{r}_p\left(u\right)\mathbf{r}_s\left(u\right)^\prime K(u)f(uh+c)\,\mathrm{d} u=:\widetilde{\mathbf{G}}_{+,2}(h),
\end{align*}
The generic element of $\widehat{\mathbf{G}}_{+,2}(h)$ is
\[\frac{1}{nh}\sum_{i=1}^n\mathbbm{1}(X_i\geq c) W_{i\ell}\cdot\left(\frac{X_i-c}{h}\right)^j K\left(\frac{X_i-c}{h}\right), \qquad \ell\in\{1,\cdots,d\},\:\:j\in\{0,\cdots,p+s\}.\]
Fix $\ell\in\{1,\cdots,d\} $ and $j\in\{0,\cdots,p+s\}$. Then
\begin{align*}
&\operatorname{\mathbb{V}}\left[\frac{1}{nh}\sum_{i=1}^n\mathbbm{1}(X_i\geq c) W_{i\ell}\cdot\left(\frac{X_i-c}{h}\right)^j K\left(\frac{X_i-c}{h}\right)\right] \\
&\qquad\qquad =\frac{1}{nh^2}\operatorname{\mathbb{V}}\left[\mathbbm{1}(X_i\geq c) W_{i\ell}\cdot\left(\frac{X_i-c}{h}\right)^j K\left(\frac{X_i-c}{h}\right)\right]\tag{Assumption \ref{ass: sampling}}\\
&\qquad\qquad \leq \frac{1}{nh^2}\operatorname{\mathbb{E}}\left[\mathbbm{1}(X_i\geq c) W_{i\ell}^2\cdot\left(\frac{X_i-c}{h}\right)^{2j} K\left(\frac{X_i-c}{h}\right)^2\right]\\
&\qquad\qquad =\frac{1}{nh}\int\limits_0^\infty \mu_{W_\ell^2}(uh+c) u^{2j} K(u)^2 f(uh+c)\,\mathrm{d} u \tag{$u=(x-c)/h$} = O(n^{-1}h^{-1}),
\end{align*}
where the last equality follows because of Assumptions \ref{ass: density running variable}, \ref{ass: kernel}, and the fact that $\mu_{W_\ell^2}(\cdot)$ is bounded in a neighborhood of $c$, which is granted by Assumption \ref{ass: SRD continuous W} and $\mathsf{k}\geq h$. Hence, by Lemma \ref{lemma: markov}
\[\widehat{\mathbf{G}}_{+,2}(h) = \widetilde{\mathbf{G}}_{+,2}(h) + O_{\mathbb{P}}(1/\sqrt{nh}).\]
Consider $\widehat{\mathbf{G}}_{+,3}(h)$. Using a similar logic as above
\begin{align*}
\operatorname{\mathbb{E}}[\widehat{\mathbf{G}}_{+,3}(h)]&=\frac{1}{h}\operatorname{\mathbb{E}}\left[ \boldsymbol{\mu}_{WW}(X_i)\otimes\mathbbm{1} ( X_i\geq c)K\left(\frac{X_i-c}{h}\right)\mathbf{r}_s\left(\frac{X_i-c}{h}\right)\mathbf{r}_s\left(\frac{X_i-c}{h}\right)^\prime\right]\\
&=\int\limits_{0}^\infty \boldsymbol{\mu}_{WW}(uh+c)\otimes K\left(u\right)\mathbf{r}_s\left(u\right)\mathbf{r}_s\left(u\right)^\prime f(uh+c)\,\mathrm{d} u =: \widetilde{\mathbf{G}}_{+,3}(h)
\end{align*}
and for $\ell, \ell' \in\{1,\cdots,d\} $ and $j\in\{0,\cdots,2s\}$. Then
\begin{align*}
&\operatorname{\mathbb{V}}\left[\frac{1}{nh}\sum_{i=1}^n\mathbbm{1}(X_i\geq c) W_{i\ell}W_{i\ell'}\cdot\left(\frac{X_i-c}{h}\right)^j K\left(\frac{X_i-c}{h}\right)\right] \\
&\qquad\qquad =\frac{1}{nh^2}\operatorname{\mathbb{V}}\left[\mathbbm{1}(X_i\geq c) W_{i\ell}W_{i\ell'}\cdot\left(\frac{X_i-c}{h}\right)^j K\left(\frac{X_i-c}{h}\right)\right]\tag{Assumption \ref{ass: sampling}}\\
&\qquad\qquad \leq \frac{1}{nh^2}\operatorname{\mathbb{E}}\left[\mathbbm{1}(X_i\geq c) W_{i\ell}^2W_{i\ell'}^2\cdot\left(\frac{X_i-c}{h}\right)^{2j} K\left(\frac{X_i-c}{h}\right)^2\right]\\
&\qquad\qquad =\frac{1}{nh}\int\limits_0^\infty \operatorname{\mathbb{E}}[W_{i\ell}^2W_{i\ell'}^2\mid X_i=uh+c] u^{2j} K(u)^2 f(uh+c)\,\mathrm{d} u \tag{$u=(x-c)/h$} = O(n^{-1}h^{-1}),
\end{align*}
where the last equality follows because of Assumptions \ref{ass: density running variable}, \ref{ass: kernel}, and the fact that $\operatorname{\mathbb{E}}[W_{i\ell}^2W_{i\ell'}^2\mid X_i=uh+c]$ is bounded in $[0,c+h]$, which is granted by Assumption \ref{ass: SRD continuous fourth moment of W} and $\mathsf{k}\geq h$.
Last, we prove that $\widehat{\boldsymbol{\Gamma}}_{+,p,s}(h)$ is asymptotically invertible. To do so, note that by taking the limit as $h\to 0$, the continuity of $f(\cdot)$ (Assumption \ref{ass: density running variable}), $\boldsymbol{\mu}_W(\cdot)$, and $\boldsymbol{\mu}_{WW}(\cdot)$ (Assumption \ref{ass: SRD continuous W}) give us that
\begin{align*}
\widetilde{\boldsymbol{\Lambda}}_{+,p,p}(h) &= f(c)\int\limits_0^\infty K(u) \mathbf{r}_p(u)\mathbf{r}_p(u)^\prime \,\mathrm{d} u + o(1) = \boldsymbol{\Lambda}_{+,p,p} + o(1), \\
\widetilde{\mathbf{G}}_{+,2}(h) &= \boldsymbol{\mu}_{W}^\prime \otimes f(c)\int\limits_0^\infty K(u) \mathbf{r}_p(u)\mathbf{r}_s(u)^\prime \,\mathrm{d} u + o(1)=\boldsymbol{\mu}_{W}^\prime \otimes \boldsymbol{\Lambda}_{+,p,s} + o(1) = \mathbf{G}_{+,3}+o(1), \\
\widetilde{\mathbf{G}}_{+,3}(h) &= \boldsymbol{\mu}_{WW} \otimes f(c)\int\limits_0^\infty K(u) \mathbf{r}_s(u)\mathbf{r}_s(u)^\prime \,\mathrm{d} u + o(1)= \boldsymbol{\mu}_{WW}\otimes \boldsymbol{\Lambda}_{+,s,s} + o(1) = \mathbf{G}_{+,3} + o(1).
\end{align*}
Therefore,
\[\widehat{\boldsymbol{\Gamma}}_{+,p,s}(h) = \boldsymbol{\Gamma}_{+,p,s} + o(1) + O_{\mathbb{P}}(1/\sqrt{nh}),\quad\text{where}\quad \boldsymbol{\Gamma}_{+,p,s} = \begin{bmatrix}
\boldsymbol{\Lambda}_{+,p,p} & \mathbf{G}_{+,2} \\
\mathbf{G}_{+,2}' & \mathbf{G}_{+,3}
\end{bmatrix},\]
which is non-singular by Assumption \ref{ass: pd gram matrix}.
\end{proof}
\subsection{Proof of Lemma \ref{lemma: vartheta hat order}}
\begin{proof}\label{proof: vartheta hat order}
First,
\begin{align*}
\operatorname{\mathbb{E}}[\widehat{\boldsymbol{\zeta}}_{+,p,s,a}(h)] &= \operatorname{\mathbb{E}}\left[\frac{1}{nh}\sum_{i=1}^n\mathbf{r}_{p,s}\left(\frac{X_i-c}{h},\mathbf{W}_i\right)K\left(\frac{X_i-c}{h}\right)\left(\frac{X_i-c}{h}\right)^{a+1}\right] \\
&= \int\limits_0^\infty K(u)\mathbf{r}_{p,s}(u,\boldsymbol{\mu}_{W}(uh+c))u^{a+1} f(c + hu)\,\mathrm{d} u \\
&= f(c)\int\limits_0^\infty K(u)\mathbf{r}_{p,s}(u,\boldsymbol{\mu}_{W})u^{a+1} \,\mathrm{d} u + O(h).
\end{align*}
Moreover, the generic element of $\widehat{\boldsymbol{\zeta}}_{+,p,s,a}(h)$ is either
\begin{align*}
\frac{1}{nh}\sum_{i=1}^n\mathbbm{1}(X_i\geq c)K\left(\frac{X_i-c}{h}\right)\left(\frac{X_i-c}{h}\right)^{j+a+1}, \qquad j\in\{0,1,\cdots, p\},
\end{align*}
or
\begin{align*}
\frac{1}{nh}\sum_{i=1}^n\mathbbm{1}(X_i\geq c)K\left(\frac{X_i-c}{h}\right)W_{i\ell}\cdot\left(\frac{X_i-c}{h}\right)^{j+a+1}, \qquad j\in\{0,1,\cdots, s\},\:\:\ell\in\{1,\cdots, d\}.
\end{align*}
The variances of these elements are of order $O(n^{-1}h^{-1})$. Indeed, the variance of the latter is
\begin{align*}
\operatorname{\mathbb{V}}&\left[\frac{1}{nh}\sum_{i=1}^n\mathbbm{1}(X_i\geq c)K\left(\frac{X_i-c}{h}\right)W_{i\ell}\cdot\left(\frac{X_i-c}{h}\right)^{j+a+1}\right]\\
&\qquad\leq \frac{1}{nh^2}\operatorname{\mathbb{E}}\left[\mathbbm{1}(X_i\geq c)K^2\left(\frac{X_i-c}{h}\right)W_{i\ell}^2\cdot\left(\frac{X_i-c}{h}\right)^{2(j+a+1)}\right] \\
&\qquad=\frac{1}{nh}\int\limits_0^\infty \mu_{W_\ell}(uh+c)u^{2(j+a+1)}K^2(u)f(uh+c) = O\left(\frac{1}{nh}\right),
\end{align*}
where for the last equality we used Assumptions \ref{ass: sampling} and \ref{ass: kernel} and $\mathsf{k}\geq h$. Then, $\widehat{\boldsymbol{\zeta}}_{+,p,s,a}(h) = \operatorname{\mathbb{E}}[\widehat{\boldsymbol{\zeta}}_{+,p,s,a}(h)] + O_{\mathbb{P}}(1/\sqrt{nh})$ by Lemma \ref{lemma: markov}. A similar argument goes through for $\widehat{\boldsymbol{\varphi}}_{+,p,s,a}(h)$.
The last part of the lemma follows by taking the limit as $h\to 0$ and using the continuity of $f(\cdot)$ (Assumption \ref{ass: density running variable}), $\boldsymbol{\mu}_W(\cdot)$, and $\boldsymbol{\mu}_{WW}(\cdot)$ (Assumption \ref{ass: SRD continuous W}).
\end{proof}
\subsection{Proof of Lemma \ref{lemma: first order variance}}
\begin{proof}\label{proof: first order variance}
Fix $h>0$ and define
\begin{align*}
\mathbf{w}_{+,i}(h):=\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),
\end{align*}
so that we can write
\[\mathbf{L}_{+,p,s}(h)= \frac{1}{nh}\sum_{i=1}^n\mathbf{w}_{+,i}(h)u_{+,i}(h).\]
By the properties of $\ell_2$-residuals, we get that $\operatorname{\mathbb{E}}[\mathbf{w}_{+,i}(h)u_{+,i}(h)]=0$ and so $\operatorname{\mathbb{E}}[\mathbf{L}_{+,p,s}(h)]=0$. The variance of $\mathbf{L}_{+,p,s}(h)$ is
\begin{align*}
\operatorname{\mathbb{V}}\left[\frac{1}{nh}\sum_{i=1}^n\mathbf{w}_{+,i}(h)u_{+,i}(h)\right] &= \frac{1}{n^2h^2}\sum_{i=1}^n\sum_{j=1}^n\operatorname{\mathbb{C}ov}(\mathbf{w}_{+,i}(h)u_{+,i}(h), \mathbf{w}_{+,j}(h)u_{+,j}(h)) \\
&= \frac{1}{nh^2}\operatorname{\mathbb{V}}[\mathbf{w}_{+,i}(h)u_{+,i}(h)] \tag{Assumption \ref{ass: sampling}} \\
&= \frac{1}{nh}\underbrace{\frac{1}{h}\operatorname{\mathbb{E}}[\mathbf{w}_{+,i}(h)\mathbf{w}_{+,i}(h)^\prime u^2_{+,i}(h)]}_{:=\mathbf{V}_{+,p,s}(h)} \tag{$\operatorname{\mathbb{E}}[\mathbf{w}_{+,i}(h)u_{+,i}(h)]=0$}.
\end{align*}
By the law of iterated expectations and changing variables, we get\footnotesize
\begin{align*}
&\mathbf{V}_{+,p,s}(h) =\\
&\quad \int\limits_0^\infty f(uh+c) K^2(u)\begin{bmatrix}
\mathbf{r}_p\left(u\right)\mathbf{r}_p\left(u\right)^\prime \operatorname{\mathbb{E}}[u_{+,i}^2(h)\mid X_i=uh+c] & \operatorname{\mathbb{E}}[\mathbf{W}_i^\prime u_{+,i}^2(h)\mid X_i=uh+c] \otimes \mathbf{r}_p\left(u\right)\mathbf{r}_s\left(u\right)^\prime \\
\operatorname{\mathbb{E}}[\mathbf{W}_i u_{+,i}^2(h)\mid X_i=uh+c] \otimes \mathbf{r}_s\left(u\right)\mathbf{r}_p\left(u\right)^\prime & \operatorname{\mathbb{E}}[\mathbf{W}_i\mathbf{W}_i^\prime u_{+,i}^2(h)\mid X_i=uh+c] \otimes \mathbf{r}_s\left(u\right)\mathbf{r}_s\left(u\right)^\prime
\end{bmatrix} \,\mathrm{d} u,
\end{align*}\normalsize
Finally, the last part of the lemma follows by taking the limit as $h\to 0$, using the properties of compact kernels (Assumption \ref{ass: kernel}), the continuity of $f(\cdot)$ (Assumption \ref{ass: density running variable}), and the bounded conditional fourth moments of $(Y_i,\mathbf{W}_i')'$ (Assumptions \ref{ass: continuous fourth moment of Y}-\ref{ass: SRD continuous fourth moment of W}), together with the Cauchy-Schwarz inequality and $\mathsf{k}\geq h$.
\end{proof}
\subsection{Proof of Theorem \ref{thm: asymptotic distribution}}
\begin{proof}\label{proof: asymptotic distribution}
Fix an arbitrary $\mathbf{c}\in\operatorname{\mathbb{R}}^{1+p+d}\setminus \{\mathbf{0}_{1+p+d}\}$. Consider
$$\phi_{+,p,s}(h) := \sqrt{nh}\mathbf{c}^\prime \mathbf{H}_{p,s}(h)(\widehat{\boldsymbol{\vartheta}}_{+,p,s}(h) - \boldsymbol{\vartheta}_{+,p,s}^\star(h))= \sqrt{nh}\mathbf{c}^\prime\widehat{\boldsymbol{\Gamma}}_{+,p,s}^{-1}(h)\mathbf{L}_{+,p,s}(h).$$
By Lemma \ref{lemma: grand hessian} and Lemma \ref{lemma: first order variance}, we get
\[\phi_{+,p,s}(h) = \dot{\phi}_{+,p,s}(h) +o_{\mathbb{P}}(1), \qquad \dot{\phi}_{+,p,s}(h)= \sqrt{nh}\mathbf{c}^\prime\widetilde{\boldsymbol{\Gamma}}_{+,p,s}^{-1}(h)\mathbf{L}_{+,p,s}(h).\]
Then, by Lemma \ref{lemma: first order variance}
\begin{align*}
\operatorname{\mathbb{V}}[\sqrt{nh}\mathbf{c}^\prime\widetilde{\boldsymbol{\Gamma}}_{+,p,s}^{-1}(h)\mathbf{L}_{+,p,s}(h)] &= \mathbf{c}^\prime\boldsymbol{\Omega}_{+,p,s}\mathbf{c}+ o(1), \qquad\boldsymbol{\Omega}_{+,p,s}:={\boldsymbol{\Gamma}}_{+,p,s}^{-1}\mathbf{V}_{+,p,s}
{{\boldsymbol{\Gamma}}_{+,p,s}^{-1\prime}}.
\end{align*}
Positive definiteness of $\boldsymbol{\Omega}_{+,p,s}$ follows from the fact that it is a product of symmetric positive definite matrices (Lemma \ref{lemma: grand hessian} and Lemma \ref{lemma: first order variance}). Finally, by appropriately scaling $\dot{\phi}_{+,p,s}(h)$ by its variance, we can rewrite it as
\[\sum_{i=1}^nw_{+,i}(h)u_{+,i}(h),\]
with
\[w_{+,i}(h) := \left(\mathbf{c}^\prime\boldsymbol{\Omega}_{+,p,s}\mathbf{c}\right)^{-1/2}\mathbf{c}^\prime\widetilde{\boldsymbol{\Gamma}}_{+,p,s}^{-1}(h)\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)/\sqrt{nh}.\]
Therefore, $\{{w}_{+,i}(h)u_{+,i}(h)\}_{i=1}^n$ is a triangular array of mean-zero row-wise independent summands with variance converging to 1 to which we can apply a Lindeberg-Feller CLT. To do so, we verify the Lindeberg condition by showing that the stronger Lyapunov condition holds for the fourth moment. To see this, by Assumption \ref{ass: DGP}
\begin{align*}
\sum_{i=1}^n\operatorname{\mathbb{E}}[|{w}_{+,i}(h)u_i(h)|^4]&\lesssim \frac{1}{nh^2}\cdot \int\limits_c^\infty \left|\breve{\mathbf{e}}_{\nu}^\prime\widetilde{\boldsymbol{\Gamma}}_{+,p,s}^{-1}(h)K\left(\frac{x-c}{h}\right)\mathbf{r}_{p,s}\left(\frac{x-c}{h},\boldsymbol{\mu}_{W}(x)\right)\right|^4f(x)\,\mathrm{d} x = O\left(\frac{1}{nh}\right),
\end{align*}
where the final order comes from the change of variables $u=\frac{x-c}{h}$ and the fact that the resulting integral is bounded by Assumptions \ref{ass: density running variable}, \ref{ass: DGP}, and \ref{ass: kernel}. Therefore, we can conclude that
\[\sqrt{nh}\mathbf{c}^\prime\mathbf{H}_{p,s}(h)\widehat{\boldsymbol{\Gamma}}_{+,p,s}^{-1}(h)\mathbf{L}_{+,p,s}(h) \rightsquigarrow \mathsf{N}(0, \mathbf{c}^\prime\boldsymbol{\Omega}_{+,p,s}\mathbf{c}).\]
As $\mathbf{c}$ was chosen arbitrarily, by the Cram\'{e}r-Wold device we get
\[\sqrt{nh}\mathbf{H}_{p,s}(h)\widehat{\boldsymbol{\Gamma}}_{+,p,s}^{-1}(h)\mathbf{L}_{+,p,s}(h) \rightsquigarrow \mathsf{N}(\mathbf{0}_{1+p+d}, \boldsymbol{\Omega}_{+,p,s}),\]
which was to be shown.
\end{proof}
\subsection{Proof of Theorem \ref{theorem: variance estimation plug-in}}
\begin{proof}\label{proof: variance estimation plug-in}
First, note that under the assumptions of Theorem \ref{thm: asymptotic distribution}, we get that
\begin{align}\label{eq: vanishing leverage}
\max_{i\in[n]}L_{-,i}=o_{\mathbb{P}}(1)\qquad\text{and}\qquad \max_{i\in[n]}L_{+,i}=o_{\mathbb{P}}(1).
\end{align}
The result above follows from Assumptions \ref{ass: sampling}, \ref{ass: DGP}, \ref{ass: kernel}, and Lemma \ref{lemma: grand hessian}.
Define $\omega_{+,i}(h):= \mathbbm{1}(X_i\geq 0)K\left(\frac{X_i-c}{h}\right)\sqrt{w_{+,i}(h)}$. First, note that
\[\widehat{u}_{+,i}(h) = u_{+,i}(h) + \widehat{v}_{+,i}(h),\]
where $u_{+,i}(h):= Y_i - \mathbf{r}_{p,s}(X_i-c,\mathbf{W}_i)^\prime\boldsymbol{\vartheta}^\star_{+,p,s}(h)$ and $\widehat{v}_{+,i}(h):= \mathbf{r}_{p,s}(X_i-c,\mathbf{W}_i)^\prime({\boldsymbol{\vartheta}}^\star_{+,p,s}(h)- \widehat{\boldsymbol{\vartheta}}_{+,p,s}(h))$. Then, focusing on the ``meat'' part of the estimator, we get
\begin{align*}
\widehat{\mathbf{V}}_{+,p,s}(h) &= \frac{1}{nh}\sum_{i=1}^n\omega_{+,i}(h)\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) \\
&= \frac{1}{nh}\sum_{i=1}^n\omega_{+,i}(h)\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 u_{+,i}^2(h) \tag{$:=\mathbf{V}_{+,1}(h)$} \\
&\quad + \frac{1}{nh}\sum_{i=1}^n\omega_{+,i}(h)\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{v}_{+,i}^2(h) \tag{$:=\mathbf{V}_{+,2}(h)$}, \\
&\quad + \frac{2}{nh}\sum_{i=1}^n\omega_{+,i}(h)\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{v}_{+,i}(h){u}_{+,i}(h) \tag{$:=\mathbf{V}_{+,3}(h)$}.
\end{align*}
For convenience, let $\mathbf{V}_{+,1}^{\mathtt{HC0}}(h)$ denote $\mathbf{V}_{+,1}(h)$ when $w_{+,i}(h) = 1$. Note that $\mathbf{V}^{\mathtt{HC0}}_{+,1}(h) \overset{\mathbb{P}}{\longrightarrow} \mathbf{V}_{+,p,s}$ by the weak law of large numbers. Then, it is immediate to see that $\mathbf{V}_{+,1}^{\mathtt{HC1}}(h) = \frac{n-k}{n}\mathbf{V}_{+,1}^{\mathtt{HC0}}(h)$ and so $\mathbf{V}_{+,1}^{\mathtt{HC1}}(h)\overset{\mathbb{P}}{\longrightarrow} \mathbf{V}_{+,p,s}$. Regarding the HC2-type estimator, it suffices to note that
\[|\mathbf{V}_{+,1}^{\mathtt{HC2}}(h)-\mathbf{V}_{+,1}^{\mathtt{HC0}}(h)| \leq \left(\frac{1}{nh}\sum_{i=1}^n\left|\mathbf{r}_{p,s}\left(\frac{X_i-c}{h},\mathbf{W}_i\right)\right|^2 u_{+,i}^2(h)\right) \cdot \left|\frac{1}{1-\max_{i\in[n]}L_i} - 1\right|,\]
where the first term is $O_{\mathbb{P}}(1)$ by Assumptions \ref{ass: sampling}-\ref{ass: DGP} and the second term is $o_{\mathbb{P}}(1)$ due to \eqref{eq: vanishing leverage}. A similar argument works for the HC3-type estimator.
At this point, if we show that $\mathbf{V}_{+,j}(h)\overset{\mathbb{P}}{\longrightarrow} 0, j=2,3,$ then we have established the consistency of our estimator for the variance of the RD estimator. Let's start with $\mathbf{V}_{+,2}(h)$:
\begin{align*}
\left|\mathbf{V}_{+,2}(h)\right| &\leq
{\frac{1}{nh}\sum_{i=1}^n|\omega_{+,i}(h)| \left|\mathbf{r}_{p,s}\left(\frac{X_i-c}{h},\mathbf{W}_i\right)\right|^2\left|\mathbf{r}_{p,s}\left(\frac{X_i-c}{h},\mathbf{W}_i\right)\right|^2\left|\mathbf{H}_{p,s}(h)\left(\widehat{\boldsymbol{\vartheta}}_{+,p,s}(h)-\boldsymbol{\vartheta}^\star_{+,p,s}(h)\right)\right|^2} \\
&\leq \left|\mathbf{H}_{p,s}(h)\left(\widehat{\boldsymbol{\vartheta}}_{+,p,s}(h)-\boldsymbol{\vartheta}^\star_{+,p,s}(h)\right)\right|^2\cdot\max_{i\in[n]}\left|\omega_{+,i}(h)\right|\cdot\frac{1}{nh}\sum_{i=1}^n\left|\mathbf{r}_{p,s}\left(\frac{X_i-c}{h},\mathbf{W}_i\right)\right|^4.
\end{align*}
The first term on the right-hand side is $O_{\mathbb{P}}(n^{-1}h^{-1})$ by Theorem \ref{thm: asymptotic distribution}, the second term is $O_{\mathbb{P}}(1)$ by \eqref{eq: vanishing leverage}, whereas the third term is $O_{\mathbb{P}}(1)$ as it converges in probability to its expectation which is
\[\int\limits_0^\infty f(uh+c)K(u) \left|\mathbf{r}_{p,s}(u,\boldsymbol{\mu}_{W}(uh+c))\right|^4\,\mathrm{d} u,\]
which is finite because of Assumptions \ref{ass: density running variable}, \ref{ass: DGP}, and \ref{ass: kernel}. Therefore, $\mathbf{V}_{+,2}(h) = o_{\mathbb{P}}(1)$.
Then
\begin{align*}
|\mathbf{V}_{+,3}(h)|
&\leq \left|\mathbf{H}_{p,s}(h)\left(\widehat{\boldsymbol{\vartheta}}_{+,p,s}(h)-\boldsymbol{\vartheta}_{+,p,s}^\star(h)\right)\right| \cdot\max_{i\in[n]}\left|\omega_{+,i}(h)\right|\cdot \frac{1}{nh}\sum_{i=1}^n\left|\mathbf{r}_{p,s}\left(\frac{X_i-c}{h},\mathbf{W}_i\right)\right|^3 |u_{+,i}(h)|.
\end{align*}
The first term on the right-hand side is $O_{\mathbb{P}}(1/\sqrt{nh})$ by Theorem \ref{thm: asymptotic distribution}. The second term is $O_{\mathbb{P}}(1)$ by\eqref{eq: vanishing leverage}. To see that the third term is $O_{\mathbb{P}}(1)$, note that by the H{\"o}lder's inequality
\[\operatorname{\mathbb{E}}\left[\left|\mathbf{r}_{p,s}\left(\frac{X_i-c}{h},\mathbf{W}_i\right)\right|^3 |u_{+,i}(h)|\right] \leq \left(\int\limits_0^\infty f(uh+c)K(u)\left|\mathbf{r}_{p,s}\left(u,\boldsymbol{\mu}_{W}(uh+c)\right)\right|^3 \,\mathrm{d} u\right)^{3/4} \operatorname{\mathbb{E}}[u_{+,i}^4(h)]^{1/4},\]
which is finite because of Assumptions \ref{ass: density running variable}, \ref{ass: DGP}, and \ref{ass: kernel}. Therefore, $\mathbf{V}_{+,3}(h) = O_{\mathbb{P}}(1/\sqrt{nh})=o_{\mathbb{P}}(1)$.
In conclusion, we showed that
\[\widehat{\mathbf{V}}_{+,p,s}(h) = \sum_{\ell=1}^3\mathbf{V}_{+,\ell}(h),\]
where $\mathbf{V}_{+,1}(h)$ converges in probability to $\mathbf{V}_{+,p,s}(h)$ and the other terms are all $o_{\mathbb{P}}(1)$. Moreover, Lemma \ref{lemma: grand hessian} gives us $\widehat{\boldsymbol{\Gamma}}_{+,p,s}^{-1}(h)\overset{\mathbb{P}}{\longrightarrow} {\boldsymbol{\Gamma}}_{+,p,s}^{-1}.$ Hence, by Slutsky's theorem
\[\widehat{\mathcal{V}}_{+,\nu,p,s}(h) \overset{\mathbb{P}}{\longrightarrow} {\mathcal{V}}_{+,\nu,p,s}. \qedhere \]
\end{proof}
\subsection{Proof of Lemma \ref{lemma: smoothing bias}}
\begin{proof}\label{proof: smoothing bias}
Let $q:=p\vee s$. First of all, if we take a $(q +2)$-th order Taylor expansion of $\mu_+(\cdot,\mathbf{w})$ around $c$, we get
\begin{align*}
{\boldsymbol{\vartheta}}^\star_{+,p,s}(h)&= \mathbf{H}_{p,s}^{-1}(h)\widetilde{\boldsymbol{\Gamma}}_{+,p,s}^{-1}(h)\operatorname{\mathbb{E}}\left[\frac{1}{h}\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)\mu(X_i,\mathbf{W}_i)\right]\\
&=\boldsymbol{\vartheta}_{+,p,s} +\mathbf{H}_{p,s}^{-1}(h)\widetilde{\boldsymbol{\Gamma}}_{+,p,s}(h)^{-1} \left(h^{p+1}\widetilde{\boldsymbol{\zeta}}_{+,p,s,p}(h) \frac{\alpha^{(p+1)}_{+}(c)}{(p+1)!} + h^{p+2}\widetilde{\boldsymbol{\zeta}}_{+,p,s,p+1}(h) \frac{\alpha^{(p+2)}_{+}(c)}{(p+2)!} \right) \\
&\:\:+ \mathbf{H}_{p,s}^{-1}(h)\widetilde{\boldsymbol{\Gamma}}_{+,p,s}(h)^{-1} \left(h^{s+1}\widetilde{\boldsymbol{\varphi}}_{+,p,s,s}(h) \frac{\boldsymbol{\lambda}^{(s+1)}_{+}(c)}{(s+1)!} + h^{s+2}\widetilde{\boldsymbol{\varphi}}_{+,p,s,s+1}(h) \frac{\boldsymbol{\lambda}^{(s+2)}_{+}(c)}{(s+2)!} \right) \\
&\:\:+ o(h^{2+p\wedge s}),
\end{align*}
where
\begin{align*}
\widetilde{\boldsymbol{\zeta}}_{+,p,s,a}(h)&=\int\limits_0^\infty K(u)\mathbf{r}_{p,s}(u,\boldsymbol{\mu}_{W}(uh+c))u^{a+1} f(uh+c)\,\mathrm{d} u,\\
\widetilde{\boldsymbol{\varphi}}_{+,p,s,a}(h)&=\int\limits_0^\infty K(u) \begin{bmatrix}
\mathbf{r}_p(u)\boldsymbol{\mu}_{W}(uh+c)^\prime \\
\boldsymbol{\mu}_{WW}(uh+c)\otimes \mathbf{r}_s(u)
\end{bmatrix} u^{a+1} f(uh+c)\,\mathrm{d} u.
\end{align*}
The $o(h^{2+p\wedge s})$ term comes from the integral form of the remainder of the Taylor expansion. Suppose $p=s$. The remainder is of the form
\begin{align*}
\mathbf{v}_{+}(h)=\int\limits_0^\infty\int\limits_{\operatorname{\mathbb{R}}^d} K\left(u\right)\mathbf{r}_{p,s}\left(u,\mathbf{w}\right)\int\limits_c^{uh+c}\Bigg(\frac{\mu_+^{(p+3)}(t, \mathbf{w})}{(p+2)!}(uh+c-t)^{p+2}\Bigg)\,\mathrm{d} tf_{X,W}(uh+c,\mathbf{w})\,\mathrm{d} \mathbf{w}\,\mathrm{d} u.
\end{align*}
Being $\mu_+^{(p+3)}(x,\mathbf{w})$ continuous by Assumption \ref{ass: SRD continuous Y}, $\mathsf{k}\geq h$, and using Assumptions \ref{ass: density running variable} and \ref{ass: kernel}, we get $\mathbf{v}_+(h) = o(h^{2+p})$.
Whenever $s\neq p$, the remainder $\mathbf{v}_+(h)$ includes additional terms we ignored in the linearization above. These terms will always include either $\alpha^{(\ell)}(c), \ell \geq p + 3$ or $\lambda^{(\ell)}(c), \ell \geq s + 3$ and be of order $o(h^{2+p})$ and $o(h^{2+s})$, respectively.
The last part of the lemma follows by taking the limit as $nh\to\infty$ and $h\to 0$ and relying on Lemma \ref{lemma: grand hessian} and \ref{lemma: vartheta hat order}.
\end{proof}
\bibliography{CCFPT_2025_HTERD--bib}
\bibliographystyle{jasa}