# Inference in Regression Discontinuity Designs with High-Dimensional
Covariates

Alexander Kreiß  Christoph Rothe

###### Abstract

We study regression discontinuity designs in which many predetermined
covariates, possibly much more than the number of observations, can be used to
increase the precision of treatment effect estimates. We consider a two-step
estimator which first selects a small number of “important” covariates through
a localized Lasso-type procedure, and then, in a second step, estimates the
treatment effect by including the selected covariates linearly into the usual
local linear estimator. We provide an in-depth analysis of the algorithm’s
theoretical properties, showing that, under an approximate sparsity condition,
the resulting estimator is asymptotically normal, with asymptotic bias and
variance that are conceptually similar to those obtained in low-dimensional
settings. Bandwidth selection and inference can be carried out using standard
methods. We also provide simulations and an empirical application.

##  1 Introduction

Regression discontinuity (RD) designs are widely used for estimating causal
effects from observational data in economics and other social sciences. These
designs exploit institutional settings in which a unit’s treatment assignment
is determined by whether its realization of a running variable falls above or
below some known cutoff value. Estimates of causal effects are then obtained
by comparing the outcomes of units that are close to, but on different sides
of the cutoff. Methods based on local linear regression are widely used in
this context, and their theoretical properties have been studied extensively
in the literature (e.g. Hahn et al., 2001; Imbens and Kalyanaraman, 2012;
Calonico et al., 2014; Armstrong and Kolesár, 2018).

While an empirical RD study can be carried out using only data on the outcome
and the running variable, in practice researchers often want to incorporate
additional covariates into their analysis to improve the precision of their
estimates. This is commonly done by including the covariates linearly and
without localization in a standard local linear RD regression (Calonico et
al., 2019). Such linear adjustment estimators generally work well if the
number of available covariates is small relative to the sample size. However,
they might yield distorted inference even with a rather moderate number of
covariates; and they are generally expected to break down in high-dimensional
settings where the number of covariates is large and possibly even exceeds the
number of observations. Such settings can occur, for example, when working
with rich administrative data sets, but also if a large number of
transformations, like interactions or polynomials, is applied to a low-
dimensional set of underlying covariates.

In this paper, we study a two-step approach that addresses this problem. In
the first step, we select a small subset of the covariates by adding an
$\ell_{1}$ or Lasso penalty (cf. Tibshirani, 1996) to the local least squares
problem that defines the linear adjustment estimator, and collect those
variables with non-zero coefficient estimates. By construction, the selected
covariates are strongly related to the outcome, and thus have the greatest
potential to “absorb” some of its variance. In the second step, we then
compute a standard linear adjustment estimator, but use only the selected
covariates. We show that the resulting “post-Lasso” estimator is
asymptotically normal under an approximate sparsity condition, with asymptotic
bias and variance that are conceptually similar to those obtained in low-
dimensional settings (here “approximate sparsity” means that only a small
number of the covariates is particularly relevant for the empirical analysis,
in the sense that including further covariates would not lead to meaningful
improvements of estimation accuracy). We also argue that one can use standard
methods for bandwidth selection and inference with the selected variables,
making the method very easy to implement in practice.

Our estimator has many parallels with the well-known post-Lasso approach for
treatment effect estimation under unconfoundedness with high-dimensional
controls in Belloni et al. (2013), including the use of a similar notion of
approximate sparsity. One important difference is that, in contrast to
unconfoundedness, it is not necessary to select the “right” covariates in our
RD framework in order to obtain a consistent estimator. This is because in our
setting the purpose of controlling for covariates is only to increase
efficiency, and not to address issues of selection bias.111This is a
conceptual parallel between our setting and the use of covariates in
randomized controlled experiments with a constant propensity score (e.g. Wager
et al., 2016). Note, however, that there is no explicit notion of random
assignment in RD designs, and thus results from the literature on experiments
do not simply carry over to our setting.  Our method therefore only requires a
single selection step that collects variables which are strongly related to
the outcome, and not a “double selection” as in Belloni et al. (2013) that
also selects variables related to treatment status. Our variable selection
step is also not based on the standard Lasso, but on a “Lasso-penalized” local
linear regression problem; and in contrast to unconfoundedness, one cannot
make use of (conditionally) random treatment assignment in RD settings, but
only exploit continuity conditions. The theoretical results that we derive in
this paper therefore do not follow from existing arguments.

Our paper contributes to a growing literature that considers covariates in RD
designs, including Armstrong and Kolesár (2018), Calonico et al. (2019),
Frölich and Huber (2019), Noack et al. (2021) and Arai et al. (2022). In
particular, Arai et al. (2022) study an estimation approach that is very
similar to the one we consider in this paper. However, their analysis relies
on strong conditions, like exact sparsity and a “$\beta$-min” condition that
puts a large lower bound on the coefficients of relevant covariates, which
imply perfect model selection.222Roughly speaking, the conditions in Arai et
al. (2022) describe a setup in which the covariates can be partitioned into
groups of “very important” and “completely irrelevant” ones, irrespective of
the chosen bandwidth. Moreover, the influence of the “very important” ones is
assumed to be large enough that they are selected by a localized Lasso
procedure with near-certainty. Such conditions seem unlikely to be satisfied
in practice. We also note that the estimator studied by Arai et al. (2022)
differs from ours in that it does not use all the covariates with non-zero
coefficients in the first stage for the “post-Lasso” step, but only those
whose estimated coefficients exceed some positive bound, which depends on an
additional tuning parameter.  Our paper uses an arguably more realistic
framework, does not require perfect model selection, and develops a complete
asymptotic theory for the final RD estimator. Our paper’s technical arguments
are also related to those in Su et al. (2019), who use a localized Lasso to
handle high-dimensional covariates in a nonparametric setup, namely a
continuous treatment model.

The remainder of this paper is structured as follows. In Section 2 we
introduce our model and our proposed estimator, and give an informal
description of its theoretical properties. In Section 3 we give the
assumptions for our theoretical analysis, and state and discuss our main
result. Section 4 explains some implementation details of our procedure, and
gives the results of a simulation study and an empirical illustration. Section
5 concludes. All proofs are collected in the Online Appendix.

##  2 Setup and Method

###  2.1 Setup and Preliminaries

Consider a sharp RD design to determine the causal effect of a binary
treatment on some outcome variable of interest. The data are an independent
sample $\\{(Y_{i},X_{i},Z_{i}),i=1,\ldots,n\\}$ of size $n$ from some large
population. Here $Y_{i}\in\mathbb{R}$ is the outcome variable,
$X_{i}\in\mathbb{R}$ is the running variable, and $Z_{i}\in\mathbb{R}^{p}$ is
a vector of pre-treatment covariates. We particularly consider high-
dimensional settings in which the covariate dimension $p$ can be large
relative to, or indeed significantly larger than, the sample size $n$. We
account for this in our framework by allowing $p=p_{n}$ to increase with the
number of observations. High-dimensional covariates occur of course if the
researcher observes a large number of conceptually distinct variables for each
unit, but also if the researcher applies a large number of transformations
from a dictionary of basis functions, that might create interactions or
polynomials, to an underlying low-dimensional vector of covariates.

Units receive the treatment if and only if the running variable exceeds some
known cutoff, which we normalize to zero without loss of generality. We denote
the resulting treatment indicator by $T_{i}$, so that
$T_{i}=\mathbbm{1}(X_{i}\geq 0)$. Units also have potential outcomes
$Y_{i}(t)$, for $t\in\\{0,1\\}$, corresponding to the outcome unit $i$ would
have experienced had it received treatment $t$, so that $Y_{i}=Y_{i}(T_{i})$.
The parameter of interest is the average treatment effect among units at the
cutoff:

$\tau_{Y}=\mathbb{E}(Y_{i}(1)-Y_{i}(0)|X_{i}=0).$

If $\mathbb{E}(Y_{i}(t)|X_{i}=x)$ is continuous around the cutoff for
$t\in\\{0,1\\}$, this parameter is identified by the jump in the conditional
expectation function $\mathbb{E}(Y_{i}|X_{i}=x)$ of the observed outcome given
the running variable at the threshold:

$\tau_{Y}=\lim_{x\downarrow 0}\mathbb{E}(Y_{i}|X_{i}=x)-\lim_{x\uparrow
0}\mathbb{E}(Y_{i}|X_{i}=x).$ (2.1)

Local linear regression (Fan and Gijbels, 1996) is arguably the most popular
framework for estimation and inference in RD designs. In the absence of
covariates, the jump $\tau_{Y}$ is estimated by fitting a linear regression of
$Y_{i}$ on $X_{i}$ locally around the cutoff, allowing for different
intercepts and slopes on each side. This estimator is the baseline procedure
for our analysis:

$\hat{\tau}_{h,\textnormal{Base}}=e_{2}^{\top}\underset{\theta\in\mathbb{R}^{4}}{\operatorname{argmin}}\sum_{i}^{n}K_{h}\left(X_{i}\right)\left(Y_{i}-V_{i}^{\top}\theta\right)^{2},$
(2.2)

with $K$ a non-negative kernel function, $h>0$ a bandwidth,
$K_{h}(x)=K(x/h)/h$, $V_{i}=\left(1,T_{i},X_{i}/h,T_{i}X_{i}/h\right)^{\top}$
a vector of appropriate transformations of the running variable, and
$e_{2}=(0,1,0,\ldots,0)^{\top}$ a unit vector of appropriate length. As
discussed in Calonico et al. (2019), practitioners often augment the local
regression in (2.2) with additional covariates in a simple linear fashion,
which yields the linear adjustment estimator

$\hat{\tau}_{h,\textnormal{CCFT}}=e_{2}^{\top}\underset{(\theta,\gamma)\in\mathbb{R}^{4+p}}{\operatorname{argmin}}\sum_{i}^{n}K_{h}\left(X_{i}\right)\left(Y_{i}-V_{i}^{\top}\theta-
Z_{i}^{\top}\gamma\right)^{2}.$ (2.3)

This estimator is consistent under standard regularity conditions if the
dimension of the covariates is fixed and if their conditional distribution
given the running variable changes smoothly around the cutoff, in the sense
that the conditional expectation of the covariates given the running variable
does not jump:

$\tau_{Z}=\lim_{x\downarrow 0}\mathbb{E}(Z_{i}|X_{i}=x)-\lim_{x\uparrow
0}\mathbb{E}(Z_{i}|X_{i}=x)=0.$ (2.4)

The linear adjustment estimator is typically more efficient than the baseline
“no covariates” estimator. It is not uniquely defined, however, if the number
of local parameters exceeds the number of observations that receive positive
kernel weights in (2.3). Moreover, due to overfitting, asymptotic
approximations based on a “fixed $p$” analysis might not provide adequate
descriptions of the estimator’s finite sample properties even in settings
where the number of covariates is moderate relative to the effective sample
size. For instance, in our simulations below we illustrate that conventional
standard errors might severely underestimate the true variability of the
linear adjustment estimator in a setting with 10–50 covariates and 1,000 data
points. Linear adjustment estimators are therefore only appropriate for very
low-dimensional settings.

###  2.2 Proposed Method

A natural way to extend linear adjustment estimators to high-dimensional
settings is to consider versions that only use a “small” active subset of the
available covariates. Formally, with
$J=\\{j_{1},\ldots,j_{s}\\}\subset\\{1,\ldots,p\\}$ a generic subset of the
covariates’ indices of size $s\equiv|J|\ll p$, and
$Z_{i}(J)=(Z^{(j_{1})}_{i},\ldots,Z^{(j_{s})}_{i})^{\top}$ the $s$-dimensional
vector of components of $Z_{i}$ whose indices are collected in $J$, such
estimators are given by

$\hat{\tau}_{h}(J)=e_{2}^{\top}\underset{(\theta,\gamma)\in\mathbb{R}^{4+s}}{\operatorname{argmin}}\sum_{i}^{n}K_{h}\left(X_{i}\right)\left(Y_{i}-V_{i}^{\top}\theta-
Z_{i}(J)^{\top}\gamma\right)^{2}.$ (2.5)

Using arguments from Calonico et al. (2019), it is easily seen that such
estimators are consistent for any fixed covariate subset $J$ under the
appropriate regularity conditions. The choice of $J$ does affect the
asymptotic variance, however, and using covariates that have high correlation
with the outcome (locally at the cutoff) can generally be expected to yield
more efficient estimates of $\tau_{Y}$. In practice, the identity of these
“most useful” covariates is typically not known _a priori_ , but can
potentially be inferred in a data driven way. We therefore consider estimators
of the form $\hat{\tau}_{h}(\hat{J}_{n})$, with $\hat{J}_{n}$ a data-dependent
subset of the covariates’ indices that is intended to contain the most
relevant ones.

Our proposed RD estimator for settings with high-dimensional covariates
determines the set $\hat{J}_{n}$ of active covariates through a “localized”
version of a Lasso regression in a preliminary model selection step. We
consider a version of the minimization problem in (2.3) that includes an
additional penalty on the sum of the absolute values of the coefficients
associated with the (appropriately standardized) covariates, and define
$\hat{J}_{n}$ as the set of covariate indices for which the corresponding
coefficient estimate is non-zero. Specifically, our procedure is as follows.

  1. 1.

Using a preliminary bandwidth $b$ and a penalty parameter $\lambda$, solve the
following “Lasso version” of the weighted least squares problem in (2.5):

$\left(\widetilde{\theta}_{n},\widetilde{\gamma}_{n}\right)=\underset{(\theta,\gamma)\in\mathbb{R}^{4+p_{n}}}{\operatorname{argmin}}\sum_{i=1}^{n}K_{b}(X_{i})\left(Y_{i}-V_{i}^{\top}\theta-\left(Z_{i}-\hat{\mu}_{Z,n}\right)^{\top}\gamma\right)^{2}+\lambda\sum_{k=1}^{p_{n}}\hat{w}_{n,k}|\gamma_{k}|,$

where

$\displaystyle\hat{\mu}_{Z,n}=\frac{1}{n}\sum_{i=1}^{n}Z_{i}K_{b}\left(X_{i}\right)\textnormal{
and
}\hat{w}_{n,k}^{2}=\frac{b}{n}\sum_{i=1}^{n}\left(K_{b}\left(X_{i}\right)Z_{i}^{(k)}-\mu_{Z,n}^{(k)}\right)^{2}$

are the local sample mean and variance, respectively, of the covariates. Note
that standardizing the covariates allows the penalty parameter $\lambda$ to be
reasonably tuned to all coefficients simultaneously.

  2. 2.

Using a final bandwidth $h$, compute the restricted post-Lasso estimate of
$\tau_{Y}$ as $\hat{\tau}_{h}(\hat{J}_{n})$ as in (2.5), where
$\hat{J}_{n}=\\{k\in\\{1,...,p_{n}\\}:\widetilde{\gamma}_{n}^{(k)}\neq 0\\}$
is the set of the indices of those covariates selected in the first step.

The tuning parameters $b$ and $\lambda$ do not appear in the asymptotic
distribution of $\hat{\tau}_{h}(\hat{J}_{n})$, and hence some _ad hoc_ choices
are needed. To choose $b$, we recommend using a method for bandwidth choice
designed for the baseline estimator without covariates, like the ones proposed
in Imbens and Kalyanaraman (2012), Calonico et al. (2014) or Armstrong and
Kolesár (2018). For $\lambda$, we compare different approaches in our
simulations, and the results suggest that the plug-in procedure of Belloni et
al. (2013) works well in practice. The final bandwidth $h$ can be chosen via
any method suitable for linear adjustment estimators with low-dimensional
covariates, such as those discussed in Calonico et al. (2019) or Armstrong and
Kolesár (2018).

###  2.3 Overview of Main Result

We now give an informal overview of the main theoretical result in this paper.
For generic random vectors $A$ and $B$, we use the notation that
$\mu_{A}(x)=\mathbb{E}(A|X=x)$, $\mu_{AB}(x)=\mathbb{E}(AB^{\top}|X=x)$,
$\sigma^{2}_{AB}(x)=\mu_{AB}(x)-\mu_{A}(x)\mu_{B}(x)^{\top}$; and write
$\sigma_{A}^{2}(x)=\sigma_{AA}^{2}(x)$ for simplicity. For a generic function
$f$, we also write $f_{+}\ =\lim_{x\downarrow 0}f(x)$ and $f_{-}\
=\lim_{x\uparrow 0}f(x)$ for its right and left limit at zero, respectively,
so that $\tau_{Y}=\mu_{Y+}-\mu_{Y-}$, for example.

A key assumption for our analysis is that the covariates satisfy an
_approximate sparsity_ condition, which intuitively means that only a small
subset of the covariates is particularly relevant for the empirical analysis,
and that including any further covariates would not lead to meaningful
improvements of estimation accuracy. To state this notion more formally, we
define the following population regression coefficients and corresponding
residuals for any $J\subset\\{1,\ldots,p_{n}\\}$ and bandwidth $h$:

$\begin{array}[]{l}(\theta_{0}(J,h),\gamma_{0}(J,h))=\underset{(\theta,\gamma)}{\operatorname{argmin}}\quad\mathbb{E}\left(K_{h}(X_{i})\left(Y_{i}-V_{i}^{\top}\theta-
Z_{i}(J)^{\top}\gamma\right)^{2}\right),\\\
r_{i}(J,h)=Y_{i}-V_{i}^{\top}\theta_{0}(J,h)-Z_{i}(J)^{\top}\gamma_{0}(J,h).\end{array}$
(2.6)

Approximate sparsity then means that there exist deterministic _target
covariate sets_ $J_{n}\subset\\{1,\ldots,p_{n}\\}$ that contain a “small”
number $s_{n}\equiv|J_{n}|\ll p_{n}$ of elements, and are such that the local
correlation between the corresponding regression errors $r_{i}(J_{n},h)$ and
_each_ component of $Z_{i}$ is small relative to the estimation error:

$\displaystyle\max_{j=1,...,p_{n}}\left|\mathbb{E}\left(K_{h}(X_{i})Z_{i}^{(j)}r_{i}(J_{n},h)\right)\right|=O\left(\sqrt{\frac{\log
p_{n}}{nh}}\right).$

Moreover, this condition needs to be satisfied for an appropriate range of
bandwidths, so that the sequence $J_{n}$ does not depend on the exact choice
of $h$.

Under this and other regularity conditions discussed below, one can show that
the post-Lasso estimator $\hat{\tau}_{h}(\hat{J}_{n})$ has the same first-
order asymptotic properties as an infeasible estimator $\hat{\tau}_{h}(J_{n})$
that uses the true target set, and then prove an asymptotic normality result
for the latter. Taken together, this yields the main result of our paper,
which is that the post-Lasso estimator $\hat{\tau}_{h}(\hat{J}_{n})$ of
$\tau_{Y}$ satisfies

$\displaystyle\frac{\sqrt{nh}\left(\hat{\tau}_{h}(\hat{J}_{n})-\tau_{Y}-h^{2}\mathcal{B}_{n}\right)}{\mathcal{S}_{n}}\overset{d}{\to}\mathcal{N}\left(0,1\right),$
(2.7)

with asymptotic bias and variance, respectively, such that

$\displaystyle\mathcal{B}_{n}$
$\displaystyle\approx\frac{C_{\mathcal{B}}}{2}\left(\mu_{\widetilde{Y}+}^{\prime\prime}-\mu_{\widetilde{Y}-}^{\prime\prime}\right)\quad\textnormal{and}\quad\mathcal{S}^{2}_{n}\approx\frac{C_{\mathcal{S}}}{f_{X}(0)}\left(\sigma^{2}_{\widetilde{Y}+}+\sigma^{2}_{\widetilde{Y}-}\right)$
(2.8)

in a sense made precise below. Here $C_{\mathcal{B}}$ and $C_{\mathcal{S}}$
are constants that depend on the kernel function $K$ only, and

$\widetilde{Y}_{i}=Y_{i}-Z_{i}(J_{n})^{\top}\gamma_{n},\textnormal{ with
}\gamma_{n}=\left(\sigma_{Z(J_{n})-}^{2}+\sigma_{Z(J_{n})+}^{2}\right)^{-1}\left(\sigma_{YZ(J_{n})-}^{2}+\sigma_{YZ(J_{n})+}^{2}\right),$

is a “covariate-adjusted” version of the outcome variable that uses a vector
$\gamma_{n}$ that can be thought of as an approximation of
$\gamma_{0}(J_{n},h)$ that is independent of the bandwidth. Our proposed
estimator is thus first-order asymptotically equivalent to a “baseline” sharp
RD estimator as in (2.2) with the covariate-adjusted outcome
$\widetilde{Y}_{i}$ replacing the original outcome $Y_{i}$.333Note, however,
that in contrast to the $Y_{i}$ the distribution of the $\widetilde{Y}_{i}$
depends on the sample size $n$, and hence the properties of this estimator do
not directly follow from existing ones for “baseline” sharp RD estimators.

Note that, as in Calonico et al. (2019), the continuity of $\mu_{Z}$ from
(2.4) is necessary to establish this result. If the components of $\mu_{Z}$
could potentially have a jump at the cutoff, the estimator
$\hat{\tau}_{h}(\hat{J}_{n})$ would generally not be consistent for
$\tau_{Y}$, but satisfy
$\hat{\tau}_{h}(\hat{J}_{n})=(\tau_{Y}-\tau_{Z}^{\top}\gamma_{n})(1+o_{P}(1)).$
In practice, researchers may want to investigate the plausibility of assuming
(2.4) by running a series of sharp RD regressions in which the covariates take
the role of the dependent variable, and testing whether the estimated jump at
the cutoff is significantly different from zero via some approach that is
appropriate for large-scale hypothesis testing (e.g. Benjamini and Hochberg,
1995). Alternatively, one could carry out this exercise only for the selected
covariates $Z_{i}(\hat{J}_{n})$.

The formulas for the bias and variance of $\hat{\tau}_{h}(\hat{J}_{n})$ in
(2.8) are analogous to those obtained in Calonico et al. (2019) for the case
that $J_{n}\equiv J$ contains only a fixed number of covariates. This suggests
that one can select the final bandwidth $h$ and create a confidence interval
for $\tau_{Y}$ by applying their proposed methods for low-dimensional setups
to the generated data set
$\\{(Y_{i},X_{i},Z_{i}(\hat{J}_{n})),i=1,\ldots,n\\}$ that only contains the
covariates selected by our algorithm. Similarly, given a bound on the second
derivative of the function $\mathbb{E}(\widetilde{Y}_{i}|X_{i}=x)$, one can
select the bandwidth $h$ and construct confidence intervals for $\tau_{Y}$ by
using the methods proposed by Armstrong and Kolesár (2018) with the generated
data set. See Section 4.1 for further discussion of implementation details.

##  3 Theoretical Analysis

###  3.1 Assumptions

We impose the following assumptions in our theoretical analysis.

Assumption (BW): (Bandwidth). _There are positive constants $c_{g,1},c_{g,2}$
such that $h,b\in[c_{g,1}g,c_{g,2}g]$, where $g\to 0$ is a reference sequence
such that $|J_{n}|\log p_{n}/\sqrt{ng}\to 0$ and $|J_{n}|g^{2}\sqrt{\log
p_{n}}\to 0$._

Assumption (BW) means that the bandwidths in the first and second stage of our
procedure are such that $b\asymp h$, i.e., $b$ and $h$ converge to zero with
the same speed. The exact role of this assumption is related to regularity
conditions in Assumption (MS) and is discussed below. We emphasize here that
(BW) can be achieved simply by selecting bandwidths from the range
$[c_{g,1}g,c_{g,2}g]$. The condition that $|J_{n}|\log p_{n}/\sqrt{ng}\to 0$
is a version of a standard assumption in the Lasso literature (cf. Chapter 6
in van de Geer and Bühlmann (2011)), adapted to our locally penalized setup.
The requirement that $|J_{n}|g^{2}\sqrt{\log p_{n}}\to 0$ is needed to control
for the bias. If the rate of the reference bandwidth $g$ is considered to be
given, Assumption (BW) can be seen as imposing restrictions on the maximal
number of covariates in the target set $J_{n}$. If, on the other hand, the
growth of $|J_{n}|$ is considered to be given, this assumption can be
interpreted as imposing limitations on the rate at which localization occurs.

Assumption (AS): (Approximate Sparsity). _It holds that $p_{n}\to\infty$ and,
with $r_{i}(J_{n},h)$ as in (2.6), that_

$\displaystyle\max_{k=1,...,p_{n}}\left|\mathbb{E}\left(Z_{n,i}^{(k)}K_{h}(X_{i})r_{i}(J_{n},h)\right)\right|=O\left(\sqrt{\frac{\log
p_{n}}{nh}}\right).$ (3.1)

_In addition, equation ( 3.1) remains true with $h$ replaced by $b$._

Assumption (AS) is similar in nature to the notion of approximate sparsity in,
for example, Belloni et al. (2013).444Note that (AS) is required to make the
Lasso work. Thus, it should not be read as a restriction, but it rather allows
for data-dependent model selection (cf. our Remark 3.8 below for alternatives
to (AS)). Online Appendix C.2 provides a numerical example for the behavior of
the Lasso in a non-sparse setting. Note that it follows from the definition of
$r_{i}(J_{n},h)$ that $\mathbb{E}(Z_{n,i}^{(k)}K_{h}(X_{i})r_{i}(J_{n},h))=0$
for all $k\in J_{n}$, and thus (3.1) only restricts the properties of
covariates that are not part of the target set. Intuitively, (AS) means that
the set $J_{n}$ contains “essentially” all relevant covariates, in the sense
that any covariate which is not contained in $J_{n}$ is almost locally
uncorrelated with the regression error $r_{i}(J_{n},h)$. Note that in a
setting with exact rather than approximate sparsity, condition (3.1) follows
automatically as, by definition, $Z_{n,i}^{(k)}$ is uncorrelated with
$(X_{i},r_{i}(J_{n},h))$ for $k\notin J_{n}$ in this case.

Assumption (D): (Differentiability). _The density of $X_{i}$, $f_{X}$, is
three times continuously differentiable in a neighborhood around zero and
$f_{X}(0)>0$. Moreover, $\mu_{Z}$ is continuous and uniformly bounded in a
neighborhood around zero. $\mu_{Z}$ and $\mu_{Y}$ are three times one-sided
differentiable at $0$, i.e., $\mu_{Z}^{\prime}$, $\mu_{Z}^{\prime\prime}$ and
$\mu_{Z}^{\prime\prime\prime}$ exist on $(-\infty,0)\cup(0,\infty)$ and the
left- and right sided limits at zero exist as well (and the same for
$\mu_{Y}$). The functions $\mu_{ZZ}$ and $\mu_{ZY}$ are one-sided
differentiable, and the derivatives fulfill_

$\displaystyle\sup_{n\in\mathbb{N}}\sup_{k\in\\{1,...,p_{n}\\}}\sup_{u\in[0,1]}\left|\mu_{Z^{(k)}}^{\prime}(uh)\right|+\left|\mu_{Z^{(k)}}^{\prime}(-uh)\right|<\infty,$
$\displaystyle\sup_{n\in\mathbb{N}}\sup_{k,l\in\\{1,...,p_{n}\\}}\sup_{u\in[0,1]}\left|\mu_{Z^{(k)}Z^{(l)}}^{\prime}(uh)\right|+\left|\mu_{Z^{(k)}Z^{(l)}}^{\prime}(-uh)\right|<\infty,$
$\displaystyle\sup_{n\in\mathbb{N}}\sup_{k\in
J_{n}}\sup_{u\in[0,1]}\left|\mu_{Z^{(k)}Y}^{\prime}(uh)\right|+\left|\mu_{Z^{(k)}Y}^{\prime}(-uh)\right|<\infty,$
$\displaystyle\sup_{n\in\mathbb{N}}\sup_{k\in\\{1,...,p_{n}\\}}\sup_{u\in[0,1]}\left|\mu_{Z^{(k)}}^{\prime\prime}(uh)\right|+\left|\mu_{Z^{(k)}}^{\prime\prime}(-uh)\right|<\infty,$
$\displaystyle\sup_{n\in\mathbb{N}}\sup_{k\in
J_{n}}\sup_{u\in[0,1]}\left|\mu_{Z^{(k)}}^{\prime\prime\prime}(uh)\right|+\left|\mu_{Z^{(k)}}^{\prime\prime\prime}(-uh)\right|<\infty.$

In a finite dimensional setting, like in Calonico et al. (2019), the above
conditions are implied by assuming existence and continuity of the one-sided
derivatives. In the high-dimensional setting, the uniformity assumption is
required in order to avoid pathological cases such as $\mu_{Z^{(k)}}$ getting
increasingly steep as $k\to\infty$. Note that the conditions on the third
derivative are only required on the target set $J_{n}$.

To state the next assumptions, we define the matrix

$M_{n}=\begin{pmatrix}\mu_{Z}(0)&\mathbf{0}&h\mu_{Z-}^{\prime}&h\left(\mu_{Z+}^{\prime}-\mu_{Z-}^{\prime}\right)\end{pmatrix}^{\top}\in\mathbb{R}^{4\times
p_{n}},$

where $\mathbf{0}$ denotes a vector of zeros and put
$\widetilde{Z}_{i}=Z_{i}-M_{n}^{\top}V_{i}$.

Assumption (TCS): (Target Covariate Set). _It holds that_

$\left\|\mathbb{E}\left(K_{h}(X_{i})\widetilde{Z}_{i}(J_{n})\widetilde{Z}_{i}^{\top}(J_{n})\right)^{-1}\right\|_{2}=O(1),$

_and there are finite numbers $\delta,\sigma_{l},\sigma_{r},C>0$ such that_

$\displaystyle\lim_{n\to\infty}\sup_{u\in[0,1]}\left|\mathbb{E}(r_{i}(J_{n},h)^{2}|X_{i}=uh)-\sigma_{r}^{2}\right|=0,$
(3.2)
$\displaystyle\lim_{n\to\infty}\sup_{u\in[0,1]}\left|\mathbb{E}(r_{i}(J_{n},h)^{2}|X_{i}=-uh)-\sigma_{l}^{2}\right|=0,$
(3.3)
$\displaystyle\sup_{x\in[-h,h]}\sup_{n\in\mathbb{N}}\left|\mathbb{E}(|r_{i}(J_{n},h)|^{2+\delta}|X_{i}=x)\right|<C,$
(3.4) $\displaystyle\sup_{n\in\mathbb{N}}\sup_{k\in
J_{n}}\sup_{u\in[0,1]}\left|\mathbb{E}\left(Z_{i}^{(k)}Y_{i}\big{|}X_{i}=uh\right)\right|+\left|\mathbb{E}\left(Z_{i}^{(k)}Y_{i}\big{|}X_{i}=-uh\right)\right|<\infty.$
(3.5)

_In addition, ( 3.2) and (3.3) also hold when $h$ is replaced by $b$._

We call the requirements (3.2) and (3.3) in (TCS) equi-continuity from the
right and left, respectively, and (3.4) and (3.5) are called uniform
boundedness. In the proof the asymptotic normality of our final RD estimator,
we use uniform boundedness to show a Lyapunov condition for the central limit
theorem, and equi-continuity to ensure that the respective asymptotic variance
converges to a finite and positive constant. While boundedness seems to be
unavoidable, the equi-continuity is assumed for convenience in our proofs.
Removing it would potentially lead to a different convergence rate for our
estimator by allowing settings where “almost” all variance is explained
through the covariates in the limit. Given that this is not a realistic
assumption, we do not consider adding this extra generality. Note that it is
necessary to distinguish the limits from left and right in (TCS) because the
conditional distribution of $r_{i}(J_{n},h)$ given the running variable may
experience a jump at zero.

Assumption (K): (Kernel). _The kernel $K:\mathbb{R}\to[0,\infty)$ integrates
to one, is continuous, symmetric and is supported on $[-1,1]$._

Such conditions on the kernel are standard in the literature, and satisfied by
the commonly used triangular and Epanechnikov kernels, for example. Kernels
with unbounded support, like the Gaussian kernel, could be accommodated at the
cost of slightly more involved theoretical arguments. Note that (K) implies
that the following quantities are finite:

$K^{(a)}=\int_{-\infty}^{\infty}u^{a}K(u)du,\quad
K_{+}^{(a)}=\int_{0}^{\infty}u^{a}K(u)du,\quad a\in\\{0,1,2,3,4\\}.$

For the following assumptions, we need further preliminary definitions. Let

$\displaystyle\mathbf{Y}=(Y_{1},...,Y_{n})^{\prime},\quad\mathbf{K}_{h}=\textrm{diag}\left(K(X_{1}/h)/h,\ldots,K(X_{n}/h)/h\right),$
$\displaystyle\mathbf{V}=\begin{pmatrix}1&T_{1}&X_{1}/h&T_{1}X_{1}/h\\\
\vdots&\vdots&\vdots&\vdots\\\
1&T_{n}&X_{n}/h&T_{n}X_{n}/h\end{pmatrix},\quad\mathbf{Z}(J)=\begin{pmatrix}Z_{1}(J)^{\top}\\\
\vdots\\\ Z_{n}(J)^{\top}\end{pmatrix},$

and for simplicity write $\mathbf{Z}=\mathbf{Z}(\\{1,...,p_{n}\\})$.

######  Definition 3.1.

Let $c>0$ and $J\subseteq\\{1,...,p_{n}\\}$, define

$k(c,J)=\inf\frac{|J|\frac{1}{n}\left\|\mathbf{K}_{b}^{\frac{1}{2}}\begin{pmatrix}\mathbf{V}&\mathbf{Z}\end{pmatrix}\begin{pmatrix}\theta\\\
\gamma\end{pmatrix}\right\|_{2}^{2}}{\left\|\begin{pmatrix}\theta&\gamma_{J}\end{pmatrix}^{\prime}\right\|_{1}^{2}},$
(3.6)

where the infimum is taken over all vectors
$(\theta\quad\gamma)^{\prime}\in\mathbb{R}^{p_{n}+4}$ for which

$\left\|\gamma_{J^{c}}\right\|_{1}\leq c\left\|\begin{pmatrix}\theta\\\
\gamma_{J}\end{pmatrix}\right\|_{1}$ (3.7)

and, for any $\gamma\in\mathbb{R}^{p_{n}}$, $\gamma_{J}^{(j)}=\gamma^{(j)}$
for $j\in J$ and $\gamma_{J}^{(j)}=0$ for $j\in J^{c}$. We say that the
compatibility condition $\textrm{CC}(c,J_{n})$ holds for a possibly random
sequence $J_{n}\subseteq\\{1,...,p_{n}\\}$ if $k(c,J_{n})^{-1}=O_{P}(1)$.

The constant $k(c,J)$ differs from the _compatibility constants_ known from
the classical Lasso literature (cf. Chapter 6.13 in van de Geer and Bühlmann
(2011)) only in the additional kernel weight and in the fact that the vector
$\theta$ is not penalized in our setup. In order to give some intuition, we
rewrite (3.6) as follows:

$k(c,J)=\frac{|J|}{nh}\inf\sum_{i=1}^{n}K\left(\frac{X_{i}}{h}\right)\left(V_{i}^{\top}\theta+Z_{i}(J)^{\top}\gamma_{J}-\left(-Z_{i}(J^{c})^{\top}\gamma_{J^{c}}\right)\right)^{2},$
(3.8)

where the infimum is taken over all pairs $(\theta,\gamma)$ for which (3.7)
and additionally $\|\theta\|_{1}+\|\gamma_{J}\|_{1}=1$ hold. Thus, $k(c,J)$ is
bounded away from zero if the covariates in $J^{c}$ with small coefficients
are unable to linearly represent the RD design vectors $V_{i}$ or the active
covariates $Z_{i}(J)$.

######  Definition 3.2.

Define for a sequence $m_{n}\in\mathbb{N}$ and set $J_{n}$

$\varphi(m_{n},J_{n})=\inf\frac{\frac{1}{n}\left\|\mathbf{K}_{h}^{\frac{1}{2}}\begin{pmatrix}\mathbf{V}&\mathbf{Z}\end{pmatrix}\alpha\right\|_{2}^{2}}{\|\alpha\|^{2}}\leq\sup\frac{\frac{1}{n}\left\|\mathbf{K}_{h}^{\frac{1}{2}}\begin{pmatrix}\mathbf{V}&\mathbf{Z}\end{pmatrix}\alpha\right\|_{2}^{2}}{\|\alpha\|^{2}}=:\Phi(m_{n},J_{n}),$
(3.9)

where $\inf$ and $\sup$ are taken over all vectors
$\alpha=(\theta,\gamma)^{\top}\in\mathbb{R}^{p_{n}+4}\setminus\\{0\\}$ such
that $|\\{i\in J_{n}^{c}:\,\gamma_{i}\neq 0\\}|\leq m_{n}$. We say that the
restricted sparse eigenvalue condition $\textrm{RSE}(m_{n},J_{n},h)$ holds for
a (random) sequence $m_{n}$ and a sequence of index sets $J_{n}$ if
$\varphi(m_{n},J_{n})^{-1}=O_{P}(1)$ and $\Phi(m_{n},J_{n})=O_{P}(1)$.

Following the pattern of the compatibility constant in Definition 3.1, we
extend the concept of restricted sparse eigenvalues to localized problems.
Continuing with the analogy, we can write down an equivalent formulation of
(3.9) in the fashion of (3.8) to see that CC and RSE are similar in terms of
their interpretation. The restricted sparse eigenvalue assumption is often
required when it comes to Lasso estimators. See for example Belloni et al.
(2013) for a discussion for non-localized estimators (Comment 3.2 therein) or
Lemma 1 in Belloni and Chernozhukov (2013). The localized case with the
additional kernel changes the problem to a conditional instead of an
unconditional variance.

Assumption (RSE & CC): (Restricted Sparse Eigenvalues and Compatibility). _The
matrix
$\widetilde{\mathbf{Z}}(\hat{J}_{n})^{\top}\mathbf{K}_{h}\widetilde{\mathbf{Z}}(J_{n})$
is almost surely invertible. The conditions RSE$(|J_{n}|\log n,J_{n},h)$ for
$\widetilde{\mathbf{Z}}$, RSE$(|J_{n}|\log n,J_{n},b)$, RSE$(0,J_{n},h)$ and
CC$(\bar{w},J_{n})$ hold true for $\bar{w}=3w^{(u)}/w^{(l)}$, where $w^{(l)}$
and $w^{(u)}$ are lower and upper bounds, respectively, on the weights
$\hat{w}_{n,k}$ given in Lemma A.16._

Recalling the discussion after Definitions 3.1 and 3.2, the assumption above
means that $Z_{i}(J_{n}^{c})$ cannot be used to represent $V_{i}$ or
$Z_{i}(J_{n})$. It can therefore be understood as excluding collinearity
between the covariates. In order to establish standard consistency results for
the Lasso (like Lemma A.18 in the Online Appendix) we only require the
compatibility condition. The restricted eigenvalue assumptions are required
for our asymptotic normality results, where they guarantee that the number of
selected covariates is growing slowly and that results about the model
selection step carry over to the RD step.

Recall the notation from (2.6) for the following assumption.

Assumption (MS): (Model Smoothness). _There is a sub-sequence
$J_{0,n}\subseteq J_{n}$, a sequence $\eta_{n}\to\infty$ and a constant $C>0$
such that for any bandwidth $\mathcal{H}$ which fulfills
$c_{g,1}g\leq\mathcal{H}\leq c_{g,2}g$, with $g,c_{g,1},c_{g,2}$ are as in
(BW), we have that for any $k\in J_{0,n}$_

$\left|\gamma_{0}^{(k)}(J_{n},\mathcal{H})\right|\geq\eta_{n}\sqrt{\frac{|J_{n}|\log
p_{n}}{ng}}$ (3.10)

_and for any $k\in J_{n}\setminus J_{0,n}$_

$\left|\gamma_{0}^{(k)}(J_{n},\mathcal{H})\right|\leq C\sqrt{\frac{|J_{n}|\log
p_{n}}{ng}}.$ (3.11)

_Moreover,_

$\sqrt{|J_{n}\setminus J_{0,n}|}\cdot\frac{|J_{n}|\log p_{n}}{\sqrt{ng}}\to
0\quad\textrm{and}\quad\sqrt{|J_{n}\setminus
J_{0,n}|}\cdot|J_{n}|g^{2}\sqrt{\log p_{n}}\to 0.$

This assumption rules out pathological settings in which a covariate’s
relevance within the target set is strongly affected by minor changes of the
bandwidth.555Note that (3.10) and (3.11) are mutually exclusive but not
exhaustive. It would be possible to formulate mutually exclusive conditions by
introducing a constant $C_{0}>0$ which depends on many unknown quantities.
This extra generality would thus bring no meaningful practical benefit. To see
this, fix $\mathcal{H}$ and consider sets $J_{0,n}(\mathcal{H})\subseteq
J_{n}$ such that (3.10) holds for $k\in J_{0,n}(\mathcal{H})$ while for $k\in
J_{n}\setminus J_{0,n}(\mathcal{H})$, (3.11) is true. Assumption (MS) reads
then as: The mapping $\mathcal{H}\mapsto J_{0,n}(\mathcal{H})$ is constant for
each $n\in\mathbb{N}$. In other words, the identity of the covariates with
large population regression coefficients in $J_{n}$ remains the same when the
bandwidth is slightly altered. There might also be covariates in the target
set $J_{n}$ with relatively small coefficients, namely those for which (3.11)
holds. Therefore, assumption (MS) is different from a $\beta$-min assumption
(cf. van de Geer and Bühlmann, 2011).

For the next assumption, we define for $k\in\\{1,...,p_{n}\\}$ and
$m\in\mathbb{N}$:

$\mu_{k,m}(x)=\mathbb{E}\left(\left|Z_{i}^{(k)}\right|^{m}\Big{|}X_{i}=x\right)\textrm{
and
}\mu_{k,m}^{(r)}(x)=\mathbb{E}\left(\left|Z_{i}^{(k)}r_{i}(J_{n},b)\right|^{m}\Big{|}X_{i}=x\right)$

Assumption (CTB): (Covariate Tail Behavior). _The functions
$\mu_{k,1},\mu_{k,2}$ and $\mu_{k,1}^{(r)}$ are uniformly bounded in a
neighborhood around zero. There are finite numbers
$\sigma_{a}^{2},c_{a},c_{a}^{*}>0$ for $a=0,1,2$ such that for all
$m\in\mathbb{N}$_

$\displaystyle\int_{\mathbb{R}}\left(1+|u|^{m}\right)K(u)^{m}\mu_{k,m}(ub)f_{X}(ub)du\leq$
$\displaystyle\frac{m!}{2}\sigma_{0}^{2}c_{0}^{m-2},$ (3.12)
$\displaystyle\int_{\mathbb{R}}K(u)^{2m}\mu_{k,2m}(ub)f_{X}(ub)du\leq$
$\displaystyle\frac{m!}{2}\sigma_{1}^{2}c_{1}^{m-2},$ (3.13)
$\displaystyle\int_{\mathbb{R}}(1+|u|^{m})K(u)\mu^{(r)}_{k,m}(ub)f_{X}(ub)du\leq$
$\displaystyle\frac{m!}{2}\sigma_{2}^{2}c_{2}^{m-2}.$ (3.14)

_Equations ( 3.12) and (3.14) hold also when $b$ is replaced by $h$._

Conditions (3.12) and (3.13) hold, for example, for bounded covariates, or
covariates which fulfill a local sub-Gaussianity condition. Similarly, (3.14)
is implied by a sub-gaussianity or boundedness condition on the covariates and
the residuals. We use these tail-constraints to prove a Bernstein-type
concentration result (cf. Proposition A.4 and the reference there for a
general statement and Lemma A.5 for a formulation of the statement which is
tailored to our setting). Condition (3.13) is a specific requirement for the
model selection step and is thus formulated in terms of the model selection
bandwidth $b$.

Assumption (CV): (Covariate Variance). _It holds that_

$\displaystyle\min_{n\in\mathbb{N}}\min_{k\in\\{1,...,p_{n}\\}}\mathbb{E}\left(\frac{1}{b}K\left(\frac{X_{i}}{b}\right)^{2}\left(Z_{i}^{(k)}\right)^{2}\right)>0,$
$\displaystyle\max_{n\in\mathbb{N}}\max_{k\in\\{1,...,p_{n}\\}}\mathbb{E}\left(\frac{1}{b}K\left(\frac{X_{i}}{b}\right)^{2}\left(Z_{i}^{(k)}\right)^{2}\right)<\infty.$

This assumption ensures that no single covariate has, asymptotically, either a
negligible or dominating variance, and thus that all covariates have a roughly
similar scale.

###  3.2 Main Result and Discussion

The main result of this paper shows asymptotic normality of our estimator
under the conditions stated in Section 3.1. In order to state it, we define
the following constants, which depend on the kernel function only:

$\displaystyle C_{\mathcal{B}}=$
$\displaystyle\frac{K_{+}^{(3)}-2K_{+}^{(1)}K_{+}^{(2)}}{K_{+}^{(2)}-2\left(K_{+}^{(1)}\right)^{2}},\quad
C_{\mathcal{S}}=\frac{(K^{2})^{(0)}(K_{+}^{(2)})^{2}+(K^{2})_{+}^{(2)}(K_{+}^{(1)})^{2}-2(K^{2})_{+}^{(1)}K_{+}^{(1)}}{\left[(K_{+}^{(1)}-\frac{1}{2}K_{+}^{(2)}\right]^{2}}.$

######  Theorem 3.1.

Suppose the assumptions from Section 3.1 hold, and that the penalty parameter
$\lambda$ is chosen such that $\lambda=O(\sqrt{\log p_{n}/(ng)})$. Then there
are sequences $\mathcal{B}_{n}$ and $\mathcal{S}_{n}$ satisfying

$\displaystyle\mathcal{B}_{n}=$
$\displaystyle\frac{C_{\mathcal{B}}}{2}\left(\mu_{\widetilde{Y}+}^{\prime\prime}-\mu_{\widetilde{Y}-}^{\prime\prime}\right)+o\left(|J_{n}|^{1/2}\right)\textnormal{
and}$ (3.15) $\displaystyle\mathcal{S}^{2}_{n}=$
$\displaystyle\frac{C_{\mathcal{S}}}{f_{X}(0)}\left(\sigma^{2}_{\widetilde{Y}+}+\sigma^{2}_{\widetilde{Y}-}\right)+o(1)$
(3.16)

such that (2.7) holds:

$\frac{\sqrt{nh}\left(\hat{\tau}_{h}(\hat{J}_{n})-\tau_{Y}-h^{2}\mathcal{B}_{n}\right)}{\mathcal{S}_{n}}\overset{d}{\to}\mathcal{N}\left(0,1\right).$

The theorem’s proof is given in Online Appendix A.1. The following remarks
discuss its implications and possible extensions.

######  Remark 3.3 (Asymptotic Bias).

The first term in the expansion of the asymptotic bias $\mathcal{B}_{n}$ in
(3.15) is proportional to

$\mu_{\widetilde{Y}+}^{\prime\prime}-\mu_{\widetilde{Y}-}^{\prime\prime}=(\mu_{Y+}^{\prime\prime}-\mu_{Y-}^{\prime\prime})-(\mu_{Z(J_{n})+}^{\prime\prime}-\mu_{Z(J_{n})-}^{\prime\prime})^{\top}\gamma_{n}\equiv
A-B_{n}.$

The term $A=\mu_{Y+}^{\prime\prime}-\mu_{Y-}^{\prime\prime}$ is the one we
would obtain for the baseline estimator, and the term
$B_{n}=(\mu_{Z(J_{n})+}^{\prime\prime}-\mu_{Z(J_{n})-}^{\prime\prime})^{\top}\gamma_{n}$
captures the covariates’ contribution to the bias. Depending on the curvature
of the components of $\mu_{Z(J_{n})}$, in theory our estimator’s overall
asymptotic bias could thus be larger or smaller than that of the baseline
estimator; and since $B_{n}$ could potentially be of order $O(|J_{n}|^{1/2})$
our estimator’s overall asymptotic bias could also vanish at the same or a
slightly slower rate than that of the baseline estimator.

We would argue that one should not be too concerned though that including
covariates could increase the bias in an empirical context. In particular, the
general notion of covariates being predetermined prior to treatment assignment
is plausibly compatible with a strengthening of (2.4) that assumes that not
only the levels but also the second derivatives of $\mu_{Z}$ are continuous
around the cutoff. In this case, we have $B_{n}=0$, and the leading bias in
(3.15) simplifies to that of the baseline RD estimator.

We note that there are some pathological cases in which the first term in
(3.15) would actually be of smaller order than the remainder term, and thus
not be leading. The following Lemma provides conditions that rule this out.

######  Lemma 3.4.

Suppose that, in addition to the assumptions of Theorem 3.1, there is a
constant $\eta>0$ such that

$\left|\left(\mu_{Z(J_{n})+}^{\prime\prime}-\mu_{Z(J_{n})-}^{\prime\prime}\right)^{\top}\gamma_{n}-\left(\mu_{Y+}^{\prime\prime}-\mu_{Y-}^{\prime\prime}\right)\right|\geq\eta$
(3.17)

and, if
$\|\mu_{Z(J_{n})+}^{\prime\prime}-\mu_{Z(J_{n})-}^{\prime\prime}\|_{2}\to\infty$,
there is a constant $c>0$ such that

$\|\mu_{Z(J_{n})+}^{\prime\prime}-\mu_{Z(J_{n})-}^{\prime\prime}\|_{2}\leq
c\left|\left(\mu_{Z(J_{n})+}^{\prime\prime}-\mu_{Z(J_{n})-}^{\prime\prime}\right)^{\top}\gamma_{n}\right|.$
(3.18)

Then

$\mathcal{B}_{n}=\frac{C_{\mathcal{B}}}{2}\left(\mu_{\widetilde{Y}+}^{\prime\prime}-\mu_{\widetilde{Y}-}^{\prime\prime}\right)(1+o\left(1\right)).$
(3.19)

Requirements (3.17) and (3.18) are the natural extensions of the standard
assumption that $\mu_{Y+}^{\prime\prime}-\mu_{Y-}^{\prime\prime}\neq 0$ to our
setting. They exclude pathological cases in which the covariates’ contribution
to the bias $B_{n}$ and the “no covariates” component $A$ happen to cancel
each other asymptotically, and cases in which some of the components of the
vector $\mu_{Z(J_{n})+}^{\prime\prime}-\mu_{Z(J_{n})-}^{\prime\prime}$ are
large, but happen to exactly offset each other such that their linear
combination $B_{n}$ vanishes.

######  Remark 3.5 (Asymptotic Variance).

The leading term of our estimator’s asymptotic variance $\mathcal{S}_{n}^{2}$
converges to a positive constant under our assumptions, and is guaranteed not
to exceed that of the baseline estimator from (2.2), or that of any estimator
of the form in (2.5) that uses a strict subset of the the target set $J_{n}$.
This can be seen by noting first that it depends on the covariates only
through the term

$\sigma^{2}_{\widetilde{Y}+}+\sigma^{2}_{\widetilde{Y}-}=\lim_{x\downarrow
0}\textrm{Var}(Y_{i}-Z_{i}(J_{n})^{\top}\gamma_{n}|X_{i}=x)+\lim_{x\uparrow
0}\textrm{Var}(Y_{i}-Z_{i}(J_{n})^{\top}\gamma_{n}|X_{i}=x);$

and second that $\gamma_{n}$ minimizes the function

$\gamma\mapsto\lim_{x\downarrow
0}\textrm{Var}(Y_{i}-Z_{i}(J_{n})^{\top}\gamma|X_{i}=x)+\lim_{x\uparrow
0}\textnormal{Var}(Y_{i}-Z_{i}(J_{n})^{\top}\gamma|X_{i}=x).$

This holds because the function is quadratic in $\gamma$, has a positive
definite Hessian matrix, and its Jacobian is set to zero by
$\gamma_{n}=(\lim_{x\downarrow
0}\textrm{Var}(Z_{i}(J_{n})|X_{i}=x)+\lim_{x\uparrow
0}\textnormal{Var}(Z_{i}(J_{n}|X_{i}=x))^{-1}(\lim_{x\downarrow
0}\textrm{Cov}(Y_{i},Z_{i}(J_{n})|X_{i}=x)+\lim_{x\uparrow
0}\textnormal{Cov}(Y_{i},Z_{i}(J_{n})|X_{i}=x))$, as defined above.666The
remarks in Calonico et al. (2019, Section IV.B) seem to suggest that this
result should only hold under additional conditions, but this does not appear
to be the case here.

At first glance, this result might seem to suggest that a sequence of larger
target sets must always lead to a smaller asymptotic variance. The following
Lemma shows that this is not the case. Indeed, it shows the stronger statement
that if there are two covariate sets
$J_{1,n},J_{2,n}\subseteq\\{1,...,p_{n}\\}$ that both satisfy our assumptions,
the corresponding estimators of the form in (2.5) must have the same
asymptotic variance.

######  Lemma 3.6.

Let $J_{1,n},J_{2,n}\subseteq\\{1,...,p_{n}\\}$ be two sequences of covariate
indices such that

$\left\|\mathbb{E}\left(K_{h}(X)\left(V_{i}^{\top}\quad Z_{i}(J_{1,n}\cup
J_{2,n})^{\top}\right)^{\top}\left(V_{i}^{\top}\quad Z_{i}(J_{1,n}\cup
J_{2,n})^{\top}\right)\right)^{-1}\right\|_{2}=O(1)$

and let $|J_{j,n}|\log p_{n}/nh\to 0$ for $j=1,2$. Suppose Assumptions (AS)
and (TCS) hold with $J_{n}$ replaced by either $J_{1,n}$ or $J_{2,n}$, and let
$\mathcal{S}_{j,n}^{2}$ be the asymptotic variance of
$\hat{\tau}_{h}(J_{j,n})$ as in Theorem 3.1, for $j=1,2$. Then
$\mathcal{S}_{1,n}^{2}-\mathcal{S}_{2,n}^{2}=o(1)$.

We conjecture that even in a setting in which (AS) does not hold for any
target set, our estimator continues to have the smallest asymptotic variance
among all linear adjustment estimators that only use a moderate (in some
appropriate sense) number of the available covariates, given some suitable
choice of $\lambda$. However, proving such a result would require developing a
“non-sparse” theory for the Lasso, which is beyond the scope of this paper.

######  Remark 3.7 (Double Selection).

As an alternative to our proposed estimator, one could also consider a “double
selection” procedure which, as in Belloni et al. (2013), additionally includes
covariates that are predictive for treatment status into the active set. That
is, one could redefine the set $\hat{J}_{n}$ as
$\hat{J}_{n}=\\{k\in\\{1,...,p_{n}\\}:\widetilde{\gamma}_{n}^{(k)}\neq
0\textnormal{ or }\bar{\gamma}_{n}^{(k)}\neq 0\\}$, where

$\bar{\gamma}_{n}=\underset{\gamma\in\mathbb{R}^{p_{n}}}{\operatorname{argmin}}\sum_{i=1}^{n}K_{b}(X_{i})\left(T_{i}-\left(Z_{i}-\hat{\mu}_{Z,n}\right)^{\top}\gamma\right)^{2}+\lambda\sum_{j=k}^{p_{n}}\hat{w}_{n,k}|\gamma_{k}|,$

and $\widetilde{\gamma}_{n}$ is as defined above. This change would not affect
the properties of the final RD estimator, however, as under (2.4) the
covariates are not informative about treatment status among units local to the
cutoff (a conceptually similar phenomenon appears in randomized experiments).
For analogous reasons, the large sample properties of our final RD estimator
would also be the same as that of an alternative procedure that excludes the
predictors $V_{i}$ from the model selection step (we choose to include them
because that substantially simplifies the algebra in some steps of our
proofs).

######  Remark 3.8 (Role of the Lasso).

Our assumptions do not imply that the Lasso recovers the target set with very
high probability, in the sense that $\mathbb{P}(\hat{J}_{n}=J_{n})\to 1$.
Existing results suggest that such a property could only be established under
substantially stronger conditions, including exact rather than approximate
sparsity and a so-called “$\beta$-min” condition that imposes a substantial
lower bound on the values of all non-zero coefficients (cf. Section 2.6 in van
de Geer and Bühlmann (2011) and the references therein). Consistent estimation
of the target set is also not required for our results. In fact, similarly as
in Belloni and Chernozhukov (2013), we only require that the covariates
obtained via the model selection step lead to regression residuals that are
“similar” to those corresponding to the target set. Precise requirements are
given in Assumption (CMS) in the Online Appendix, and we show in the proof of
Theorem 3.1 that these are satisfied by the Lasso. Any model selection
procedure that satisfies Assumption (CMS) could be used instead of the Lasso
in our procedure as well. The simplest example being the _ad-hoc_ selection of
a pre-defined set of covariates in which case $p_{n}\equiv p_{0}$ and
$\hat{J}_{n}\equiv\\{1,...,p_{0}\\}$.

##  4 Numerical Results

###  4.1 Implementation

In order to implement our proposed method in practice, one has to choose the
initial bandwidths $b$, the Lasso penalty $\lambda$, and the final bandwidth
$h$. As discussed in Section 2, the latter can in principle be chosen by
applying any approach deemed suitable for settings with low-dimensional
covariates to the generated data set
$\\{(Y_{i},X_{i},Z_{i}(\hat{J}_{n})),i=1,\ldots,n\\}$, such as those proposed
by Calonico et al. (2019) or Armstrong and Kolesár (2018). We conduct our
simulations using both frameworks, with results based on the methods in
Armstrong and Kolesár (2018) reported in this section, and results based on
the methods in Calonico et al. (2019) reported in Online Appendix C.1.

The choice of $b$ and $\lambda$ is complicated by the fact that these
quantities do not appear in the limiting distribution of our final RD
estimator. Our heuristic recommendation is to use a method for bandwidth
choice designed for settings without covariates, like the ones proposed in
Imbens and Kalyanaraman (2012), Calonico et al. (2014) or Armstrong and
Kolesár (2018) to select $b$; and we focus on the method proposed in Armstrong
and Kolesár (2018) in this section. We also consider choosing $\lambda$ via
adaptions of three methods for non-localized Lasso estimators to our RD
setting: standard cross-validation, the plug-in procedure of Belloni et al.
(2013), and a recently proposed bootstrap-based method by Lederer and Vogt
(2021). The three procedures are described formally in Online Appendix B, and
we refer to them by the acronyms (CV), (BCH) and (LV), respectively, below.
Our computations in this section use the R packages glmnet for implementing
Lasso-based covariate selection, and RDHonest for bandwidth selection,
standard errors, and confidence intervals.

###  4.2 Simulations

For the simulations, we consider the following DGP, which is a variation of
“Model 2” in Calonico et al. (2019) and corresponds to an RD setting with
$p=200$ covariates and parameter of interest $\tau_{Y}=0.02$:

$\displaystyle X$ $\displaystyle\sim 2\cdot\textrm{beta}(2,4)-1,\quad
T=\mathbbm{1}(X\geq
0),\quad(\varepsilon,Z^{\top})^{\top}\sim\mathcal{N}\left(\mathbf{0},\Sigma\right),\quad\Sigma=\begin{pmatrix}\sigma_{\varepsilon}^{2}&v^{\top}\\\
v&\sigma_{Z}^{2}I_{200}\end{pmatrix},$ $\displaystyle Y$
$\displaystyle=\varepsilon+\left\\{\begin{array}[]{lll}0.36+0.96\cdot
X+5.47\cdot X^{2}&&\\\ \quad\quad\quad\quad+15.28\cdot x^{3}+15.87\cdot
X^{4}+5.14\cdot X^{5}+0.22\cdot Z^{\top}\alpha,&\textrm{if }T=0,\\\
0.38+0.62\cdot X-2.84\cdot X^{2}&&\\\ \quad\quad\quad\quad+8.42\cdot
X^{3}-10.24\cdot X^{4}+4.31\cdot X^{5}+0.28\cdot Z^{\top}\alpha,&\textrm{if
}T=1,\end{array}\right.$

with $\sigma_{\varepsilon}^{2}=0.1295^{2}$, $\sigma_{Z}^{2}=0.1353^{2}$,
$I_{200}$ denoting the $200\times 200$ identity matrix, $v\in\mathbb{R}^{200}$
a vector whose $k$th component is equal to
$v_{k}=0.8\sqrt{6}\sigma_{\varepsilon}^{2}/\pi k$, and
$\alpha\in\mathbb{R}^{200}$ a vector whose $k$th component is equal to
$\alpha_{k}=2/k^{2}$. This choice of $\alpha$ implies that the $k$th covariate
$Z^{(k)}$ becomes less important for variance reduction as $k$ increases. We
consider the sample size $n=1,000$ and set the number of Monte Carlo
replications to $10,000$.

We report results for our post-Lasso procedure with the penalty parameter
$\lambda$ selected via either of (CV), (BCH) and (VL). For comparison, we also
consider linear adjustment estimators that use different fixed subsets of the
covariates, namely either no covariates, only the first covariate, only the
first 10 covariates, only the first 30 covariates, only the first 50
covariates, or only the “optimal” linear combination of covariates
$Z^{\top}\alpha$. Note that the covariates are ordered by their “importance”
in our DGP, the procedures that use a fixed non-zero number of covariates are
infeasible, as in practice the econometrician would generally not know which
covariates are the most important ones. The “optimal covariate” estimator is
not feasible as well, as the vector $\alpha$ is generally unknown in
applications. These estimators serve as (oracle) performance benchmarks in our
simulation study.

Table 4.1: Simulation Results

Covariate Selection | #Cov. | Bias | SD | Avg. SE | CI Length | Coverage
---|---|---|---|---|---|---
Lasso (CV) | 9.5 | 0.0054 | 0.0464 | 0.0325 | 0.1559 | 87.7
Lasso (BCH) | 1.2 | 0.0051 | 0.0512 | 0.0488 | 0.2189 | 96.1
Lasso (LV) | 1.9 | 0.0047 | 0.0482 | 0.0445 | 0.2005 | 95.6
Fixed: No Covariates | – | 0.0070 | 0.0744 | 0.0735 | 0.3278 | 96.7
Fixed: Most Important Covariate Only | – | 0.0050 | 0.0516 | 0.0499 | 0.2236 | 96.3
Fixed: 10 Most Important Covariates | – | 0.0044 | 0.0427 | 0.0364 | 0.1670 | 94.6
Fixed: 30 Most Important Covariates | – | 0.0055 | 0.0447 | 0.0282 | 0.1390 | 87.5
Fixed: 50 Most Important Covariates | – | 0.0063 | 0.0487 | 0.0213 | 0.1176 | 77.0
Fixed: Optimal Covariate | – | 0.0042 | 0.0441 | 0.0424 | 0.1900 | 96.2

Results based on $10,000$ Monte Carlo replications. For each estimator, the
table shows shows average number of selected covariates (#Cov.), the bias
(Bias), the standard deviation (SD), the average value of the final
estimator’s standard error (SE), the average length of the corresponding
confidence interval for the parameter of interest (CI Length), and the share
of simulation runs in which the respective confidence interval covered the
true parameter value (Coverage).

Our simulation results are summarized in Table 4.1. Regarding our proposed
procedures, we see that choosing the Lasso penalty via (CV) leads to
substantially more covariates being selected relative to (BCH) or (LV), with
the latter being roughly similar. All Lasso-based estimators have similarly
low bias and similar empirical standard deviations; and the latter are both
lower than that of the “no covariates” baseline, and close to that of the
“optimal covariate” oracle estimator. Standard errors and CI coverage are
accurate with (BCH) and (LV). However, for (CV) the standard errors notably
underestimate the true standard deviation, which leads to slight CI
undercoverage. This phenomenon seems to occur because the slightly larger
number of covariates that (CV) tends to select already leads to overfitting in
the post-Lasso stage.

A similar effect occurs with the linear adjustment estimator that uses fixed
sets of the covariates. While the performance of this estimator is good in our
simulations if only the single most important covariate is used, with the 10,
30 or 50 most important covariates we see a progressively severe downward bias
in the standard error, with corresponding CI undercoverage.

Overall, the simulation results are in line with our asymptotic theory, and
show that our procedures can obtain near-oracle performance in practice. They
also highlight the need for working with a small number of covariates to
obtain reliable inference, and thus the need for covariate selection even if
the number of available covariates is only moderate relative to the sample
size.

###  4.3 Empirical Application

In this section, we apply our methodology to data on Austrian workers from
Card et al. (2007), to whom we refer for an extensive description of its
construction. During the sample period, workers are eligible for severance
payments when losing their job if they have at least 36 months of job tenure
at the time of separation. One part of the analysis in Card et al. (2007)
concerns the question whether severance payments lead to higher wages in
future jobs (by enabling workers to search longer for a new position, and thus
find better matches). We use our method to reanalyze this question, taking
previous job tenure as the running variable, with a cutoff at 36 months, and
the difference in log wages between old and new jobs as the outcome. The data
include a large number of covariates containing information about workers’
socio-demographic characteristics and the nature of their employment. We
select 60 of these covariates, and split them into a basic and an extended set
as follows:

Basic Covariates:

Gender, marital status, Austrian nationality, “blue collar” occupation, age
and its square, $\log$ of previous wage and its square, indicators for month
and year of job termination (38 covariates)

Additional Covariates:

Work experience and its square, number of employees in firm at job just lost,
indicator of having a job before the one just lost, “blue collar” status at
job prior to the one lost, indicator of having a prior spell of nonemployment,
duration of last nonemployment, total number of spells of nonemployment in
career, indicator of being recalled to the job before the one just lost,
indicator for higher education, indicators for industry sector and region (22
covariates)

We also create further covariates by including all non-trivial interaction
terms and trigonometric series transformations of all non-dummy variables,
which are of the form $\sin(2\pi k\times\textrm{variable})$ and $\cos(2\pi
k\times\textrm{variable})$ for $k=1,...,5$. This results in a total of 1,958
covariates. After removing all observations with at least one missing
covariate value, data on 288,175 workers is available for the empirical
analysis. We then compute an estimate of the RD parameter, with associated
standard error and confidence intervals. In view of our simulation results, we
only consider (BCH) for selecting the penalty parameter. We compare the result
to those based on the baseline estimator or linear adjustment estimators that
either use only the basic set of covariates, or both the basic and the
additional set of covariates.

The results are shown in Table 4.2. All four methods produce similar point
estimates close to zero, which is in line with the results in Card et al.
(2007). Our method’s standard error is about $11\%$ lower than that of the
baseline estimator without covariates, showing that the use of a small number
of carefully selected covariates can meaningfully improve estimation accuracy.
The corresponding confidence interval is also shorter by a similar factor.
Importantly, (BCH) only selects three of the 1,958 covariates, all of which
are interactions of the squared logarithm of the previous wage with some other
variable. The results for the two remaining linear adjustment estimators show
that controlling for a larger number of covariates does not yield meaningfully
smaller standard errors, which suggest that approximate sparsity is a
reasonable assumption in this context. In view of our simulation results
above, there also is concern that using 38 or 60 covariates could lead to
downward biased standard errors for the corresponding linear adjustment
estimators.

Table 4.2: Estimation Results

|  |  |  | Confidence Interval
---|---|---|---|---
Covariate Selection | #Cov. | Estimator | SE | Lower | Upper | Length
Lasso (BCH) | 3 | 0.0387 | 0.0207 | -0.0066 | 0.0839 | 0.0905
Fixed: None | 0 | -0.0063 | 0.0233 | -0.0570 | 0.0443 | 0.1013
Fixed: Basic only | 38 | 0.0116 | 0.0200 | -0.0322 | 0.0554 | 0.0876
Fixed: Basic + Additional | 60 | 0.0168 | 0.0203 | -0.0277 | 0.0612 | 0.0890

Results based on 288,175 observations. Column #Cov. shows the number of
selected covariates. Based on our simulation results above, there are concerns
that SEs for estimators “Fixed: Basic only” and “Fixed: Basic + Additional”
could be downward biased.

##  5 Concluding Remarks

Our results on sharp RD estimation with a potentially large number of
covariates can be extended to other settings, such as fuzzy RD or regression
kink (RK) designs. In fuzzy RD designs, for instance, units are assigned to
treatment if their realization of the running variable falls above the
threshold value, but they do not necessarily comply with this assignment. The
conditional treatment probability hence jumps at the cutoff, but in contrast
to sharp RD designs it generally does not jump from zero to one. The parameter
of interest in fuzzy RD designs is

$\tau_{\textnormal{fuzzy}}=\frac{\tau_{Y}}{\tau_{T}}\equiv\frac{\mu_{Y+}-\mu_{Y-}}{\mu_{T+}-\mu_{T-}},$

which is the ratio of two sharp RD estimands. It can be estimated by running
our proposed procedure twice, once with $Y_{i}$ and once with $T_{i}$ as the
dependent variable, and taking the ratio of the two estimates. If the
assumptions from Section 3.1 also hold with $T_{i}$ replacing $Y_{i}$, and
$\tau_{T}$ is bounded away from zero, the asymptotic normality of the
resulting estimator of $\tau_{\textnormal{fuzzy}}$ simply follows from the
delta method.

## References

  * Arai et al. (2022) Arai, Y., T. Otsu, and M. H. Seo (2022): “Regression Discontinuity Design with Potentially Many Covariates,” _Working Paper_.
  * Armstrong and Kolesár (2018) Armstrong, T. B. and M. Kolesár (2018): “Optimal inference in a class of regression models,” _Econometrica_ , 86, 655–683.
  * Belloni and Chernozhukov (2013) Belloni, A. and V. Chernozhukov (2013): “Least squares after model selection in high-dimensional sparse models,” _Bernoulli_ , 19, 521–547.
  * Belloni et al. (2013) Belloni, A., V. Chernozhukov, and C. Hansen (2013): “Inference on treatment effects after selection among high-dimensional controls,” _The Review of Economic Studies_ , 81, 608–650.
  * Benjamini and Hochberg (1995) Benjamini, Y. and Y. Hochberg (1995): “Controlling the false discovery rate: a practical and powerful approach to multiple testing,” _Journal of the Royal statistical society: series B (Methodological)_ , 57, 289–300.
  * Bickel et al. (2009) Bickel, P. J., Y. Ritov, and A. B. Tsybakov (2009): “Simultaneous analysis of Lasso and Dantzig selector,” _Annals of Statistics_ , 37, 1705–1732.
  * Calonico et al. (2019) Calonico, S., M. D. Cattaneo, M. H. Farrell, and R. Titiunik (2019): “Regression discontinuity designs using covariates,” _The Review of Economics and Statistics_ , 101, 442–451.
  * Calonico et al. (2014) Calonico, S., M. D. Cattaneo, and R. Titiunik (2014): “Robust nonparametric confidence intervals for regression-discontinuity designs,” _Econometrica_ , 82, 2295–2326.
  * Card et al. (2007) Card, D., R. Chetty, and A. Weber (2007): “Cash-on-Hand and Competing Models of Intertemporal Behavior: New Evidence from the Labor Market,” _The Quarterly Journal of Economics_ , 122.
  * Fan and Gijbels (1996) Fan, J. and I. Gijbels (1996): _Local polynomial modelling and its applications_ , Chapman & Hall/CRC.
  * Frölich and Huber (2019) Frölich, M. and M. Huber (2019): “Including Covariates in the Regression Discontinuity Design,” _Journal of Business & Economic Statistics_, 37, 736–748.
  * Giné and Nickl (2016) Giné, E. and R. Nickl (2016): _Mathematical foundations of infinite-dimensional statistical models_ , Cambridge series in statistical and probabilistic mathematics, New York: Cambridge University Press.
  * Hahn et al. (2001) Hahn, J., P. Todd, and W. Van der Klaauw (2001): “Identification and Estimation of Treatment Effects with a Regression-Discontinuity Design,” _Econometrica_ , 69, 201–209.
  * Imbens and Kalyanaraman (2012) Imbens, G. and K. Kalyanaraman (2012): “Optimal bandwidth choice for the regression discontinuity estimator,” _The Review of Economic Studies_ , 79, 933–959.
  * Klenke (2008) Klenke, A. (2008): _Probability Theory_ , Springer, London.
  * Lederer and Vogt (2021) Lederer, J. and M. Vogt (2021): “Estimating the Lasso’s Effective Noise,” _Journal of Machine Learning Research_ , 22, 1–32.
  * Noack et al. (2021) Noack, C., T. Olma, and C. Rothe (2021): “Flexible Covariate Adjustments in Regression Discontinuity Designs,” _Working Paper_.
  * Su et al. (2019) Su, L., T. Ura, and Y. Zhang (2019): “Non-separable models with high-dimensional data,” _Journal of Econometrics_ , 212, 646–677.
  * Tibshirani (1996) Tibshirani, R. (1996): “Regression shrinkage and selection via the Lasso,” _Journal of the Royal Statistical Society. Series B (Methodological)_ , 58, 267–288.
  * van de Geer and Bühlmann (2011) van de Geer, S. and P. Bühlmann (2011): _Statistics for high-dimensional data_ , Heidelberg: Springer.
  * Wager et al. (2016) Wager, S., W. Du, J. Taylor, and R. J. Tibshirani (2016): “High-dimensional regression adjustments in randomized experiments,” _Proceedings of the National Academy of Sciences_ , 113, 12673–12678.

##  Appendix A Proofs

###  A.1 Notation and Overview

For subsets $J\subseteq\\{1,...,p_{n}\\}$ we denote by $|J|$ its size and for
an arbitrary vector $a\in\mathbb{R}^{p_{n}}$ we denote its restricted version
by $a_{J}$, i.e., $a_{J}^{(i)}=a^{(i)}$ for $i\in J$ and $a_{J}^{(i)}=0$ for
$i\notin J$, where $a^{(i)}$ refers to the $i$-th entry of $a$. The restricted
version of $Z_{i}$ will be denoted by $Z_{i}(J)$. Furthermore,

$\displaystyle\mathbf{Y}=$
$\displaystyle(Y_{1},...,Y_{n})^{\top},\quad\mathbf{K}_{h}=\textrm{diag}\left(h^{-1}K(X_{1}h^{-1}),...,h^{-1}K(X_{n}h^{-1})\right),$
$\displaystyle\mathbf{V}=$
$\displaystyle\begin{pmatrix}1&T_{1}&X_{1}/h&T_{1}X_{1}/h\\\
\vdots&\vdots&\vdots&\vdots\\\
1&T_{n}&X_{n}/h&T_{n}X_{n}/h\end{pmatrix},\quad\mathbf{Z}(J)=\begin{pmatrix}Z_{1}(J)^{\top}\\\
\vdots\\\ Z_{n}(J)^{\top}\end{pmatrix},$
$\displaystyle\mathbf{r}_{n}(J_{n},h)=$
$\displaystyle(r_{1}(J_{n},h),...,r_{n}(J_{n},h))^{\top},$

where $\textrm{diag}(v)$ denotes a diagonal matrix which diagonal is given by
the vector $v$. In this notation, the (constrained to $J$) minimizer in the
$\underset{}{\operatorname{argmin}}$ of (2.3) is given by

$\left(\hat{\theta}_{n}(J),\hat{\gamma}_{n}(J)\right)=\underset{\gamma_{J^{c}=0}}{\underset{(\theta,\gamma)\in\mathbb{R}^{4+p_{n}}}{\operatorname{argmin}}}\frac{1}{n}\left\|\mathbf{K}_{h}^{\frac{1}{2}}\left(\mathbf{Y}-\mathbf{V}\theta_{0}-\mathbf{Z}\gamma\right)\right\|^{2}$

and $\hat{\tau}_{n}(J)$ is the second entry of $\hat{\theta}_{n}(J)$. Here
$\mathbf{K}_{h}^{\frac{1}{2}}$ denotes a diagonal matrix whose diagonal equals
the square roots of the corresponding diagonal elements of $\mathbf{K}_{h}$.
Note that $\hat{\tau}_{n}=\hat{\tau}_{n}(\hat{J}_{n})$ and
$\hat{\theta}_{n}=\hat{\theta}_{n}(\hat{J}_{n})$. Moreover, we can find an
explicit formula for $\hat{\theta}_{n}(J)$:

$\displaystyle\hat{\theta}_{n}(J)$
$\displaystyle=\left[\left(\mathbf{V}-\mathbf{Z}(J)\hat{\gamma}_{V,n}(J)\right)^{\top}\mathbf{K}_{h}\left(\mathbf{V}-\mathbf{Z}(J)\hat{\gamma}_{V,n}(J)\right)\right]^{-1}$
$\displaystyle\quad\quad\quad\quad\times\left(\mathbf{V}-\mathbf{Z}(J)\hat{\gamma}_{V,n}(J)\right)^{\top}\mathbf{K}_{h}\left(\mathbf{Y}-\mathbf{Z}(J)\hat{\gamma}_{n}^{*}(J)\right),$
$\displaystyle\hat{\gamma}_{n}^{*}(J)$
$\displaystyle=\left(\mathbf{Z}(J)^{\top}\mathbf{K}_{h}\mathbf{Z}(J)\right)^{-1}\mathbf{Z}(J)^{\top}\mathbf{K}_{h}\mathbf{Y},\quad\hat{\gamma}_{V,n}(J)=\left(\mathbf{Z}(J)^{\top}\mathbf{K}_{h}\mathbf{Z}(J)\right)^{-1}\mathbf{Z}(J)^{\top}\mathbf{K}_{h}\mathbf{V}.$

Note that $\hat{\gamma}_{n}^{*}(J)$ is the regression coefficient in a
regression of $Y_{i}$ on $Z_{i}(J)$, and $\hat{\gamma}_{V,n}(J)$ is a matrix
of regression coefficients from a regression of each component of the vector
$V_{i}=(1,T_{i},X_{i}/h,T_{i}X_{i}/h)^{\top}$ on $Z_{i}(J)$.

In the proof of Theorem 3.1 we will begin by generalizing the results from
Calonico et al. (2019) to a setting where there is a large (possibly larger
than $n$) set of covariates available of which we choose a subset which can be
reasonably handled with the given number of observations in a data driven way.
Therefore, we have to extend the considerations from Calonico et al. (2019) to
the case where the covariates are a random subset of a large covariate set
which grows as the number of observations grows. While doing this, we do not
specify the exact model selection algorithm but we study the asymptotic
behavior of a general model selection procedure. In the proof of Theorem 3.1
we will formulate conditions that a general model selection procedure has to
fulfill and we will show that the Lasso procedure which we introduced in the
main text has these properties.

###### Proof of Theorem 3.1.

We first show how we can reduce the situation to the case where $\mu_{Z}$ is
continuously differentiable. Recall that
$\widetilde{Z}_{i}=Z_{i}-M_{n}^{\top}V_{i}$ and define
$\widetilde{\mathbf{Z}}$ like $\mathbf{Z}$ but where $Z_{i}$ is replaced by
$\widetilde{Z}_{i}$. We can directly see that

$\displaystyle\mu_{\widetilde{Z}}(x)=\mathbb{E}(\widetilde{Z}_{i}|X_{i}=x)=$
$\displaystyle\mu_{Z}(x)-\mu_{Z-}-x\mu_{Z-}^{\prime}-x\mathbbm{1}(x\geq
0)(\mu_{Z+}^{\prime}-\mu_{Z-}^{\prime})$

is differentiable with
$\mu_{\widetilde{Z}}(0)=\mu_{\widetilde{Z}}^{\prime}(0)=0$. We compare an
estimator based on $V_{i}$ and $\widetilde{Z}_{i}$,

$\left(\check{\theta}_{n},\check{\gamma}_{n}\right)=\underset{(\theta,\gamma)}{\operatorname{argmin}}\sum_{i=1}^{n}K_{h}(X_{i})\left(Y_{i}-V_{i}^{\top}\theta-\widetilde{Z}_{i}^{\top}\gamma\right)^{2},$

with our estimator (note that it is no problem that we have here $V_{i}$
rather than $V_{i}$ in (2.5) because only the second component of
$\hat{\theta}_{n}$ will be of interest later):

$\displaystyle\left(\hat{\theta}_{n},\hat{\gamma}_{n}\right)=$
$\displaystyle\underset{(\theta,\gamma)}{\operatorname{argmin}}\sum_{i=1}^{n}K_{h}(X_{i})\left(Y_{i}-V_{i}^{\top}\theta-
Z_{i}^{\top}\gamma\right)^{2}$ $\displaystyle=$
$\displaystyle\underset{(\theta,\gamma)}{\operatorname{argmin}}\sum_{i=1}^{n}K_{h}(X_{i})\left(Y_{i}-V_{i}^{\top}\left(\theta+M_{n}\gamma\right)-\widetilde{Z}_{i}^{\top}\gamma\right)^{2},$

where all argmin above are over the set of all $(\theta,\gamma)$ with
$\gamma_{\hat{J}_{n}^{c}}=0$. Since both estimators are otherwise
unconstrained, it is true that the optimal values of the two objective
functions above are identical. From this we conclude that

$\check{\gamma}_{n}=\hat{\gamma}_{n}\quad\textrm{and}\quad\check{\theta}_{n}=\hat{\theta}_{n}+M_{n}\hat{\gamma}_{n}.$

If the optima are not unique, such solutions exist which we keep for the
remainder of the proof.

In a completely analogous fashion we also compare the population quantities

$\left(\check{\theta}_{0}(J_{n},h),\check{\gamma}_{0}(J_{n},h)\right)=\underset{(\theta,\gamma)}{\operatorname{argmin}}\,\mathbb{E}\left(K_{h}(X_{i})\left(Y_{i}-V_{i}^{\top}\theta-\widetilde{Z}_{i}^{\top}\gamma\right)^{2}\right)$

and

$\left(\theta_{0}(J_{n},h),\gamma_{0}(J_{n},h)\right)=\underset{(\theta,\gamma)}{\operatorname{argmin}}\,\mathbb{E}\left(K_{h}(X_{i})\left(Y_{i}-V_{i}^{\top}\left(\theta+M_{n}\gamma\right)-\widetilde{Z}_{i}^{\top}\gamma\right)^{2}\right),$

where the argmin are over the set of all $(\theta,\gamma)$ for which
$\gamma_{J_{n}^{c}}=0$. Then, we obtain

$\check{\gamma}_{0}(J_{n},h)=\gamma_{0}(J_{n},h)\quad\textrm{and}\quad\check{\theta}_{0}(J_{n},h)=\theta_{0}(J_{n},h)+M_{n}\gamma_{0}(J_{n},h).$

In particular, we read from these formulas

$\check{\theta}_{n}^{(2)}(J_{n},h)=\hat{\theta}_{n}^{(2)}(J_{n},h)+\left(M_{n}\right)_{2\cdot}\hat{\gamma}_{n}(J_{n},h)=\hat{\theta}_{n}^{(2)}(J_{n},h).$

We formulate now the general condition on a model selection procedure. In
principle, we can distinguish two different types of such selection
procedures: ones that (like the Lasso) generate parameter estimates
$\widetilde{\theta}_{n}$ and $\widetilde{\gamma}_{n}$ and then put
$\hat{J}_{n}=\\{j:\widetilde{\gamma}^{(j)}\neq 0\\}$; and ones that generate a
set $\hat{J}_{n}$ directly. Depending on the type of the procedure at least
one of the following quantities is well defined:

$\displaystyle B_{n}=$
$\displaystyle\frac{1}{n}\left\|\mathbf{K}_{h}^{\frac{1}{2}}\left(\mathbf{Y}-\mathbf{V}\check{\theta}_{0}(J_{n},h)-\widetilde{\mathbf{Z}}\widetilde{\gamma}_{n}\right)\right\|_{2}^{2}$
$\displaystyle\quad\quad\quad-\frac{1}{n}\left\|\mathbf{K}_{h}^{\frac{1}{2}}\left(\mathbf{Y}-\mathbf{V}\check{\theta}_{0}(J_{n},h)-\widetilde{\mathbf{Z}}\check{\gamma}_{0}(J_{n},h)\right)\right\|_{2}^{2},$
(A.1) $\displaystyle C_{n}(I_{n})=$
$\displaystyle\frac{1}{n}\left\|\mathbf{K}_{h}^{\frac{1}{2}}\left(\mathbf{Y}-\mathbf{V}\check{\theta}_{0}(J_{n},h)-\widetilde{\mathbf{Z}}\left(\check{\gamma}_{0}(J_{n},h)\right)_{\hat{J}_{n}\cap
I_{n}}\right)\right\|_{2}^{2}$
$\displaystyle\quad\quad\quad-\frac{1}{n}\left\|\mathbf{K}_{h}^{\frac{1}{2}}\left(\mathbf{Y}-\mathbf{V}\check{\theta}_{0}(J_{n},h)-\widetilde{\mathbf{Z}}\check{\gamma}_{0}(J_{n},h)\right)\right\|_{2}^{2},$
(A.2)

where $I_{n}$ is an arbitrary sequence of subsets of indices. If we use a
procedure which does not output $\widetilde{\gamma}_{n}$, we just put
$B_{n}=\infty$. Both of the quantities $B_{n}$ and $C_{n}(I_{n})$ can function
as performance measures for model selection. Let $\delta_{n}$ be such that
(recall the definition of $r_{i}(J_{n},h)$ from (2.6)):

$\sup_{k=1,...,p_{n}}\frac{1}{nh}\sum_{i=1}^{n}K\left(\frac{X_{i}}{h}\right)\widetilde{Z}_{i}^{(k)}r_{i}(J_{n},h)=O_{P}(\delta_{n}(h)).$
(A.3)

The assumption (CMS) below formulates precisely what is required of the model
selection $\hat{J}_{n}$ in order to obtain asymptotic normality of the final
RD estimator.

Assumption (CMS): (Conditions on Model Selection)
_$\widetilde{\mathbf{Z}}(\hat{J}_{n})\mathbf{K}_{h}\widetilde{\mathbf{Z}}(\hat{J}_{n})$
is almost surely invertible and RSE$(|\hat{J}_{n}|,J_{n},h)$ holds for
$\widetilde{Z}_{i}$. Denote
$\zeta_{n}=\sqrt{\max(0,\min(B_{n},C_{n}(I_{n})))}$ (with $Z_{i}$ replaced by
$\widetilde{Z}_{i}$ in the definition) for some subset
$I_{n}\subseteq\\{1,...,p_{n}\\}$ and suppose that_

$\displaystyle\beta_{1,n}$ $\displaystyle=\frac{|\hat{J}_{n}|\log
p_{n}}{nh}=o_{P}(1),\quad\beta_{2,n}=|\hat{J}_{n}|h^{4}=o_{P}(1),$
$\displaystyle\beta_{3,n}$ $\displaystyle=\left(\sqrt{\log
p_{n}}+\sqrt{nh^{5}}\right)\left(|J_{n}|+|\hat{J}_{n}|\right)^{\frac{1}{2}}\zeta_{n}=o_{P}(1),$
$\displaystyle\alpha_{n}$ $\displaystyle=\left(\sqrt{\log
p_{n}}+\sqrt{nh^{5}}\right)\left(|J_{n}|+|\hat{J}_{n}|\right)\left(\frac{1}{\sqrt{nh}}+\delta_{n}\right)=o_{P}(1),$
$\displaystyle\sqrt{|\hat{J}_{n}|}\left(\sqrt{\beta_{1,n}}+\sqrt{\beta_{2,n}}\right)\left(\alpha_{n}+\beta_{3,n}\right)=o_{P}(1)$

Suppose for the moment that (CMS) holds. Our notation from above yields (where
$\mathcal{S}_{n}^{2}$ is defined in Theorem A.1):

$\displaystyle\sqrt{\frac{nh}{\mathcal{S}_{n}^{2}}}\left(\hat{\theta}_{n}^{(2)}-\check{\theta}_{0}^{(2)}(J_{n},h)\right)=$
$\displaystyle\sqrt{\frac{nh}{\mathcal{S}_{n}^{2}}}\left(\check{\theta}_{n}^{(2)}-\check{\theta}_{0}^{(2)}(J_{n},h)\right).$
(A.4)

Hence, the asymptotics of the left hand side are determined by the asymptotics
of the right hand side. Note that

$\check{r}_{i}(J_{n},h)=Y_{i}-V_{i}^{\top}\check{\theta}_{n}-\widetilde{Z}_{i}^{\top}\check{\gamma}_{n}=Y_{i}-V_{i}^{\top}\hat{\theta}_{n}-Z_{i}^{\top}\hat{\gamma}_{n}=r_{i}(J_{n},h).$

Hence, all assumptions which we make on $r_{i}(J_{n},h)$ carry over to
$\check{r}_{i}(J_{n},h)$. Moreover,
$\mu_{Z}^{\prime\prime}=\mu_{\widetilde{Z}}^{\prime\prime}$. Thus, all
assumptions of Theorem A.1 hold true and we may thus apply this theorem with
$Z_{i}=\widetilde{Z}_{i}$ to show that

$\sqrt{\frac{nh}{\mathcal{S}_{n}^{2}}}\left(\check{\theta}_{n}^{(2)}-\check{\theta}_{0}(J_{n},h)^{(2)}\right)\to\mathcal{N}(0,1).$

We also see that all assumptions of Lemma A.10 hold for
$Z_{i}=\widetilde{Z}_{i}$ and can hence invoke this lemma to show that
$\check{\theta}_{0}(J_{n},h)^{(2)}=\tau+h^{2}\mathcal{B}_{n}$. This completes
the asymptotic normality. We show now that (3.15), (3.16) and (CMS) hold.

Discussion of Bias: Define for ease of notation

$\check{\beta}_{n}=\mathbb{E}\left(K_{h}(X_{i})\widetilde{Z}_{i}(J_{n})\widetilde{Z}_{i}(J_{n})^{\top}\right)^{-1}\mathbb{E}\left(K_{h}(X_{i})\widetilde{Z}_{i}(J_{n})Y_{i}\right),$

and use the formula for $\mathcal{B}_{n}$ from Lemma A.10 below. In order to
prove (3.15), we have to show that

$\displaystyle\frac{C_{\mathcal{B}}}{2}\left(\mu_{Y+}^{\prime\prime}-\mu_{Y-}^{\prime\prime}-\sum_{k\in
J_{n}}\left(\mu_{Z^{(k)}+}^{\prime\prime}-\mu_{Z^{(k)}-}^{\prime\prime}\right)\gamma_{n}^{(k)}\right)+o\left(|J_{n}|^{1/2}\right)$
$\displaystyle\quad=\frac{C_{\mathcal{B}}}{2}\left(\mu_{Y+}^{\prime\prime}-\mu_{Y-}^{\prime\prime}-\sum_{k\in
J_{n}}\left(\mu_{Z^{(k)}+}^{\prime\prime}-\mu_{Z^{(k)}-}^{\prime\prime}\right)\check{\beta}_{n}^{(k)}\right)+o(1)+O(|J_{n}|h^{2})+O(|J_{n}|^{\frac{1}{2}}h)$
$\displaystyle\Leftrightarrow$ $\displaystyle\sum_{k\in
J_{n}}\left(\mu_{Z^{(k)}+}^{\prime\prime}-\mu_{Z^{(k)}-}^{\prime\prime}\right)\left(\check{\beta}_{n}^{(k)}-\gamma_{n}^{(k)}\right)=o(1)+O(|J_{n}|h^{2})+O(|J_{n}|^{\frac{1}{2}}h)+o\left(|J_{n}|^{1/2}\right).$

Since we assume $|J_{n}|h^{2}\to 0$ in (BW), the dominating term on the right
hand side is $o(|J_{n}|^{1/2})$. By the Cauchy-Schwarz Inequality we have

$\left|\sum_{k\in
J_{n}}\left(\mu_{Z^{(k)}+}^{\prime\prime}-\mu_{Z^{(k)}-}^{\prime\prime}\right)\left(\check{\beta}_{n}^{(k)}-\gamma_{n}^{(k)}\right)\right|\leq\left(\sum_{k\in
J_{n}}\left(\mu_{Z^{(k)}+}^{\prime\prime}-\mu_{Z^{(k)}-}^{\prime\prime}\right)^{2}\right)^{\frac{1}{2}}\left\|\check{\beta}_{n}-\gamma_{n}\right\|_{2}.$

From the boundedness assumptions in (D) we conclude that

$\left(\sum_{k\in
J_{n}}\left(\mu_{Z^{(k)}+}^{\prime\prime}-\mu_{Z^{(k)}-}^{\prime\prime}\right)^{2}\right)^{\frac{1}{2}}=O\left(|J_{n}|^{1/2}\right)$

and hence we have left to prove that
$\left\|\check{\beta}_{n}-\gamma_{n}\right\|_{2}\to 0$. We have

$\displaystyle\left\|\check{\beta}_{n}-\gamma_{n}\right\|_{2}=$
$\displaystyle\left\|\mathbb{E}\left(K_{h}(X_{i})\widetilde{Z}_{i}(J_{n})\widetilde{Z}_{i}(J_{n})^{\top}\right)^{-1}\mathbb{E}\left(K_{h}(X_{i})\widetilde{Z}_{i}(J_{n})\left(Y_{i}-\widetilde{Z}_{i}(J_{n})^{\top}\gamma_{n}\right)\right)\right\|_{2}$
$\displaystyle=$ $\displaystyle
O(1)\left\|\mathbb{E}\left(K_{h}(X_{i})\widetilde{Z}_{i}(J_{n})\left(Y_{i}-\widetilde{Z}_{i}(J_{n})^{\top}\gamma_{n}\right)\right)\right\|_{2}$

by assumption. Define

$\mu(x)=\mathbb{E}\left((\widetilde{Z}_{i}(J_{n})\left(Y_{i}-\widetilde{Z}_{i}(J_{n})^{\top}\gamma_{n}\right)\big{|}X_{i}=x\right).$

It holds that

$\mu_{+}=\sigma_{Z(J_{n})Y+}^{2}-\sigma_{Z(J_{n})+}^{2}\gamma_{n}\quad\textrm{and}\quad\mu_{+}=\sigma_{Z(J_{n})Y+}^{2}-\sigma_{Z(J_{n})-}^{2},$

which in turn implies $\mu_{+}+\mu_{-}=0$. Using an argument as in (A.7) and
the bounded derivatives we obtain that each entry of $\mu(x)$ can be bounded
in absolute value by $Ch$, where $C$ is a suitable constant. This yields

$\left\|\mathbb{E}\left(K_{h}(X_{i})\widetilde{Z}_{i}(J_{n})\left(Y_{i}-\widetilde{Z}_{i}(J_{n})^{\top}\gamma_{n}\right)\right)\right\|_{2}^{2}\leq
C^{2}|J_{n}|h^{2}\to 0$

by the assumptions on $h$ and $|J_{n}|$. This completes the proof of (3.15).

Discussion of Variance: Recall that $\sigma_{l}^{2},\sigma_{r}^{2}$ denote the
left and right limits of $\mathbb{E}(r_{i}(J_{n},h)^{2}|X_{i}=x)$ at $x=0$,
respectively, in (TCS). It can be computed that for $h\to 0$

$\mathcal{S}_{n}^{2}\to\frac{C_{\mathcal{S}}}{f_{X}(0)}\left(\sigma_{l}^{2}+\sigma_{r}^{2}\right).$

Thus, in order to show that (3.16) holds, we have to prove that
$\sigma_{\widetilde{Y}+}^{2}\to\sigma_{r}^{2}$ and
$\sigma_{\widetilde{Y}-}^{2}\to\sigma_{l}^{2}$. In the beginning of the proof
of Lemma A.10, cf. (A.40), we provide a formula for
$\check{\theta}_{0}(J_{n},h)$ (note that each $Z_{i}$ in Lemma A.10 has to be
replaced by $\widetilde{Z}_{i}$). By combining this with (A.9), Lemma A.12 and
the boundedness assumptions on the derivatives we obtain (recall that
$\mu_{\widetilde{Z}(J_{n})}$ is continuous with
$\mu_{\widetilde{Z}(J_{n})}(0)=0$)

$\check{\theta}_{0}(J_{n},h)\to\begin{pmatrix}\mu_{Y-}&\tau&0&0\end{pmatrix}^{\top}.$

We can obtain a formula for $\gamma_{0}(J_{n},h)=\check{\gamma}_{0}(J_{n},h)$
in the fashion of (A.40) simply by interchanging the roles of $V_{i}$ and
$\widetilde{Z}_{i}$. Doing this, we obtain
$\gamma_{0}(J_{n},h)=\gamma_{n}+O(h)$ uniformly. Using this and the relation
between $\theta_{0}(J_{n},h)$ and $\check{\theta}_{0}(J_{n},h)$ we obtain
(recall that $|J_{n}|h\to 0$):

$\displaystyle\theta_{0}(J_{n},h)-\begin{pmatrix}\mu_{\widetilde{Y}-}&\mu_{\widetilde{Y}+}-\mu_{\widetilde{Y}-}&0&0\end{pmatrix}^{\top}$
$\displaystyle=$
$\displaystyle\check{\theta}_{0}(J_{n},h)-M_{n}\gamma_{0}(J_{n},h)-\begin{pmatrix}\mu_{Y-}-\mu_{Z-}^{\top}\gamma_{n}&\tau&0&0\end{pmatrix}^{\top}\to
0.$

Using this convergence we obtain

$\displaystyle\mathbb{E}\left(r(J_{n},h)^{2}|X_{i}=0+\right)-\textrm{Var}(\widetilde{Y}_{i}|X_{i}=0+)$
$\displaystyle=$
$\displaystyle\theta_{0}(J_{n},h)^{\top}\mathbb{E}(V_{i}V_{i}^{\top}|X_{i}=0+)\theta_{0}(J_{n},h)-\left(\mu_{Y+}-\mu_{Z(J_{n})+}^{\top}\gamma_{n}\right)^{2}$
$\displaystyle+\gamma_{0}(J_{n},h)^{\top}\mu_{Z(J_{n})Z(J_{n})+}\gamma_{0}(J_{n},h)-\gamma_{n}^{\top}\mu_{Z(J_{n})Z(J_{n})+}\gamma_{n}$
$\displaystyle-2\mathbb{E}(Y_{i}V_{i}^{\top}|X_{i}=0+)\theta_{0}(J_{n},h)+2\mu_{Y+}^{2}-2\mu_{Z(J_{n})+}^{\top}\mu_{Y+}\gamma_{n}$
$\displaystyle-2\mathbb{E}(Y_{i}Z_{i}(J_{n})^{\top}|X_{i}=0+)\gamma_{0}(J_{n},h)+2\mu_{Y(J_{n})Z+}^{\top}\gamma_{n}$
$\displaystyle+2\theta_{0}(J_{n},h)^{\top}\mathbb{E}(V_{i}Z_{i}(J_{n})^{\top}|X_{i}=0+)\gamma_{0}(J_{n},h)+2\gamma_{n}^{\top}\mu_{Z(J_{n})+}\mu_{Z(J_{n})+}^{\top}\gamma_{n}-2\mu_{Y+}\mu_{Z(J_{n})+}^{\top}\gamma_{n}$
$\displaystyle\to$ $\displaystyle 0.$

Thus, by definition of $\sigma_{r}^{2}$ in (TCS) we conclude that
$\textrm{Var}(\widetilde{Y}_{i}|X_{i}=0+)\to\sigma_{r}^{2}$. We can prove in a
similar fashion that
$\textrm{Var}(\widetilde{Y}_{i}|X_{i}=0-)\to\sigma_{l}^{2}$ which concludes
the proof of (3.16).

Discussion of (CMS): We show finally that the Lasso model selection fulfils
the Assumption (CMS). We assume that
$\widetilde{\mathbf{Z}}(\hat{J}_{n})\mathbf{K}_{h}\widetilde{\mathbf{Z}}(\hat{J}_{n})$
is almost surely invertible. By Theorem A.2 we have that
$|\hat{J}_{n}|=O_{P}(|J_{n}|)$ which means in particular
$\mathbb{P}(|\hat{J}_{n}|\leq\log n\cdot|J_{n}|)\to 1$. From this and
$\textrm{RSE}(|J_{n}|\log n,J_{n},h)$ for $\widetilde{Z}_{i}$, we conclude
that $\textrm{RSE}(|\hat{J}_{n}|,J_{n},h)$ holds as well. From
$|\hat{J}_{n}|=O(|J_{n}|)$ and $\delta_{n}=O(\sqrt{\log p_{n}/nh})$ which we
get from Lemma A.1 for $h$ we conclude that
$\alpha_{n},\beta_{1,n},\beta_{2,n}\to 0$ by the assumptions on the rates. By
employing Theorem A.3 we find that also $\beta_{3,n}\to 0$ and the last
condition of (CMS) holds too. ∎

In the following we discuss the use of the Lasso as model selector in Theorem
3.1, that is, we need to show that (CMS) is true. To this end, we use many
standard arguments for the Lasso and the post Lasso. See, for example, van de
Geer and Bühlmann (2011); Belloni and Chernozhukov (2013); Bickel et al.
(2009). Note that in the following our interest will lie solely in $\gamma$.
Therefore, similarly to the discussion in the proof of Theorem 3.1, we note
that we may shift $Z_{i}$ by $\alpha^{\prime}V_{i}$ for any $\alpha$ without
changing the value of $\widetilde{\gamma}_{n}$ (this is because $\theta$ is
not penalized). In contrast to Theorem 3.1 we shift for the theoretical
analysis in a way such that $\mathbb{E}(K_{b}(X_{i})Z_{i})=0$. This makes some
notation simpler. However, for computational stability it might be useful to
centralize the covariates by their empirical mean. This practice is hence not
problematic. In order to understand well what the Lasso estimator is doing we
recall in particular the notation in (2.6). Note moreover that the bandwidth
$b$ is different from $h$. Recall the definition of $\delta_{n}(h)$ in (A.3).
Under approximate sparsity, we can prove an exact rate.

######  Lemma A.1.

Let (CTB, 3.14) and (AS) hold. Let $p_{n}\to\infty$, $b\to 0$ and $\log
p_{n}/nb\to 0$. Then, for $C>0$ large enough

$\displaystyle\mathbb{P}\left(\sup_{k=1,...,p_{n}}\left|\frac{1}{nb}\sum_{i=1}^{n}Z_{i}^{(k)}K\left(\frac{X_{i}}{b}\right)r_{i}(J_{n},b)\right|>C\sqrt{\frac{\log
p_{n}}{nb}}\right)\to 0.$

Particularly, if in addition, (TCS (3.2), (3.3)) and (D conditions on
$\mu_{Z}$ and $\mu_{Z}^{\prime}$) hold for $h$ we have

$\delta_{n}(b)=O\left(\sqrt{\frac{\log p_{n}}{nb}}\right).$ (A.5)

If, moreover, all of (CTB) and $\mathbb{E}(K_{b}(X_{i})Z_{i})=0$ hold (but
possibly not (TCS)), then

$\displaystyle\mathbb{P}\left(\sup_{k=1,...,p_{n}}\left|\frac{1}{nb}\sum_{i=1}^{n}\hat{\omega}_{n,k}^{-1}Z_{i}^{(k)}K\left(\frac{X_{i}}{b}\right)r_{i}(J_{n},b)\right|>C\sqrt{\frac{\log
p_{n}}{nb}}\right)\to 0.$

###### Proof.

We only show the proof for the case with weights. Otherwise, put
$\hat{\omega}_{n,k}=1$ below and carry out the steps analogously. Since the
conditions of Lemma A.16 hold we have that $\mathbb{P}(\mathcal{A}_{n}^{c})\to
0$ for

$\mathcal{A}_{n}=\left\\{\forall k\in\\{1,...,p_{n}\\}:\hat{\omega}_{n,k}\geq
w^{(l)}\right\\}.$

Thus, we get for any $C>0$

$\displaystyle\mathbb{P}\left(\sup_{k=1,...,p_{n}}\left|\frac{1}{nb}\sum_{i=1}^{n}\hat{\omega}_{n,k}^{-1}Z_{i}^{(k)}K\left(\frac{X_{i}}{b}\right)r_{i}(J_{n},b)\right|>C\sqrt{\frac{\log
p_{n}}{nb}}\right)$ $\displaystyle=$ $\displaystyle\mathbb{P}\left(\exists
k\in\\{1,...,p_{n}\\}:\left|\frac{1}{nb}\sum_{i=1}^{n}Z_{i}^{(k)}K\left(\frac{X_{i}}{b}\right)r_{i}(J_{n},b)\right|>C\sqrt{\frac{\log
p_{n}}{nb}}\hat{\omega}_{n,k}\right)$ $\displaystyle\leq$
$\displaystyle\mathbb{P}\left(\exists
k\in\\{1,...,p_{n}\\}:\left|\frac{1}{nb}\sum_{i=1}^{n}Z_{i}^{(k)}K\left(\frac{X_{i}}{b}\right)r_{i}(J_{n},b)\right|>C\sqrt{\frac{\log
p_{n}}{nb}}w^{(l)}\right)+\mathbb{P}(\mathcal{A}_{n}^{c})$ $\displaystyle\leq$
$\displaystyle
p_{n}\max_{k=1,...,p_{n}}\mathbb{P}\left(\left|\frac{1}{nb}\sum_{i=1}^{n}Z_{i}^{(k)}K\left(\frac{X_{i}}{b}\right)r_{i}(J_{n},b)\right|>C\sqrt{\frac{\log
p_{n}}{nb}}w^{(l)}\right)+\mathbb{P}(\mathcal{A}_{n}^{c})$ $\displaystyle\leq$
$\displaystyle
p_{n}\max_{k=1,...,p_{n}}\mathbb{P}\Bigg{(}\left|\frac{1}{nb}\sum_{i=1}^{n}Z_{i}^{(k)}K\left(\frac{X_{i}}{b}\right)r_{i}(J_{n},b)-\mathbb{E}\left(Z_{i}^{(k)}\frac{1}{b}K\left(\frac{X_{i}}{b}\right)r_{i}(J_{n},b)\right)\right|$
$\displaystyle\quad\quad\quad\quad\quad>C\sqrt{\frac{\log
p_{n}}{nb}}w^{(l)}\Bigg{)}$ (A.6)
$\displaystyle+p_{n}\max_{k=1,...,p_{n}}\mathbb{P}\left(\left|\mathbb{E}\left(Z_{i}^{(k)}\frac{1}{b}K\left(\frac{X_{i}}{b}\right)r_{i}(J_{n},b)\right)\right|>C\sqrt{\frac{\log
p_{n}}{nb}}w^{(l)}\right)+\mathbb{P}(\mathcal{A}_{n}^{c}).$

We have just argued that $\mathbb{P}(\mathcal{A}_{n}^{c})\to 0$ and we assume
that (cf. (AS)) the expectations are smaller than $C\sqrt{\frac{\log
p_{n}}{nb}}$ for a large enough choice of $C>0$. Thus, we get that the second
line converges to zero. The first line converges to zero by Lemma A.5 which we
may apply because we assume the moment conditions in (CTB, (3.14)).

Note that when proving the statement without weights, we do not need to rely
on Lemma A.16 and hence we do not need to require
$\mathbb{E}(K_{b}(X_{i})Z_{i})=0$.

In order to see what $\delta_{n}(b)$ is, we note that

$\frac{1}{nh}\sum_{i=1}^{n}K\left(\frac{X_{i}}{h}\right)\widetilde{Z}_{i}^{(k)}r_{i}(J_{n},h)=\frac{1}{nh}\sum_{i=1}^{n}K\left(\frac{X_{i}}{h}\right)Z_{i}^{(k)}r_{i}(J_{n},h)-\left[M_{n}\right]_{\cdot
k}^{\top}\frac{1}{nh}\sum_{i=1}^{n}V_{i}r_{i}(J_{n},h).$

We have just shown that the first part is of order $\sqrt{\log p_{n}/nb}$
after taking the $\sup$. The second part will be shown in Lemma A.8 to be of
order $1/\sqrt{nb}$ because $M_{n}$ remains bounded by assumption. ∎

We finish this section with the proofs of Lemmas 3.4 and 3.6 which we stated
in the main text.

###### Proof of Lemma 3.4.

Let $a_{n}$ be such that the following is true ($\check{\beta}_{n}$ was
defined in the proof of Theorem 3.1 and the big-O and small-o terms are the
same as in the definition of $\mathcal{B}_{n}$ in Lemma A.10):

$\displaystyle\frac{C_{\mathcal{B}}}{2}\left(\mu_{\widetilde{Y}+}^{\prime\prime}-\mu_{\widetilde{Y}-}^{\prime\prime}\right)(1+a_{n})=\frac{C_{\mathcal{B}}}{2}\Bigg{[}\mu_{Y+}^{\prime\prime}-\mu_{Y-}^{\prime\prime}-\sum_{k\in
J_{n}}\left(\mu_{Z^{(k)}+}^{\prime\prime}-\mu_{Z^{(k)}-}^{\prime\prime}\right)\check{\beta}_{n}^{(k)}\Bigg{]}$
$\displaystyle\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad+o(1)+O(|J_{n}|h^{2})+O(|J_{n}|^{\frac{1}{2}}h)$
$\displaystyle\Leftrightarrow$
$\displaystyle\left(\mu_{\widetilde{Y}+}^{\prime\prime}-\mu_{\widetilde{Y}-}^{\prime\prime}\right)a_{n}=\Bigg{[}-\sum_{k\in
J_{n}}\left(\mu_{Z^{(k)}+}^{\prime\prime}-\mu_{Z^{(k)}-}^{\prime\prime}\right)\left(\check{\beta}_{n}^{(k)}-\gamma_{n}^{(k)}\right)\Bigg{]}$
$\displaystyle\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad+o(1)+O(|J_{n}|h^{2})+O(|J_{n}|^{\frac{1}{2}}h)$

In order to prove (3.19), we need to show that $a_{n}\to 0$. Note next, that
we assume $|J_{n}|h^{2}\to 0$ in (BW) and hence all big-O-terms above are
$o(1)$. Since (3.17) holds by assumption, we have left to show that

$\displaystyle\left|\frac{\sum_{k\in
J_{n}}\left(\mu_{Z^{(k)}+}^{\prime\prime}-\mu_{Z^{(k)}-}^{\prime\prime}\right)\left(\check{\beta}_{n}^{(k)}-\gamma_{n}^{(k)}\right)}{\mu_{\widetilde{Y}+}^{\prime\prime}-\mu_{\widetilde{Y}-}^{\prime\prime}}\right|$
$\displaystyle\leq$ $\displaystyle\left|\frac{\sum_{k\in
J_{n}}\left(\mu_{Z^{(k)}+}^{\prime\prime}-\mu_{Z^{(k)}-}^{\prime\prime}\right)^{2}}{\left(\mu_{\widetilde{Y}+}^{\prime\prime}-\mu_{\widetilde{Y}-}^{\prime\prime}\right)^{2}}\right|^{\frac{1}{2}}\left\|\check{\beta}_{n}-\gamma_{n}\right\|_{2}\to
0$

(use the Cauchy-Schwarz Inequality). The first part is $O(1)$ by (3.17) and
(3.18) and $\|\check{\beta}_{n}-\gamma_{n}\|_{2}\to 0$ converges to zero which
was proven in the proof of Theorem 3.1. ∎

###### Proof of Lemma 3.6.

Suppose firstly that $J_{1,n}\supseteq J_{2,n}$. By definition (cf. A.1)

$\mathcal{S}_{n}^{2}(J_{1,n})=\frac{1}{h}\mathbb{E}\left(K\left(\frac{X_{i}}{h}\right)^{2}\xi\left(\frac{X_{i}}{h}\right)r_{i}(J_{1,n},h)^{2}\right),$

where $\xi$ is a given function. Hence, we need to study the limit of the
following difference (below $\gamma_{0}(J_{2,n},h)$ as a vector in
$\mathbb{R}^{|J_{1,n}|}$ which has zeros at the places $J_{1,n}\setminus
J_{2,n}$)

$\displaystyle\frac{1}{h}\mathbb{E}\left(K\left(\frac{X_{i}}{h}\right)^{2}\xi\left(\frac{X_{i}}{h}\right)r_{i}(J_{1,n},h)^{2}\right)-\frac{1}{h}\mathbb{E}\left(K\left(\frac{X_{i}}{h}\right)^{2}\xi\left(\frac{X_{i}}{h}\right)r_{i}(J_{2,n},h)^{2}\right)$
$\displaystyle=$
$\displaystyle\frac{1}{h}\mathbb{E}\left(K\left(\frac{X_{i}}{h}\right)^{2}\xi\left(\frac{X_{i}}{h}\right)\left(\begin{pmatrix}\theta_{0}(J_{2,n},h)-\theta_{0}(J_{1,n},h)\\\
\gamma_{0}(J_{2,n},h)-\gamma_{0}(J_{1,n},h)\end{pmatrix}^{\top}\begin{pmatrix}V_{i}\\\
Z_{i}(J_{1,n})\end{pmatrix}\right)^{2}\right)$
$\displaystyle+\frac{2}{h}\mathbb{E}\left(K\left(\frac{X_{i}}{h}\right)^{2}\xi\left(\frac{X_{i}}{h}\right)r_{i}(J_{2,n},h)\begin{pmatrix}\theta_{0}(J_{2,n},h)-\theta_{0}(J_{1,n},h)\\\
\gamma_{0}(J_{2,n},h)-\gamma_{0}(J_{1,n},h)\end{pmatrix}^{\top}\begin{pmatrix}V_{i}\\\
Z_{i}(J_{1,n})\end{pmatrix}\right).$

We suppose in (TCS) that $\mathbb{E}(r_{i}(J_{a,n},h)^{2}|X_{i}=x)$ behaves
nicely around $x=0$ for $a=1,2$ and hence the qualitative behavior of the
above is determined by
$\left\|\theta_{0}(J_{1,n},h)-\theta_{0}(J_{2,n},h)\right\|_{2}^{2}$ and
$\left\|\gamma_{0}(J_{1,n},h)-\gamma_{0}(J_{2,n},h)\right\|_{2}^{2}$. These
can be controlled by applying least squares algebra as follows:

$\displaystyle\begin{pmatrix}\theta_{0}(J_{1,n},h)\\\
\gamma_{0}(J_{1,n},h)\end{pmatrix}-\begin{pmatrix}\theta_{0}(J_{2,n},h)\\\
\gamma_{0}(J_{2,n},h)\end{pmatrix}$ $\displaystyle=$
$\displaystyle\mathbb{E}\left(K_{h}(X_{i})\begin{pmatrix}V_{i}\\\
Z_{i}(J_{1,n})\end{pmatrix}\begin{pmatrix}V_{i}\\\
Z_{i}(J_{1,n})\end{pmatrix}^{\top}\right)^{-1}$
$\displaystyle\quad\quad\quad\quad\times\mathbb{E}\left(K_{h}(X_{i})\begin{pmatrix}V_{i}\\\
Z_{i}(J_{1,n})\end{pmatrix}\left(Y_{i}-\begin{pmatrix}V_{i}\\\
Z_{i}(J_{1,n})\end{pmatrix}^{\top}\begin{pmatrix}\theta_{0}(J_{2,n},h)\\\
\gamma_{0}(J_{2,n},h)\end{pmatrix}\right)\right)$ $\displaystyle=$
$\displaystyle\mathbb{E}\left(K_{h}(X_{i})\begin{pmatrix}V_{i}\\\
Z_{i}(J_{1,n})\end{pmatrix}\begin{pmatrix}V_{i}\\\
Z_{i}(J_{1,n})\end{pmatrix}^{\top}\right)^{-1}\mathbb{E}\left(K_{h}(X_{i})\begin{pmatrix}V_{i}\\\
Z_{i}(J_{1,n})\end{pmatrix}r_{i}(J_{n,2},h)\right).$

Since

$\left\|\mathbb{E}\left(K_{h}(X_{i})\begin{pmatrix}V_{i}\\\
Z_{i}(J_{1,n})\end{pmatrix}\begin{pmatrix}V_{i}\\\
Z_{i}(J_{1,n})\end{pmatrix}^{\top}\right)^{-1}\right\|_{2}=O(1)$

and by Assumption (AS) and noting that
$\mathbb{E}(K_{h}(X_{i})V_{i}r_{i}(J_{2,n},h))=0$ we obtain

$\displaystyle\left\|\begin{pmatrix}\theta_{0}(J_{1,n},h)\\\
\gamma_{0}(J_{1,n},h)\end{pmatrix}-\begin{pmatrix}\theta_{0}(J_{2,n},h)\\\
\gamma_{0}(J_{2,n},h)\end{pmatrix}\right\|_{2}^{2}=O(1)\left\|\mathbb{E}\left(K_{h}(X_{i})\begin{pmatrix}V_{i}\\\
Z_{i}(J_{1,n})\end{pmatrix}r_{i}(J_{n,2},h)\right)\right\|_{2}^{2}$
$\displaystyle=$ $\displaystyle O\left(\frac{\log p_{n}}{nh}|J_{1,n}|\right)$

which converges to zero. If $J_{2,n}\supseteq J_{1,n}$ we can apply the same
arguments with the roles of $J_{1,n}$ and $J_{2,n}$ interchanged. If $J_{1,n}$
and $J_{2,n}$ are not nested, define $J_{3,n}=J_{1,n}\cup J_{2,n}$. (AS) and
(TCS) continue to hold for $J_{3,n}$ and by the above argument,
$\mathcal{S}_{n}^{2}(J_{3,n})-\mathcal{S}_{n}^{2}(J_{1,n})\to 0$ and
$\mathcal{S}_{n}^{2}(J_{3,n})-\mathcal{S}_{n}^{2}(J_{2,n})\to 0$ which implies
$\mathcal{S}_{n}^{2}(J_{1,n})-\mathcal{S}_{n}^{2}(J_{2,n})\to 0$ and the proof
is complete. ∎

###  A.2 Preliminary Results

Let $L:\mathbb{R}\to\mathbb{R}$ be an arbitrary function and recall the
definitions

$L_{-}^{(\alpha)}=\int_{-\infty}^{0}L(u)u^{\alpha}du,\quad
L_{+}^{(\alpha)}=\int_{0}^{\infty}L(u)u^{\alpha}du,\quad
L^{(\alpha)}=\int_{-\infty}^{\infty}L(u)u^{\alpha}du$

for $\alpha\in\\{0,1,2,3,4\\}$. If $L$ is a symmetric second order kernel, we
have $L^{(0)}=1$, $L_{-}^{(0)}=L_{+}^{(0)}=\frac{1}{2}$,
$L_{-}^{(1)}=-L_{+}^{(1)}$ and $L_{-}^{(2)}=L_{+}^{(2)}$. This proves the
following lemma.

######  Lemma A.2.

For every symmetric kernel $K$ with $K^{(2)}<\infty$ the matrix

$\kappa(K)=\begin{pmatrix}K^{(0)}&K_{+}^{(0)}&K^{(1)}&K_{+}^{(1)}\\\
K_{+}^{(0)}&K_{+}^{(0)}&K_{+}^{(1)}&K_{+}^{(1)}\\\
K^{(1)}&K_{+}^{(1)}&K^{(2)}&K_{+}^{(2)}\\\
K_{+}^{(1)}&K_{+}^{(1)}&K_{+}^{(2)}&K_{+}^{(2)}\end{pmatrix}$

is invertible.

###### Proof.

By using the above relations, we find that the determinant of $\kappa(K)$ is
given by

$\left(\left(K_{+}^{(1)}\right)^{2}-\frac{1}{2}K^{(2)}\right)^{2}.$

By using Jensen’s Inequality for integrals we get that
$\left(K_{+}^{(1)}\right)^{2}<\frac{1}{2}K^{(2)}$ and we conclude that the
determinant is strictly positive which completes the proof. ∎

The matrix $\kappa(K)$ will play a role in Lemma A.3 below. We will also often
need computations of the following type and therefore we do it here once as a
reference: Suppose that $f:\mathbb{R}\to\mathbb{R}$ is twice one-sided
differentiable at zero, that is, the first two derivatives of
$f_{m}:(-\infty,0)\to\mathbb{R},x\mapsto f(x)$ and
$f_{p}:(0,\infty)\to\mathbb{R},x\mapsto f(x)$ exist and can be continuously
extended to zero. We write $f_{-}$ for the extension of $f_{m}$ to zero and
$f_{+}$ for the extension of $f_{p}$ at zero. Define
$f_{-}^{\prime},f_{+}^{\prime},f_{-}^{\prime\prime}$, and
$f_{+}^{\prime\prime}$ in a similar way. Then,

$\displaystyle\mathbb{E}\left(\frac{1}{h}L\left(\frac{X_{i}}{h}\right)f(X_{i})\right)=\int_{-\infty}^{0}L(u)f(uh)f_{X}(uh)du+\int_{0}^{-\infty}L(u)f(uh)f_{X}(uh)du$
$\displaystyle=$ $\displaystyle
L_{-}^{(0)}f_{-}f_{X}(0)+L_{+}^{(0)}f_{+}f_{X}(0)+h\left[L_{-}^{(1)}\left(f\cdot
f_{X}\right)_{-}^{\prime}+L_{+}^{(1)}\left(f\cdot
f_{X}\right)_{+}^{\prime}\right]$
$\displaystyle\quad+\frac{1}{2}h^{2}\left[L_{-}^{(2)}\left(f\cdot
f_{X}\right)_{-}^{\prime\prime}+L_{+}^{(2)}\left(f\cdot
f_{X}\right)_{+}^{\prime\prime}\right]+o(h^{2}).$ (A.7)

A simple consequence of the above is the following Lemma:

######  Lemma A.3.

Let $f_{X}$ be twice continuously differentiable and let $K^{(\alpha)}$ for
$\alpha\in\\{0,...,4\\}$ and $\left(K^{2}\right)^{(\alpha)}$ for
$\alpha\in\\{0,1,2\\}$ be finite. If $h\to 0$ and $nh\to\infty$, then

$\displaystyle\frac{1}{n}\sum_{i=1}^{n}K_{h}(X_{i})V_{i}V_{i}^{\top}=\mathbb{E}\left(K_{h}(X_{i})V_{i}V_{i}^{\top}\right)+O_{P}\left(\frac{1}{nh}\right),$
(A.8)
$\displaystyle\mathbb{E}(K_{h}(X_{i})V_{i}V_{i}^{\top})=f_{X}(0)\kappa(K)+f_{X}^{\prime}(0)h\begin{pmatrix}K^{(1)}&K^{(1)}_{+}&K^{(2)}&K_{+}^{(2)}\\\
K_{+}^{(1)}&K_{+}^{(1)}&K_{+}^{(2)}&K_{+}^{(2)}\\\
K^{(2)}&K_{+}^{(2)}&K^{(3)}&K_{+}^{(3)}\\\
K_{+}^{(2)}&K_{+}^{(2)}&K_{+}^{(3)}&K_{+}^{(3)}\end{pmatrix}$
$\displaystyle\quad+\frac{h^{2}}{2}f_{X}^{\prime\prime}(0)\begin{pmatrix}K^{(2)}&K_{+}^{(2)}&K^{(3)}&K_{+}^{(3)}\\\
K_{+}^{(2)}&K_{+}^{(2)}&K_{+}^{(3)}&K_{+}^{(3)}\\\
K^{(3)}&K_{+}^{(3)}&K^{(4)}&K_{+}^{(4)}\\\
K_{+}^{(3)}&K_{+}^{(3)}&K_{+}^{(4)}&K_{+}^{(4)}\end{pmatrix}+o(h^{2})$ (A.9)

###### Proof.

For the proof we have to compute the expectation and the variance of the
average by means of (A.7). ∎

Lastly, we consider in this preliminary discussion a local version of the
Bernstein Inequality. For the convenience of the reader, we state the regular
Bernstein Inequality as it can be found e.g. in Giné and Nickl (2016).

######  Proposition A.4.

Let $A_{i}$, $i=1,...,n$ be a sequence of independent, centered random
variables such that there are numbers $c$ and $\sigma_{i}$ such that for all
$m\in\mathbb{N}$
$\mathbb{E}(|A_{i}|^{m}|)\leq\frac{m!}{2}\sigma_{i}^{2}c^{m-2}$. Set
$\sigma^{2}=\sum_{i=1}^{n}\sigma_{i}^{2}$, $S_{n}=\sum_{i=1}^{n}A_{i}$. Then,
for all $t\geq 0$, $\mathbb{P}(S_{n}\geq
t)\leq\exp\left(-\frac{t^{2}}{2(\sigma^{2}+ct)}\right)$.

In our setting, the following local version will be relevant.

######  Lemma A.5.

Let $B_{1},...,B_{n}$ be iid and denote
$\mu_{m}(x)=\mathbb{E}\left(\left|B_{i}\right|^{m}\big{|}X_{i}=x\right)$
(which is independent of $i=1,...n$). Suppose that for all $m\geq 2$ and
$n\in\mathbb{N}$

$\displaystyle\int_{\mathbb{R}}K(u)^{m}\mu_{m}(uh)f_{X}(uh)du\leq\frac{m!}{2}\sigma_{0}^{2}c^{m-2},$
$\displaystyle\left|\int_{\mathbb{R}}K(u)\mathbb{E}(B_{i}|X_{i}=uh)f_{X}(uh)du\right|\leq
c^{*}$

for some constants $\sigma_{0}^{2},c,c^{*}>0$. Let furthermore $p_{n}$ and $h$
be such that (after possibly increasing $\sigma_{0}$ and $c$) for all $m\geq
2$ and all $n\in\mathbb{N}$

$\frac{m!}{2}\sigma_{0}^{2}c^{m-2}\left(c^{*}\right)^{-m}\geq h^{m-1}\textrm{
and }\sqrt{\frac{\log p_{n}}{nh}}\leq\frac{1}{8c}.$

It holds for all $x\geq 16\sigma_{0}^{2}$ that

$\mathbb{P}\left(\frac{1}{n}\sum_{i=1}^{n}\left(K_{h}(X_{i})B_{i}-\mathbb{E}\left(K_{h}(X_{i})B_{i}\right)\right)>x\sqrt{\frac{\log
p_{n}}{nh}}\right)\leq\left(\frac{1}{p_{n}}\right)^{x}.$

###### Proof.

We begin by rewriting the term of interest as follows. Let
$\varepsilon_{n}=x\sqrt{\frac{\log p_{n}}{nh}}$. Then,

$\displaystyle\mathbb{P}\left(\frac{1}{n}\sum_{i=1}^{n}\left(K_{h}(X_{i})B_{i}-\mathbb{E}\left(K_{h}(X_{i})B_{i}\right)\right)>\varepsilon_{n}\right)=\mathbb{P}\left(\sum_{i=1}^{n}A_{n,i}>nh\varepsilon_{n}\right),$

where

$A_{n,i}=K\left(\frac{X_{i}}{h}\right)B_{i}-\mathbb{E}\left(K\left(\frac{X_{i}}{h}\right)B_{i}\right).$

We apply now Bernstein’s Inequality (cf. Proposition A.4) to $A_{n,i}$. Since
for $a,b\geq 0$, we have that $(a+b)^{m}\leq 2^{m-1}(a^{m}+b^{m})$, we have
for any $i=1,...,n$ and any $m\in\mathbb{N}$ by assumption

$\displaystyle\mathbb{E}\left(\left|A_{i,n}\right|^{m}\right)\leq
2^{m-1}\left(\mathbb{E}\left(K\left(\frac{X_{i}}{h}\right)^{m}\left|B_{i}\right|^{m}\right)+\left|\mathbb{E}\left(K\left(\frac{X_{i}}{h}\right)B_{i}\right)\right|^{m}\right)$
$\displaystyle\leq$ $\displaystyle
2^{m-1}\left(h\int_{\mathbb{R}}K(u)^{m}\mu_{m}(uh)f_{X}(uh)du+\left|h\int_{\mathbb{R}}K(u)\mathbb{E}(B_{i}|X_{i}=uh)f_{X}(uh)du\right|^{m}\right)$
$\displaystyle\leq$ $\displaystyle
2^{m-1}\left(\frac{m!}{2}h\sigma_{0}^{2}c^{m-2}+h^{m}\left(c^{*}\right)^{m}\right)\leq\frac{m!}{2}\cdot
4h\sigma_{0}^{2}\cdot(2c)^{m-2}.$

We may thus apply Proposition A.4 with $\sigma^{2}=4nh\sigma_{0}^{2}$ and
”$c=2c$”. We conclude

$\displaystyle\mathbb{P}\left(\frac{1}{n}\sum_{i=1}^{n}\left(K_{h}(X_{i})B_{i}-\mathbb{E}\left(K_{h}(X_{i})B_{i}\right)\right)>\varepsilon_{n}\right)$
$\displaystyle\leq$
$\displaystyle\exp\left(-\frac{n^{2}h^{2}\varepsilon_{n}^{2}}{2\left(4nh\sigma_{0}^{2}+2cnh\varepsilon_{n}\right)}\right)=\exp\left(-\frac{x^{2}nh\log
p_{n}}{2\left(4nh\sigma_{0}^{2}+2cx\sqrt{nh\log p_{n}}\right)}\right)$
$\displaystyle\leq$ $\displaystyle\exp\left(-\frac{x^{2}\log
p_{n}}{2\left(4\sigma_{0}^{2}+2cx\sqrt{\frac{\log
p_{n}}{nh}}\right)}\right)\leq\exp\left(-x\log p_{n}\right)$

by the assumptions on $x$ and $\log p_{n}/nh$. ∎

###  A.3 The Result for a General Model Selection Algorithm

####  A.3.1 Statement of the Result and Proof Structure

######  Theorem A.1.

Let $p_{n}\to\infty$, $h\to 0$, $nh\to\infty$ and $\log p_{n}/nh\to 0$ and let
$K$ be symmetric, compactly supported with $K^{(4)},(K^{2})^{(2)}<\infty$.
Suppose that (CMS) holds and let $f_{X}$ and $\mu_{Z^{(k)}}$ be twice
differentiable in a neighbourhood around zero with $f_{X}(0)>0$ and
$f_{X}^{\prime\prime}$ being continuous at zero and
$\mu_{Z^{(k)}}(0)=\mu_{Z^{(k)}}^{\prime}(0)=0$ as well as

$\sup_{n\in\mathbb{N}}\sup_{k\in\\{1,...,p_{n}\\}}\sup_{u\in[0,1]}\left|\mu_{Z^{(k)}}^{\prime\prime}(uh)\right|+\left|\mu_{Z^{(k)}}^{\prime\prime}(-uh)\right|<\infty.$
(A.10)

Moreover, suppose that for
$\mu_{k,m}(x)=\mathbb{E}\left(\left|Z_{i}^{(k)}\right|^{m}\Big{|}X_{i}=x\right)$
there are finite numbers $\sigma_{0}^{2},c,c^{*}$ such that for all natural
$m\geq 2$ and all $k$

$\displaystyle\int_{\mathbb{R}}\left(1+|u|^{m}\right)K(u)^{m}\mu_{k,m}(uh)f_{X}(uh)du\leq$
$\displaystyle\frac{m!}{2}\sigma_{0}^{2}c^{m-2},$ (A.11)
$\displaystyle\int_{\mathbb{R}}\left(1+|u|\right)K(u)^{m}\mu_{k,1}(uh)f_{X}(uh)du\leq$
$\displaystyle c^{*}.$ (A.12)

Suppose that for the target set $J_{n}$, there are $\delta>0$ and finite
numbers $\sigma_{l},\sigma_{r},C>0$ such that

$\displaystyle\lim_{n\to\infty}\sup_{u\in[0,1]}\left|\mathbb{E}(r_{i}(J_{n},h)^{2}|X_{i}=uh)-\sigma_{r}^{2}\right|=0,$
$\displaystyle\quad\quad\quad\quad\lim_{n\to\infty}\sup_{u\in[0,1]}\left|\mathbb{E}(r_{i}(J_{n},h)^{2}|X_{i}=-uh)-\sigma_{l}^{2}\right|=0,$
(A.13)
$\displaystyle\sup_{n\in\mathbb{N}}\sup_{x\in[-h,h]}\left|\mathbb{E}(|r_{i}(J_{n},h)|^{2+\delta}|X_{i}=x)\right|<C.$
(A.14)

Set $w=\left(\left(f_{X}(0)\kappa(K)\right)^{-1}\right)_{2\cdot}^{\top}$ to be
the scaled second row of the inverse of $\kappa(K)$. Define

$\displaystyle\mathcal{S}_{n}^{2}=\frac{1}{h}\mathbb{E}\left(K\left(\frac{X_{i}}{h}\right)^{2}\left(w^{\top}V_{i}\right)^{2}r_{i}(J_{n},h)^{2}\right).$

Then,

$\displaystyle\sqrt{\frac{nh}{\mathcal{S}^{2}_{n}}}\left(\hat{\tau}_{n}(\hat{J}_{n})-\theta_{0,n}^{(2)}\right)\overset{d}{\to}\mathcal{N}\left(0,1\right).$

###### Proof of Theorem A.1.

Note that the conditions of this theorem contain all conditions of the
supporting results (from the following section) or imply them (e.g. (A.13) is
stronger than (A.35)). Therefore we can use all results from Section A.3.2.
Let

$\mathbf{M}_{n}(\hat{J}_{n})=\mathbf{I}_{n}-\mathbf{K}_{h}^{\frac{1}{2}}\mathbf{Z}(\hat{J}_{n})\left(\mathbf{Z}(\hat{J}_{n})^{\top}\mathbf{K}_{h}\mathbf{Z}(\hat{J}_{n})\right)^{-1}\mathbf{Z}(\hat{J}_{n})^{\top}\mathbf{K}_{h}^{\frac{1}{2}}$
(A.15)

denote the projection matrix on the chosen covariates ($\mathbf{I}_{n}$
denotes $n$-dimensional identity matrix). We write here
$\mathbf{r}_{n}=\mathbf{r}(J_{n},h)$ and $\gamma_{0,n}=\gamma_{0,n}(J_{n},h)$
because $J_{n}$ and $h$ will be the same sequences throughout the proof.
Moreover, $\gamma_{0,n}$ will be understood as element of $\mathbb{R}^{p_{n}}$
with $\gamma_{0,n}^{(k)}=0$ for $k\notin J_{n}$. Note that
$\mathbf{M}_{n}\mathbf{K}_{h}^{\frac{1}{2}}\mathbf{Z}(\hat{J}_{n})=0$. We thus
obtain by calculation (or the Frisch-Waugh-Lovell Theorem for weighted
regression) the following representation of our estimator

$\displaystyle\hat{\theta}_{n}=$
$\displaystyle\left(\mathbf{V}^{\top}\mathbf{K}_{h}^{\frac{1}{2}}\mathbf{M}_{n}(\hat{J}_{n})\mathbf{K}_{h}^{\frac{1}{2}}\mathbf{V}\right)^{-1}\mathbf{V}^{\top}\mathbf{K}_{h}^{\frac{1}{2}}\mathbf{M}_{n}(\hat{J}_{n})\mathbf{K}_{h}^{\frac{1}{2}}\mathbf{Y}$
$\displaystyle=$
$\displaystyle\theta_{0,n}+\left(\mathbf{V}^{\top}\mathbf{K}_{h}^{\frac{1}{2}}\mathbf{M}_{n}(\hat{J}_{n})\mathbf{K}_{h}^{\frac{1}{2}}\mathbf{V}\right)^{-1}\mathbf{V}^{\top}\mathbf{K}_{h}^{\frac{1}{2}}\mathbf{M}_{n}(\hat{J}_{n})\mathbf{K}_{h}^{\frac{1}{2}}\left(\mathbf{Z}\gamma_{0,n}+\mathbf{r}_{n}\right).$

Suppose for the moment that we know that

$\frac{1}{n}\mathbf{V}^{\top}\mathbf{K}^{\frac{1}{2}}M_{n}(\hat{J}_{n})\mathbf{K}_{h}^{\frac{1}{2}}\mathbf{Z}\gamma_{0,n}=o_{P}\left(\frac{1}{\sqrt{nh}}\right).$
(A.16)

If the above is true we obtain together with Proposition A.7, Lemma A.8, the
assumptions on $\delta_{n}$ and the fact that $\kappa(K)$ is invertible by
Lemma A.2 that

$\displaystyle\left(\frac{1}{n}\mathbf{V}^{\top}\mathbf{K}_{h}^{\frac{1}{2}}\mathbf{M}_{n}(\hat{J}_{n})\mathbf{K}_{h}^{\frac{1}{2}}\mathbf{V}\right)^{-1}=\left(f_{X}(0)\kappa(K)\right)^{-1}+o_{P}(1),$
$\displaystyle\frac{1}{n}\mathbf{V}^{\top}\mathbf{K}_{h}^{\frac{1}{2}}\mathbf{M}_{n}(\hat{J}_{n})\mathbf{K}_{h}^{\frac{1}{2}}\left(\mathbf{Z}\gamma_{0,n}+\mathbf{r}_{n}\right)=\frac{1}{n}\mathbf{V}^{\top}\mathbf{K}_{h}\mathbf{r}_{n}+o_{P}\left(\frac{1}{\sqrt{nh}}\right),$
$\displaystyle\frac{1}{n}\mathbf{V}^{\top}\mathbf{K}_{h}\mathbf{r}_{n}=O_{P}\left(\frac{1}{\sqrt{nh}}\right).$

Hence,

$\sqrt{nh}\left(\hat{\theta}_{n}-\theta_{0,n}\right)=\left(f_{X}(0)\kappa(K)\right)^{-1}\sqrt{nh}\frac{1}{n}\mathbf{V}^{\top}\mathbf{K}_{h}\mathbf{r}_{n}+o_{P}(1).$

Thus, to find the asymptotics of the estimator for the treatment effect we
study the second entry of the above vector. Recall to this end that $w$
denotes the second row of the kernel matrix written as a column. Thus we have

$\sqrt{nh}\left(\hat{\tau}_{n}-\theta_{0,n}^{(2)}\right)=\sqrt{nh}\frac{1}{n}w^{\top}\mathbf{V}^{\top}\mathbf{K}_{h}\mathbf{r}_{n}+o_{P}(1)=\frac{1}{\sqrt{nh}}\sum_{i=1}^{n}K\left(\frac{X_{i}}{h}\right)w^{\top}V_{i}r_{i}(J_{n},h)+o_{P}(1).$

For simplicity of notation we write $\nu(X_{i}/h)=w^{\top}V_{i}$. Now we can
employ Lyapunov’s central limit theorem (cf. Lemma 15.41 and Theorem 15.43 in
Klenke (2008)). We have by definition

$\displaystyle\textrm{Var}\left(\frac{1}{\sqrt{nh}}\sum_{i=1}^{n}K\left(\frac{X_{i}}{h}\right)\nu\left(\frac{X_{i}}{h}\right)r_{i}(J_{n},h)\right)$
$\displaystyle=$
$\displaystyle\frac{1}{nh}\sum_{i=1}^{n}\mathbb{E}\left(K\left(\frac{X_{i}}{h}\right)^{2}\nu\left(\frac{X_{i}}{h}\right)^{2}r_{i}(J_{n},h)^{2}\right)$
$\displaystyle=$
$\displaystyle\frac{1}{h}\mathbb{E}\left(K\left(\frac{X_{i}}{h}\right)^{2}\nu\left(\frac{X_{i}}{h}\right)^{2}r_{i}(J_{n},h)^{2}\right)=\mathcal{S}_{n}^{2}.$

By the continuity and boundedness assumptions (A.13) and (A.14) on
$\mathbb{E}(r_{i}(J_{n},h)^{2}|X_{i})$ and
$\mathbb{E}(|r_{n,ij}|^{2+\delta}|X_{i})$, we conclude that there are
constants $\alpha_{0},\alpha_{1}>0$ such that for some $\delta>0$ and
$n\to\infty$

$\displaystyle\frac{1}{h}\mathbb{E}\left(K\left(\frac{X_{i}}{h}\right)^{2}\nu\left(\frac{X_{i}}{h}\right)^{2}r_{i}(J_{n},h)^{2}\right)\to\alpha_{0},$
$\displaystyle\frac{1}{h}\mathbb{E}\left(K\left(\frac{X_{i}}{h}\right)^{2+\delta}\left|\nu\left(\frac{X_{i}}{h}\right)\right|^{2+\delta}|r_{i}(J_{n},h)|^{2+\delta}\right)\leq\alpha_{1}.$

Thus the Lyapunov condition is fulfilled: For $n\to\infty$

$\displaystyle\sum_{i=1}^{n}\frac{n^{-\frac{2+\delta}{2}}h^{-\frac{2+\delta}{2}}\mathbb{E}\left(K\left(\frac{X_{i}}{h}\right)^{2+\delta}\left|\nu\left(\frac{X_{i}}{h}\right)\right|^{2+\delta}|r_{i}(J_{n},h)|^{2+\delta}\right)}{\mathcal{S}_{n}^{2+\delta}}\leq\frac{\alpha_{1}(nh)^{-\frac{\delta}{2}}}{\left(\frac{\alpha_{0}}{2}\right)^{\frac{2+\delta}{2}}}\to
0$

and we conclude

$\frac{1}{\sqrt{nh\mathcal{S}_{n}^{2}}}\sum_{i=1}^{n}K\left(\frac{X_{i}}{h}\right)\nu\left(\frac{X_{i}}{h}\right)r_{i}(J_{n},h)\to\mathcal{N}(0,1).$

In order to finish the proof, we have to show (A.16) which we will do next.
Recall to this end the post selection estimators

$\left(\hat{\theta}_{n},\hat{\gamma}_{n}\right)=\underset{\theta,\gamma:\gamma_{\hat{J}_{n}^{c}}=0}{\operatorname{argmin}}\left\|\mathbf{K}_{h}^{\frac{1}{2}}\left(\mathbf{Y}-\mathbf{V}\theta-\mathbf{Z}\gamma\right)\right\|_{2}^{2}.$

Note that $\hat{\gamma}_{n}$ is constrained to be zero for covariates not
included in $\hat{J}_{n}$. Therefore, we have
$M_{n}(\hat{J}_{n})\mathbf{K}^{\frac{1}{2}}_{h}\mathbf{Z}\hat{\gamma}_{n}=0$.
We can hence write

$\displaystyle\frac{1}{n}\mathbf{V}^{\top}\mathbf{K}_{h}^{\frac{1}{2}}M_{n}(\hat{J}_{n})\mathbf{K}_{h}^{\frac{1}{2}}\mathbf{Z}\gamma_{0,n}=\frac{1}{n}\mathbf{V}^{\top}\mathbf{K}_{h}^{\frac{1}{2}}M_{n}(\hat{J}_{n})\mathbf{K}_{h}^{\frac{1}{2}}\mathbf{Z}\left(\gamma_{0,n}-\hat{\gamma}_{n}\right)$
$\displaystyle=$
$\displaystyle\frac{1}{n}\mathbf{V}^{\top}\mathbf{K}_{h}\mathbf{Z}\left(\gamma_{0,n}-\hat{\gamma}_{n}\right)+\frac{1}{n}\mathbf{V}^{\top}\mathbf{K}_{h}\mathbf{Z}(\hat{J}_{n})\left(\mathbf{Z}(\hat{J}_{n})^{\top}\mathbf{K}_{h}\mathbf{Z}(\hat{J}_{n})\right)^{-1}\mathbf{Z}(\hat{J}_{n})^{\top}\mathbf{K}_{h}\left(\mathbf{Z}\gamma_{0,n}-\mathbf{Z}\hat{\gamma}_{n}\right)$
$\displaystyle=$
$\displaystyle\frac{1}{n}\mathbf{V}^{\top}\mathbf{K}_{h}\mathbf{Z}\left(\gamma_{0,n}-\hat{\gamma}_{n}\right)$
(A.17)
$\displaystyle+\frac{1}{n}\mathbf{V}^{\top}\mathbf{K}_{h}\mathbf{Z}(\hat{J}_{n})\left(\mathbf{Z}(\hat{J}_{n})^{\top}\mathbf{K}_{h}\mathbf{Z}(\hat{J}_{n})\right)^{-1}\mathbf{Z}(\hat{J}_{n})^{\top}\mathbf{K}_{h}\mathbf{V}\left(\hat{\theta}_{n}-\theta_{0,n}\right)$
(A.18)
$\displaystyle-\frac{1}{n}\mathbf{V}^{\top}\mathbf{K}_{h}\mathbf{Z}(\hat{J}_{n})\left(\mathbf{Z}(\hat{J}_{n})^{\top}\mathbf{K}_{h}\mathbf{Z}(\hat{J}_{n})\right)^{-1}\mathbf{Z}(\hat{J}_{n})^{\top}\mathbf{K}_{h}\mathbf{r}_{n}$
(A.19)
$\displaystyle+\frac{1}{n}\mathbf{V}^{\top}\mathbf{K}_{h}\mathbf{Z}(\hat{J}_{n})\left(\mathbf{Z}(\hat{J}_{n})^{\top}\mathbf{K}_{h}\mathbf{Z}(\hat{J}_{n})\right)^{-1}\mathbf{Z}(\hat{J}_{n})^{\top}\mathbf{K}_{h}\left(\mathbf{Y}-\mathbf{Z}\hat{\gamma}_{n}-\mathbf{V}\hat{\theta}_{n}\right).$
(A.20)

Note that (A.20) equals zero because it contains the empirical correlation of
a covariate with the empirical residuals (which is zero). From the definition
of $M_{n}(\hat{J}_{n})$ it follows from (A.33) of Proposition A.7 that
$\eqref{eq:l3}=o_{P}(1/\sqrt{nh})$ by the conditions on $\delta_{n}$. Hence,
in order to prove (A.16), we have to prove that (A.17) and (A.18) are both of
order $o_{P}(1/\sqrt{nh})$. We do this by studying the rate of convergence of
$\hat{\theta}_{n}$ and $\hat{\gamma}_{n}$. We note firstly that

$\frac{1}{n}\begin{pmatrix}\theta_{0,n}-\hat{\theta}_{n}\\\
\gamma_{0,n}-\hat{\gamma}_{n}\end{pmatrix}^{\top}\begin{pmatrix}\mathbf{V}^{\top}\\\
\mathbf{Z}^{\top}\end{pmatrix}\mathbf{K}_{h}\begin{pmatrix}\mathbf{V}&\mathbf{Z}\end{pmatrix}\begin{pmatrix}\theta_{0,n}-\hat{\theta}_{n}\\\
\gamma_{0,n}-\hat{\gamma}_{n}\end{pmatrix}\geq\varphi(|\hat{J}_{n}|,J_{n})\left\|\begin{pmatrix}\theta_{0,n}-\hat{\theta}_{n}\\\
\gamma_{0,n}-\hat{\gamma}_{n}\end{pmatrix}\right\|_{2}^{2},$

where $\varphi$ is defined as in Definition 3.2. Since we assume the
restricted sparse eigenvalue condition $\textrm{RSE}(|\hat{J}_{n}|,J_{n},h)$
we have that $\varphi(|\hat{J}_{n}|,J_{n})^{-1}=O_{P}(1)$ and we conclude that

$\left\|\begin{pmatrix}\theta_{0,n}-\hat{\theta}_{n}\\\
\gamma_{0,n}-\hat{\gamma}_{n}\end{pmatrix}\right\|_{2}=O_{P}\left(\frac{1}{\sqrt{n}}\left\|\mathbf{K}_{h}^{\frac{1}{2}}\begin{pmatrix}\mathbf{V}&\mathbf{Z}\end{pmatrix}\begin{pmatrix}\theta_{0,n}-\hat{\theta}_{n}\\\
\gamma_{0,n}-\hat{\gamma}_{n}\end{pmatrix}\right\|_{2}\right).$ (A.21)

The proof strategy is now similar to Belloni and Chernozhukov (2013). At first
we note that, by definition for any index set $I_{n}$

$\displaystyle\frac{1}{n}\left\|\mathbf{K}^{\frac{1}{2}}\left(\mathbf{Y}-\mathbf{V}\hat{\theta}_{n}-\mathbf{Z}\hat{\gamma}_{n}\right)\right\|_{2}^{2}-\frac{1}{n}\left\|\mathbf{K}^{\frac{1}{2}}\left(\mathbf{Y}-\mathbf{V}\theta_{0,n}-\mathbf{Z}\gamma_{0,n}\right)\right\|_{2}^{2}$
$\displaystyle\leq$
$\displaystyle\left\\{\begin{array}[]{l}\frac{1}{n}\left\|\mathbf{K}^{\frac{1}{2}}\left(\mathbf{Y}-\mathbf{V}\theta_{0,n}-\mathbf{Z}\widetilde{\gamma}_{n}\right)\right\|_{2}^{2}-\frac{1}{n}\left\|\mathbf{K}^{\frac{1}{2}}\left(\mathbf{Y}-\mathbf{V}\theta_{0,n}-\mathbf{Z}\gamma_{0,n}\right)\right\|_{2}^{2}=:B_{n}\\\
\frac{1}{n}\left\|\mathbf{K}^{\frac{1}{2}}\left(\mathbf{Y}-\mathbf{V}\theta_{0,n}-\mathbf{Z}\left(\gamma_{0,n}\right)_{\hat{J}_{n}\cap
I_{n}}\right)\right\|_{2}^{2}\\\
\quad\quad\quad\quad\quad\quad\quad-\frac{1}{n}\left\|\mathbf{K}^{\frac{1}{2}}\left(\mathbf{Y}-\mathbf{V}\theta_{0,n}-\mathbf{Z}\gamma_{0,n}\right)\right\|_{2}^{2}=C_{n}(I_{n})\end{array}\right..$
(A.25)

Since the index set $I_{n}$ will not change during this proof, we write in the
following $C_{n}=C_{n}(I_{n})$. Let now,
$\hat{\alpha}_{n}=(\hat{\theta}_{n}^{\top},\hat{\gamma}_{n}^{\top})^{\top}-(\theta_{0,n}^{\top},\gamma_{0,n}^{\top})^{\top}$.
Then,
$\left\|\left(\hat{\theta}_{n}-\theta_{0,n}\right)_{J_{n}^{c}}\right\|_{0}\leq|\hat{J}_{n}|$
and hence we can apply Lemma A.9 with $m_{n}=|\hat{J}_{n}|$ to obtain that

$\displaystyle\frac{1}{n}\Bigg{|}\left\|\mathbf{K}_{h}^{\frac{1}{2}}\left(\mathbf{Y}-\begin{pmatrix}\mathbf{V}&\mathbf{Z}\end{pmatrix}\begin{pmatrix}\hat{\theta}_{n}\\\
\hat{\gamma}_{n}\end{pmatrix}\right)\right\|_{2}^{2}-\left\|\mathbf{K}_{h}^{\frac{1}{2}}\left(\mathbf{Y}-\begin{pmatrix}\mathbf{V}&\mathbf{Z}\end{pmatrix}\begin{pmatrix}\theta_{0,n}\\\
\gamma_{0,n}\end{pmatrix}\right)\right\|_{2}^{2}$
$\displaystyle\quad\quad\quad\quad\quad-\left\|\mathbf{K}_{h}^{\frac{1}{2}}\begin{pmatrix}\mathbf{V}&\mathbf{Z}\end{pmatrix}\hat{\alpha}_{n}\right\|_{2}^{2}\Bigg{|}$
$\displaystyle\leq$
$\displaystyle\rho_{n}(|\hat{J}_{n}|)\frac{1}{\sqrt{n}}\left\|\mathbf{K}_{h}^{\frac{1}{2}}\begin{pmatrix}\mathbf{V}&\mathbf{Z}\end{pmatrix}\hat{\alpha}_{n}\right\|_{2}$

which in turn implies together with (A.25)

$\displaystyle\frac{1}{n}\left\|\mathbf{K}_{h}^{\frac{1}{2}}\begin{pmatrix}\mathbf{V}&\mathbf{Z}\end{pmatrix}\hat{\alpha}_{n}\right\|_{2}^{2}-\rho_{n}(|\hat{J}_{n}|)\frac{1}{\sqrt{n}}\left\|\mathbf{K}_{h}^{\frac{1}{2}}\begin{pmatrix}\mathbf{V}&\mathbf{Z}\end{pmatrix}\hat{\alpha}_{n}\right\|_{2}$
$\displaystyle\leq$
$\displaystyle\frac{1}{n}\left\|\mathbf{K}_{h}^{\frac{1}{2}}\left(\mathbf{Y}-\begin{pmatrix}\mathbf{V}&\mathbf{Z}\end{pmatrix}\begin{pmatrix}\hat{\theta}_{n}\\\
\hat{\gamma}_{n}\end{pmatrix}\right)\right\|_{2}^{2}-\frac{1}{n}\left\|\mathbf{K}_{h}^{\frac{1}{2}}\left(\mathbf{Y}-\begin{pmatrix}\mathbf{V}&\mathbf{Z}\end{pmatrix}\begin{pmatrix}\theta_{0,n}\\\
\gamma_{0,n}\end{pmatrix}\right)\right\|_{2}^{2}\leq\min(B_{n},C_{n}).$

It is elementary to prove that $x_{n}^{2}-p_{n}x_{n}\leq q_{n}$ for non-
negative sequences $x_{n},p_{n}$ and an arbitrary sequence $q_{n}$ implies
$x_{n}\leq p_{n}+\sqrt{\max(q_{n},0)}$. Hence, the above gives us

$\frac{1}{\sqrt{n}}\left\|\mathbf{K}_{h}^{\frac{1}{2}}\begin{pmatrix}\mathbf{V}&\mathbf{Z}\end{pmatrix}\hat{\alpha}_{n}\right\|_{2}\leq\rho_{n}(|\hat{J}_{n}|)+\sqrt{\max(0,\min(B_{n},C_{n}))}.$

Thus, we obtain from (A.21) that

$\left\|\begin{pmatrix}\theta_{0,n}-\hat{\theta}_{n}\\\
\gamma_{0,n}-\hat{\gamma}_{n}\end{pmatrix}\right\|_{2}=O_{P}\left(\rho_{n}(|\hat{J}_{n}|)+\sqrt{\max(0,\min(B_{n},C_{n}))}\right).$
(A.26)

We use this to show that (A.17) and (A.18) are both $o_{P}(1/\sqrt{nh})$. Note
that (A.17) and (A.18) are both vectors of length $4$ and hence we may use any
norm to study the asymptotics. Moreover, for any matrix $M$, we denote by
$|M|_{\infty}$ the supremum over the absolute values of the entries of $M$.
Set moreover $\zeta_{n}=\sqrt{\max(0,\min(B_{n},C_{n}))}$. We begin with
(A.17). Note therefore that
$\left\|\gamma_{0,n}-\hat{\gamma}_{n}\right\|_{1}\leq(|J_{n}|+|\hat{J}_{n}|)^{\frac{1}{2}}\left\|\gamma_{0,n}-\hat{\gamma}_{n}\right\|_{2}$.
Then we obtain for some constant $C>0$ from Lemma A.6 and by using (A.26) in
which we substitute the rate of $\rho(|\hat{J}_{n}|)$ from Lemma A.9:

$\displaystyle\sqrt{nh}\|\eqref{eq:l1}\|_{1}\leq\sqrt{nh}\left|\frac{1}{n}\mathbf{V}^{\top}\mathbf{K}_{h}\mathbf{Z}\right|_{\infty}\left\|\gamma_{0,n}-\hat{\gamma}_{n}\right\|_{1}$
$\displaystyle\leq$ $\displaystyle C\sqrt{nh}\left(\sqrt{\frac{\log
p_{n}}{nh}}+h^{2}\right)\left(|J_{n}|+|\hat{J}_{n}|\right)^{\frac{1}{2}}\left(\left(|J_{n}|+|\hat{J}_{n}|\right)^{\frac{1}{2}}\left(\frac{1}{\sqrt{nh}}+\delta_{n}\right)+\zeta_{n}\right)$
$\displaystyle=$ $\displaystyle C\left(\sqrt{\log
p_{n}}+\sqrt{nh^{5}}\right)\left(\left(|J_{n}|+|\hat{J}_{n}|\right)\left(\frac{1}{\sqrt{nh}}+\delta_{n}\right)+\left(|J_{n}|+|\hat{J}_{n}|\right)^{\frac{1}{2}}\zeta_{n}\right).$

The above converges to zero by the assumptions on the rates in (CMS). Hence,
$\eqref{eq:l1}=o_{P}(1/\sqrt{nh})$. For (A.18) we obtain, for a different
constant $C>0$, by applying the bound (A.34) together with Lemma A.6 and
(A.26) where we again replace $\rho(|\hat{J}_{n}|)$ as in Lemma A.9)

$\displaystyle\sqrt{nh}\|\eqref{eq:l2}\|_{1}\leq\sqrt{nh}\left|\frac{1}{n}\mathbf{V}^{\top}\mathbf{K}_{h}\mathbf{Z}(\hat{J}_{n})\left(\mathbf{Z}(\hat{J}_{n})^{\top}\mathbf{K}_{h}\mathbf{Z}(\hat{J}_{n})\right)^{-1}\mathbf{Z}(\hat{J}_{n})^{\top}\mathbf{K}_{h}\mathbf{V}\right|_{\infty}\left\|\hat{\theta}_{n}-\theta_{0,n}\right\|_{1}$
$\displaystyle\leq$ $\displaystyle
C|\hat{J}_{n}|\sqrt{nh}\sqrt{\Phi\left(\left|\hat{J}_{n}\right|,J_{n}\right)}\left(\sqrt{\frac{\log
p_{n}}{nh}}+h^{2}\right)^{2}$
$\displaystyle\quad\quad\quad\quad\quad\times\left(\left(|J_{n}|+|\hat{J}_{n}|\right)\left(\frac{1}{\sqrt{nh}}+\delta_{n}\right)+\left(|J_{n}|+|\hat{J}_{n}|\right)^{\frac{1}{2}}\zeta_{n}\right)$
$\displaystyle\leq$ $\displaystyle
C|\hat{J}_{n}|\sqrt{\Phi\left(\left|\hat{J}_{n}\right|,J_{n}\right)}\left(\sqrt{\frac{\log
p_{n}}{nh}}+h^{2}\right)$ $\displaystyle\quad\quad\quad\times\left(\sqrt{\log
p_{n}}+\sqrt{nh^{5}}\right)\left(\left(|J_{n}|+|\hat{J}_{n}|\right)\left(\frac{1}{\sqrt{nh}}+\delta_{n}\right)+\left(|\hat{J}_{n}|+|J_{n}|\right)^{\frac{1}{2}}\zeta_{n}\right).$

Since we assume $\textrm{RSE}(|\hat{J}_{n}|,J_{n},h)$ and the rates in (CMS),
we have that the above converges to zero and the proof is complete. ∎

####  A.3.2 Supporting Results

In the following Lemma the difference between regression discontinuity design
and and regular treatment effects becomes visible: We assume only
asymptotically that the covariates and the treatment indicator are
uncorrelated, moreover we do not model explicitly the relation between them.
Thus we can only obtain that the correlation converges to zero.

######  Lemma A.6.

Let $K$ be second order and compactly supported and let $h\to 0$,
$p_{n}\to\infty$ and $\log p_{n}/nh\to 0$. Suppose that $f_{X}$ and
$\mu_{Z^{(k)}}(x)$ are twice differentiable on a neighborhood around $0$ with
$\mu_{Z^{(k)}}(0)=\mu_{Z^{(k)}}^{\prime}(0)=0$ and $f_{X}^{\prime\prime}$
continuous and $\mu_{Z^{(k)}}^{\prime\prime}$ fulfill (A.10). Let furthermore
(A.11) and (A.12) hold. Then, we have

$\displaystyle\sup_{a\in\\{1,...,4\\}}\left\|\mathbb{E}\left(\frac{1}{n}\mathbf{Z}^{\top}\mathbf{K}_{h}\mathbf{V}_{\cdot
a}\right)\right\|_{\infty}=O(h^{2}),$ (A.27)
$\displaystyle\sup_{a\in\\{1,...,4\\}}\left\|\frac{1}{n}\mathbf{Z}^{\top}\mathbf{K}_{h}\mathbf{V}_{\cdot
a}\right\|_{\infty}=O_{P}\left(\sqrt{\frac{\log p_{n}}{nh}}\right)+O(h^{2}).$
(A.28)

###### Proof.

Since $a=1,...,4$ can only take finitely many values, it suffices to prove the
statements (A.27) and (A.28) for an arbitrary $a\in\\{1,...,4\\}$. Let $a$ be
thus fixed. Note that

$\frac{1}{n}\mathbf{Z}^{\top}\mathbf{K}_{h}\mathbf{V}_{\cdot
a}=\frac{1}{n}\sum_{i=1}^{n}K_{h}(X_{i})Z_{i}V_{i}^{(a)}.$

Note that the bounds on $\mu_{Z^{(k)}}^{\prime\prime}$ and the
differentiability of $f_{X}$ imply a similar bound for
$\left(\mu_{Z^{(k)}}f_{X}\right)^{\prime\prime}$, i.e.,

$C=\sup_{n\in\mathbb{N}}\sup_{k\in\\{1,...,p_{n}\\}}\sup_{u\in[0,1]}\left|\left(\mu_{Z^{(k)}}f_{X}\right)^{\prime\prime}(uh)\right|+\left|\left(\mu_{Z^{(k)}}f_{X}\right)^{\prime\prime}(-uh)\right|<\infty.$

We begin by computing the expectation in (A.27). By using that
$\mu_{Z^{(k)}}(x)$ is differentiable and that
$\mu_{Z^{(k)}}(0)=\mu_{Z^{(k)}}^{\prime}(0)=0$ we obtain for all
$k=1,...,p_{n}$ by an argument as in (A.7) that:

$\displaystyle\left|\mathbb{E}\left(K_{h}(X_{i})Z_{i}^{(k)}V_{i}^{(a)}\right)\right|\leq$
$\displaystyle\frac{1}{2}h^{2}\cdot C.$ (A.29)

This implies (A.27) and that

$\left\|\frac{1}{n}\mathbf{Z}^{\top}\mathbf{K}_{h}\mathbf{V}_{\cdot
a}\right\|_{\infty}\leq\left\|\frac{1}{n}\mathbf{Z}^{\top}\mathbf{K}_{h}\mathbf{V}_{\cdot
a}-\mathbb{E}\left(\frac{1}{n}\mathbf{Z}^{\top}\mathbf{K}_{h}\mathbf{V}_{\cdot
a}\right)\right\|_{\infty}+O(h^{2}).$ (A.30)

For the first part above we will apply Lemma A.5 with
$B_{i}=Z_{i}^{(k)}V_{i}^{(a)}$. The integral conditions hold true by
assumption because
$\mathbb{E}\left(\left|Z_{i}^{(k)}V_{i}^{(a)}\right|^{m}\big{|}X_{i}=x\right)\leq(1+(|x|/h)^{m})\mu_{k,m}(x)$.
The conditions on $n,p_{n},h$ follow because we assume that $h\to 0$ and $\log
p_{n}/nh\to 0$. Note furthermore that all constants in the assumptions do not
depend on $k$. We may thus apply Lemma A.5 simultaneously for all
$k=1,...,p_{n}$ and obtain for $x$ large enough as in Lemma A.5 (but $x>1$)
for $\varepsilon_{n}=x\sqrt{\log p_{n}/nh}$ by the union bound that

$\displaystyle\mathbb{P}\left(\max_{k\in\\{1,...,p_{n}\\}}\frac{1}{n}\sum_{i=1}^{n}\left(K_{h}(X_{i})Z_{i}^{(k)}V_{i}^{(a)}-\mathbb{E}\left(K_{h}(X_{i})Z_{i}^{(k)}V_{i}^{(a)}\right)\right)>\varepsilon_{n}\right)$
$\displaystyle\leq$ $\displaystyle
p_{n}\max_{k\in\\{1,...,p_{n}\\}}\mathbb{P}\left(\frac{1}{n}\sum_{i=1}^{n}\left(K_{h}(X_{i})Z_{i}^{(k)}V_{i}^{(a)}-\mathbb{E}\left(K_{h}(X_{i})Z_{i}^{(k)}V_{i}^{(a)}\right)\right)>\varepsilon_{n}\right)$
$\displaystyle\leq$ $\displaystyle p_{n}^{1-x}\to 0$ (A.31)

which implies

$\left\|\frac{1}{n}\mathbf{Z}^{\top}\mathbf{K}_{h}\mathbf{V}_{\cdot
a}-\mathbb{E}\left(\frac{1}{n}\mathbf{Z}^{\top}\mathbf{K}_{h}\mathbf{V}_{\cdot
a}\right)\right\|_{\infty}=O_{P}\left(\sqrt{\frac{\log p_{n}}{nh}}\right).$

This, together with (A.30), completes the proof of (A.28). ∎

The previous results are used in the following way to understand the behaviour
of the projection matrix $\mathbf{M}_{n}(\hat{J}_{n})$.

######  Proposition A.7.

Let all notation be as above and let the assumptions of Lemmas A.3 and A.6 be
true. Suppose that $\textrm{RSE}(|\hat{J}_{n}|,J_{n},h)$ holds and
$\frac{|\hat{J}_{n}|\log p_{n}}{nh}=o_{P}(1)$ as well as
$|\hat{J}_{n}|h^{4}=o_{P}(1)$. We have then that

$\displaystyle\frac{1}{n}\mathbf{V}^{\top}\mathbf{K}_{h}^{\frac{1}{2}}\mathbf{M}_{n}(\hat{J}_{n})\mathbf{K}_{h}^{\frac{1}{2}}\mathbf{V}\overset{\mathbb{P}}{\to}f_{X}(0)\kappa(K)$
(A.32)
$\displaystyle\frac{1}{n}\mathbf{V}^{\top}\mathbf{K}_{h}^{\frac{1}{2}}\mathbf{M}_{n}(\hat{J}_{n})\mathbf{K}_{h}^{\frac{1}{2}}\mathbf{r}_{n}=\frac{1}{n}\mathbf{V}^{\top}\mathbf{K}_{h}\mathbf{r}_{n}+O_{P}\left(\left\|\frac{1}{n}\mathbf{Z}^{\top}\mathbf{K}_{h}\mathbf{r}_{n}\right\|_{\infty}\cdot|\hat{J}_{n}|\left(\sqrt{\frac{\log
p_{n}}{nh}}+h^{2}\right)\right).$ (A.33)

###### Proof.

By definition (A.15) we have

$\frac{1}{n}\mathbf{V}^{\top}\mathbf{K}_{h}^{\frac{1}{2}}\mathbf{M}_{n}(\hat{J}_{n})\mathbf{K}_{h}^{\frac{1}{2}}\mathbf{V}=\frac{1}{n}\mathbf{V}^{\top}\mathbf{K}_{h}\mathbf{V}-\frac{1}{n}\mathbf{V}^{\top}\mathbf{K}_{h}\mathbf{Z}(\hat{J}_{n})\left(\mathbf{Z}(\hat{J}_{n})^{\top}\mathbf{K}_{h}\mathbf{Z}(\hat{J}_{n})\right)^{-1}\mathbf{Z}(\hat{J}_{n})^{\top}\mathbf{K}_{h}\mathbf{V}.$

The first term above converges by Lemma A.3 to the quantity we claim in the
lemma. It remains to show that the second part converges to zero in
probability. Recall to this end the definition of $\Phi(m_{n},J_{n})$ in
Definition 3.2. We note that for all $a,b\in\\{1,...,4\\}$ (below $\|.\|$
denotes the Euclidean norm for vectors and the spectral norm for matrices)

$\displaystyle\left[\frac{1}{n}\mathbf{V}^{\top}\mathbf{K}_{h}\mathbf{Z}(\hat{J}_{n})\left(\mathbf{Z}(\hat{J}_{n})^{\top}\mathbf{K}_{h}\mathbf{Z}(\hat{J}_{n})\right)^{-1}\mathbf{Z}(\hat{J}_{n})^{\top}\mathbf{K}_{h}\mathbf{V}\right]_{a,b}$
$\displaystyle\leq$
$\displaystyle\left\|\left(\frac{1}{n}\mathbf{Z}(\hat{J}_{n})^{\top}\mathbf{K}_{h}\mathbf{Z}(\hat{J}_{n})\right)^{-1}\frac{1}{n}\mathbf{Z}(\hat{J}_{n})^{\top}\mathbf{K}_{h}\mathbf{V}_{\cdot
a}\right\|_{1}\cdot\left\|\frac{1}{n}\mathbf{Z}^{\top}\mathbf{K}_{h}\mathbf{V}_{\cdot
b}\right\|_{\infty}$ $\displaystyle\leq$
$\displaystyle|\hat{J}_{n}|^{\frac{1}{2}}\sqrt{\Phi(|\hat{J}_{n}|,J_{n})}\cdot\left\|\frac{1}{n}\mathbf{Z}(\hat{J}_{n})^{\top}\mathbf{K}_{h}\mathbf{V}_{\cdot
a}\right\|\cdot\left\|\frac{1}{n}\mathbf{Z}(\hat{J}_{n})^{\top}\mathbf{K}_{h}\mathbf{V}_{\cdot
b}\right\|_{\infty}$ $\displaystyle\leq$
$\displaystyle|\hat{J}_{n}|\sqrt{\Phi(|\hat{J}_{n}|,J_{n})}\cdot\left\|\frac{1}{n}\mathbf{Z}^{\top}\mathbf{K}_{h}\mathbf{V}_{\cdot
a}\right\|_{\infty}\cdot\left\|\frac{1}{n}\mathbf{Z}^{\top}\mathbf{K}_{h}\mathbf{V}_{\cdot
b}\right\|_{\infty}.$ (A.34)

By assumption we have $\Phi(|\hat{J}_{n}|,J_{n})=O_{P}(1)$ and we just have to
deal with the two infinity-norms. From Lemma A.6 we find that

$\displaystyle|\hat{J}_{n}|\left(\sup_{a=1,...,4}\left\|\frac{1}{n}\mathbf{Z}^{\top}\mathbf{K}_{h}\mathbf{V}_{\cdot
a}\right\|_{\infty}\right)^{2}=O_{P}\left(\frac{|\hat{J}_{n}|\log
p_{n}}{nh}\right)+O_{P}(|\hat{J}_{n}|h^{4})=o_{P}(1)$

by the assumptions on $|\hat{J}_{n}|$ and $h$. Thus, we have shown (A.32). In
order to show (A.33) we use very similar arguments: Firstly,

$\frac{1}{n}\mathbf{V}^{\top}\mathbf{K}_{h}^{\frac{1}{2}}\mathbf{M}_{n}(\hat{J}_{n})\mathbf{K}_{h}^{\frac{1}{2}}\mathbf{r}_{n}=\frac{1}{n}\mathbf{V}^{\top}\mathbf{K}_{h}\mathbf{r}_{n}-\frac{1}{n}\mathbf{V}^{\top}\mathbf{K}_{h}\mathbf{Z}(\hat{J}_{n})\left(\mathbf{Z}(\hat{J}_{n})^{\top}\mathbf{K}_{h}\mathbf{Z}(\hat{J}_{n})\right)^{-1}\mathbf{Z}(\hat{J}_{n})^{\top}\mathbf{K}_{h}\mathbf{r}_{n}.$

The first part equals exactly the first part of (A.33). Thus we have left to
prove that the second part above has the rate which appears in (A.33). This
can be seen by the same reasoning as in (A.34) and statement (A.28) from Lemma