Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
64,761 characters · 13 sections · 53 citation commands
Inference in Regression Discontinuity Designs with High-Dimensional CovariatesFirst version: October 26, 2021. This version: . The authors gratefully acknowledge financial support by the European Research Council (ERC) through grant SH1-77202. Computations for this work were done in part using resources of the Leipzig University Computing Centre. Author contact information: Alexander Kreiß, Mathematical Institute, Leipzig University and Department of Statistics, London School of Economics. E-Mail: [email removed]. Christoph Rothe, Department of Economics, University of Mannheim. E-Mail: [email removed]. Website: http://www.christophrothe.net.
\onehalfspacing
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 hahn2001identification, IK12, calonico2014robust, AK18.
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 CCFT19. 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 T96 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 BCH14, 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.\footnote{This is a conceptual parallel between our setting and the use of covariates in randomized controlled experiments with a constant propensity score wager2016high. 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 BCH14 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 AK18, CCFT19, froelich2019including, noack2021flexible and arai2021regression. In particular, arai2021regression 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.\footnote{Roughly speaking, the conditions in arai2021regression 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 arai2021regression 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 su2019non, 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 (ref) we introduce our model and our proposed estimator, and give an informal description of its theoretical properties. In Section (ref) we give the assumptions for our theoretical analysis, and state and discuss our main result. Section (ref) 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.
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:
Local linear regression fan1996local 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:
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_iX_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 CCFT19, practitioners often augment the local regression in (ref) with additional covariates in a simple linear fashion, which yields the linear adjustment estimator
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:
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 (ref). 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.
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
Using arguments from CCFT19, 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 (ref) 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.
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 IK12, calonico2014robust or AK18. For $\lambda$, we compare different approaches in our simulations, and the results suggest that the plug-in procedure of BCH14 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 CCFT19 or AK18.
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$:
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:
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
with asymptotic bias and variance, respectively, such that
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 (ref) with the covariate-adjusted outcome $\widetilde Y_i$ replacing the original outcome $Y_i$.\footnote{Note, 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 CCFT19, the continuity of $\mu_Z$ from (ref) 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 (ref) 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 benjamini1995controlling. 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 (ref) are analogous to those obtained in CCFT19 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 AK18 with the generated data set. See Section (ref) for further discussion of implementation details.
We impose the following assumptions in our theoretical analysis.
\noindentAssumption (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\to0$ is a reference sequence such that $|J_n|\log p_n/\sqrt{ng}\to0$ and $|J_n|g^2\sqrt{\log p_n}\to0$.
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}\to0$ is a version of a standard assumption in the Lasso literature (cf. Chapter 6 in vdGB11), adapted to our locally penalized setup. The requirement that $|J_n|g^2\sqrt{\log p_n}\to0$ 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.
\noindentAssumption (AS): (Approximate Sparsity). It holds that $p_n\to\infty$ and, with $r_i(J_n,h)$ as in (ref), that
In addition, equation (ref) remains true with $h$ replaced by $b$.}
Assumption (AS) is similar in nature to the notion of approximate sparsity in, for example, BCH14.\footnote{Note 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 (ref) below for alternatives to (AS)). Online Appendix (ref) 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 (ref) 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 (ref) 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.
\noindentAssumption (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'$, $\mu_Z''$ and $\mu_Z'''$ 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
}
In a finite dimensional setting, like in CCFT19, 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=
^\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$.
\noindentAssumption (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
In addition, (ref) and (ref) also hold when $h$ is replaced by $b$.}
We call the requirements (ref) and (ref) in (TCS) equi-continuity from the right and left, respectively, and (ref) and (ref) 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.
\noindentAssumption (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^aK(u)du,\quad K_+^{(a)}=\int_0^{\infty}u^aK(u)du, \quad a\in\{0,1,2,3,4\}.$$
For the following assumptions, we need further preliminary definitions. Let
and for simplicity write $\mathbf{Z}=\mathbf{Z}(\{1,...,p_n\})$.
The constant $k(c,J)$ differs from the compatibility constants known from the classical Lasso literature (cf.\ Chapter 6.13 in vdGB11) 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 (ref) as follows:
where the infimum is taken over all pairs $(\theta,\gamma)$ for which (ref) 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)$.
Following the pattern of the compatibility constant in Definition (ref), we extend the concept of restricted sparse eigenvalues to localized problems. Continuing with the analogy, we can write down an equivalent formulation of (ref) in the fashion of (ref) 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 BCH14 for a discussion for non-localized estimators (Comment 3.2 therein) or Lemma 1 in BC13. The localized case with the additional kernel changes the problem to a conditional instead of an unconditional variance.
\noindentAssumption (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 (ref).
Recalling the discussion after Definitions (ref) and (ref), 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 (ref) 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 (ref) for the following assumption.
\noindentAssumption (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}$
and for any $k\in J_n\setminus J_{0,n}$
Moreover, $$\sqrt{|J_n\setminus J_{0,n}|}\cdot\frac{|J_n|\log p_n}{\sqrt{ng}}\to0\quad\textrm{and}\quad\sqrt{|J_n\setminus J_{0,n}|}\cdot|J_n|g^2\sqrt{\log p_n}\to0.$$}
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.\footnote{Note that (ref) and (ref) 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 (ref) holds for $k\in J_{0,n}(\mathcal{H})$ while for $k\in J_n\setminus J_{0,n}(\mathcal{H})$, (ref) 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 (ref) holds. Therefore, assumption (MS) is different from a $\beta$-min assumption vdGB11.
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)$$
\noindentAssumption (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}$
Equations (ref) and (ref) hold also when $b$ is replaced by $h$.}
Conditions (ref) and (ref) hold, for example, for bounded covariates, or covariates which fulfill a local sub-Gaussianity condition. Similarly, (ref) 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 (ref) and the reference there for a general statement and Lemma (ref) for a formulation of the statement which is tailored to our setting). Condition (ref) is a specific requirement for the model selection step and is thus formulated in terms of the model selection bandwidth $b$.
\noindentAssumption (CV): (Covariate Variance). It holds that
}
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.
The main result of this paper shows asymptotic normality of our estimator under the conditions stated in Section (ref). In order to state it, we define the following constants, which depend on the kernel function only:
The theorem's proof is given in Online Appendix (ref). The following remarks discuss its implications and possible extensions.
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 CCFT19 or AK18. We conduct our simulations using both frameworks, with results based on the methods in AK18 reported in this section, and results based on the methods in CCFT19 reported in Online Appendix (ref).
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 IK12, calonico2014robust or AK18 to select $b$; and we focus on the method proposed in AK18 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 BCH14, and a recently proposed bootstrap-based method by VL20. The three procedures are described formally in Online Appendix (ref), and we refer to them by the acronyms (CV), (BCH) and (LV), respectively, below. Our computations in this section use the {\tt R} packages {\tt glmnet} for implementing Lasso-based covariate selection, and {\tt RDHonest} for bandwidth selection, standard errors, and confidence intervals.
For the simulations, we consider the following DGP, which is a variation of “Model 2” in CCFT19 and corresponds to an RD setting with $p=200$ covariates and parameter of interest $\tau_Y=0.02$:
with $\sigma_{\varepsilon}^2=0.1295^2$, $\sigma_Z^2=0.1353^2$, $I_{200}$ denoting the $200\times200$ 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.
Our simulation results are summarized in Table (ref). 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.
In this section, we apply our methodology to data on Austrian workers from CCW07, 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 CCW07 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:
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 (ref). All four methods produce similar point estimates close to zero, which is in line with the results in CCW07. 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.
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 (ref) 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.