Regression Discontinuity and Heteroskedasticity Robust Standard Errors: Evidence from a Fixed-Bandwidth Approximation

1 Introduction

Regression discontinuity (RD) designs have been propelled to the spotlight of economic analysis in recent years,^[1] especially in the policy and treatment evaluation literatures, as a form of estimating treatment effects in a non-experimental setting. In this design, treatment is assigned based on values of an observed characteristic, with the probability of receiving treatment jumping discontinuously at a known threshold. RD’s appeal stems from the relatively weak assumptions necessary for the nonparametric identification of treatment effects, specifically smoothness of the conditional expectation of outcome for treated and untreated individuals at the cutoff. This paper examines the effect of the additional assumptions that are explicitly and implicitly made to justify estimation and inference and shows that those assumptions can be weakened by using a heteroskedasticity-robust variance estimator, which is a common approach in empirical research but has not been supported by previous theoretical research.

Local polynomial estimators are the most common choice in empirical and theoretical work due to ease of implementation and flexibility. These kernel-based estimators rely on fitting a polynomial function to a range of the data, the size of which is determined by a bandwidth, h, just around the threshold. To perform inference in this setting, additional assumptions are usually imposed. In particular, asymptotic approximations assume smoothness of the conditional variance of the outcome, σ²(x). That assumption is somewhat restrictive and not necessary as it will be discussed below. Implementation of the local polynomial estimators also depends on the choice of the bandwidth, which can be benchmarked against several bandwidth selectors available in the literature.

The traditional nonparametric asymptotic approximations (“small-h”) provide formulas for the variance of such estimators as the bandwidth shrinks towards zero asymptotically. Those formulas reflect the fact that the distribution of the kernel-based estimator is asymptotically equivalent to the distribution with homoskedastic data when the smoothing parameter shrinks, which may have led researchers to overlook the need to develop robust inference in small samples. A similar point has been recently made by Kim, Sun, and Yang (2017) when discussing time series dependence. Moreover, in practice variance estimators that rely on homoskedasticity are discouraged even though they are supported by the theory. Hence, most researchers often use a different “pre-asymptotics” variance formula to compute standard errors, based on the White-Huber-Eicker robust standard errors or some nearest neighbor alternative.^[2] This paper provides an analysis of when and why the common practice of using White-Huber-Eicker robust standard errors in the RD context is appropriate.

I focus on presenting an alternative asymptotic approximation for the standard RD treatment effects estimator using a fixed bandwidth that relaxes the smoothness conditions imposed on σ²(x). This “fixed-h” approximation provides expressions for the estimator’s asymptotic variance that incorporate the bandwidth used by the researcher, and lead naturally to the standard error formulas used in the applied literature, clarifying the assumptions behind its use. The results show that the heteroskedasticity robust standard errors continue to provide valid inference under these relaxed assumptions. Furthermore, these standard errors are appropriate for any specific bandwidth chosen.

This paper shows that this form of the variance estimator is valid both under an asymptotic nesting where the bandwidth parameter, h, is fixed, and as under conventional RD asymptotics where h → 0. However, the traditional asymptotic variance estimator is not consistent for fixed bandwidths. Consequently the common empirical strategy of using robust standard errors in RD studies is appropriate especially when heteroskedasticity is a concern.

While I focus on a setting where the bias is “small” and does not affect the asymptotic distribution of the RD estimator, an alternative approximation for the estimator’s asymptotic bias that differs from the usual approximations in the literature is presented and provides additional intuition on the robustness-precision trade-off facing the researcher in RD designs. A discussion on how to address bias in practice is presented.

We first consider the procedure proposed by Calonico, Cattaneo, and Titiunik (2014) and implement an analytical bias correction reduction based on a first order expansion of the bias, reducing the order of the leading bias term. The framework proposed here is compatible with such bias reduction procedures. Furthermore, I explore a novel approach to bias mitigation based on a higher order correction inspired by the fixed bandwidth asymptotic bias approximations. That strategy is easily implemented by modifying the wild bootstrap RD procedure in Bartalotti, Calhoun, and He (2017). Simulation exercises attest that such approach can significantly improve test performance relative to first order based corrections.

This work contributes to the emerging literature on inference for treatment effects in the context of RD designs. Hahn, Todd, and van der Klaauw (1999 and 2001) and Lee (2008) presented the conditions for identification and estimation in RD designs. Porter (2003) provided results on the asymptotic properties of the estimators for the treatment effect of interest, obtaining limiting distributions for estimators based on local polynomial regression and partially linear estimation. Calonico, Cattaneo, and Titiunik (2014) studied asymptotic approximations for the bias-corrected local polynomial RD estimator as described above. Other studies about estimation, inference, and bandwidth choice in RD designs include Bartalotti and Brummet (2017), Calonico, Cattaneo, and Farrell (Forthcoming), Cattaneo, Frandsen, and Titiunik (2015), and Imbens and Kalyanaraman (2012) among others. McCrary (2008) studied specification testing. Finally, a broad review of the theoretical and applied literature, with emphasis on the identification of the parameter of interest and its potential interpretations, can be found in Imbens and Lemieux (2008) and Lee and Lemieux (2009).

The asymptotic framework using fixed bandwidths follows a growing literature that recognizes the potential for practical improvements in inference procedures in nonparametric methods. Notably, Neave (1970), in the context of spectral density estimation, obtained more accurate approximations to the variance of nonparametric spectral estimates. Neave’s work was later extended by Hashimzade and Vogelsang (2008). Similarly, Fan (1998) provided an alternative approximation for goodness-of-fit tests for density function estimates, obtaining improved approximations to the asymptotic behavior of the test and critical values for inference. More recent studies, on the context of time-series or spatial dependence, analyze and justify the use asymptotic variance formulas based on fixed-bandwidth approximations when pursuing heteroskedasticity, autocorrelation and spatial dependence robust inference (Bester et al. 2016; Chen, Liao, and Sun 2014; Kim, Sun, and Yang 2017; Sun 2014). This paper contributes to that literature in the context of local polynomial estimators and RD designs.

Section 5 provides simulation evidence that standard errors used by practitioners work well relative to the ones based in formulas suggested by the traditional approach (h → 0). The robust standard errors based on the fixed bandwidth formulas perform well for larger bandwidths, especially when heteroskedasticity is present. The intuition is that, for larger bandwidths, the homoskedasticity condition becomes less likely to hold as it would need to be valid over a larger support of the running variable. Finally, I provide an empirical application using Lee (2008), exemplifying with actual data the improvements obtained.

2 Model and Estimator

The interest lies in estimating the average treatment effect, τ, of a certain policy affecting part of a population of interest.^[3] There are two types of RD designs: sharp and fuzzy. They differ in regard to the assignment of treatment and to the impact of the discontinuity on the assignment process. At first, I focus on sharp RD designs. The discussion and extensions for the fuzzy design are presented afterwards.

In the sharp design, the treatment status, d, is a deterministic function of a so-called “running” variable, x, such that,

where x¯ is the known cut-off point. Then let Y₁ and Y₀ be the potential outcomes corresponding to the two possible treatment assignments. As usual, we cannot observe both potential outcomes, having access only to Y = dY₁ + (1 − d)Y₀. As described by Hahn, Todd, and van der Klaauw (2001) and Porter (2003), under some weak smoothness assumptions, the average treatment effect can be estimated by comparing points just above and just below the discontinuity. If the running variable is continuous, as well as the conditional expectations of Y₀ and Y₁, the discontinuity in treatment assignment at x¯ provides the opportunity to identify the average treatment effect at the cutoff without additional parametric functional form restrictions on the conditional expectations of the outcome variable. The average causal effect of the treatment at the discontinuity is

The sharp RD design uses the discontinuity in the conditional expectation of Y given X to uncover the average treatment effect at the cutoff. For a comprehensive review of RD designs and their applications and interpretation, see Lee and Lemieux (2009).

I focus on estimates of τ obtained using local polynomial estimation, which is the most common in applied work. The order p local polynomial estimator is defined as follows. In the sharp design case, given data (yi,xi)i=1,2,…,n, let di=1[xi≥x¯], k(⋅) be a kernel function, and h denote a bandwidth that controls the size of the local neighborhood to be averaged over. Also, define the p + 1 × 1 vector Z(x)=(1,(x−x¯h),(x−x¯h)2,…,(x−x¯h)p)′ and let(α^p+,β^p+)′ be the solution to the minimization problem:

while, similarly, (α^p−,β^p−) minimizes the same objective function, but with 1 − d_i replacing d_i. The estimator of the parameter of interest is given by

In the fuzzy design the probability of receiving treatment still changes discontinuously at the threshold, but is not required to go from 0 to 1,

This framework allows for a greater range of applications since it includes cases in which the incentives to receive treatment change discontinuously at the threshold, but are not strong enough to induce all individuals above it to be treated. The average treatment effect at the cutoff can be identified by the ratio of the change in the conditional expectation for the outcome variable to the change in the conditional probability of receiving treatment (Imbens and Lemieux 2008):

The estimator of the parameter of interest is given by the ratio

where α^_p is the estimator described above, and θ^ is the estimator for the change in the probability of being in the treated group at the cutoff which is obtained by using the same estimators with the treatment assignment variable, D_i, as the dependent variable.

In practice, the implementation of RD designs is not without challenges, even in simple cases, that require implicit assumptions or choices by the researcher not fully reflected in this basic set up. Importantly for the aims of this paper, is that several of these situations can be better understood by relying on an alternative asymptotic approximation based on a fixed bandwidth. That not only sheds light on the assumptions and issues dealt implicitly by the practitioners, but also naturally allows for the careful use tools in parametric analysis to approach these challenges.

To develop this framework, I consider the following assumptions, which are more flexible than the ones used in the literature on RD designs, e.g., Hahn, Todd, and van der Klaauw (2001) and Porter (2003), etc. Let F_o(x) denote the density of x, m(x)=E[y∣x], and x¯ is the cutoff that determines treatment. Finally, define ε=y−E[y∣X=x]=y−m(x) and σ2(x)=E[ε2∣X=x]., Importantly, we allow for the conditional variance, σ²(x), to be not continuous and non-smooth in a region around the cutoff.

Assumption 1 is very standard in the RD literature and does not restricts the support of the kernel being used in estimation significantly for most applications. Assumption 2 and Assumption 3 relax some of the usual smoothness conditions present in the literature. A relevant feature of the analysis performed in the next section is that fixed bandwidth asymptotics will allow the researcher to directly incorporate heteroskedasticity into the asymptotic distribution of the estimates of interest, including more complex dependence problems that are generally intractable in a shrinking bandwidth framework, e.g., clustering (Bartalotti and Brummet 2017). A similar point has been made by Kim, Sun, and Yang (2017) in the time series context. Most importantly, the conditional variance is allowed to not be continuous within the bandwidth around the cutoff. Hence, the framework developed here can be used to analyze cases in which the running variable might exhibit discreteness, which could be due to the nature of the running variable, as described by Lee and Card (2008), heaping (Barreca, Lindo, and Waddell 2016), measurement error (Bartalotti, Brummet, and Dieterle 2017 and Dong 2015) or some other characteristic of the application at hand.

3 Asymptotic Distributions

To derive the asymptotic properties of the estimator for τ, the usual regularity and smoothness conditions in the literature are sufficient for both the traditional and fixed-h approximations, e.g., Hahn, Todd, and van der Klaauw (2001) and Porter (2003), etc. As discussed in Section 2, we are able to relax some of the smoothness conditions by using a fixed bandwidth approach, which is particularly useful when dealing with scenarios in which the usual nonparametric approach is not feasible. In the following, let F_o(x) denote the density of x and m(x) denote the conditional expectation of y given x, i.e., m(x)=E[y∣x], and x¯ is the value of the running variable in which the discontinuity occurs. Finally, define ε=y−E[y∣X=x]=y−m(x) and σ2(x)=E[ε2∣X=x].

Using the traditional nonparametric asymptotic approximations, small-h, Hahn, Todd, and van der Klaauw (2001), Porter (2003), and Imbens and Lemieux (2008) among others obtain the following asymptotic variance formula for the (scaled) τ^,

where e1′Γ−1ΔΓ−1e1 is a constant scalar which depends only on the order of the local polynomial, p, and the kernel used.^[4] Crucially for practical purposes, it does not depend on the bandwidth (h) by construction. That is natural since the approximation is based on a shrinking bandwidth around the cutoff. However, that does not capture heteroskedasticity or changes in the values for the density of x in the range of data actually used in practice. That issue has been recognized in the literature, e.g. Imbens and Lemieux (2008) and Calonico, Cattaneo, and Titiunik (2014) and by practitioners who commonly implement the usual heteroskedasticity robust White-Huber-Eicker standard errors (or weighted analogues depending on the choice of kernel) when implementing RD based on intuitive or “pre-asymptotics” arguments.

Interestingly these heteroskedasticity robust standard errors can also be justified by an approximation under the asymptotic nesting where the bandwidth parameter, h, is fixed. The fixed-h asymptotic approximation can be summarized in the following result.

The fixed-h asymptotic variance formula is given by V_fixed-h, and depends directly on the behavior of σ²(x), f_o(x) inside the bandwidth around the cutoff, p, and the kernel function through (Γ+∗)−1Δ+∗(Γ+∗)−1 and (Γ−∗)−1Δ−∗(Γ−∗)−1. Hence, it captures the impact of h on the asymptotic variance accounting for potential local heteroskedasticity. This asymptotic approximation provides a natural and intuitive way to understand the sources of improvement in inference observed in practice when using these standard errors, as discussed below.

There are two noteworthy cases in which the formulas for the fixed-h asymptotic variance (and bias) simplify to the small-h formulas. First, if in the bandwidth around the cutoff, f_o(x), σ²(x) and m(x) are continuous, then when h → 0, ^[5]

Hence, if h is small, fixed-h and small-h provide similar approximations to the asymptotic behavior of α^_p. This is intuitive since the fixed bandwidth approximation is valid for any value of h, and the difference between the formulas arising from the behavior of σ²(x), f_o(x) inside the bandwidth will likely be small as the support of the running variable on which they are being evaluated becomes smaller.

Second, if f_o(x) and σ²(x) are constant around the cutoff, the fixed-h asymptotic variance formula simplifies to the small-h formula, i.e., V_fixed-h = V_small-h.^[6] This fact makes clear that the differences between the fixed bandwidth and traditional approximations are due to incorporating the behavior of f_o(x) and σ²(x) in the ranges around the cutoff in the first case, while considering only its values at the cutoff, f_o(x¯) and σ²(x¯) in the second.

4 Variance Estimators

To perform inference about α, appropriate estimates for the asymptotic variance formulas from Theorem 1 are necessary. The components of V_fixed-h, Δ+∗ and Γ+∗, have typical elements given by

and similarly for γj− and δj−. A natural estimator for the asymptotic variance is given by their sample analogues, using the estimated residuals ε^=y−m^(x),

which are consistent by standard arguments in both asymptotic frameworks.

Then, the natural plug-in estimator of the fixed-h variance-covariance matrix is given by

These are simple averages of the data and kernel weights, and they have the familiar “sandwich form” (Fan and Gijbels 1996). This estimator is analogous to the usual heteroskedasticity robust standard errors in a general weighted least squares framework, and it comes naturally from the fixed-h framework developed above. The fixed-h approach provides a framework that justifies the use of such estimators by practitioners. Intuitively, the variance estimators are “robust to the choice of bandwidths,” since they are valid for any finite h, take into consideration the impact of higher order polynomials on the estimator’s variance and are flexible regarding the conditional variance and density of X around the cutoff. Note that,[(Γ^+∗)−1Δ^+∗(Γ^+∗)−1+(Γ^−∗)−1Δ^−∗(Γ^−∗)−1]→p[(Γ+∗)−1Δ+∗(Γ+∗)−1+(Γ−∗)−1Δ−∗(Γ−∗)−1]=Vfixed-h.

In the rectangular kernel case, the variance estimator in equation (6) simplifies to the usual heteroskedastic robust variance estimator when using the data just above and below the cutoff.

Variance estimators based on small-h asymptotic variance formulas, as proposed by Porter (2003) and Lee and Lemieux (2009), are not fully robust to local heteroskedasticity. For example, Porter (2003) suggested an estimator for the variance of α^ based on the small-h approximation that requires only the estimation of the density of x andconditional variance of the errors at the cutoff. Let

and

is the estimator for the traditional asymptotic variance matrix. The matrix Γ−1ΔΓ−1 can be calculated directly because it is a deterministic function of the kernel and local polynomial order.

An additional drawback of the variance estimator in formula (10) is the need to estimate f_o(x¯), which is sidestepped if formula (6) is used. To obtain f^o(x¯), we need to choose a kernel and a bandwidth for the density estimator, increasing the number of tuning parameters to be chosen.^[9]

5 Simulations

This section presents simulation evidence displaying the empirical coverage of a standard t-statistic used to perform inference about the treatment effect of interest. The objective of the simulations is to determine the extent to which the heteroskedasticity robust variance estimator discussed in Section 4 improves on the small-h, homoskedastic variance estimator when the bias introduced by local misspecification is small. As shown previously, it is expected that both approaches yield similar test performance when the bandwidths are small, and differences in empirical coverage should be of greater importance when local heteroskedasticity is present around the cutoff.

To evaluate the relative performance of tests based on the fixed-h and small-h asymptotic variance approximations and their respective estimators, in this section I present evidence from two data generating processes for which the local polynomial estimator has negligible or mild asymptotic bias (in the sense that the bias does not overwhelm inference completely), in line with the condition imposed on Theorem 1. Obviously, if the bias in the local linear estimator is important, the inference on both approaches would suffer equally, since they use the same estimator.

Evidence from the simulations presented below indicates that, both in the theoretical (unfeasible) and feasible cases, inference using the heteroskedasticity robust standard errors performs well, especially for larger bandwidths and when local heteroskedasticity is present.

The simulations are based on a sharp RD design and consist of 2,000 replications with sample size n equal to 1,000 observations; the effective sample size included depends on the bandwidth used.^[10] For the DGPs used the running variable is drawn from X∼2Beta(2,4)−1, and σ(x) = 0.1295 for the homoskedastic case, σ(x)=0.1295+(x)2 for the first (mild) heteroskedastic case, and σ(x)=0.1295+(5x)2 for the second (acute) case. DGP 1 is linear in X and will serve as a benchmark of the refinements offered by the robust asymptotic approximation developed in Section 3. DGP 2 follows Calonico, Cattaneo, and Titiunik (2014) and is based on an empirical RD problem, corresponding to the regression function fitted to Lee’s (2008) data both above and below the cutoff. Since this DGP introduces bias for some bandwidths, it will allow us to compare the respective coverage obtained by both approaches in the presence of some bias.

The bandwidths used range from 0.1 to 1. While bandwidth selection algorithms have been proposed in the literature, in practice implementation requires the researcher to use a particular set of bandwidths. For example, for discrete running variables, a discretely positive bandwidth is necessary. Even with a continuous running variable, sample sizes often are small enough that concerns about precision impel researchers to use a relatively large bandwidth. As pointed out by Calonico, Cattaneo, and Titiunik (2014), bandwidth selectors in the literature “typically lead to bandwidth choices that are too large for the usual distributional assumptions to be valid.” Furthermore, as pointed out in several recent studies, these optimal bandwidth selectors could change significantly from setting to setting in a way not captured by the bandwidth selection algorithm [see Calonico et al. (2017) or Bartalotti, Calhoun, and He (2017) for cluster cases]. Hence, most researchers use data-driven bandwidth selectors to benchmark their results while reporting results for a set of chosen ad hoc bandwidths, which is adequate. The simulations presented here cover a wide range of fixed bandwidth values including those suggested by the selectors, being representative of the practice in the applied literature, highlighting the robustness of the fixed bandwidth approximation regardless of bandwidth choice.

In Appendix A.1 I present two additional DGPs which have high bias and compare two potential bias mitigation approaches, one based on small-h bias approximation in equation 4 (Calonico, Cattaneo, and Titiunik 2014) and a novel approach inspired by the fixed-h bias approximation in equation 3, which is based on a modified version of the bias correction bootstrap described in Bartalotti, Calhoun, and He (2017).

The empirical coverages presented are the fraction of rejections in the 2,000 repetitions for a two-sided test of nominal size 5%. More specifically, the DGPs are:

• DGP 1: yi=0.48+1.27xi+uiif x<0,0.52+0.84xi+uiif x≥0.

• DGP 2: yi=0.48+1.27xi+7.18xi2+20.21xi3+21.54xi4+7.33xi5+uiif x<0,0.52+0.84xi−3.00xi2+7.99xi3−9.01xi4+3.56xi5+uiif x≥0.

  The simulations use the local linear estimator (p = 1), since it is the preferred choice in applied work.^[11] The next subsection compares the test coverages obtained by the theoretical fixed-h and small-h asymptotic distributions derived in Section 3. Subsection 5.2 compares the empirical coverages obtained with (feasible) estimated standard errors.

6 Empirical Example

To illustrate the potential differences in the small-h and fixed-h approximations discussed above, this section uses data from Lee’s (2008) study of the electoral advantage of incumbency in the United States.^[13]Lee (2008) argues that the U.S. Congressional electoral system has a “built-in” RD design. That is, being the incumbent party in a congressional district is a deterministic function of the candidate-party’s vote share in the district during the last electoral cycle. This feature can be described in the following model:

where v_it is the democratic candidate’s vote share in congressional district i in election year t, w_it is a vector of characteristics or agents’ choices (potentially unobserved) as of election day on period t, and d_it indicates if the Democratic party is the incumbent in district i at election period t. We also assume that fi1(v∣w), the density of v_i1 conditional on w_i1, is continuous in v. The main issue in the analysis, as discussed in detail by Lee (2008), is that w_it is potentially unobserved and would likely be correlated with being incumbent in a certain district. For example, w_i1 would include party resources, demographic characteristics, and political leaning of districts, all of which could affect both the vote share in periods 1 and 2, thus biasing the estimates for the causal effect of incumbency.

The thought experiment to find the causal effect of incumbency would be to randomly allocate incumbency in districts to Democrats and Republicans while keeping all other characteristics constant. This clearly cannot be done, but by looking at closely contested elections, we can consider that the incumbents in those districts, whether Democrats or Republicans, were decided randomly, given the inherent uncertainty regarding the outcome of such events.^[14] This can provide reasonable estimates of the causal effect of incumbency in closely contested elections.

The data includes U.S. congressional election returns from 1946 to 1998, excluding the years ending in ‘0’ and ‘2’ due to the decennial redistricting which characterizes the U.S. congressional electoral system. The running variable is defined as the difference in vote share between the Democratic candidate and the strongest opponent. Hence, the Democrat wins the election when this variable crosses the 0 threshold – i.e. there is a positive difference in vote share, indicating that the Democrat has received more votes.

Table 3 shows the estimated advantage of incumbency and the estimated standard errors given by formulas based on a fixed bandwidth approximation (6), which in the rectangular kernel case is the usual heteroskedasticity robust s.e., and a small-h approximation (10), respectively.

Panel A presents estimates of the Nadaraya–Watson estimator (p = 0) for different bandwidths. Column 1 uses all the data available, column 2 looks only at elections for which the margin of victory in period t − 1 was within 50% of the total votes, and column 3 uses only elections with margins lower than 5%. Panel B presents similar estimates and standard errors obtained by a local linear estimator (p = 1), which is the preferred specification in several RD applications in the literature and is expected to significantly reduce bias in the estimates of the incumbency advantage effect.

Panel C presents the results Lee (2008) called the “parametric fit” in the first column, which uses a polynomial of order 4 to fit the whole data. The other two columns use smaller bandwidths to emphasize how the order of the polynomial chosen to fit the data can significantly impact estimates and standard errors, especially in small samples, by over-fitting the data.

The point estimates are exactly the same as presented by Lee (2008) for panel A and the first column in panel C. The remaining estimates are new. The results indicate a significant incumbent advantage in U.S. congressional races, even when comparing districts that had close elections in the previous electoral cycle for which the determination of incumbent status can be considered “as good as randomized.”^[15] See Table 5 in the Appendix A for similar results that include the pre-determined variables that Lee (2008) uses and the comparisons for fixed-h and small-h standard errors. The observed differences in pre-determined variables between incumbents and challengers vanish as we compare districts that previously had competitive races, lending credibility to the RD as an identification strategy for the incumbency effect. More relevant to this paper are the differences between the competing standard error estimates.

The estimated standard errors, shown in Table 3, differ significantly, with the fixed-h standard errors being smaller in most of the cases – as one would expect, given the simulations in Section 5. Also as expected, using smaller bandwidths comes at a large cost in terms of precision of the estimates due to the smaller amount of data available, negatively affecting both standard error estimators. The two last columns in all three panels show an increase of the estimated standard errors when the data is restricted to districts that had victory margins smaller than 5% in the previous election.

Interestingly, it is usual for the relative gap between the two standard errors estimates to decline as the bandwidth shrinks, as predicted in Section 3. For panels A and B and for the first two columns of panel C, this pattern is confirmed as we compare the percent difference between the standard errors within panels.^[16] For the third entry in panel C, note that both standard errors become larger than those for the wider bandwidth in column 2, but the relative gap increases. That can be due to the fact that the fixed-h standard error estimate requires the calculation of 4(2p + 1) terms (see equation 6). Hence, it is more susceptible to the combination of large polynomial orders and small sample sizes induced by the smaller bandwidth. Nevertheless, it still provides tighter confidence intervals than the small-h estimated standard error.

Finally, note that the differences between standard errors increase as the order of the polynomial used to fit the data increases for the same bandwidth. This is due to the fact that the small-h standard error estimator is based on a fixed scaling term for a given kernel and polynomial order choice, e1′Γ−1ΔΓ−1e1. Intuitively, as the polynomial order increases, the greater the distortion is likely to be between using this approximation and the more refined formula implied by the fixed-h approach.

7 Conclusion

This paper focuses on presenting an alternative asymptotic approximation for the standard RD treatment effects estimator using a fixed bandwidth. This approximation provides expressions for the estimator’s asymptotic variance that incorporate the bandwidth used by the researcher, and leads naturally to the standard error formulas used in the applied literature, i.e., White-Huber-Eicker robust standard errors when a rectangular kernel is implemented.

The proposed approximation is valid under an asymptotic framework where the bandwidth parameter, h, is fixed, as well as under conventional RD asymptotics where h → 0, while the traditional asymptotic approximations based on a shrinking bandwidth suggested, for example, in Imbens and Lemieux (2008) and Porter (2003) are not adequate for relatively larger bandwidths used in practice, especially in the presence of heteroskedasticity. Furthermore, the approximations provided in this study relax the assumption of continuity of the conditional variance of the outcome at the cutoff, providing a theoretical framework for that case. The common empirical strategy of using robust standard errors in RD studies is appropriate when the misspecification is small and heteroskedasticity is a concern.

Simulations document the theoretical refinements provided by relying on asymptotic variances based on fixed bandwidth asymptotics, and illustrate that tests using heteroskedasticity robust standard errors produces significantly more reliable test coverage than its counterparts using standard errors based on asymptotic approximations where h → 0.

Furthermore, an alternative, fixed-h, approximation for the estimator’s potential asymptotic bias is obtained and provides additional intuition on the robustness-precision trade-off facing the researcher in RD designs. I investigate two alternatives to address the bias in practice. The first is the widely used analytical approach based on the traditional asymptotics suggested by Calonico, Cattaneo, and Titiunik (2014). The second follows more naturally from the fixed bandwidth asymptotics developed here and is based on a iterative bootstrap approach that provides higher order bias mitigation by building upon the wild bootstrap algorithm in Bartalotti, Calhoun, and He (2017).