# Difference-in-Discontinuities: Estimation, Inference and Validity Tests

Pedro Picchetti
<EMAIL_ADDRESS>Department of Economics, Insper Institute of
Education and Research, Rua Quatá, 300, São Paulo, SP, Brazil 04546-042.
Cristine C. X. Pinto 11footnotemark: 1
<EMAIL_ADDRESS>Stéphanie T Shinoki11footnotemark: 1
<EMAIL_ADDRESS>
###### Abstract

This paper investigates the econometric theory behind the newly developed
difference-in-discontinuities design (DiDC). Despite its increasing use in
applied research, there are currently limited studies of its properties. The
method combines elements of regression discontinuity (RDD) and difference-in-
differences (DiD) designs, allowing researchers to eliminate the effects of
potential confounders at the discontinuity. We formalize the difference-in-
discontinuity theory by stating the identification assumptions and proposing a
nonparametric estimator, deriving its asymptotic properties and examining the
scenarios in which the DiDC has desirable bias properties when compared to the
standard RDD. We also provide comprehensive tests for one of the
identification assumption of the DiDC. Monte Carlo simulation studies show
that the estimators have good performance in finite samples. Finally, we
revisit Grembi et al. (2016), that studies the effects of relaxing fiscal
rules on public finance outcomes in Italian municipalities. The results show
that the proposed estimator exhibits substantially smaller confidence
intervals for the estimated effects.

Keywords: Difference-in-discontinuities, regression discontinuity design,
difference-in-difference

##  1 Introduction

The difference-in-discontinuities design (DiDC) aims to solve the limitations
of both regression discontinuity designs (RDD) and difference-in-difference
designs (DiD) by combining both the temporal and discontinuity-based sources
of variation from the data-generating processes (DGP). Grembi et al. (2016)
and Eggers et al. (2018) proposed this quasi-experimental approach, examining
the difference between pre- and post-treatment periods around a threshold for
both treated and not treated groups. Although its use in applied
microeconomics works is increasing (Azuaga and Sampaio, 2017; Chicoine, 2017;
García, 2020; Albright, 2023), research on the econometric theory of the DiDC,
including identification assumptions, estimation, and asymptotic properties,
remains limited.

The method offers more flexibility over the standard RDD in some settings and
can often be used in contexts where neither the RDD nor the DiD is applicable.
In particular, the DiDC can handle cases where control and treatment groups
differ significantly so as not to satisfy the parallel trends assumption of
the DiD, or when multiple time-invariant confounders exist at the threshold in
an RDD setting. Additionally, by incorporating more information into the
estimation, the DiDC eliminates the bias in the RDD estimates providing more
accurate and reliable estimates of treatment effects under certain
assumptions, including ones on the DGP’s functional form.

In this paper, we develop the econometric theory for the design with a
discontinuity-based approach to the differences. While Galindo-Silva et al.
(2021) develop an identification theory for the fuzzy difference-in-
discontinuity design based on the difference of RDD estimations, this work
will focus on the complementary approach of studying the discontinuity of the
differences.

Drawing from standard assumptions in cross-sectional RDD and assumptions
specific to the DiDC framework, we establish conditions for reliable
identification. These new assumptions, specific to the DiDC design, restrict
how the confounding effects behave around the threshold and through time.

We show that the parameter of interest can be recovered using local polynomial
estimation of the differences in the outcomes employing the methodology
proposed by Calonico et al. (2014). We derive the asymptotic properties of the
estimator, highlighting scenarios where the asymptotic bias of the DiDC can be
smaller than that of the RDD. Additionally, we introduce straightforward tests
to evaluate the testable implications of necessary assumptions for
identification. We explore the conditions under which DiDC mitigates bias more
effectively than the RDD offering practical guidance for unbiased treatment
effect estimation.

Through Monte Carlo simulations, we evaluate the finite-sample properties of
the estimator comparing its performance to that of the local linear RDD
estimator proposed by Calonico et al. (2014) and the nonparametric DiD
regression estimator proposed by Sant’Anna and Zhao (2020). By comparing
DiDC’s performance with DiD and RDD methods, we identify the scenarios where
DiDC outperforms RDDs and we determine the appropriate application contexts
for each method. Finally, we apply our estimator to the study by Grembi et al.
(2016) on the impact of fiscal rules on municipal deficits in Italy.

This work is closely connected to Calonico et al. (2014), whose contributions
have been instrumental to robust nonparametric RDD estimation and whose
estimation methods are used for our DiDC approach. Frölich and Sperlich (2019)
has a section exploring the possibilities of an intersection between RDD and
DiD however, the absence of a formal derivation of econometric properties
leaves a gap in understanding, with only preliminary discussions on
identification and estimation procedures. This work also aligns with the
broader literature on the intersection of RDD and panel data. Pettersson-
Lidbom (2012) explores settings where RDD is combined with fixed effects to
address small sample issues and violations of the continuous support
assumption. Lemieux and Milligan (2008), use a first-difference RD approach to
eliminate individual-specific fixed effects by capitalizing on the
longitudinal nature of the Finnish Census data. Lastly, Cellini et al. (2010)
introduce dynamic RD models, accommodating scenarios with multiple treatment
opportunities and looking into the dynamics of treatment effects.

The advantages of the DiDC have been explored across several applied
microeconomic studies. For instance, Azuaga and Sampaio (2017) uses the DiDC
to analyze the impact of introducing legislation on domestic violence in
Brazil. Butts (2021) explore the use of Diff-in-Disc in geographic settings,
Chicoine (2017) investigates the expiration of the Assault Weapon Ban (AWB) by
comparing municipalities in which the incumbent mayoral party wins a close
election with those where the incumbent is defeated and Albright (2023)
isolates the causal effects of algorithmic recommendations on decision-makers
using DiDC. These real-life examples demonstrate how effective this approach
can be in various research settings, highlighting its value as a practical
analysis method.

The remainder of this article is organized as follows. Section 2 presents the
DiDC as the discontinuity of differences in a potential outcomes model and
shows the main identification results, along with the necessary assumptions.
Estimation procedures for the treatment effect in the sharp setting are
presented in Section 3, along with the derivation of large-sample properties,
optimal bandwidths and robust confidence intervals. In Section 4 we present
tests for 2 of the identifying assumptions. Monte Carlo simulations are
conducted in Section 5 to examine the properties of the estimators. Section 6
provides an empirical illustration and Section 7 concludes. In the Appendix we
provide detailed notation, proofs and other methodological results.

##  2 Difference-in-Discontinuities as the discontinuity of differences

###  2.1 Setup

We define the potential outcomes as a function of both treatment assignment
and time. There are two time periods, $t=0$ and $t=1$, with treatment
determined by a threshold $z_{0}$ on the running variable $Z_{i}$ at $t=0$.
Treatment is introduced sometime between $t=0$ and $t=1$ so that no units are
treated by the treatment of interest in $t=0$. If there’s a pre-existing
confounder due to the discontinuity at $z_{0}$ (a confounding and time-
invariant discontinuity in the outcomes at $t=0$), we refer to it as the
confounding treatment. We indicate whether units were treated in period $t$
with $D_{i,t}=1$ and units that were not treated with $D_{i,t}=0$, where
$D_{i,0}=1$ implies treatment by the confounding policy and $D_{i,1}=1$
implies treatment by the actual treatment of interest. We assume that the
assignment rule and the running variables remain constant throughout both
periods.

We define four potential outcomes ($Y_{i,t}(D_{1},D_{0})$) for each unit at
each $t$, which depend on the treatment of interest and the confounding
treatment:

$\displaystyle Y_{i,t}(0,0)\quad Y_{i,t}(0,1)\quad Y_{i,t}(1,0)\quad
Y_{i,t}(1,1)$

We are interested in the local causal effect of the treatment of interest,
$D_{1}=1$, at $Z=z_{0}$. To estimate this we compare the difference in
outcomes between pre- and post-treatment periods for units near the threshold,
both for treated and untreated groups. This comparison can be approached
through two methods: the differences of RDDs and the RDD of the differences.
While prior studies have predominantly employed the former, this research opts
for the latter approach. This choice is because the RDD of the differences
approach offers a closer analogy to the well-established and robust RDD
literature. However, it is essential to highlight that both methods lead to
the same estimand, ultimately identifying the same underlying treatment
effects, as demonstrated in Claim 1 shown in the next section.

The observed outcome under the perspective of the difference of RDDs is:

$\displaystyle Y_{i,1}$
$\displaystyle=\left[Y_{i,1}\left(1,1\right)D_{0}+Y_{i,1}\left(0,1\right)\left(1-D_{0}\right)\right]D_{1}+\left[Y_{i,1}\left(1,0\right)D_{0}+Y_{i,1}\left(0,0\right)\left(1-D_{0}\right)\right]\left(1-D_{1}\right)$
(1) $\displaystyle Y_{i,0}$
$\displaystyle=Y_{i,0}\left(1,0\right)D_{0}+Y_{i,0}\left(0,0\right)\left(1-D_{0}\right)$
(2)

Where $Y_{i,0}$ takes this form because there is no treatment $D_{1}$ at time
$t=0$.

Meanwhile, the observed outcome under the approach of the RDD of the
differences between outcomes at times $t=0$ and $t=1$ is:

$\displaystyle\Delta Y_{i}$ $\displaystyle=\Delta Y_{i}(0)+\left[\Delta
Y_{i}(1)-\Delta Y_{i}(0)\right]D_{1}$ (3) $\displaystyle\Delta
Y_{i}\left(D_{1}\right)$ $\displaystyle=Y_{i,1}\left(D_{1}\right)-Y_{i,0}$
$\displaystyle Y_{i,1}\left(D_{1}\right)$
$\displaystyle=Y_{i,1}\left(1,D_{1}\right)+\left[Y_{i,1}\left(1,D_{1}\right)-Y_{i,1}\left(0,D_{1}\right)\right]D_{0}$
$\displaystyle Y_{i,0}$
$\displaystyle=Y_{i,0}\left(0,0\right)+\left[Y_{i,0}\left(1,0\right)-Y_{i,0}\left(0,0\right)\right]D_{0}$

###  2.2 Identification

Identification in difference-in-discontinuities design relies on the standard
assumptions from the cross-sectional RDD given by Hahn et al. (2001)
(Assumptions 1 and 2, adapted to include two time-periods), and assumptions
concerning the behavior of confounding effects over time (Assumptions 3 and
4). We state the assumptions here.

######  Assumption 1 (Continuity).

All potential outcomes are continuous in $Z=z_{0}$, for any
$D_{0},D_{1}\in\\{0,1\\}$ and $t\in\\{0,1\\}$:

$\displaystyle\lim_{\varepsilon\rightarrow
0}E[Y_{it}(D_{0},D_{1})|Z_{i}=z_{0}+\varepsilon]=\lim_{\varepsilon\rightarrow
0}E[Y_{it}(D_{0},D_{1})|Z_{i}=z_{0}-\varepsilon]$

######  Assumption 2 (Discontinuity of treatment probability).

The limits $D_{t}^{+}=\lim_{\epsilon\rightarrow
0}E(D_{it}|Z_{it}=z_{0}+\epsilon)$ and $D_{t}^{-}=\lim_{\epsilon\rightarrow
0}E(D_{it}|Z_{it}=z_{0}-\epsilon)$ exist and $D_{t}^{+}\neq D_{t}^{-}$.

Assumption 2 is a relevance condition, which needs to hold for identification
in both sharp and fuzzy settings. In the sharp settings, however, it
simplifies to $D_{1}^{+}=1$ and $D_{1}^{-}=0$, while $D_{0}^{+}=D_{0}^{-}=0$.

######  Claim 1.

Under Assumptions 1 and 2, the difference of RDDs, and the RDD of the
differences are the same

$\displaystyle\tau^{DiDC}=\tau_{1}^{RDD}-\tau_{0}^{RDD}=\lim_{\varepsilon\rightarrow
0}E[\Delta Y_{t}|Z_{it}=z_{0}+\epsilon]-\lim_{\varepsilon\rightarrow
0}E[\Delta Y_{t}|Z_{it}=z_{0}-\epsilon]$

###### Proof.

See Appendix C. ∎

To establish the identification of the treatment effect,
$E[Y_{i}(1)-Y_{i}(0)|Z_{i}=z_{0}]$ we invoke Assumptions 3 and 4, which are
specific to the DiDC design and concern the behavior of the confounding
effects around the threshold.

######  Assumption 3 (Time-invariance of confounding effects).

The treatment of interest is the only time-variant effect at the threshold
$z_{0}$:

$\displaystyle
E\left[Y_{i,0}(1,0)-Y_{i,0}(0,0)|Z_{i}=z_{0}\right]=E\left[Y_{i,1}(1,0)-Y_{i,1}(0,0)|Z_{i}=z_{0}\right]$

Assumption 3 assumption asserts that the effect of confounding policy $D_{0}$
when there is no treatment remains constant over time. That is, any
differences between the treatment and control groups at the threshold in
$t=1$, not caused by the treatment, should have existed in $t=0$ before the
treatment was introduced. The importance of this assumption cannot be
understated. Its validity directly affects the credibility of all subsequent
estimates and without it, the estimates from DiDC analyses are unreliable and
biased. In Section 4.1, we introduce a comprehensive test to assess its
validity.

######  Assumption 4 (Treatment effect is independent from the confounding
policy).

$\displaystyle
E\left[Y_{i,t}(1,1)-Y_{i,t}(1,0)|Z_{i}=z_{0}\right]=E\left[Y_{i,t}(0,1)-Y_{i,t}(0,0)|Z_{i}=z_{0}\right]=E\left[Y_{i,t}(1)-Y_{i,t}(0)|Z_{i}=z_{0}\right]$

This final assumption needed for identification states that the treatment
effect should not be affected by the confounding policy. This condition
implies that the total effect captured by an RDD at time $t=1$ would be the
sum of two effects: the effect of the confounding treatment and the effect of
the treatment of interest. That is, there is no complementarity or interaction
between the treatment of interest and the confounding treatment.

Under Assumptions 1 \- 4, we derive the Average Treatment Effect of treatment
$D_{1}$:

$\displaystyle\tau^{DiDC}$ $\displaystyle=\lim_{\varepsilon\rightarrow
0}E[\Delta Y_{t}|Z_{it}=z_{0}+\epsilon]-\lim_{\varepsilon\rightarrow
0}E[\Delta Y_{t}|Z_{it}=z_{0}-\epsilon]$
$\displaystyle=E\left[Y_{i}\left(1\right)-Y_{i}\left(0\right)|Z_{i}=z_{0}\right]\left(D_{1}^{+}-D_{1}^{-}\right)$

Therefore the estimand for the DiDC is:

$\displaystyle\hat{\tau}^{DiDC}$ $\displaystyle=\frac{\Delta Y^{+}-\Delta
Y^{-}}{\left(D_{1}^{+}-D_{1}^{-}\right)}$ (4)

####  2.2.1 Relaxing the Independence Assumption

Assumption 4 can be overly restrictive and sometimes there could be scenarios
where the confounding policy and the treatment of interest might interact.
Rather than assuming additive effects at time $t=1$, we can relax this
assumption with one that allows the inclusion of multiplicative effects:

######  Assumption 4’ (Multiplicative Effects).

$\displaystyle
E\left[Y_{i,t}(1,1)-Y_{i,t}(1,0)-Y_{i,t}(0,1)+Y_{i,t}(0,0)|Z_{i}=z_{0}\right]\left(D_{0}^{+}D_{1}^{+}-D_{0}^{-}D_{1}^{-}\right)$
$\displaystyle\quad=E\left[Y_{i,t}(0,1)-Y_{i,t}(0,0)|Z_{i}=z_{0}\right]\left(D_{1}^{+}-D_{1}^{-}\right)E\left[Y_{i,t}(1,0)-Y_{i,t}(0,0)|Z_{i}=z_{0}\right]\left(D_{0}^{+}-D_{0}^{-}\right)$

This assumption allows for the possibility that the combined effect of the
confounding policy and the treatment of interest is not simply the sum of
their individual effects but it also includes an additional component that
reflects how both treatments might interact. This additional component is
determined by multiplying the effects together. By incorporating this concept,
we can derive a new estimand that uses both the RDD at time 0 and the DiDC to
identify the causal effect of the treatment of interest:

$\displaystyle\frac{\tau^{DiDC}}{1+\tau_{0}^{RDD}}=E\left[Y_{i,t}(0,1)-Y_{i,t}(0,0)|Z_{i}=z_{0}\right]\left(D_{1}^{+}-D_{1}^{-}\right)$

This estimand can capture more complex relationships between the treatment and
the confounding factors. Our future work will focus on the properties and
robustness of this estimator.

##  3 Estimation and Inference for the Sharp DiDC

Following standard practice in the regression discontinuity literature, we
adopt local polynomial estimation to recover the parameter of interest. This
nonparametric method involves fitting a polynomial function to the data near
the threshold and using the estimated function to calculate the differences in
outcomes between the treatment and control groups at the threshold.

To implement the estimation at $z_{0}$, we lean on the methodologies proposed
by Calonico et al. (2014) for local polynomial estimation of the RDD
estimator. In our case, we estimate the local polynomial regression on the
differences in the outcomes with respect to time ($\Delta Y_{i}$).

###  3.1 Notation

Given a known threshold value $z_{0}$, which can be set to $z_{0}=0$ without
loss of generality, the observed value $Z_{i}$ determines treatment
assignment, with individuals assigned to the treatment group if $Z_{i}\geq 0$
and to the control group if $Z_{i}<0$ in the sharp setting.

The random variables $Y_{i,t}(1)$ and $Y_{i,t}(0)$ denote the potential
outcomes of unit $i$ with and without treatment, respectively, at period $t$.
The parameter of interest is
$E\left[Y_{i}\left(1\right)-Y_{i}\left(0\right)|Z_{i}=z_{0}\right]$, the
average treatment effect of treatment $D_{1}$ at the threshold. Under a mild
continuity condition, Hahn et al. (2001) showed that
$\tau^{RDD}=E[Y_{i}(1)-Y_{i}(0)|Z=z_{0}]$ is nonparametrically identifiable as
the difference of two conditional expectations evaluated at the (induced)
boundary point $z_{0}=0$. Similarly, the sharp DiDC parameter can be
identified as the difference of the difference (in time) of two conditional
expectations evaluated at $z_{0}=0$ at each side of $z_{0}$:

$\displaystyle\tau^{DIDC}$ $\displaystyle=\Delta\mu_{+}-\Delta\mu_{-},$
$\displaystyle\Delta\mu_{+}$
$\displaystyle=\raisebox{2.15277pt}{\scalebox{0.8}{$\displaystyle\lim_{z\rightarrow
0^{+}}\;$}}\Delta\mu(z),\quad\Delta\mu_{-}=\raisebox{2.15277pt}{\scalebox{0.8}{$\displaystyle\lim_{z\rightarrow
0^{-}}\;$}}\Delta\mu(z)$ $\displaystyle\Delta\mu(z)$
$\displaystyle=E\left[\Delta Y_{i}|Z_{i}=z\right]$

We employ Assumptions A.1 and A.2 in Appendix A for the non-parametric local
polynomial regression estimation. These assumptions impose restrictions on the
kernel function as well as impose the existence of moments, enforce continuity
of the running variable in the region, impose smoothness conditions on the
regression functions, and restrict the conditional variance of the observed
outcome.

###  3.2 Local Polynomial Estimator

For a given $\nu\leq p\in\mathbf{N}$, the general estimand of interest is
$\tau^{DIDC}_{\nu}=\Delta\mu_{+}-\Delta\mu_{-}$ with
$\Delta\mu^{(\nu)}_{+}=\nu!e^{\prime}_{\nu}\delta_{+,p}$,
$\Delta\mu^{(\nu)}_{-}=\nu!e^{\prime}_{\nu}\delta_{-,p}$ being the $\nu$th-
order derivatives of the $p$th-order local polynomial of the difference. The
$p$th-order local polynomial estimators of the $\nu$th-order derivatives
$\Delta\mu^{(\nu)}_{+,p}$ and $\Delta\mu^{(\nu)}_{-,p}$ are:

$\displaystyle\Delta\hat{\mu}^{(\nu)}_{+,p}(h_{n})$
$\displaystyle=\nu!e^{\prime}_{\nu}\hat{\delta}_{+,p}(h_{n})$
$\displaystyle\Delta\hat{\mu}^{(\nu)}_{-,p}(h_{n})$
$\displaystyle=\nu!e^{\prime}_{\nu}\hat{\delta}_{-,p}(h_{n})$
$\displaystyle\hat{\delta}_{\Delta Y+,p}(h_{n})$
$\displaystyle=arg\min_{\delta\in\mathbf{R}^{p+1}}\sum_{i=1}^{n}\mathbbm{1}(Z_{i}\geq
0)(\Delta Y_{i}-r_{p}(Z_{i})^{\prime}\delta)^{2}K_{h_{n}}(Z_{i})$
$\displaystyle\hat{\delta}_{\Delta Y-,p}(h_{n})$
$\displaystyle=arg\min_{\delta\in\mathbf{R}^{p+1}}\sum_{i=1}^{n}\mathbbm{1}(Z_{i}<0)(\Delta
Y_{i}-r_{p}(Z_{i})^{\prime}\delta)^{2}K_{h_{n}}(Z_{i})$

where $e_{\nu}$ is a conformable $(\nu+1)$ unit vector, $K_{h}(u)=K(u/h)/h$,
$h_{n}$ is a positive bandwidth sequence,
$r_{p}(x)=\begin{bmatrix}1&x&\ldots&x^{p}\end{bmatrix}^{\prime}$ and $\Delta
Y=\begin{bmatrix}\Delta Y_{1}&\Delta Y_{2}&\ldots&\Delta
Y_{n}\end{bmatrix}^{\prime}$. Therefore, for a positive bandwidth $h_{n}$, the
nonparametric estimator of $\tau_{\nu,p}$ is

$\displaystyle\hat{\tau}^{DiDC}_{\nu,p}(h_{n})=\Delta\hat{\mu}^{(\nu)}_{+,p}(h_{n})-\Delta\hat{\mu}^{(\nu)}_{-,p}(h_{n})$
(5)

####  3.2.1 Bandwidth choice

In order to perform the local polynomial estimation, it is necessary to choose
an appropriate bandwidth $h_{n}$. This parameter determines the range of
observations used for the estimation and impacts the trade-off between bias
and variance in the estimated treatment effect $(\tau^{DiDC})$. When
performing point estimation, the standard approach is to find the bandwidth
that minimizes the asymptotic Mean Squared Error (MSE) of the estimator. We
define the MSE as:

$\displaystyle
MSE_{\nu,p,s}(h_{n})=E\left[\left(\hat{\tau}_{\nu,p,s}(h_{n})-\tau_{\nu,p}\right)^{2}|\chi_{n}\right]$

######  Lemma 1.

Under Assumptions 1 and 2 with $S\geq p+1$, $\nu\leq p$, $h_{n}\rightarrow 0$
and $nh_{n}\rightarrow\infty$, the asymptotic MSE-optimal bandwidth is given
by:

$\displaystyle
h_{n}^{MSE}{}_{\nu,p}=\left(\frac{\left(1+2\nu\right)\mathbf{V}_{\nu,p}}{n2\left(1+p-\nu\right)\textbf{B}^{2}_{\nu,p,p+1,s}}\right)^{\frac{1}{2p+3}}$

where
$\mathbb{V_{\nu,p}}=\nu!\frac{\sigma_{+}^{2}-\sigma_{-}^{2}}{f}e^{\prime}_{\nu}\Gamma_{p}^{-1}\Psi_{p}\Gamma_{p}^{-1}$
and
$\mathbb{B_{\nu,p,p+1,s}}=\frac{\Delta\mu^{(p+1)}_{+}-(-1)^{\nu+p+s}\Delta\mu^{(p+1)}_{-}}{(p+1)!}\nu!e^{\prime}_{\nu}\Gamma_{p}^{-1}\varphi_{p,p+1}$,
provided that $\textbf{B}_{\nu,p,p+1,s}\neq 0$.

###### Proof.

in Appendix C. ∎

###  3.3 Inference

We now discuss the asymptotic properties of the estimator. Using Lemma B.4
found in the Appendix B, it is possible to recover the leading asymptotic
bias, expressed as:

$\displaystyle
E\left[\hat{\tau}^{DiDC}_{\nu,p}(h_{n})|\chi_{n}\right]-\tau_{\nu}$
$\displaystyle=h_{n}^{p+1-\nu}\mathbf{B}_{\nu,p,p+1}(h_{n})+h_{n}^{p+2-\nu}\mathbf{B}_{\nu,p,p+2}(h_{n})$
(6) $\displaystyle\mathbf{B}_{\nu,p,r}(h_{n})$
$\displaystyle=\frac{\Delta\mu_{+}^{(r)}\mathcal{B}_{+,\nu,p,r}(h_{n})-\Delta\mu_{-}^{(r)}\mathcal{B}_{-,\nu,p,r}(h_{n})}{r!}$

where
$\mathcal{B}_{+,\nu,p,r}=\nu!e_{\nu}^{\prime}\Gamma_{+,p}^{-1}(h_{n})\vartheta_{+,p,r}(h_{n})$
and
$\mathcal{B}_{-,\nu,p,r}=\nu!e_{\nu}^{\prime}\Gamma_{-,p}^{-1}(h_{n})\vartheta_{-,p,r}(h_{n})$
are asymptotically bounded. Further notation is available in Appendix A and a
detailed proof for the statement above can be found in Appendix B under Lemma
B.4.

This is where we believe is one of the main contributions of this research:
the asymptotic bias of the DiDC can be zero if we include an assumption
similar to that of parallel trends, or if the shapes of the data-generating
processes for both groups are time-invariant. In the first case, imposition of
parallel trends restricts the functional form in such a way that
$\Delta\mu_{+}^{(r)}\mathcal{B}_{+,\nu,p,r}(h_{n})=\Delta\mu_{-}^{(r)}\mathcal{B}_{-,\nu,p,r}(h_{n})$.
It’s important to note that the symmetry of the kernel function, as imposed in
Assumption A.2, plays a significant role in this result. In the case of time-
invariant data-generating processes, both $\Delta\mu_{+}^{(\nu)}$ and
$\Delta\mu_{-}^{(\nu)}$ equate to zero.

We present Claim 2, demonstrating the bias of difference-in-discontinuities in
relation to RDDs:

######  Claim 2.

The bias of $\hat{\tau}_{\nu,p}^{DiDC}$ can be decomposed as

$\displaystyle B\left[\hat{\tau}^{DiDC}_{\nu,p}(h_{n})\right]$
$\displaystyle=B\left[\hat{\tau}^{RDD}_{1,\nu,p}(h_{n})\right]-B\left[\hat{\tau}^{RDD}_{0,\nu,p}(h_{n})\right]$
(7)

where $B\left[\hat{\tau}^{RDD}_{1,\nu,p}(h_{n})\right]$ is the bias of the RDD
estimated at time $t=1$, after the intervention happened, and
$B\left[\hat{\tau}^{RDD}_{0,\nu,p}(h_{n})\right]$ is the bias of the RDD
estimated at time $t=0$, before the intervention happened.

###### Proof.

in Appendix C. ∎

Equation 7 demonstrates how incorporating additional data can help reduce the
bias in RDD analysis. Traditional RDD analysis of interventions typically uses
only solely on a cross-section of post-intervention data. However, by
incorporating pre-intervention data into the analysis, bias can be
substantially reduced and, under specific conditions, eliminated. Even in
cases where the bias is not completely eliminated, it can still be improved if
the RDD estimation of pre-treatment data exhibits bias in the same direction
as that of post-treatment RDD and if its magnitude is not too large to
overcome the original bias.

####  3.3.1 Bias Correction

Using the MSE-optimal bandwidth for point inference can result in bandwidths
that are too large, which may seem attractive for reducing variance but can
introduce first-order asymptotic bias111Calonico et al. (2014). To address
this issue, we employ robust bias-corrected confidence intervals (CIs)
proposed by Calonico et al. (2014). These intervals take into account the
asymptotic bias of the point estimate by 1) estimating the bias and
recentering the CI, and 2) incorporating the additional variance from
estimating the bias for bias correction into the CI. This process requires
estimating a separate local polynomial of order $q$, with $q>p\geq\nu$. The
bias-corrected estimator is defined as:

$\displaystyle\hat{\tau}^{bc}_{\nu,p,q}\left(h_{n},b_{n}\right)=$
$\displaystyle\tau_{\nu,p}\left(h_{n}\right)-h_{n}^{p+1-\nu}\hat{\mathbf{B}}_{\nu,p,q}\left(h_{n},b_{n}\right),$
$\displaystyle\hat{\mathbf{B}}_{\nu,p,q,s}\left(h_{n},b_{n}\right)=$
$\displaystyle\frac{\Delta\hat{\mu}^{(p+1)}_{+,q}(b_{n})}{(p+1)!}\mathcal{B}_{+,\nu,p,p+1}-\frac{(-1)^{\nu+p+s}\Delta\hat{\mu}^{(p+1)}_{-,q}(b_{n})}{(p+1)!}\mathcal{B}_{-,\nu,p,p+1}$

with $\Delta\hat{\mu}^{(p+1)}_{+,q}(b_{n})$ and
$\Delta\hat{\mu}^{(p+1)}_{-,q}(b_{n})$ being local polynomial estimations as
described in Appendix A. $\hat{\mathbf{B}}_{\nu,p,q}\left(h_{n},b_{n}\right)$
is the estimationg we get of the bias from the $q$th-order local polynomial.

To determine the MSE-optimal bandwidth for estimating the bias, we need to
conduct a separate local polynomial estimation of order $q$ with $q>p\geq\nu$.
Once again, we aim to minimize the MSE. Following Lemma 1 and applying it to
the bias estimate, we can find the MSE-optimal bandwidth for the bias
estimation.

$\displaystyle
b_{n}=h_{n}^{MSE}{}_{p+1,q,\nu+p+1}=\left(\frac{\left(3+2p\right)\mathbf{V}_{p+1,q}}{n2\left(q-p\right)\textbf{B}^{2}_{p+1,q,q+1,\nu+p+1}}\right)^{\frac{1}{2q+3}}$

####  3.3.2 Robust Confidence Interval

Following Calonico et al. (2014), we derive a large-sample distributional
approximation that accounts for the added variability introduced by the bias
estimate. The large-sample approximation for the standardized t-statistic is,
if $S\geq q+1$,
$n\min\\{h_{n}^{2p+3},b_{n}^{2p+3}\\}\times\max\\{h_{n}^{2},b_{n}^{2(q-p)}\\}\rightarrow\infty$,
then

$\displaystyle
T^{rbc}_{\nu,p,q}\left(h_{n},b_{n}\right)=\frac{\hat{\tau}^{bc}_{\nu,p,q}\left(h_{n},b_{n}\right)-\tau_{\nu}}{\sqrt{\mathbf{V}^{bc}_{\nu,p,q}\left(h_{n},b_{n}\right)}}\xrightarrow{d}\mathcal{N}\left(0,1\right)$

where $\mathbf{V}^{bc}_{\nu,p,q}\left(h_{n},b_{n}\right)$ is described in the
appendix C. This justifies the following confidence intervals:

$\displaystyle
CI_{\nu,p,q}\left(h_{n},b_{n}\right)=\left[\hat{\tau}^{bc}_{\nu,p,q}\left(h_{n},b_{n}\right)\pm\Phi^{-1}_{1-\alpha/2}\sqrt{\mathbf{V}^{bc}_{\nu,p,q}\left(h_{n},b_{n}\right)}\right]$

These confidence intervals have better properties compared to conventional
bias-corrected intervals. They are more robust to bandwidth selection, feature
coverage error decays at a faster rate, and offer shorter interval lengths, as
explained by Calonico et al. (2014).

##  4 Validity Test

###  4.1 Testing if Confounding Effect is Time-Invariant

As highlighted in Section 2, the importance of Assumption 3 \- that the
confounding effect is constant over time - cannot be overstated. Without this
assumption, all estimations become inherently flawed, resulting in biased
treatment effect estimates. Therefore, an initial step when considering the
applicability of a differences-in-discontinuity approach is to assess the
validity of this assumption.

To accomplish this, we propose a simple test to check violations of this
assumption. In essence, it involves using the available pre-treatment periods
to conduct an estimation of stacked RDDs. The goal is to examine whether the
RDD coefficients remain consistent across multiple periods. This involves a
set of $k$ periods where it is known that no changes occurred at the
threshold.

It is crucial to note that this estimation method is not suitable if any
alterations occurred at the threshold between the initial period in this
sample and the period just preceding the introduction of the treatment of
interest. To ensure the validity of this approach, it is necessary to select a
period during which no events occurred at this threshold that could influence
the observed outcome.

The procedure involves considering the following stacked RDDs regression
model:

$\displaystyle
Y_{i}=\alpha_{i}+\sum_{k=1}^{K}\left[\beta_{-k}\left(Z_{i}-z_{0}\right)+\theta_{-k}D\left(Z_{i}-z_{0}\right)\right]T_{-k}+\varepsilon_{i}$
(8)

where $T_{-k}$ is a dummy indicating the period to which the data belong, that
is, the RDD of which period. The hypothesis to be tested is:

$\displaystyle H_{0}:\quad\theta_{0}=\theta_{-1}=...=\theta_{-K}$

A Wald test is sufficient to test these hypotheses. Rejection of the null
hypothesis indicates that the differences-in-discontinuity design may not be
suitable for estimating the effect of treatment in this setting.

A question arises about concerning bandwidth to use when conducting this test.
There are several options available, and Appendix D.1 provides details on the
simulations used to assess the better bandwidth. In summary, in practice, any
bandwidth derived from the data that is optimal for the specific RDD should
work well.

###  4.2 Testing the Time-Invariance of Data-Generating Functional Form

In Section 3.3, we briefly discussed an important aspect of the DiDC method:
its potential to achieve zero asymptotic bias when choosing DiDC over the
standard RDD, provided that the functional forms of the data-generating
processes for both groups are time-invariant. Essentially, this would mean
that the derivatives of the functions on each side of the threshold at each
$t=\\{0,1\\}$ ($\mu_{+,1}(z)$ and $\mu_{+,0}(z)$, $\mu_{-,1}(z)$ and
$\mu_{-,0}(z)$) are equal at each point of the running variable $z$, resulting
in zero bias in the local polynomial estimation.

######  Claim 3.

If the data-generating processes for both groups (above and below the
threshold) are time-invariant, then $\Delta\mu_{+}^{(\nu)}$ and
$\Delta\mu_{-}^{(\nu)}$ will be zero.

When $\Delta\mu_{+}^{(\nu)}=\Delta\mu_{-}^{(\nu)}=0$, the bias term
$\mathbf{B}_{\nu,p,r}(h_{n})=\frac{\Delta\mu_{+}^{(r)}\mathcal{B}_{+,\nu,p,r}(h_{n})-\Delta\mu_{-}^{(r)}\mathcal{B}_{-,\nu,p,r}(h_{n})}{r!}$
is also zero. This condition implies that $\Delta\mu(z)$ is a constant,
allowing unbiased point estimates through linear estimations on both sides of
the threshold. As a result, two simple OLS regressions would be sufficient for
accurately estimating the treatment effect.

Therefore, it can be very beneficial for researchers to know if the shapes of
the data-generating functions remain stable over time. This knowledge would
allow them to determine if they can use simpler and less biased estimation
methods. To test this, we propose comparing the derivatives of the estimated
functions on each side of the threshold. We can estimate $\mu_{+,t}(z)$ and
$\mu_{-,t}(z)$ nonparametrically, for $t=\\{0,1\\}$. This is, we locally
approximate $\mu_{t}(z)$ by a polynomial of order $p$ for z in the
neighborhood of $\bar{z}$:

$\displaystyle\mu_{t}(z)\sim\mu_{t}(\bar{z})+\mu^{(1)}_{t}(\bar{z})\left(z-\bar{z}\right)+...+\frac{\mu^{(p)}_{t}(\bar{z})}{p!}\left(z-\bar{z}\right)^{p}$

We perform a local polynomial regression fit for each side of the threshold
and each period $t={0,1}$:

$\displaystyle\hat{\beta}_{+,t}$
$\displaystyle=arg\min_{\beta\in\mathbf{R}^{p+1}}\sum_{i=1}^{n}\mathbbm{1}(Z_{i}\geq
0)(Y_{i,t}-r_{p}(Z_{i})^{\prime}\beta)^{2}K_{h_{n}}(Z_{i}-\bar{z})$
$\displaystyle\hat{\beta}_{-,t}$
$\displaystyle=arg\min_{\beta\in\mathbf{R}^{p+1}}\sum_{i=1}^{n}\mathbbm{1}(Z_{i}<0)(Y_{i,t}-r_{p}(Z_{i})^{\prime}\beta)^{2}K_{h_{n}}(Z_{i}-\bar{z})$

where
$\beta_{+,p,t}=\begin{bmatrix}\mu_{+,t}&\mu_{+,t}^{(1)}&\ldots&\mu_{+,t}^{(p)}/p!\end{bmatrix}^{\prime}$,
and
$\beta_{-,p,t}=\begin{bmatrix}\mu_{-,t}&\mu_{-,t}^{(1)}&\ldots&\mu_{-,t}^{(p)}/p!\end{bmatrix}^{\prime}$.
The estimation procedure is similar to the one described in Section 3. We are
examining whether the derivatives of the functions, rather than the levels,
are consistent between periods $t=0$ and $t=1$ on each side. To do this, we
create a test statistic inspired by the "two-sample" Kolmogorov-Smirnov test:

$\displaystyle
D_{+}=\sup(\mathbf{1}_{p+1}-e_{1})\left|\left|\beta_{+,p,1}-\beta_{+,p,0}\right|\right|$
$\displaystyle
D_{-}=\sup(\mathbf{1}_{p+1}-e_{1})\left|\left|\beta_{-,p,1}-\beta_{-,p,0}\right|\right|$

The critical values of the test statistic depend on the choice of $p$, so
determining the optimal $p$ is on the agenda.

The critical values of the test statistic depend on the choice of $p$.
Selecting the right value of $p$ can significantly affect the outcome of the
test and the conclusions drawn from the results. Therefore, determining the
optimal $p$ is on the agenda.

##  5 Monte Carlo Simulations

We now demonstrate the finite-sample properties of the estimator through Monte
Carlo simulations and compare the performance of the proposed estimator to
that of the local linear RDD estimator proposed by Calonico et al. (2014) and
the nonparametric DiD regression estimator proposed by Sant’Anna and Zhao
(2020). We consider Data Generating Processes (DGP) based on model 3 from
Calonico et al. (2014), with small modifications that will be described.

We conduct our simulation studies in four distinct settings: with time-
invariant confounding factors at the threshold, and without any confounding
factors, which is the typical scenario for RDDs, for both time-invariant DGPs
and time-changing DGPs. For each simulation, we conduct 1000 replications and
for each replication, we consider a sample size $n=500$, with
$Z_{i}\sim(2\mathcal{B}(2,4)-1)$ where $\mathcal{B}(p_{1},p_{2})$ is a beta
distribution with parameters $p_{1}$ and $p_{2}$. We also consider
$\varepsilon_{it}\sim N(0,\sigma^{2}_{\varepsilon})$,
$\sigma_{\varepsilon}=0.1295$ and the outome generated is
$Y_{i,t}=\mu_{it}(Z_{i})+\varepsilon_{i,t}$. The detailed specifications and
further functional forms of both models can be found in the Appendix D for
added clarity. In both scenarios, we observe that our DiDC estimator performs
better than the RDD estimator, with a smaller bias and improved coverage.
Additionally, we find that the DiDC estimator has a similar bias and improved
coverage when compared to the DiD estimator if the DGPs change over time.

###  5.1 Identical DGPs over time

In the first simulation, we mimic a scenario where one or more confounding
discontinuities are present at the threshold $z_{0}=0$ but the DGPs are the
same in both periods. This setting is comparable to that of Grembi et al.
(2016), where the treatment of interest was introduced at some point between
$t=0$ and $t=1$ and was given to units whose running variable values $Z_{i}$
are above the threshold $z_{0}=0$ however there were other pre-existing
treatments determined by the same threshold $z_{0}=0$ on the same running
variable.

The second model we consider is a scenario with identical DGPs over time, the
only distinction being that in period $t=1$, there is an effect of treatment
$\tau$ for units with $Z_{i}\geq 0$. Importantly, there are no other sources
of discontinuity in the outcome at the threshold $z_{0}=0$, making this what
would be considered an ideal scenario for estimating the effect of the
treatment with a regular RDD estimation.

Figure 1: DGPs for models 1 and 2

We present the results of these simulation studies in Table 1 and Table 2. We
estimate the DiD along with the RDD and DiD and compare the average bias,
median bias, root-mean-squared errors, 95% coverage probability, and 95%
confidence interval length for each estimator when the treatment effect $\tau$
is equal to 1.

Table 1: Monte Carlo Simulations - Identical DGPs - Yes Confounder

| Av. Bias | RMSE | Coverage | CIL | Av. SE | Bandwiths |
---|---|---|---|---|---|---|---
| $h_{n}$ | $b_{n}$ |
RDD (Robust) | 1.042 | 1.042 | 0 | 0.315 | 0.080 | 0.183 | 0.328 |
DiDC (RDD of $\Delta$s) | 0.001 | 0.061 | 0.928 | 0.278 | 0.070 | 0.215 | 0.348 |
DiDC ($\Delta$ of RDDs) | 0.001 | 0.068 | 0.924 | 0.301 | 0.077 | 0.183 | 0.328 |
DiD | 0.001 | 0.017 | 1 | 0.559 | 0.143 | - | - |
OLS of $\Delta$s | 0.001 | 0.023 | 1 | 0.144 | 0.258 | - | - |

Note: Simulations based on 10,000 Monte Carlo experiments with a sample size
$n=500$. RDD is the non-bias corrected RD estimator from CCT (2014), DiD is
the outcome regression DiD estimator from Sant’Anna and Zhao (2020), DiDC is
the estimator proposed in this paper. “Av. Bias”, “RMSE”, “Cover”, “CIL’ and
"Av. SE, stand for the average simulated bias, simulated root mean-squared
errors, 95% coverage probability, 95% confidence interval length and average
standard error, respectively. "Bandwidths" $h_{n}$ and $b_{n}$ report the
plug-in bandwidths for point and bias estimation, respectively.

Table 1 shows that there is a significant improvement over the RDD estimation
when there is a confounding factor at the threshold. The reason for this is
that the RDD method relies solely on post-treatment data and fails to account
for pre-existing confounding factors, making it unsuitable for this setting.
In comparison to the DiD, our estimator has closely the same bias, but the
coverage seems to be more desirable since it is more informative. While the
coverage of DiD is 1, meaning it covers the entire interval of the data, the
coverage of the DiDC for the 95% interval is very close to 95%, at 92.8%.

L[1]>m#1 Z[1]>m#1 C[1]> m#1

Table 2: Monte Carlo Simulations - Identical DGPs - No Confounder

| Av. Bias | RMSE | Coverage | CIL | Av. SE | Bandwiths |
---|---|---|---|---|---|---|---
| $h_{n}$ | $b_{n}$ |
RDD (Robust) | 0.044 | 0.073 | 0.903 | 0.314 | 0.080 | 0.183 | 0.328 |
DiDC (RDD of $\Delta$s) | 0.001 | 0.060 | 0.924 | 0.278 | 0.070 | 0.215 | 0.348 |
DiDC ($\Delta$ of RDDs) | 0.001 | 0.067 | 0.919 | 0.301 | 0.077 | 0.183 | 0.328 |
DiD | 0.001 | 0.017 | 1 | 0.559 | 0.143 | - | - |
OLS of $\Delta$s | 0.001 | 0.023 | 1 | 0.144 | 0.258 | - | - |

Note: Simulations based on 10,000 Monte Carlo experiments with a sample size
$n=500$. RDD is the non-bias corrected RD estimator from CCT (2014), DiD is
the outcome regression DiD estimator from Sant’Anna and Zhao (2020), DiDC is
the estimator proposed in this paper. “Av. Bias”, “RMSE”, “Cover”, “CIL’ and
"Av. SE, stand for the average simulated bias, simulated root mean-squared
errors, 95% coverage probability, 95% confidence interval length and average
standard error, respectively. "Bandwidths" $h_{n}$ and $b_{n}$ report the
plug-in bandwidths for point and bias estimation, respectively.

Table 2 reveals that once again our estimator has a similar bias to that of
the DiD, but with more informative coverage. However, the key takeaway is that
this table provides a valuable comparison between the DiDC and the RDD. In
this simulation, as there were no other discontinuities at the threshold
besides the treatment assignment probabilities, the scenario would be
considered for the traditional use of RDD.

The results show that difference-in-discontinuities estimation has an apparent
improvement over the RDD, yielding a smaller bias and better coverage. This is
noteworthy as it highlights that even in cases where the RDD would
traditionally be regarded as suitable, the differences-in-discontinuity
approach can give more desirable results by incorporating more data into the
estimation. The results hold when the functional format exhibits little
temporal variation, as shown in the next section, with the additional period
contributing to more reliable estimations compared to RDD, whenever the bias
from the extra period $t=0$ does not surpass the one from period $t=1$.

###  5.2 Time-varying DGPs

Next, we introduce scenarios where the Data Generating Processes (DGPs) change
between periods, contrasting to the prior section where DGPs were identical
across time. This setup allows us to assess the performance of the difference-
in-discontinuities in scenarios with small temporal DGP variations.

Figure 2: DGPs for models 3 and 4

We again do a version with confounder and a version without confounder, but
now we alter the model for the period $t=1$ so that it is a linear model
derived from the original. The DGP for time $t=0$ is identical to that of
models 1 and 2. Results are shown in Table 3.

L[1]>m#1 Z[1]>m#1 C[1]> m#1

Table 3: Monte Carlo Simulations - Different DGPs

| Av. Bias | RMSE | Coverage | CIL | Av. SE | Bandwiths |
---|---|---|---|---|---|---|---
| $h_{n}$ | $b_{n}$ |
Yes confounder |  |  |  |  |  |  |  |
RDD (Robust) | 1.040 | 1.040 | 0 | 0.289 | 0.074 | 0.215 | 0.348 |
DiDC (RDD of $\Delta$s) | 0.016 | 0.067 | 0.906 | 0.290 | 0.074 | 0.195 | 0.337 |
DiDC ($\Delta$ of RDDs) | 0.014 | 0.066 | 0.914 | 0.290 | 0.074 | 0.214 | 0.348 |
DiD | -1.568 | 1.568 | 0 | 0.384 | 0.098 | - | - |
OLS of $\Delta$s | 0.684 | 0.684 | 0 | 0.244 | 0.363 | - | - |
No confounder |  |  |  |  |  |  |  |
RDD (Robust) | 0.040 | 0.069 | 0.899 | 0.291 | 0.074 | 0.215 | 0.348 |
DiDC (RDD of $\Delta$s) | 0.016 | 0.067 | 0.905 | 0.290 | 0.074 | 0.195 | 0.336 |
DiDC ($\Delta$ of RDDs) | 0.014 | 0.066 | 0.914 | 0.290 | 0.074 | 0.214 | 0.348 |
DiD | -1.568 | 1.568 | 0 | 0.384 | 0.098 | - | - |
OLS of $\Delta$s | 0.684 | 0.684 | 0 | 0.244 | 0.363 | - | - |

Note: Simulations based on 10,000 Monte Carlo experiments with a sample size
$n=500$. RDD is the non-bias corrected RD estimator from CCT (2014), DiD is
the outcome regression DiD estimator from Sant’Anna and Zhao (2020), DiDC is
the estimator proposed in this paper. “Av. Bias”, “RMSE”, “Cover”, “CIL’ and
"Av. SE, stand for the average simulated bias, simulated root mean-squared
errors, 95% coverage probability, 95% confidence interval length and average
standard error, respectively. "Bandwidths" $h_{n}$ and $b_{n}$ report the
plug-in bandwidths for point and bias estimation, respectively.

This table shows that in the scenario without confounding factors and time-
varying DGPs, the DiD cannot be used due to the lack of parallel trends.
However, the RDD can be employed, but it shows a higher bias compared to the
differences-in-discontinuity approach because the GDPs, although different,
are similar enough.

With confounding factors, both the RDD and DiD estimators are unsuitable. RDD
remains problematic due to the confounding effects as it does not account for
pre-existing discontinuities. DiD is also unsuitable because no parallel
trends exist in this context. The differences-in-discontinuities method is the
only approach that correctly estimates the treatment effect.

##  6 Empirical Illustration

We demonstrate the application of our estimator by revisiting Grembi et al.
(2016), which uses data from Italian municipalities to analyze the impact of
fiscal rules on municipal deficit by exploiting a government fiscal rule
relaxation as a natural experiment. Our objective is to implement the DiDC
approach and compare our results with the original findings. Additionally, we
conduct the validity test from Section 4.1 to assess the suitability of using
DiDC in this context.

The dataset comprises data from Italian municipalities, focusing on the period
surrounding the government’s relaxation of fiscal rules in 2001.
Municipalities with fewer than 5,000 inhabitants experienced a relaxation of
fiscal rules and are defined as treatment, while those with more than 5,000
inhabitants served as controls.

We reframe the data into a more manageable 2x2 design, simplifying the
estimation. This transformation involves computing the differences in outcomes
before and after the policy change. Grembi et al. (2016)’s original study
utilized a large panel dataset, a rectangular kernel and a polynomial of
degree one as in the model below:

$\displaystyle Y_{it}=$
$\displaystyle\delta_{0}+\delta_{1}(X_{i}-c)+S_{i}(\gamma_{0}+\gamma_{1}(X_{i}-c))$
$\displaystyle
T_{i}[\alpha_{0}+\alpha_{1}(X_{i}-c)+S_{i}(\beta_{0}+\beta_{1}(X_{i}-c))]+\varepsilon_{it}.$

where $S_{i}$ is a dummy variable for cities below 5,000 (treatment
indicator), $T_{i}$ is a dummy variable for the post-treatment period and
$\beta_{0}$ is the parameter of interest.

For our implementation, we employ the methodology described in Section 3.
Table 4 compares the results from both our DiDC methods (RDD of the
differences and the differences of RDDs) to Grembi et al. (2016)’s
estimations. It also reports bandwidths for point estimation and number of
observations.

Table 4: Effects of relaxing a fiscal rule - CCT (2014) Bandwidths Estimators | Deficit | Fiscal Gap | Taxes
---|---|---|---
Diif-in-Disc (Grembi et al, 2016) |  | 17.495
---
(7.737)
| 56.468
---
(32.079)
| -76.083
---
(32.597)
Bandwidth | 600 | 513 | 378
Observations | 2414 | 2136 | 1536
Diff-in-Disc (RDD of $\Delta$) |  | 23.100
---
(12.763)
| 44.508
---
(16.414)
| -21.645
---
(11.690)
Bandwidth | 619 | 600 | 420
Observations | 1245 | 1028 | 1030
Diff-in-Disc ($\Delta$ of RDDs) |  | 26.864
---
(4.971)
| 56.043
---
(10.850)
| -18.876
---
(10.850)
Bandwidth | 995 | 809 | 819
Observations | 1456 | 1161 | 1164

We conduct the test from Section 4.1 to evaluate the assumption of time-
invariant confounding effects. We use data from 1997 to 2000 to perform
"stacked" RDD regressions. The results are reported in Table 5.

Table 5:  | Dependent variable:
---|---
| imposte
Z:$T_{1997}$ |  $-$0.040
| (0.071)
Z:$T_{1998}$ | 0.265∗∗∗
| (0.069)
Z:$T_{1999}$ | 0.176∗∗
| (0.069)
Z:$T_{2000}$ | 0.054
| (0.068)
Z:$T_{1997}$:$D_{1997}$ | 0.280∗∗∗
| (0.108)
Z:$T_{1998}$:$D_{1998}$ |  $-$0.369∗∗∗
| (0.103)
Z:$T_{1999}$:$D_{1999}$ |  $-$0.213∗∗
| (0.103)
Z:$T_{2000}$:$D_{2000}$ | 0.038
| (0.102)
Constant | 121.487∗∗∗
| (3.420)
Observations | 4,810
R2 | 0.059
Adjusted R2 | 0.048
Residual Std. Error | 2.137 (df = 707)
F Statistic | 5.518∗∗∗ (df = 8; 707)
Note: |  ∗p$<$0.1; ∗∗p$<$0.05; ∗∗∗p$<$0.01

The hypothesis test for
$Z:T_{1997}:D_{1997}=Z:T_{1998}:D_{1998}=Z:T_{1999}:D_{1999}=Z:T_{2000}:D_{2000}$
indicates a rejection of the null hypothesis, suggesting that the coefficients
are not constant over time, as reported in Table 6:

Table 6:  Statistic | N | Mean | St. Dev. | Min | Max
---|---|---|---|---|---
Res.Df | 2 | 708.500 | 2.121 | 707 | 710
RSS | 2 | 3,298.776 | 98.196 | 3,229.341 | 3,368.211
Df | 1 | 3.000 |  | 3 | 3
Sum of Sq | 1 | 138.871 |  | 138.871 | 138.871
F | 1 | 10.134 |  | 10.134 | 10.134
Pr(>F) | 1 | 0.00000 |  | 0.00000 | 0.00000

In summary, our estimated treatment effects are similar in magnitude to the
difference-in-discontinuities approach utilized in Grembi et al. (2016) study.
Notably, our difference of RDDs approach yields smaller confidence intervals,
suggesting it to be the most powerful approach. However, our validity test
reveals that confounding effects vary over time, so the use of DiDC in this
setting appears to be leading to biased estimates. We also replicate our
estimates using Ludwig and Miller (2007)’s bandwidths with similar results.
These estimates are shown in Table E.1 in Appendix E.

##  7 Conclusion

The difference-in-discontinuities (DiDC) design is emerging as a promising
method for estimating causal inference, addressing the limitations of both
regression discontinuity (RDD) and difference-in-difference (DiD) approaches.
This paper lays the theoretical groundwork for DiDC, examines its
identification assumptions, estimation procedures, and asymptotic properties.
We showcase its advantages through Monte Carlo simulations and an empirical
application.

DiDC can handle scenarios where control and treatment groups differ
significantly, violating the parallel trends assumption of DiD, or when RDD
encounters confounding factors at the threshold. By incorporating more
information, DiDC eliminates bias in RDD estimates under specific assumptions
about the data-generating processes.

However, it is important to give the necessary importance to the
identification assumptions, particularly the time-invariance of confounding
effects. We propose a test based on stacked RDDs to assess its validity in
practic. Additionally, DiDC requires the treatment effect to be independent of
confounding policy, though we introduce a possible relaxation for potential
interaction effects.

We find that the DiDC method can eliminate bias completely if the shape of the
functions remains stable over time at each side of the threshold. This
suggests that it could offer significant advantages over standard RDD
estimations, even in settings where no other confounding variables exist at
the threshold. We also propose test to compare the derivatives of estimated
functions on either side of the threshold, allowing researchers to evaluate
the stability of data-generating processes over time.

Monte Carlo simulations demonstrate DiDC’s potential for improvement over RDD,
yielding lower bias and better coverage. It can provide more desirable results
by incorporating more data, especially when the functional form exhibits
minimal temporal variation. Notably, DiDC is the only viable approach when
confounding factors render both RDD and DiD unsuitable. The empirical
application highlights the importance of the time-invariance assumption.

Future research directions include the development of alternative estimators
that are robust to violations of identification assumptions, as well as the
exploration of other confidence interval methods specifically tailored for the
DiDC design. Overall, the DiDC method offers a valuable addition to the causal
inference toolkit. It is applicable in settings where no other methods were
previously available and shows potential to reduce bias in estimation in other
settings.

## References

  * Albright (2023) A. Albright.  The hidden effects of algorithmic recommendations.  _Preprint. Last accessed March_ , 28:2023, 2023.
  * Azuaga and Sampaio (2017) F. Azuaga and B. Sampaio.  Violência contra mulher: o impacto da lei maria da penha sobre o feminicídio no brasil.  _Anais do 45o Encontro ANPEC_ , 2017.
  * Butts (2021) K. Butts.  Geographic difference-in-discontinuities.  _Applied Economics Letters_ , pages 1–5, 2021.
  * Calonico et al. (2014) S. Calonico, M. D. Cattaneo, and R. Titiunik.  Robust nonparametric confidence intervals for regression-discontinuity designs.  _Econometrica_ , 82(6):2295–2326, 2014.
  * Calonico et al. (2022) S. Calonico, M. D. Cattaneo, M. H. Farrell, and R. Titiunik.  _rdrobust: Robust Data-Driven Statistical Inference in Regression-Discontinuity Designs_ , 2022.  URL https://CRAN.R-project.org/package=rdrobust.  R package version 2.1.1.
  * Card and Krueger (2000) D. Card and A. B. Krueger.  Minimum wages and employment: a case study of the fast-food industry in new jersey and pennsylvania: reply.  _American Economic Review_ , 90(5):1397–1420, 2000.
  * Cellini et al. (2010) S. R. Cellini, F. Ferreira, and J. Rothstein.  The value of school facility investments: Evidence from a dynamic regression discontinuity design.  _The Quarterly Journal of Economics_ , 125(1):215–261, 2010.
  * Chicoine (2017) L. E. Chicoine.  Homicides in mexico and the expiration of the us federal assault weapons ban: a difference-in-discontinuities approach.  _Journal of economic geography_ , 17(4):825–856, 2017.
  * Eggers et al. (2018) A. C. Eggers, R. Freier, V. Grembi, and T. Nannicini.  Regression discontinuity designs based on population thresholds: Pitfalls and solutions.  _American Journal of Political Science_ , 62(1):210–229, 2018.
  * Frölich and Sperlich (2019) M. Frölich and S. Sperlich.  _Impact evaluation_.  Cambridge University Press, 2019.
  * Galindo-Silva et al. (2021) H. Galindo-Silva, N. H. Some, and G. Tchuente.  Fuzzy difference-in-discontinuities: Identification theory and application to the affordable care act.  _arXiv preprint arXiv:1812.06537_ , 2021.
  * García (2020) R. E. García.  Plant behavior & equal pay: The effect on female employment, capital investments, & productivity—a difference-in-discontinuities design.  _Capital Investments, & Productivity—A Difference-in-Discontinuities Design (August 14, 2020)_, 2020.
  * Grembi et al. (2016) V. Grembi, T. Nannicini, and U. Troiano.  Do fiscal rules matter?  _American Economic Journal: Applied Economics_ , pages 1–30, 2016.
  * Hahn et al. (2001) J. Hahn, P. Todd, and W. Van der Klaauw.  Identification and estimation of treatment effects with a regression-discontinuity design.  _Econometrica_ , 69(1):201–209, 2001.
  * Imbens and Lemieux (2008) G. W. Imbens and T. Lemieux.  Regression discontinuity designs: A guide to practice.  _Journal of econometrics_ , 142(2):615–635, 2008.
  * Lemieux and Milligan (2008) T. Lemieux and K. Milligan.  Incentive effects of social assistance: A regression discontinuity approach.  _Journal of Econometrics_ , 142(2):807–828, 2008.
  * Ludwig and Miller (2007) J. Ludwig and D. L. Miller.  Does head start improve children’s life chances? evidence from a regression discontinuity design.  _The Quarterly journal of economics_ , 122(1):159–208, 2007.
  * Pettersson-Lidbom (2012) P. Pettersson-Lidbom.  Does the size of the legislature affect the size of government? evidence from two natural experiments.  _Journal of Public Economics_ , 96(3-4):269–278, 2012.
  * Sant’Anna and Zhao (2020) P. H. Sant’Anna and J. Zhao.  Doubly robust difference-in-differences estimators.  _Journal of Econometrics_ , 219(1):101–122, 2020.

##  Appendix A Setup, assumptions and notation for estimation

We construct the local polynomial estimator following Calonico et al. (2014).
For a given $\nu\leq p\in\mathbf{N}$, the general estimand of interest is
$\tau_{\nu}=\Delta\mu_{+}-\Delta\mu_{-}$ with
$\Delta\mu^{(\nu)}_{+}=\nu!e^{\prime}_{\nu}\delta_{+,p}$,
$\Delta\mu^{(\nu)}_{-}=\nu!e^{\prime}_{\nu}\delta_{-,p}$ being the $\nu$th-
order derivatives of the $p$th-order local polynomial of the difference. The
$p$th-order local polynomial estimators of the $\nu$th-order derivatives
$\Delta\mu^{(\nu)}_{+,p}$ and $\Delta\mu^{(\nu)}_{-,p}$ are:

$\displaystyle\Delta\hat{\mu}^{(\nu)}_{+,p}(h_{n})$
$\displaystyle=\nu!e^{\prime}_{\nu}\hat{\delta}_{+,p}(h_{n})$
$\displaystyle\Delta\hat{\mu}^{(\nu)}_{-,p}(h_{n})$
$\displaystyle=\nu!e^{\prime}_{\nu}\hat{\delta}_{-,p}(h_{n})$
$\displaystyle\hat{\delta}_{\Delta Y+,p}(h_{n})$
$\displaystyle=arg\min_{\delta\in\mathbf{R}^{p+1}}\sum_{i=1}^{n}\mathbbm{1}(Z_{i}\geq
0)(\Delta Y_{i}-r_{p}(Z_{i})^{\prime}\delta)^{2}K_{h_{n}}(Z_{i})$
$\displaystyle\hat{\delta}_{\Delta Y-,p}(h_{n})$
$\displaystyle=arg\min_{\delta\in\mathbf{R}^{p+1}}\sum_{i=1}^{n}\mathbbm{1}(Z_{i}<0)(\Delta
Y_{i}-r_{p}(Z_{i})^{\prime}\delta)^{2}K_{h_{n}}(Z_{i})$

where $e_{\nu}$ is a conformable $(\nu+1)$ unit vector, $K_{h}(u)=K(u/h)/h$,
$h_{n}$ is a positive bandwidth sequence,
$r_{p}(x)=\begin{bmatrix}1&x&\ldots&x^{p}\end{bmatrix}^{\prime}$, $\Delta
Y=\begin{bmatrix}\Delta Y_{1}&\Delta Y_{2}&\ldots&\Delta
Y_{n}\end{bmatrix}^{\prime}$. We set
$\chi_{n}=\begin{bmatrix}Z_{1}&\ldots&Z_{n}\end{bmatrix}^{\prime}$,
$\varepsilon_{\Delta Y}=\begin{bmatrix}\varepsilon_{\Delta
Y,1}&\ldots&\varepsilon_{\Delta Y,n}\end{bmatrix}^{\prime}$ with
$\varepsilon_{\Delta Y,i}=\Delta Y_{i}-\mu_{\Delta Y}(Z_{i})$, $\mu_{\Delta
Y}(Z)=E\left(\Delta Y|Z\right)$ and $\sigma^{2}_{\Delta W,\Delta
Y}(Z)=Cov\left(\Delta W,\Delta Y|Z\right)$. We also set

$\displaystyle S_{p}(h)$
$\displaystyle=\begin{bmatrix}\left(Z_{1}/h\right)^{p}&\ldots&\left(Z_{n}/h\right)^{p}\end{bmatrix}^{\prime}$
$\displaystyle Z_{p}(h)$
$\displaystyle=\begin{bmatrix}r_{p}(Z_{1}/h)&\ldots&r_{p}(Z_{n}/h)\end{bmatrix}^{\prime}$
$\displaystyle W_{+}(h)$ $\displaystyle=diag\left(\mathbbm{1}(Z_{i}\geq
0)K_{h}(Z_{1}),\ldots,\mathbbm{1}(Z_{i}\geq 0)K_{h}(Z_{n})\right)$
$\displaystyle W_{-}(h)$
$\displaystyle=diag\left(\mathbbm{1}(Z_{i}<0)K_{h}(Z_{1}),\ldots,\mathbbm{1}(Z_{i}<0)K_{h}(Z_{n})\right)$
$\displaystyle\Sigma_{\Delta W\Delta Y}$
$\displaystyle=diag\left(\sigma_{\Delta W\Delta
Y}^{2}(Z_{i}),\ldots,\sigma_{\Delta W\Delta Y}^{2}(Z_{n})\right)$
$\displaystyle\Gamma_{+,p}(h)$
$\displaystyle=Z_{p}(h)^{\prime}W_{+}(h)Z_{p}(h)/n$
$\displaystyle\Gamma_{-,p}(h)$
$\displaystyle=Z_{p}(h)^{\prime}W_{-}(h)Z_{p}(h)/n$
$\displaystyle\vartheta_{+,p,q}(h)$
$\displaystyle=Z_{p}(h)^{\prime}W_{+}(h)S_{q}(h)/n$
$\displaystyle\vartheta_{-,p,q}(h)$
$\displaystyle=Z_{p}(h)^{\prime}W_{-}(h)S_{q}(h)/n$ $\displaystyle\Psi_{\Delta
W\Delta Y+,p,q}(h,b)$ $\displaystyle=Z_{p}(h)^{\prime}W_{+}(h)\Sigma_{\Delta
W\Delta Y}W_{+}(b)Z_{q}(b)/n$ $\displaystyle\Psi_{\Delta W\Delta Y-,p,q}(h,b)$
$\displaystyle=Z_{p}(h)^{\prime}W_{-}(h)\Sigma_{\Delta W\Delta
Y}W_{-}(b)Z_{q}(b)/n$

It follows that with $H_{p}(h)=diag\left(1,h^{-1},\ldots,h^{-p}\right)$:

$\displaystyle\hat{\delta}_{\Delta Y+,p}(h_{n})$
$\displaystyle=H_{p}(h_{n})\Gamma^{-1}_{+,p}(h_{n})Z_{p}(h_{n})^{\prime}W_{+}(h_{n})\Delta
Y/n$ $\displaystyle\hat{\delta}_{\Delta Y-,p}(h_{n})$
$\displaystyle=H_{p}(h_{n})\Gamma^{-1}_{-,p}(h_{n})Z_{p}(h_{n})^{\prime}W_{-}(h_{n})\Delta
Y/n$

The estimand and estimators are

$\displaystyle\tau^{DiDC}_{\nu}=\Delta\mu_{+}^{(\nu)}-\Delta\mu_{-}^{(\nu)},\quad\Delta\mu_{+}^{(\nu)}=\nu!e_{\nu}^{\prime}\delta_{+,p},\quad\Delta\mu_{-}^{(\nu)}=\nu!e_{\nu}^{\prime}\delta_{-,p},$
$\displaystyle\hat{\tau}^{DiDC}_{\nu,p}\left(h_{n}\right)=\Delta\hat{\mu}_{+,p}^{(\nu)}\left(h_{n}\right)-\Delta\hat{\mu}_{-,p}^{(\nu)}\left(h_{n}\right),$
$\displaystyle\Delta\hat{\mu}_{+,p}^{(\nu)}\left(h_{n}\right)=\nu!e_{\nu}^{\prime}\hat{\delta}_{+,p}\left(h_{n}\right),\quad\Delta\hat{\mu}_{-,p}^{(\nu)}\left(h_{n}\right)=\nu!e_{\nu}^{\prime}\hat{\delta}_{-,p}\left(h_{n}\right),$

where, for any random variables $W$ and $X$, and $s\in\mathbb{N}$,

$\displaystyle\Delta\mu_{X+}^{(s)}=\lim_{x\rightarrow
0^{+}}\frac{\partial^{s}}{\partial
z^{s}}\Delta\mu_{X}(z),\quad\Delta\mu_{X-}^{(s)}=\lim_{x\rightarrow
0^{-}}\frac{\partial^{s}}{\partial z^{s}}\Delta\mu_{X}(z),$
$\displaystyle\Delta\mu_{X}(z)=\mathbb{E}[X\mid Z=z],$
$\displaystyle\sigma_{X+}^{2}=\lim_{z\rightarrow
0^{+}}\sigma_{X}^{2}(z),\quad\sigma_{X-}^{2}=\lim_{z\rightarrow
0^{-}}\sigma_{X}^{2}(z),$ $\displaystyle\sigma_{X}^{2}(z)=\mathbb{V}[X\mid
Z=z],$ $\displaystyle\sigma_{WX+}^{2}=\lim_{z\rightarrow
0^{+}}\sigma_{WX}^{2}(z),\quad\sigma_{WX_{-}}^{2}=\lim_{z\rightarrow
0^{-}}\sigma_{WX}^{2}(z),$ $\displaystyle\sigma_{WX}^{2}(z)=\mathbb{C}[W,X\mid
X=z].$

We employ the following assumptions on the sharp model for the non-parametric
local polynomial regression estimation:

######  Assumption A.1.

For some $\mathcal{K}_{0}>0$, the following holds in the neighborhood
$\left(-\mathcal{K}_{0},\mathcal{K}_{0}\right)$ around the cutoff $z_{0}=0$:

  1. (a)

$E\left[\Delta Y_{i}^{4}|Z_{i}=z\right]$ is bounded, and the density $f(z)$ of
the random sample $Z_{i}$ is continuous and bounded away from zero.

  2. (b)

$\Delta\mu_{-}(z)=E\left[\Delta Y_{i}(0)|Z_{i}=z\right]$ and
$\Delta\mu_{+}(z)=E\left[\Delta Y_{i}(1)|Z_{i}=z\right]$ are S times
continuously differentiable.

  3. (c)

$\sigma^{2}_{-}(z)=V\left[\Delta Y_{i}(0)|Z_{i}=z\right]$ and
$\sigma^{2}_{+}(z)=V\left[\Delta Y_{i}(1)|Z_{i}=z\right]$ are continuous and
bounded away from zero.

We also impose the following assumption on the kernel function to be employed
in the estimator. This assumption allows for most of the commonly used
kernels.

######  Assumption A.2.

For some $\mathcal{K}>0$, the kernel function
$k(\cdot):[0,\mathcal{K}]\rightarrow\mathbbm{R}$ is bounded and nonnegative,
zero outside its support, symmetric around $z_{0}$ and positive and continuous
on $\left(0,\mathcal{K}\right)$.

Assumptions A.1 and A.2 limit the behavior of $E\left(\Delta
Y_{i}|Z_{i}=z_{0}\right)$ in the vicinity of the cutoff $z_{0}=0$.

##  Appendix B Preliminary lemmas and results

Before proceeding, please refer to Appendix A for notation. This appendix
restates, with minor adaptations, several lemmas, results and proofs from
Calonico et al. (2014) that are necessary for deriving the asymptotic results.

The following lemma establishes convergence in probability of the sample
matrices $\Gamma_{-,p}\left(h_{n}\right)$,
$\boldsymbol{\vartheta}_{-,p,q}\left(h_{n}\right)$,
$\Psi_{-,p}\left(h_{n}\right)$ and $\Gamma_{+,p}\left(h_{n}\right)$,
$\boldsymbol{\vartheta}_{+,p,q}\left(h_{n}\right)$,
$\Psi_{+,p}\left(h_{n}\right)$ to their expectation counterparts, and
characterizes those limits.

######  Lemma B.1.

Suppose Assumptions $1-2$ hold, and $nh_{n}\rightarrow\infty$.

  1. (a)

If $\kappa h_{n}<\kappa_{0}$, then:

    1. (a.1)

$\Gamma_{+,p}\left(h_{n}\right)=\tilde{\Gamma}_{p}\left(h_{n}\right)+o_{p}(1)$
with $\tilde{\Gamma}_{+,p}\left(h_{n}\right)=\int_{0}^{\infty}K(u)\allowbreak
r_{p}(u)\allowbreak r_{p}(u)^{\prime}\allowbreak
f\left(uh_{n}\right)\allowbreak\mathrm{d}u\allowbreak\asymp\Gamma_{p}$

    2. (a.2)

$\Gamma_{-,p}\left(h_{n}\right)=H_{p}(-1)\tilde{\Gamma}_{p}\left(h_{n}\right)H_{p}(-1)+o_{p}(1)$
with $\tilde{\Gamma}_{-,p}\left(h_{n}\right)=\int_{0}^{\infty}\allowbreak
K(u)\allowbreak r_{p}(u)\allowbreak r_{p}(u)^{\prime}\allowbreak
f\left(-uh_{n}\right)\allowbreak\mathrm{d}u\asymp\Gamma_{p}$,

    3. (a.3)

$\boldsymbol{\vartheta}_{+,p,q}\left(h_{n}\right)=\tilde{\boldsymbol{\vartheta}}_{+,p,q}\left(h_{n}\right)+o_{p}(1)\text{
with
}\tilde{\boldsymbol{\vartheta}}_{+,p,q}\left(h_{n}\right)=\int_{0}^{\infty}K(u)r_{p}(u)\allowbreak
u^{q}\allowbreak
f\left(uh_{n}\right)\allowbreak\mathrm{d}u\allowbreak\asymp\boldsymbol{\vartheta}_{p,q}$,

    4. (a.4)

$\boldsymbol{\vartheta}_{-,p,q}\left(h_{n}\right)=(-1)^{q}H_{p}(-1)\tilde{\boldsymbol{\vartheta}}_{-,p,q}\left(h_{n}\right)+o_{p}(1)\text{
with
}\tilde{\boldsymbol{\vartheta}}_{-,p,q}\left(h_{n}\right)=\int_{0}^{\infty}\allowbreak
K(u)\allowbreak r_{p}(u)\allowbreak u^{q}\allowbreak
f\left(-uh_{n}\right)\allowbreak\mathrm{d}u\allowbreak\asymp\boldsymbol{\vartheta}_{p,q}$,

    5. (a.5)

$h_{n}\Psi_{+,p}\left(h_{n}\right)=\tilde{\Psi}_{+,p}\left(h_{n}\right)+o_{p}(1)\text{
with }\tilde{\Psi}_{+,p}\left(h_{n}\right)=\int_{0}^{\infty}\allowbreak
K(u)^{2}\allowbreak r_{p}(u)\allowbreak
r_{p}(u)^{\prime}\allowbreak\boldsymbol{\sigma}_{+}^{2}\left(uh_{n}\right)\allowbreak
f\left(uh_{n}\right)\allowbreak\mathrm{d}u\asymp\Psi_{p}$,

    6. (a.6)

$h_{n}\Psi_{-,p}\left(h_{n}\right)=H_{p}(-1)\tilde{\Psi}_{-,p}\left(h_{n}\right)H_{p}(-1)+o_{p}(1)\text{
with }\tilde{\Psi}_{-,p}\left(h_{n}\right)=\int_{0}^{\infty}\allowbreak
K(u)^{2}\allowbreak r_{p}(u)\allowbreak
r_{p}(u)^{\prime}\allowbreak\sigma_{-}^{2}\left(-uh_{n}\right)\allowbreak
f\left(-uh_{n}\right)\allowbreak\mathrm{d}u\asymp\Psi_{p}$.

  2. (b)

If $h_{n}\rightarrow 0$, then

    1. (b.1)

$\tilde{\Gamma}_{+,p}\left(h_{n}\right)=f\Gamma_{p}+o(1)$ and
$\tilde{\Gamma}_{-,p}\left(h_{n}\right)=f\Gamma_{p}+o(1)$,

    2. (b.2)

$\tilde{\vartheta}_{+,p,q}\left(h_{n}\right)=f\vartheta_{p,q}+o(1)$ and
$\tilde{\vartheta}_{-,p,q}\left(h_{n}\right)=f\boldsymbol{\vartheta}_{p,q}+o(1)$,

    3. (b.3)

$\tilde{\Psi}_{+,p}^{+,p,q}\left(h_{n}\right)=\sigma_{+}^{2}f\Psi_{p}+o(1)$
and $\tilde{\Psi}_{-,p}\left(h_{n}\right)=\sigma_{-}^{2}f\Psi_{p}+o(1)$.

###### Proof.

For part (a.5), change of variable implies

$\displaystyle\mathbb{E}\left[\Psi_{+,p}\left(h_{n}\right)\right]$
$\displaystyle=\mathbb{E}\left[h_{n}Z_{p}(h)^{\prime}W_{+}(h)\Sigma_{\Delta
W\Delta Y}W_{+}(b)Z_{q}(b)/n\right]$
$\displaystyle=\int_{0}^{\infty}K\left(\frac{z}{h_{n}}\right)^{2}r_{p}\left(\frac{z}{h_{n}}\right)r_{p}\left(\frac{z}{h_{n}}\right)^{\prime}\sigma_{+}^{2}f(z)\mathrm{d}z$
$\displaystyle=\int_{0}^{\infty}K(u)^{2}r_{p}(u)r_{p}(u)^{\prime}\sigma_{+}^{2}\left(uh_{n}\right)f\left(uh_{n}\right)\mathrm{d}u$
$\displaystyle=\tilde{\Psi}_{+,p}\left(h_{n}\right),$

and
$h_{n}^{2}\mathbb{E}\left[\left|\Psi_{+,p}\left(h_{n}\right)-\mathbb{E}\left[\Psi_{+,p}\left(h_{n}\right)\right]\right|^{2}\right]=n^{-1}h_{n}^{-1}\int_{0}^{\infty}K(u)^{4}\left|r_{p}(u)\right|^{4}f\left(uh_{n}\right)\mathrm{d}u=$
$O\left(n^{-1}h_{n}^{-1}\right)$, provided $\kappa h_{n}<\kappa_{0}$. For part
(a.6),

$\displaystyle\mathbb{E}\left[h_{n}\Psi_{-,p}\left(h_{n}\right)\right]$
$\displaystyle=h_{n}^{-1}\int_{-\infty}^{0}K\left(u/h_{n}\right)^{2}r_{p}\left(u/h_{n}\right)r_{p}\left(u/h_{n}\right)^{\prime}\sigma_{-}^{2}(u)f(u)\mathrm{d}u$
$\displaystyle=H_{p}(-1)\tilde{\Psi}_{-,p}\left(h_{n}\right)H_{p}(-1),$

and the rest is proven as above. Also, note that
$\tilde{\Psi}_{+,p}^{+,p,q}\left(h_{n}\right)=\sigma_{+}^{2}f\Psi_{p}+o(1)$
and $\tilde{\Psi}_{-,p}\left(h_{n}\right)=\sigma_{-}^{2}f\Psi_{p}+o(1)$ if
$h_{n}\rightarrow 0$, by continuity of $\sigma_{+}^{2}(u)$,
$\sigma_{-}^{2}(u)$ and $f(u)$, which proves part (b.3).

Proofs for the other items follow similarly to the one above and can be found
in Calonico et al. (2014), as they are identical to those provided there. ∎

Let $s,\ell\in\mathbb{N}$ with $s\leq\ell$. The following lemma gives the
asymptotic bias, variance, and distribution for the $\ell$ th-order local
polynomial estimator of $\Delta\mu_{+}^{(s)}$ and $\Delta\mu_{-}^{(s)}:$

$\displaystyle\Delta\hat{\mu}_{+,\ell}^{(s)}\left(h_{n}\right)=s!e_{s}^{\prime}\delta_{+,\ell}\left(h_{n}\right),$
$\displaystyle\Delta\hat{\delta}_{+,\ell}\left(h_{n}\right)=H_{\ell}\left(h_{n}\right)\Gamma_{+,\ell}^{-1}\left(h_{n}\right)Z_{\ell}\left(h_{n}\right)^{\prime}W_{+}\left(h_{n}\right)Y/n,$
$\displaystyle\Delta\hat{\mu}_{-,\ell}^{(s)}\left(h_{n}\right)=s!e_{s}^{\prime}\hat{\delta}_{-,\ell}\left(h_{n}\right),$
$\displaystyle\hat{\delta}_{-,\ell}\left(h_{n}\right)=H_{\ell}\left(h_{n}\right)\Gamma_{-,\ell}^{-1}\left(h_{n}\right)Z_{\ell}\left(h_{n}\right)^{\prime}W_{-}\left(h_{n}\right)Y/n.$

######  Lemma B.2.

Suppose Assumptions 1 and 2 hold with $S\geq\ell+2$, and
$nh_{n}\rightarrow\infty$.

  1. (B)

If $h_{n}\rightarrow 0$, then

$\displaystyle\mathbb{E}\left[\Delta\hat{\mu}_{+,\ell}^{(s)}\left(h_{n}\right)\mid\mathcal{X}_{n}\right]=s!e_{s}^{\prime}\delta_{+,\ell}+h_{n}^{1+\ell-s}\frac{\Delta\mu_{+}^{(\ell+1)}}{(\ell+1)!}\mathcal{B}_{+,s,\ell,\ell+1}\left(h_{n}\right)$
$\displaystyle+h_{n}^{2+\ell-s}\frac{\Delta\mu_{+}^{(\ell+2)}}{(\ell+2)!}\mathcal{B}_{+,s,\ell,\ell+2}\left(h_{n}\right)+o_{p}\left(h_{n}^{2+\ell-s}\right),$
$\displaystyle\mathcal{B}_{+,s,\ell,r}\left(h_{n}\right)=s!e_{s}^{\prime}\Gamma_{+,\ell}^{-1}\left(h_{n}\right)\vartheta_{+,\ell,r}\left(h_{n}\right)=s!e_{s}^{\prime}\Gamma_{\ell}^{-1}\vartheta_{\ell,r}+o_{p}(1),$

and

$\displaystyle\mathbb{E}\left[\Delta\hat{\mu}_{-,\ell}^{(s)}\left(h_{n}\right)\mid\mathcal{X}_{n}\right]=s!e_{s}^{\prime}\delta_{-,\ell}+h_{n}^{1+\ell-s}\frac{\Delta\mu_{-}^{(\ell+1)}}{(\ell+1)!}\mathcal{B}_{-,s,\ell,\ell+1}\left(h_{n}\right)$
$\displaystyle+h_{n}^{2+\ell-s}\frac{\Delta\mu_{-}^{(\ell+2)}}{(\ell+2)!}\mathcal{B}_{-,s,\ell,\ell+2}\left(h_{n}\right)+o_{p}\left(h_{n}^{2+\ell-s}\right),$
$\displaystyle\mathcal{B}_{-,s,\ell,r}\left(h_{n}\right)=s!e_{s}^{\prime}\Gamma_{-,\ell}^{-1}\left(h_{n}\right)\vartheta_{-,\ell,r}\left(h_{n}\right)=(-1)^{s+r}s!e_{s}^{\prime}\Gamma_{\ell}^{-1}\boldsymbol{\vartheta}_{\ell,r}+o_{p}(1).$

  2. (V)

If $h_{n}\rightarrow 0$, then
$\mathbb{V}\left[\Delta\hat{\mu}_{+,\ell}^{(s)}\left(h_{n}\right)\mid\mathcal{X}_{n}\right]=\mathcal{V}_{+,s,\ell}\left(h_{n}\right)$
with

$\displaystyle\mathcal{V}_{+,s,\ell}\left(h_{n}\right)$
$\displaystyle=\frac{1}{nh_{n}^{2s}}s!^{2}e_{s}^{\prime}\Gamma_{+,\ell}^{-1}\left(h_{n}\right)\Psi_{+,\ell}\left(h_{n}\right)\Gamma_{+,\ell}^{-1}\left(h_{n}\right)e_{s}$
$\displaystyle=\frac{1}{nh_{n}^{1+2s}}\frac{\sigma_{+}^{2}}{f}s!^{2}e_{s}^{\prime}\Gamma_{\ell}^{-1}\Psi_{\ell}\Gamma_{\ell}^{-1}e_{s}\left[1+o_{p}(1)\right],$

and
$\mathbb{V}\left[\Delta\hat{\mu}_{-,\ell}^{(s)}\left(h_{n}\right)\mid\mathcal{X}_{n}\right]=\mathcal{V}_{-,s,\ell}\left(h_{n}\right)$
with

$\displaystyle\mathcal{V}_{-,s,\ell}\left(h_{n}\right)$
$\displaystyle=\frac{1}{nh_{n}^{2s}}s^{2}e_{s}^{\prime}\Gamma_{-\ell}^{-1}\left(h_{n}\right)\Psi_{-,\ell}\left(h_{n}\right)\Gamma_{-,\ell}^{-1}\left(h_{n}\right)e_{s}$
$\displaystyle=\frac{1}{nh_{n}^{1+2s}}\frac{\sigma_{-}^{2}}{f}s!^{2}e_{s}^{\prime}\Gamma_{\ell}^{-1}\Psi_{\ell}\Gamma_{\ell}^{-1}e_{s}\left[1+o_{p}(1)\right].$

  3. (D)

If $nh_{n}^{2\ell+5}\rightarrow 0$, then

$\displaystyle\frac{\Delta\hat{\mu}_{+,\ell}^{(s)}\left(h_{n}\right)-\Delta\mu_{+}^{(s)}-h_{n}^{1+\ell-s}\frac{\Delta\mu_{+}^{(\ell+1)}}{(\ell+1)!}\mathcal{B}_{+,s,\ell,\ell+1}\left(h_{n}\right)}{\sqrt{\mathcal{V}_{+,s,\ell}\left(h_{n}\right)}}\rightarrow{}_{d}\mathcal{N}(0,1)$

and

$\displaystyle\frac{\Delta\hat{\mu}_{-,\ell}^{(s)}\left(h_{n}\right)-\Delta\mu_{-}^{(s)}-h_{n}^{1+\ell-s}\frac{\Delta\mu_{-}^{(\ell+1)}}{(\ell+1)!}\mathcal{B}_{-,s,\ell,\ell+1}\left(h_{n}\right)}{\sqrt{\mathcal{V}_{-,s,\ell}\left(h_{n}\right)}}\rightarrow_{d}\mathcal{N}(0,1).$

###### Proof.

For part (B), a Taylor series expansion gives

$\displaystyle\mathbb{E}[s!$
$\displaystyle\left.\hat{\delta}_{+,\ell}\left(h_{n}\right)\mid\mathcal{X}_{n}\right]$
$\displaystyle=$ $\displaystyle
s!\delta_{+,\ell}+h_{n}^{\ell+1}H_{\ell}\left(h_{n}\right)\Gamma_{+,\ell}^{-1}\left(h_{n}\right)Z_{\ell}\left(h_{n}\right)W_{+}\left(h_{n}\right)S_{\ell+1}\left(h_{n}\right)s!\frac{\Delta\mu_{+}^{(\ell+1)}}{(\ell+1)!}$
$\displaystyle+h_{n}^{\ell+2}H_{\ell}\left(h_{n}\right)\Gamma_{+,\ell}^{-1}\left(h_{n}\right)Z_{\ell}\left(h_{n}\right)W_{+}\left(h_{n}\right)S_{\ell+2}\left(h_{n}\right)s!\frac{\Delta\mu_{+}^{(\ell+2)}}{(\ell+2)!}+H_{\ell}\left(h_{n}\right)o_{p}\left(h_{n}^{\ell+2}\right)$
$\displaystyle=$ $\displaystyle
s!\delta_{+,\ell}+h_{n}^{\ell+1}H_{\ell}\left(h_{n}\right)s!\frac{\Delta\mu_{+}^{(\ell+1)}}{(\ell+1)!}\Gamma_{+,\ell}^{-1}\left(h_{n}\right)\vartheta_{+,\ell,\ell+1}\left(h_{n}\right)$
$\displaystyle+h_{n}^{\ell+2}H_{\ell}\left(h_{n}\right)s!\frac{\Delta\mu_{+}^{(\ell+2)}}{(\ell+2)!}\Gamma_{+,\ell}^{-1}\left(h_{n}\right)\vartheta_{+,\ell,\ell+2}\left(h_{n}\right)+H_{\ell}\left(h_{n}\right)o_{p}\left(h_{n}^{\ell+2}\right),$

and the result for
$\mathbb{E}\left[\Delta\hat{\mu}_{+,\ell}^{(s)}\left(h_{n}\right)\mid\mathcal{X}_{n}\right]$
follows by $e_{s}^{\prime}H_{\ell}\left(h_{n}\right)=h_{n}^{-s}$ and Lemma
B.1. Next, for
$\mathbb{E}\left[\Delta\hat{\mu}_{-,\ell}^{(s)}\left(h_{n}\right)\mid\mathcal{X}_{n}\right]$
the same calculations apply, with only a modification for
$\mathcal{B}_{-,s,\ell,r}\left(h_{n}\right)$ because, by Lemma B.1,

$\displaystyle\mathcal{B}_{-,s,\ell,r}\left(h_{n}\right)=$ $\displaystyle
s!e_{s}^{\prime}\Gamma_{-,\ell}^{-1}\left(h_{n}\right)\vartheta_{-,\ell,r}\left(h_{n}\right)$
$\displaystyle=$ $\displaystyle
s!e_{s}^{\prime}\left[H_{\ell}(-1)\tilde{\Gamma}_{-,\ell}^{-1}\left(h_{n}\right)H_{\ell}(-1)\right]\left[(-1)^{r}H_{\ell}(-1)\vartheta_{-,\ell,r}\left(h_{n}\right)\right]$
$\displaystyle+o_{p}(1)$ $\displaystyle=$
$\displaystyle(-1)^{s+r}s!e_{s}^{\prime}\tilde{\Gamma}_{-,\ell}^{-1}\left(h_{n}\right)\vartheta_{-,\ell,r}\left(h_{n}\right)+o_{p}(1),$

because $e_{s}^{\prime}H_{\ell}(-1)=(-1)^{s}$ and
$H_{\ell}(-1)H_{\ell}(-1)=I_{\ell+1}$. For part (V), simply note that

$\displaystyle\mathbb{V}\left[s!e_{s}^{\prime}\hat{\delta}_{+,\ell}\left(h_{n}\right)\mid\mathcal{X}_{n}\right]=$
$\displaystyle
s!^{2}e_{s}^{\prime}H_{\ell}\left(h_{n}\right)\Gamma_{+,\ell}^{-1}\left(h_{n}\right)Z_{\ell}\left(h_{n}\right)W_{+}\left(h_{n}\right)\Sigma$
$\displaystyle\times
W_{+}\left(h_{n}\right)Z_{\ell}\left(h_{n}\right)\Gamma_{+,\ell}^{-1}\left(h_{n}\right)H_{\ell}\left(h_{n}\right)e_{s}/n$
$\displaystyle=$ $\displaystyle
h_{n}^{-2s}s!^{2}e_{s}^{\prime}\Gamma_{+,\ell}^{-1}\left(h_{n}\right)Z_{\ell}\left(h_{n}\right)W_{+}\left(h_{n}\right)\Sigma$
$\displaystyle\times
W_{+}\left(h_{n}\right)Z_{\ell}\left(h_{n}\right)\Gamma_{+,\ell}^{-1}\left(h_{n}\right)e_{s}/n$
$\displaystyle=$ $\displaystyle
n^{-1}h_{n}^{-2s}s!^{2}e_{s}^{\prime}\Gamma_{+,\ell}^{-1}\left(h_{n}\right)\Psi_{+,\ell}\left(h_{n}\right)\Gamma_{+,\ell}^{-1}\left(h_{n}\right)e_{s}$
$\displaystyle=$ $\displaystyle\mathcal{V}_{+,\ell,s}\left(h_{n}\right),$

and the result follows by Lemma B.1. The proof of
$\mathbb{V}\left[s!e_{s}^{\prime}\hat{\delta}_{-,\ell}\left(h_{n}\right)\mid\mathcal{X}_{n}\right]$
is analogous.

For part (D), using the previous results, we have

$\displaystyle=\frac{s!e_{s}^{\prime}\hat{\boldsymbol{\delta}}_{+,\ell}\left(h_{n}\right)-s!e_{s}^{\prime}\delta_{+,\ell}-h_{n}^{1+\ell-s}\frac{\Delta\mu_{+}^{(\ell+1)}}{(\ell+1)!}\mathcal{B}_{+,s,\ell,\ell+1}\left(h_{n}\right)}{\sqrt{\mathcal{V}_{+,s,\ell}\left(h_{n}\right)}}$

The result for $\Delta\hat{\mu}_{-,\ell}^{(s)}\left(h_{n}\right)$ can be
established the same way. This concludes the proof. Q.E.D. ∎

Let $\nu$, $p$, $q$ $\in\mathbb{N}$ with $\nu\leq p<q$. The nex lemma dives
the asymptotic bias, variance and distribution for the $p$th-order local
polynomial estimator of $\Delta\mu_{+}^{(\nu)}$ and $\Delta\mu_{-}^{(\nu)}$
with bias correction constructed using a $q$th-order local polynomial.

######  Lemma B.3.

Suppose Assumptions 1 and 2 hold with $S\geq q+1$, and
$n\min\left\\{h_{n},b_{n}\right\\}\rightarrow\infty$. (B) If
$\max\left\\{h_{n},b_{n}\right\\}\rightarrow 0$, then

$\displaystyle\mathbb{E}\left[\Delta\hat{\mu}_{+,p,q}^{(\nu)\text{ bc
}}\left(h_{n},b_{n}\right)\mid\mathcal{X}_{n}\right]$
$\displaystyle=\nu!e_{\nu}^{\prime}\delta_{+,p}+h_{n}^{2+p-\nu}\frac{\Delta\mu_{+}^{(p+2)}}{(p+2)!}\mathcal{B}_{+,\nu,p,p+2}\left(h_{n}\right)\left\\{1+o_{p}(1)\right\\}$
$\displaystyle\quad-
h_{n}^{p+1-\nu}b_{n}^{q-p}\frac{\Delta\mu_{+}^{(q+1)}}{(q+1)!}\mathcal{B}_{+,p+1,q,q+1}\left(b_{n}\right)\frac{\mathcal{B}_{+,\nu,p,p+1}\left(h_{n}\right)}{(p+1)!}\left\\{1+o_{p}(1)\right\\}$

and

$\displaystyle\mathbb{E}\left[\Delta\hat{\mu}_{-,p,q}^{(\nu)\mathrm{bc}}\left(h_{n},b_{n}\right)\mid\mathcal{X}_{n}\right]$
$\displaystyle=\nu!e_{\nu}^{\prime}\delta_{-,p}+h_{n}^{2+p-\nu}\frac{\Delta\mu_{-}^{(p+2)}}{(p+2)!}\mathcal{B}_{-,\nu,p,p+2}\left(h_{n}\right)\left\\{1+o_{p}(1)\right\\}$
$\displaystyle\quad-
h_{n}^{p+1-\nu}b_{n}^{q-p}\frac{\Delta\mu_{-}^{(q+1)}}{(q+1)!}\mathcal{B}_{-,p+1,q,q+1}\left(b_{n}\right)\frac{\mathcal{B}_{-,\nu,p,p+1}\left(h_{n}\right)}{(p+1)!}\left\\{1+o_{p}(1)\right\\}.$

(V) If $n\min\left\\{h_{n},b_{n}\right\\}\rightarrow\infty$, then
$\mathbb{V}\left[\Delta\hat{\mu}_{+,p,q}^{(\nu)\mathrm{bc}}\left(h_{n},b_{n}\right)\mid\mathcal{X}_{n}\right]=\mathcal{V}_{+,v,p,q}^{\mathrm{bc}}\left(h_{n},b_{n}\right)$,
where

$\displaystyle\mathcal{V}_{+,\nu,p,q}^{\mathrm{bc}}(h,b)=$
$\displaystyle\mathcal{V}_{+,\nu,p}(h)+h_{n}^{2(p+1-\nu)}\mathcal{V}_{+,p+1,q}(b)\frac{\mathcal{B}_{+,\nu,p,p+1}(h)^{2}}{(p+1)!^{2}}$
$\displaystyle-2h^{p+1-\nu}\mathcal{C}_{+,\nu,p,q}(h,b)\frac{\mathcal{B}_{+,v,p,p+1}(h)}{(p+1)!},$
$\displaystyle\mathcal{C}_{+,\nu,p,q}(h,b)=$
$\displaystyle\frac{1}{nh^{\nu}b^{p+1}}\nu!(p+1)!e_{\nu}^{\prime}\Gamma_{+,p}^{-1}(h)\Psi_{+,p,q}(h,b)\Gamma_{+,q}^{-1}(b)e_{p},$

and $\mathbb{V}\left[\Delta\hat{\mu}_{-,p,q}^{(\nu)\text{ bc
}}\left(h_{n},b_{n}\right)\mid\mathcal{X}_{n}\right]=\mathcal{V}_{-,\nu,p,q}^{\mathrm{bc}}\left(h_{n},b_{n}\right)$,
where

$\displaystyle\mathcal{V}_{-,\nu,p,q}^{\mathrm{bc}}(h,b)=$
$\displaystyle\mathcal{V}_{-,\nu,p}(h)+h_{n}^{2(p+1-\nu)}\mathcal{V}_{-,p+1,q}(b)\frac{\mathcal{B}_{-,\nu,p,p+1}(h)^{2}}{(p+1)!^{2}}$
$\displaystyle-2h^{p+1-\nu}\mathcal{C}_{-,\nu,p,q}(h,b)\frac{\mathcal{B}_{-,v,p,p+1}(h)}{(p+1)!}$
$\displaystyle\mathcal{C}_{-,\nu,p,q}(h,b)=$
$\displaystyle\frac{1}{nh^{\nu}b^{p+1}}\nu!(p+1)!e_{\nu}^{\prime}\Gamma_{-,p}^{-1}(h)\Psi_{-,p,q}(h,b)\Gamma_{-,q}^{-1}(b)e_{p}$

(D) If
$n\min\left\\{h_{n}^{2p+3},b_{n}^{2p+3}\right\\}\max\left\\{h_{n}^{2},b_{n}^{2(q-p)}\right\\}\rightarrow
0$, and $\kappa\max\left\\{h_{n},b_{n}\right\\}<\kappa_{0}$, then

$\displaystyle\frac{\Delta\hat{\mu}_{+,p,q}^{(\nu)\mathrm{bc}}\left(h_{n},b_{n}\right)-\nu!e_{\nu}^{\prime}\delta_{+,p}}{\sqrt{\mathcal{V}_{+,\nu,p,q}^{\mathrm{bc}}\left(h_{n},b_{n}\right)}}\rightarrow_{d}\mathcal{N}(0,1)$

and

$\displaystyle\frac{\Delta\hat{\mu}_{-,p,q}^{(\nu)\mathrm{bc}}\left(h_{n},b_{n}\right)-\nu!e_{\nu}^{\prime}\delta_{-,p}}{\sqrt{\mathcal{V}_{-,\nu,p,q}^{\mathrm{bc}}\left(h_{n},b_{n}\right)}}\rightarrow_{d}\mathcal{N}(0,1)$

###### Proof.

For part (B), first note that
$\mathbb{E}\left[\Delta\hat{\mu}_{+,p,q}^{(\nu)bc}\left(h_{n},b_{n}\right)\mid\mathcal{X}_{n}\right]=B_{1}-B_{2}$
with
$B_{1}=\mathbb{E}\left[\nu!e_{\nu}^{\prime}\hat{\delta}_{+,p}\left(h_{n}\right)\mid\mathcal{X}_{n}\right]$
and
$B_{2}=h_{n}^{p+1-\nu}\mathbb{E}\left[e_{p+1}^{\prime}\hat{\delta}_{+,q}\left(b_{n}\right)\mid\mathcal{X}_{n}\right]\mathcal{B}_{+,\nu,p}\left(h_{n}\right)$.
By Lemma B.2, with $s=\nu$ and $\ell=p$, we have

$\displaystyle B_{1}=$
$\displaystyle\nu!e_{\nu}^{\prime}\delta_{+,p}+h_{n}^{1+p-\nu}\frac{\Delta\mu_{+}^{(p+1)}}{(p+1)!}\mathcal{B}_{+,v,p,p+1}\left(h_{n}\right)$
$\displaystyle+h_{n}^{2+p-\nu}\frac{\Delta\mu_{+}^{(p+2)}}{(p+2)!}\mathcal{B}_{+,v,p,p+2}\left(h_{n}\right)+o_{p}\left(h_{n}^{2+p-\nu}\right)$

Similarly, by Lemma B.2, with $s=p+1$ and $\ell=q$, we have

$\displaystyle\mathbb{E}\left[(p+1)!e_{p+1}^{\prime}\hat{\delta}_{+,q}\left(b_{n}\right)\mid\mathcal{X}_{n}\right]$
$\displaystyle\quad=(p+1)!e_{p+1}^{\prime}\delta_{+,q}+b_{n}^{q-p}\frac{\Delta\mu_{+}^{(q+1)}}{(q+1)!}\mathcal{B}_{+,p+1,q,q+1}\left(b_{n}\right)+o_{p}\left(b_{n}^{q-p}\right)$

and hence

$\displaystyle B_{2}=$ $\displaystyle
h_{n}^{p+1-\nu}\mathbb{E}\left[(p+1)!e_{p+1}^{\prime}\hat{\delta}_{+,q}\left(b_{n}\right)\mid\mathcal{X}_{n}\right]\frac{\mathcal{B}_{+,v,p,p+1}\left(h_{n}\right)}{(p+1)!}$
$\displaystyle=$ $\displaystyle
h_{n}^{p+1-\nu}\left(e_{p+1}^{\prime}\delta_{+,q)}\right)\mathcal{B}_{+,\nu,p,p+1}\left(h_{n}\right)$
$\displaystyle+h_{n}^{p+1-\nu}b_{n}^{q-p}\frac{\Delta\mu_{+}^{(q+1)}}{(q+1)!}\mathcal{B}_{+,p+1,q,q+1}\left(b_{n}\right)\frac{\mathcal{B}_{+,v,p,p+1}\left(h_{n}\right)}{(p+1)!}$
$\displaystyle+h_{n}^{p+1-\nu}o_{p}\left(b_{n}^{q-p}\right)\mathcal{B}_{+,\nu,p,p+1}\left(h_{n}\right)$

Collecting terms, the result in part (B) follows:

$\displaystyle\mathbb{E}$
$\displaystyle{\left[\nu!e_{\nu}^{\prime}\hat{\delta}_{+,p,q}^{\mathrm{bc}}\left(h_{n},b_{n}\right)\mid\mathcal{X}_{n}\right]}$
$\displaystyle=$
$\displaystyle\nu!e_{\nu}^{\prime}\delta_{+,p}+h_{n}^{2+p-\nu}\frac{\Delta\mu_{+}^{(p+2)}}{(p+2)!}\mathcal{B}_{+,\nu,p,p+2}\left(h_{n}\right)\left\\{1+o_{p}(1)\right\\}$
$\displaystyle-
h_{n}^{p+1-\nu}b_{n}^{q-p}\frac{\Delta\mu_{+}^{(q+1)}}{(q+1)!}\mathcal{B}_{+,p+1,q,q+1}\left(b_{n}\right)\frac{\mathcal{B}_{+,\nu,p}\left(h_{n}\right)}{(p+1)!}\left\\{1+o_{p}(1)\right\\}$

For part $(\mathrm{V})$, first note that
$\mathbb{V}\left[\Delta\hat{\mu}_{+,p,q}^{(\nu)\text{ bc
}}\left(h_{n},b_{n}\right)\mid\mathcal{X}_{n}\right]=V_{1}+V_{2}-2C_{12}$
where, using Lemma B.2 with $s=\nu$ and $\ell=p$,

$\displaystyle
V_{1}=\mathbb{V}\left[\nu!e_{\nu}^{\prime}\hat{\delta}_{+,p}\left(h_{n}\right)\mid\mathcal{X}_{n}\right]=\mathbb{V}\left[\Delta\hat{\mu}_{+,p}^{(\nu)}\left(h_{n}\right)\mid\mathcal{X}_{n}\right]=\mathcal{V}_{+,v,p}\left(h_{n}\right)$

and, using LemmaB.2 with $s=p+1$ and $\ell=q$,

$\displaystyle V_{2}$
$\displaystyle=\mathbb{V}\left[h_{n}^{p+1-\nu}\left(e_{p+1}^{\prime}\hat{\delta}_{+,q}\left(b_{n}\right)\right)\mathcal{B}_{+,\nu,p,p+1}\left(h_{n}\right)\mid\mathcal{X}_{n}\right]$
$\displaystyle=h_{n}^{2(p+1-\nu)}\mathbb{V}\left[(p+1)!e_{p+1}^{\prime}\hat{\boldsymbol{\delta}}_{+,q}\left(b_{n}\right)\mid\mathcal{X}_{n}\right]\frac{\mathcal{B}_{+,\nu,p,p+1}\left(h_{n}\right)^{2}}{(p+1)!^{2}}$
$\displaystyle=h_{n}^{2(p+1-\nu)}\mathcal{V}_{+,p+1,q}\left(b_{n}\right)\frac{\mathcal{B}_{+,v,p,p+1}\left(h_{n}\right)^{2}}{(p+1)!^{2}}$

and

$\displaystyle C_{12}$
$\displaystyle=\mathbb{C}\left[\nu!e_{\nu}^{\prime}\hat{\delta}_{+,p}\left(h_{n}\right),h_{n}^{p+1-\nu}\left(e_{p+1}^{\prime}\hat{\delta}_{+,q}\left(b_{n}\right)\right)\mathcal{B}_{+,\nu,p,p+1}\left(h_{n}\right)\mid\mathcal{X}_{n}\right]$
$\displaystyle=h_{n}^{p+1-\nu}\mathbb{C}\left[\nu!e_{\nu}^{\prime}\hat{\delta}_{+,p}\left(h_{n}\right),(p+1)!e_{p+1}^{\prime}\hat{\delta}_{+,q}\left(b_{n}\right)\mid\mathcal{X}_{n}\right]\frac{\mathcal{B}_{+,\nu,p,p+1}\left(h_{n}\right)}{(p+1)!}$

with

$\displaystyle\mathbb{C}$
$\displaystyle{\left[e_{\nu}^{\prime}\hat{\delta}_{+,p}\left(h_{n}\right),e_{p+1}^{\prime}\hat{\delta}_{+,q}\left(b_{n}\right)\mid\mathcal{X}_{n}\right]}$
$\displaystyle=$ $\displaystyle
h_{n}^{-\nu}e_{\nu}^{\prime}\Gamma_{+,p}^{-1}\left(h_{n}\right)X_{p}\left(h_{n}\right)W_{+}\left(h_{n}\right)\mathbb{C}\left[Y,Y\mid\mathcal{X}_{n}\right]$
$\displaystyle\times
W_{+}\left(b_{n}\right)X_{q}\left(b_{n}\right)\Gamma_{+,q}^{-1}\left(b_{n}\right)e_{p+1}b_{n}^{-p-1}/n^{2}$
$\displaystyle=$
$\displaystyle\frac{1}{nh_{n}^{v}b_{n}^{p+1}}\nu!(p+1)!e_{\nu}^{\prime}\Gamma_{+,\ell}^{-1}\left(h_{n}\right)\Psi_{+,p,q}\left(h_{n},b_{n}\right)\Gamma_{+,q}^{-1}\left(b_{n}\right)e_{p+1}.$

Thus, collecting terms, we obtain the result in part (V).

Proof for the control group (subindex "-") is analogous. ∎

######  Lemma B.4.

Suppose Assumptions $1-2$ hold with $S\geq p+2$, and
$nh_{n}\rightarrow\infty$. Let $r\in\mathbb{N}$.

(B) If $h_{n}\rightarrow 0$, then

$\displaystyle\mathbb{E}\left[\hat{\tau}_{\nu,p}\left(h_{n}\right)\mid\mathcal{X}_{n}\right]=$
$\displaystyle\tau_{\nu}+h_{n}^{p+1-\nu}\mathbf{B}_{\nu,p,p+1}\left(h_{n}\right)$
$\displaystyle+h_{n}^{p+2-\nu}\mathbf{B}_{\nu,p,p+2}\left(h_{n}\right)+o_{p}\left(h_{n}^{p+2-\nu}\right)$

where

$\displaystyle\mathbf{B}_{\nu,p,r}\left(h_{n}\right)=\frac{\Delta\mu_{+}^{(r)}}{r!}\mathcal{B}_{+,\nu,p,r}\left(h_{n}\right)-\frac{\Delta\mu_{-}^{(r)}}{r!}\mathcal{B}_{-,\nu,p,r}\left(h_{n}\right)$
$\displaystyle\mathcal{B}_{+,v,p,r}\left(h_{n}\right)=\nu!e_{\nu}^{\prime}\Gamma_{+,p}^{-1}\left(h_{n}\right)\vartheta_{+,p,r}\left(h_{n}\right)=\nu!e_{\nu}^{\prime}\Gamma_{p}^{-1}\boldsymbol{\vartheta}_{p,r}+o_{p}(1)$
$\displaystyle\mathcal{B}_{-,\nu,p,r}\left(h_{n}\right)=\nu!e_{\nu}^{\prime}\Gamma_{-,p}^{-1}\left(h_{n}\right)\vartheta_{-,p,r}\left(h_{n}\right)=(-1)^{\nu+r}\nu!e_{\nu}^{\prime}\Gamma_{p}^{-1}\boldsymbol{\vartheta}_{p,r}+o_{p}(1)$

(V) If $h_{n}\rightarrow 0$, then
$\mathbf{V}_{\nu,p}\left(h_{n}\right)=\mathbb{V}\left[\hat{\tau}_{\nu,p}\left(h_{n}\right)\mid\mathcal{X}_{n}\right]=\mathcal{V}_{+,v,p}\left(h_{n}\right)+\mathcal{V}_{-,\nu,p}\left(h_{n}\right)$,
where

$\displaystyle\mathcal{V}_{+,\nu,p}\left(h_{n}\right)$
$\displaystyle=\frac{1}{nh_{n}^{2\nu}}\nu!^{2}e_{\nu}^{\prime}\Gamma_{+,p}^{-1}\left(h_{n}\right)\Psi_{+,p}\left(h_{n}\right)\Gamma_{+,p}^{-1}\left(h_{n}\right)e_{\nu}$
$\displaystyle=\frac{1}{nh_{n}^{1+2\nu}}\frac{\sigma_{+}^{2}}{f}\nu!^{2}e_{\nu}^{\prime}\Gamma_{p}^{-1}\Psi_{p}\Gamma_{p}^{-1}e_{\nu}\left\\{1+o_{p}(1)\right\\}$
$\displaystyle\mathcal{V}_{-,v,p}\left(h_{n}\right)$
$\displaystyle=\frac{1}{nh_{n}^{2\nu}}\nu!^{2}e_{\nu}^{\prime}\Gamma_{-,p}^{-1}\left(h_{n}\right)\Psi_{-,p}\left(h_{n}\right)\Gamma_{-,p}^{-1}\left(h_{n}\right)e_{\nu}$
$\displaystyle=\frac{1}{nh_{n}^{1+2\nu}}\frac{\sigma_{-}^{2}}{f}\nu!^{2}e_{\nu}^{\prime}\Gamma_{p}^{-1}\Psi_{p}\Gamma_{p}^{-1}e_{\nu}\left\\{1+o_{p}(1)\right\\}$

(D) If $nh_{n}^{2p+5}\rightarrow 0$, then

$\displaystyle\frac{\hat{\tau}_{\nu,p}\left(h_{n}\right)-\tau_{\nu}-h_{n}^{p+1-\nu}\mathbf{B}_{\nu,p,p+1}\left(h_{n}\right)}{\sqrt{\mathbf{V}_{\nu,p}\left(h_{n}\right)}}\rightarrow_{d}\mathcal{N}(0,1)$

###### Proof.

Part (B) follows immediately from Lemma B.2(B), its analogue for the left-side
estimator
$\left(s!e_{s}^{\prime}\hat{\beta}_{-,\ell}\left(h_{n}\right)\right)$, and the
linearity of conditional expectations. Part (V) also follows immediately from
Lemma B.2(V), its analogue for the left-side estimator
$\left(s!e_{s}^{\prime}\hat{\beta}_{-,\ell}\left(h_{n}\right)\right)$, and the
conditional independence of observations at either side of the threshold
$(x=0)$. Finally, part (D) follows by the same argument given in the proof of
Lemma B.2(D), but now applied to the estimator
$\hat{\tau}_{\nu,p}\left(h_{n}\right)=\Delta\hat{\mu}_{+,p}^{(\nu)}\left(h_{n}\right)-\Delta\hat{\mu}_{-,p}^{(\nu)}\left(h_{n}\right)=\nu!e_{\nu}^{\prime}\hat{\beta}_{+,p}\left(h_{n}\right)-\nu!e_{\nu}^{\prime}\hat{\beta}_{-,p}\left(h_{n}\right)$.
This completes the proof. ∎

##  Appendix C Proofs of main results

We provide the proofs of our results in this appendix.

###### Proof of Claim 1.

Given Assumptions 1 and 2, we can prove that
$\tau^{DiDC}=\tau_{1}^{RDD}-\tau_{0}^{RDD}=\Delta Y^{+}-\Delta Y^{-}$.

  1. A)

Difference of RDDs. From 1 and 2 we can derive two RDDs estimands, one for
$t=0$ and one for $t=1$.

$\displaystyle\tau^{RDD}_{1}$ $\displaystyle=Y_{1}^{+}-Y_{1}^{-}$
$\displaystyle=E\left[Y_{i,1}\left(1,1\right)-Y_{i,1}\left(0,1\right)-Y_{i,1}\left(1,0\right)+Y_{i,1}\left(0,0\right)|Z_{i}=z_{0}\right]\left(D_{0}^{+}D_{1}^{+}-D_{0}^{-}D_{1}^{-}\right)+$
$\displaystyle\quad\quad
E\left[Y_{i,1}\left(0,1\right)-Y_{i,1}\left(0,0\right)|Z_{i}=z_{0}\right]\left(D_{1}^{+}-D_{1}^{-}\right)+$
(C.1) $\displaystyle\quad\quad
E\left[Y_{i,0}\left(1,0\right)-Y_{i,1}\left(0,0\right)|Z_{i}=z_{0}\right]\left(D_{0}^{+}-D_{0}^{-}\right)$
$\displaystyle\tau^{RDD}_{0}$ $\displaystyle=Y_{0}^{+}-Y0^{-}$
$\displaystyle=E\left[Y_{i,0}\left(1,0\right)-Y_{i,0}\left(0,0\right)|Z_{i}=z_{0}\right]\left(D_{0}^{+}-D_{0}^{-}\right)$
(C.2)

where $X^{+}=\lim_{\epsilon\rightarrow 0}E(X_{i}|Z_{i}=z_{0}+\epsilon)$ and
$X^{-}=\lim_{\epsilon\rightarrow 0}E(X_{i}|Z_{i}=z_{0}-\epsilon)$. By taking
the difference between the equations C.1 and C.2 and under assumptions 1-4, we
arrive at the differences-in-discontinuity estimand:

$\displaystyle\tau^{DiDC}$ $\displaystyle=\tau^{RDD}_{1}-\tau^{RDD}_{0}$
$\displaystyle=E\left[\left(Y_{i,1}(1,1)-Y_{i,1}(0,1)\right)\left(Y_{i,1}(1,0)-Y_{i,1}(0,0)\right)|Z_{i}=z_{0}\right]\left(D_{0}^{+}D_{1}^{+}-D_{0}^{-}D_{1}^{-}\right)+$
$\displaystyle\quad
E\left[Y_{i,1}(0,1)-Y_{i,1}(0,0)|Z_{i}=z_{0}\right]\left(D_{1}^{+}-D_{1}^{-}\right)+$
$\displaystyle\quad
E\left[\left(Y_{i,1}(1,0)-Y_{i,1}(0,0)\right)\left(Y_{i,0}(1,0)-Y_{i,0}(0,0)\right)|Z_{i}=z_{0}\right]\left(D_{0}^{+}-D_{0}^{-}\right)$

  2. B)

RDD of the differences. Using equation3 and taking the difference of the
limits above and below the threshold:

$\displaystyle\tau^{DiDC}$ $\displaystyle=\Delta Y^{+}-\Delta Y^{-}$
$\displaystyle=E\left[\left(Y_{i,1}(1,1)-Y_{i,1}(0,1)\right)\left(Y_{i,1}(1,0)-Y_{i,1}(0,0)\right)|Z_{i}=z_{0}\right]\left(D_{0}^{+}D_{1}^{+}-D_{0}^{-}D_{1}^{-}\right)+$
$\displaystyle\quad
E\left[Y_{i,1}(0,1)-Y_{i,1}(0,0)|Z_{i}=z_{0}\right]\left(D_{1}^{+}-D_{1}^{-}\right)+$
$\displaystyle\quad
E\left[\left(Y_{i,1}(1,0)-Y_{i,1}(0,0)\right)\left(Y_{i,0}(1,0)-Y_{i,0}(0,0)\right)|Z_{i}=z_{0}\right]\left(D_{0}^{+}-D_{0}^{-}\right)$

where $\Delta X^{+}=\lim_{\epsilon\rightarrow 0}E(\Delta
X_{i}|Z_{i}=z_{0}+\epsilon)$ and $\Delta X^{-}=\lim_{\epsilon\rightarrow
0}E(\Delta X_{i}|Z_{i}=z_{0}-\epsilon)$.

Therefore $\tau^{DiDC}=\tau^{RDD}_{1}-\tau^{RDD}_{0}=\Delta Y^{+}-\Delta
Y^{-}$. ∎

###### Proof.

Proof of Lemma 1

Recall the definition
$MSE_{\nu,p,s}(h_{n})=E\left[\left(\hat{\tau}_{\nu,p,s}(h_{n})-\tau_{\nu,p}\right)^{2}|\chi_{n}\right]$.
We can rewrite as

$\displaystyle MSE_{\nu,p,s}(h_{n})=$ $\displaystyle
V\left[\hat{\tau}_{\nu,p,s}(h_{n})|\chi_{n}\right]+\left(E\left[\hat{\tau}_{\nu,p,s}(h_{n})-\tau_{\nu,p}|\chi_{n}\right]\right)^{2}$

Then, from Lemma B.4,

$\displaystyle V\left[\hat{\tau}_{\nu,p,s}(h_{n})|\chi_{n}\right]$
$\displaystyle=\frac{\nu!}{nh_{n}^{1+2\nu}}\frac{\sigma_{+}^{2}-\sigma_{-}^{2}}{f}e^{\prime}_{\nu}\Gamma_{p}^{-1}\Psi_{p}\Gamma_{p}^{-1}\\{1+o_{p}(1)\\}$

and

$\displaystyle
E\left[\hat{\tau}_{\nu,p,s}(h_{n})-\tau_{\nu,p}|\chi_{n}\right]=h_{n}^{1+p-\nu}\frac{\Delta\mu^{(p+1)}_{+}-(-1)^{\nu+p+s}\Delta\mu^{(p+1)}_{-}}{(p+1)!}\nu!e^{\prime}_{\nu}\Gamma_{p}^{-1}\varphi_{p,p+1}\\{1+o_{p}(1)\\}$
(C.3)

so we can rewrite

$\displaystyle
MSE_{\nu,p,s}(h_{n})=\frac{\mathbb{V_{\nu,p}}}{nh_{n}^{1+2\nu}}\\{1+o_{p}(1)\\}+\left(h_{n}^{1+p-\nu}\mathbb{B_{\nu,p,p+1,s}}\\{1+o_{p}(1)\\}\right)^{2}$

where
$\mathbb{V_{\nu,p}}=\nu!\frac{\sigma_{+}^{2}-\sigma_{-}^{2}}{f}e^{\prime}_{\nu}\Gamma_{p}^{-1}\Psi_{p}\Gamma_{p}^{-1}$
and
$\mathbb{B_{\nu,p,p+1,s}}=\frac{\Delta\mu^{(p+1)}_{+}-(-1)^{\nu+p+s}\Delta\mu^{(p+1)}_{-}}{(p+1)!}\nu!e^{\prime}_{\nu}\Gamma_{p}^{-1}\varphi_{p,p+1}$.
The MSE-optimal bandwidth is

$\displaystyle h_{n}^{MSE}{}_{\nu,p}$
$\displaystyle=arg\min_{h_{n}}MSE_{\nu,p,s}(h_{n})$
$\displaystyle=\left(\frac{\left(1+2\nu\right)\mathbf{V}_{\nu,p}}{n2\left(1+p-\nu\right)\textbf{B}^{2}_{\nu,p,p+1,s}}\right)^{\frac{1}{2p+3}}$

∎

###### Proof.

Proof of Claim 2

We define higher-order derivatives notation of the unknown regression
functions as:

$\displaystyle\Delta\mu_{+}^{(\nu)}(z)=\frac{d^{\nu}\Delta\mu_{+}(z)}{dz^{\nu}},$
$\displaystyle\Delta\mu_{-}^{(\nu)}(z)=\frac{d^{\nu}\Delta\mu_{(}z)}{dz^{\nu}}$
$\displaystyle\Delta\mu_{+}^{(\nu)}=\lim_{x\rightarrow
0^{+}}{\Delta\mu_{+}^{(\nu)}(z)},$
$\displaystyle\Delta\mu_{-}^{(\nu)}=\lim_{x\rightarrow
0^{-}}{\Delta\mu_{-}^{(\nu)}(z)}$
$\displaystyle\Delta\mu_{+}=\lim_{x\rightarrow 0^{+}}{\Delta\mu(z)},$
$\displaystyle\Delta\mu_{-}=\lim_{x\rightarrow 0^{=}}{\Delta\mu(z)}$

where $\Delta\mu(z)=E\left(\Delta
Y_{i}|Z_{i}=z\right)=E\left(Y_{i,1}-Y_{i,0}|Z_{i}=z\right)$.

The notation for the RDD at time $t$ where
$\mu_{t}(z)=E\left(Y_{i,t}|Z_{i}=z\right)$ is

$\displaystyle\mu_{t,+}^{(\nu)}(z)=\frac{d^{\nu}\mu_{t,+}(z)}{dz^{\nu}},$
$\displaystyle\mu_{t,-}^{(\nu)}(z)=\frac{d^{\nu}\mu_{(}z)}{dz^{\nu}}$
$\displaystyle\mu_{t,+}^{(\nu)}=\lim_{x\rightarrow
0^{+}}{\mu_{t,+}^{(\nu)}(z)},$
$\displaystyle\mu_{t,-}^{(\nu)}=\lim_{x\rightarrow
0^{-}}{\mu_{t,-}^{(\nu)}(z)}$ $\displaystyle\mu_{t,+}=\lim_{x\rightarrow
0^{+}}{\mu(z)},$ $\displaystyle\mu_{t,-}=\lim_{x\rightarrow 0^{=}}{\mu(z)}$

Then we derive the following:

$\displaystyle\Delta\mu_{+}^{(\nu)}=\mu_{1,+}^{(\nu)}-\mu_{0,+}^{(\nu)}$
$\displaystyle\Delta\mu_{-}^{(\nu)}=\mu_{1,}^{(\nu)}-\mu_{0,+}^{(\nu)}$ (C.4)

It follows that
$B\left[\hat{\tau}^{DiDC}_{\nu,p}(h_{n})\right]=B\left[\hat{\tau}^{RDD}_{1,\nu,p}(h_{n})\right]-B\left[\hat{\tau}^{RDD}_{0,\nu,p}(h_{n})\right]$
∎

###### Proof.

Proof of Calim 3

$\displaystyle\Delta\mu(z)=$ $\displaystyle E\left[\Delta
Y_{i}|Z_{i}=z\right]$ $\displaystyle=$ $\displaystyle
E\left[Y_{i,1}-Y_{i,0}|Z_{i}=z\right]$ $\displaystyle=$ $\displaystyle
E\left[Y_{i,1}|Z_{i}=z\right]-E\left[Y_{i,0}|Z_{i}=z\right]$ $\displaystyle=$
$\displaystyle\mu_{1}(z)-\mu_{0}(z)$

Therefore,

$\displaystyle\Delta\mu_{+}^{(\nu)}=$
$\displaystyle\mu_{+,1}^{(\nu)}-\mu_{+,0}^{(\nu)}$

Then, if the shapes of the DGPs are time-invariant for each side of the
threshold, then $\mu_{+,1}^{(\nu)}-\mu_{+,0}^{(\nu)}$ and
$\mu_{-,1}^{(\nu)}-\mu_{-,0}^{(\nu)}$, leading to
$\Delta\mu_{+}^{(\nu)}=\Delta\mu_{-}^{(\nu)}=0$ ∎

##  Appendix D Simulations

###  D.1 Bandwidth for Testing if Confounding Effects are Constant

This appendix details into the critical issue of bandwidth selection the Wald
test to test Assumption3. As mentioned before, this assumption is vital for
the difference-in-discontinuitis method, as the non validity means all
estimations would be unvalid. The procedure for testing involves considering
the following stacked RDDs regression model:

$\displaystyle
Y_{i}=\alpha_{i}+\sum_{k=1}^{K}\left[\beta_{-k}\left(Z_{i}-z_{0}\right)+\theta_{-k}D\left(Z_{i}-z_{0}\right)\right]T_{-k}+\varepsilon_{i}$
(8)

where $T_{-k}$ is a dummy indicating which period the data belongs to, that
is, the RDD of which period it’s running. The null hypothesis for the Wald
test is:

$\displaystyle H_{0}:\quad\theta_{0}=\theta_{-1}=...=\theta_{-K}$

One central question arises: how should we determine the optimal bandwidth for
Equation 8. In this section, we outline the methodology used for bandwidth
selection and present simulation results in two different scenarios: one with
identical Data Generating Processes for pre- and post-treatment periods and
another with different DGPs, both while considering confounding factors. The
first scenario aligns with the framework presented in Model 1, as detailed in
Appendix D.2.1 while the second encompasses both models 1 and 3 (Appendix
D.2.1).

To select bandwidths, we employ the `rdrobust` package and estimate bandwidths
separately for each RDD. Next, we cross the bandwidths and conduct RDD
regressions with each of them, evaluating the bias and root-mean-squared
errors (RMSE) of the estimated parameters. The chosen bandwidth are the one
that minimizes RMSE, alongside the largest and smallest bandwidths, and the
CCT (Calonico et al., 2014) bandwidth for all observations regardless of time
dummies.

The various bandwidths selected for produced similar results. This finding
suggests that the choice of bandwidth does not substantially affect the
accuracy of the test. It is essential to note that while this results says
this, the choice of bandwidth for the DiDC estimation of the treatment effect
remains a crucial consideration of the method.

###  D.2 Data Generating Processes

We provide further details on each of the data generating processes (DGPs)
employed in our simulation studies. Both GDPs employ the same simulation setup
based on model 3 from Calonico et al. (2014), with the key difference being
the presence of time-fixed confounding factors at the threshold in one of the
settings. We perform 1000 replications for each simulation and for each
replication, we use a sample size of $n=500$, with

$\displaystyle Y_{i,t}$ $\displaystyle=\mu_{it}(Z_{i})+\varepsilon_{i,t},\quad
Z_{i}\sim(2\mathcal{B}(2,4)-1),$ $\displaystyle\quad\varepsilon_{it}\sim
N(0,\sigma^{2}_{\varepsilon})$

where $\mathcal{B}(p_{1},p_{2})$ is a beta distribution with parameters
$p_{1}$ and $p_{2}$, $\sigma_{\varepsilon}=0.1295$.

####  D.2.1 Model 1: Identical DGPs with time-invariant confounders at the
threshold

This model includes a time-invariant confounder, c, to represent pre-existing
discontinuities at the threshold that affect the outcome $Y_{i}$. Here, units
with $Z_{i}\geq 0$ receive a treatment c and units with $Z_{i}<0$ do not. At
the same time, units with $Z_{i}\geq 0$ have a different data generating
function than units with $Z_{i}<0$. The data generating functions $\mu_{0}(z)$
and $\mu_{1}(z)$, for periods $t=0$ and $t=1$ respectively are as follows

$\displaystyle\mu_{0}(z)$ $\displaystyle=\begin{cases}0.48+1.27z-0.5\cdot
7.18z^{2}\vspace{-0.5cm}\\\ \quad+0.7\cdot 20.21z^{3}+1.1\cdot
21.54z^{4}+1.5\cdot 7.33z^{5}&\text{if }z<0\\\ 0.52+0.84z-0.1\cdot
3z^{2}\vspace{-0.5cm}\\\ \quad-0.3\cdot 7.99z^{3}-0.1\cdot
9.01z^{4}+3.56z^{5}+\textbf{c}&\text{if }z\geq 0\end{cases}$
$\displaystyle\mu_{1}(z)$ $\displaystyle=\begin{cases}0.48+1.27z-0.5\cdot
7.18z^{2}\vspace{-0.5cm}\\\ \quad+0.7\cdot 20.21z^{3}+1.1\cdot
21.54z^{4}+1.5\cdot 7.33z^{5}&\text{if }z<0\\\ 0.52+0.84z-0.1\cdot
3z^{2}\vspace{-0.5cm}\\\ \quad-0.3\cdot 7.99z^{3}-0.1\cdot
9.01z^{4}+3.56z^{5}+\textbf{c}+\tau&\text{if }z\geq 0\end{cases}$

####  D.2.2 Model 2: Identical DGPs

This model follows precisely model 3 from Calonico et al. (2014) for both
periods, the only difference between them being that in period $t=1$ there is
the effect of the treatment $\tau$ for units with $Z-i\geq 0$. No other
discontinuities that could affect the outcome exist at the threshold
$z_{0}=0$. The data generating functions $\mu_{0}(z)$ and $\mu_{1}(z)$ are as
follows

$\displaystyle\mu_{0}(z)$ $\displaystyle=\begin{cases}0.48+1.27z-0.5\cdot
7.18z^{2}\vspace{-0.5cm}\\\ \quad+0.7\cdot 20.21z^{3}+1.1\cdot
21.54z^{4}+1.5\cdot 7.33z^{5}&\text{if }z<0\\\ 0.52+0.84z-0.1\cdot
3z^{2}\vspace{-0.5cm}\\\ \quad-0.3\cdot 7.99z^{3}-0.1\cdot
9.01z^{4}+3.56z^{5}&\text{if }z\geq 0\end{cases}$ $\displaystyle\mu_{1}(z)$
$\displaystyle=\begin{cases}0.48+1.27z-0.5\cdot 7.18z^{2}\vspace{-0.5cm}\\\
\quad+0.7\cdot 20.21z^{3}+1.1\cdot 21.54z^{4}+1.5\cdot 7.33z^{5}&\text{if
}z<0\\\ 0.52+0.84z-0.1\cdot 3z^{2}\vspace{-0.5cm}\\\ \quad-0.3\cdot
7.99z^{3}-0.1\cdot 9.01z^{4}+3.56z^{5}+\tau&\text{if }z\geq 0\end{cases}$

####  D.2.3 Model 3: Different DGPs with time-invariant confounders at the
threshold

This model follows model 3 from Calonico et al. (2014) and includes a time-
invariant confounder, c, to represent pre-existing discontinuities at the
threshold that affect the outcome $Y_{i}$ at $t=0$, and a linear model for
time $t=1$. The data generating functions $\mu_{0}(z)$ and $\mu_{1}(z)$, for
periods $t=0$ and $t=1$ respectively are as follows

$\displaystyle\mu_{0}(z)$ $\displaystyle=\begin{cases}0.48+1.27z-0.5\cdot
7.18z^{2}\vspace{-0.5cm}\\\ \quad+0.7\cdot 20.21z^{3}+1.1\cdot
21.54z^{4}+1.5\cdot 7.33z^{5}&\text{if }z<0\\\ 0.52+0.84z-0.1\cdot
3z^{2}\vspace{-0.5cm}\\\ \quad-0.3\cdot 7.99z^{3}-0.1\cdot
9.01z^{4}+3.56z^{5}+\textbf{c}&\text{if }z\geq 0\end{cases}$
$\displaystyle\mu_{1}(z)$ $\displaystyle=\begin{cases}0.48+1.4z&\text{if
}z<0\\\ 0.52+0.1z+\textbf{c}+\tau&\text{if }z\geq 0\end{cases}$

####  D.2.4 Model 4: Different DGPs

Here the model for $t=0$ follows precisely model 3 from Calonico et al. (2014)
and the model for $t=1$ is linear. No discontinuities that could affect the
outcome exist at the threshold $z_{0}=0$ other than the treatment effect
$\tau$ introduced for $t=1$. The data generating functions $\mu_{0}(z)$ and
$\mu_{1}(z)$ are as follows

$\displaystyle\mu_{0}(z)$ $\displaystyle=\begin{cases}0.48+1.27z-0.5\cdot
7.18z^{2}\vspace{-0.5cm}\\\ \quad+0.7\cdot 20.21z^{3}+1.1\cdot
21.54z^{4}+1.5\cdot 7.33z^{5}&\text{if }z<0\\\ 0.52+0.84z-0.1\cdot
3z^{2}\vspace{-0.5cm}\\\ \quad-0.3\cdot 7.99z^{3}-0.1\cdot
9.01z^{4}+3.56z^{5}+\textbf{c}&\text{if }z\geq 0\end{cases}$
$\displaystyle\mu_{1}(z)$ $\displaystyle=\begin{cases}0.48+1.4z&\text{if
}z<0\\\ 0.52+0.1z+\tau&\text{if }z\geq 0\end{cases}$

##  Appendix E Empirical illustration

We replicate the analysis conducted by Grembi et al. (2016) once more, this
time using the bandwidth selection method proposed by Ludwig and Miller (2007)
as she did. The use of these alternative bandwidths does not significantly
alter the estimated treatment effects; they remain similar in magnitude.

Table E.1: Effects of relaxing a fiscal rule - Ludwig & Miller (2007) Bandwidths Estimators | Deficit | Fiscal Gap | Taxes
---|---|---|---
Diif-in-Disc (Grembi et al, 2016) |  | 9.454
---
(4.343)
| 48.469
---
(23.315)
| -34.748
---
(20.166)
Bandwidth | 1498 | 833 | 684
Observations | 5858 | 3438 | 1536
Diff-in-Disc (RDD of $\Delta$) |  | 16.979
---
(8.859)
| 35.302
---
(12.449)
| -16.958
---
(7.473)
Bandwidth | 887 | 908 | 1047
Observations | 1379 | 1278 | 1254
Diff-in-Disc ($\Delta$ of RDDs) |  | 17.895
---
(4.713)
| 39.585
---
(12.729)
| -27.620
---
(10.761)
Bandwidth | 961 | 647 | 695
Observations | 1163 | 1156 | 1246