Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
71,240 characters · 9 sections · 53 citation commands
Non-parametric Causal Inference in Dynamic Thresholding Designs
\allowdisplaybreaks
Dynamic threshold-based rules \begingroup \footnote{Draft version \ifcase\month\or Jan\or Feb\or Mar\or Apr\or May\or Jun\or Jul\or Aug\or Sep\or Oct\or Nov\or Dec\fi \space \number\year . This research was supported by the Office of Naval Research under grant number N00014-24-1-2091.} \addtocounter{footnote}{-1} \endgroup
govern many consequential decisions in healthcare, education, credit markets, and public policy---and present numerous opportunities for policy evaluation. For example, IIZUKA2021 study benefits of health signals using data from the Japanese healthcare system, where health signals are driven by threshold rules. Patients receive yearly health checkups at which their fasting blood sugar (FBS) is measured; they are diagnosed as pre-diabetic and recommended lifestyle interventions along with follow-up care if their FBS level crosses 110 mg/dL. The main research question in IIZUKA2021 is whether such pre-diabetic health signals are efficient from a public health perspective relative to the cost of the induced follow-up care.
The goal of this paper is to develop methods for non-parametric policy evaluation in such dynamic thresholding designs. The fact that thresholding designs open the door to non-parametric causal inference has been recognized for a long time thistlethwaite1960regression, and recent decades have seen a flurry of work on regression discontinuity designs following this insight Hahn-et-al-2001,IL2008review,CCT2014robustCI,ArmstrongKolesar2018. Existing work on regression discontinuity designs, however, are focused on cross-sectional settings where each unit only receives treatment once and then experiences an outcome---and so are amenable to causal analysis using the basic potential outcomes model imbens2015causal. In contrast, we are interested in settings where each unit is eligible for treatment multiple times (e.g., in the setting of IIZUKA2021, each patient gets a new pre-diabetes diagnosis each year), thus resulting in complex treatment dynamics that need to be modeled in order to achieve correct inferences about the overall effect of a policy robins1986new. For example, if prescribing lifestyle interventions is effective in lowering FBS, then having a patient be diagnosed as prediabetic one year may make them less likely to receive the same diagnosis in subsequent years.
Here, we propose a framework for analyzing thresholding designs all while non-parametrically accounting for dynamics as they arise in the model of robins1986new. Our main finding is that, despite the apparent complexity of accounting to flexible dynamics, we are able to use a carefully tailored local linear regression estimator to consistently estimate the marginal policy effect of the thresholding rule carneiro2010evaluating. In the setting of IIZUKA2021, this marginal policy effect corresponds to the net-present health benefit of infinitesimally lowering the FBS cutoff for prediabetes, scaled by the discounted increase in the number of prediabetes diagnoses from lowering the cutoff. Our approach draws from the literature on reinforcement learning and Markov decision processes SuttonBarto2018; our analysis is in particular motivated by the policy gradient theorem and the work of sutton1999policy on first-order optimization of reinforcement learning models.
As background for our results on causal inference in dynamic threshold designs, we first briefly review standard regression discontinuity {(RD)} designs (or, what could be called cross-sectional thresholding designs) and how the resulting estimand can be interpreted as a marginal policy effect. Following IL2008review, assume that we have data on IID-sampled pairs $(Z_i,\, Y_i)$ for units $i = 1, \, \ldots, \, n$, where $Z_i \in \mathbb{R}$ is the running variable and $Y_i \in \mathbb{R}$ is the outcome of interest.
The sharp {RD} design assumes that there exists a cutoff $c \in \mathbb{R}$ such that treatment is assigned as $A_i = \mathbf{1}\left\{Z_i\ge c\right\}$. We posit potential outcomes $\{Y_i(0), \, Y_i(1)\}$ such that $Y_i = Y_i(A_i)$, and define the conditional average treatment effect (CATE) as $\tau(z) = \mathbb{E}[Y_i(1) - Y_i(0) \, |\, Z_i = z]$. Then, the sharp RD design enables us to estimate $\tau(z)$, i.e., the CATE at the cutoff, as a discontinuity in the conditional response surface at $Z_i = c$ Hahn-et-al-2001,
provided that the $\mu_{a}(z) = \mathbb{E}[Y_i(a) \, |\, Z_i = z]$ are continuous in $z$ and that the running variable has continuous support around $c$.
The classical interpretation of the RD estimand as the CATE for units whose running variable $Z_i$ straddles the cutoff relies crucially on $Z_i$ being causally prior to any actions induced by our thresholding policy. In dynamic thresholding designs, however, past thresholding actions can affect future values of the running variable; and, as argued in frangakis2002principal, conditioning on observed variables whose value may be affected by treatment generally precludes causal interpretation of resulting estimands. To avoid this issue, we find it helpful to re-interpret the classical RD estimand as a marginal policy effect in the sense of carneiro2010evaluating, i.e., as essentially the answer to a cost-benefit analysis. Given a threshold $c$, the sharp RD design with threshold introduces a treatment policy $\pi_c(Z_i)=\mathbf{1}\left\{Z_i\ge c\right\}$ with associated policy value
Then, provided the running variable is exogenous to the policy cutoff, $\tau_{\,\mathrm{RD}}$ can be interpreted as the policy gradient of lowering the cutoff (i.e., of treating more units), divided by the corresponding increase in the number of units treated as we lower the cutoff. Further results on policy counterfactuals in thresholding designs are given in dong2015identifying.
In other words, provided the running variables $Z_i$ are exogenous to the thresholding policy, and if the cost of providing treatment to a unit is $\lambda$, then a social planner could achieve cost-adjusted welfare benefits by reducing $c$ (and thus marginally increasing the treatment rate) if and only if $\tau_{\,\mathrm{RD}} > \lambda$. As we move to a multi-period setting, we will find this alternative characterization of the RD estimand as the solution to a cost-benefit analysis to be remarkably resilient to challenges induced by treatment dynamics.
Now consider a setting where units are observed at times $t = 0, \, 1, \, 2, \, \ldots, \, T$, where the horizon $T$ may be either finite or infinite\footnote{In the infinite-horizon setting, we use $T=\infty$ to define the population quantity of interest. However, the methods we propose for estimation and inference on this estimand are designed to operate in the realistic setting where the observed trajectories are finite but long (i.e., when $T$ is sufficiently large for our asymptotic results to provide useful approximations).}. At each time period $t$ we observe a running variable $Z_{i,t} \in \mathbb{R} \, \cup \, \{-\infty\}$,\footnote{We allow for the case $Z_{i,t} = -\infty$ to account the possibility that the running variable may not be observed in every time period hsu2024dynamic. For example, in the case of yearly health checkups, it's possible a patient misses their health checkup one year and so no health measurements are taken. We assume that units are not treated in periods where the running variable is unobserved.} take a thresholding action $A_{i,t} = \mathbf{1}\left\{Z_{i,t}\ge c\right\}$ for some threshold $c \in \mathbb{R}$, and observe an outcome $Y_{i,t} \in \mathbb{R}$. The causal structure of dynamic problems is considerably richer than in cross-sections ones: In addition to affecting outcomes $Y_{i,t}$, actions $A_{i,t}$ taken at time $t$ can affect state---and thus also actions---at all times $t' > t$.
The induced potential outcomes then acquire a tree-like branching structure indexing over all possible past treatment sequences robins1986new. This branching makes a direct reduced-form approach to dynamic thresholding designs intractable---or, at the very least, subject to an exponential blow-up in dimensionality as the time horizon (and thus action space) grows. Instead, it is usually more fruitful to proceed via what Robins refers to as the $g$-formula which provides a useful factorization for the observed-data distribution under natural temporal consistency assumptions. The probability factorization in Robins' $g$-formula is equivalent to what arises in the study of Markov decision processes SuttonBarto2018; recent textbook discussions are given in HernanRobins2020 and wager2024causal. Throughout, we will assume that conditions required for the $g$-formula to hold are satisfied.
Our main question of interest is how treatment---as determined by dynamic thresholding as in Assumption (ref)---affects net-present expected welfare and treatment frequency,
where $0 < \gamma \leq 1$ is a discount rate (if $T = \infty$ then we must have $\gamma < 1$). Because of the branching structure of potential outcomes a direct analogue to (ref) does not immediately enable meaningful program evaluation. However, perhaps surprisingly, we will find that marginal policy effect characterizations of the form (ref) remain useful: An RD estimand defined as
can still be effectively estimated in a sharp RD design---and can still be used to resolve policy-relevant cost-benefit tradeoffs.
The modern literature on regression discontinuity designs goes back to Hahn-et-al-2001; influential contributions to this literature include IL2008review, imbens2012optimal, CCT2014robustCI and ArmstrongKolesar2018. Most of the existing methodological literature on regression discontinuity designs, however, is focused on the cross-sectional setting where treatment is only assigned once. And, when faced with the longitudinal setting, empirical researchers have reduced the problem to a cross-sectional setting by simply considering various reduced-form regression discontinuities, e.g., by running a standard RDD of $Y_t$ on $Z_t$ or of $Y_{t+1}$ on $Z_t$. This is, for example, the strategy taken in the original analysis of IIZUKA2021. Such reduced form analyses can be of considerable substantive interest in applications; however, they do not capture full treatment dynamics (e.g., how actions taken in one period may change the running variable---and thus actions---taken in subsequent ones), and are thus not directly interpretable as policy-relevant treatment effects heckman2016dynamic.
One notable exception is cellini2010value, who use regression discontinuities to identify a type of treatment on the treated (ATT) effect in dynamic designs. They then use their estimator to identify the effect of local school spending via public bonds on house prices in California by comparing outcomes in school districts where bond measures are just barely accepts vs. rejected by voters; and their estimator allows them to formally consider the fact that approving bonds in the past makes it less likely that additional bonds will be approved in the future. The approach of cellini2010value, however, makes crucial use of a linear parametric model whereby
i.e., treatment effects are homogeneous and decay uniformly over time. And, as shown by hsu2024dynamic, their approach no longer recovers an ATT if we allow for treatment heterogeneity. hsu2024dynamic propose an alternative analysis that avoids (ref). But they in turn require a strong conditional mean independence assumption (CIA) which, e.g., in the 2-period case requires that the time-2 control potential outcomes be independent of the time-2 running variable for all units whose time-1 running variable is near the cutoff.\footnote{See Assumption 3.1.2 of hsu2024dynamic for a precise statement. This assumption is substantive, and would not hold in generic dynamic thresholding designs. In particular, the CIA assumption will generally not hold if control potential outcomes and the running variable both vary smoothly with some time-varying latent confounder; e.g., in our motivating example, it would generally not hold if FBS and health outcomes both vary smoothly with unobserved and time-varying health-seeking behaviors.} To the best of our knowledge our paper is the first to provide results on non-parametric causal inference for dynamic thresholding designs with generality that's comparable to standard results in the cross-sectional setting following Hahn-et-al-2001.
Our flexible potential-outcomes based model for dynamic causal inference goes back to robins1986new. This model is widely used in biostatistics in the context of, e.g., marginal structural models robins2000marginal and optimal treatment regimes robins2004optimal. To the best of our knowledge, this model has not been previously used in the context of dynamic thresholding designs---the one exception being hsu2024dynamic, who pair the model of robins1986new with their potentially restrictive CIA assumption to make progress. Our approach is motivated by results from the reinforcement learning literature SuttonBarto2018, and especially the policy-gradient theorem. The policy-gradient theorem is widely used for optimizing reinforcement-learning systems via first-order algorithms sutton1999policy, and has recently been deployed for estimating global treatment effects in nonstationary Markovian A/B tests with temporal interference johari2025. However, we are not aware of previous uses of policy-gradient theorems for observational study causal inference in settings of the type we consider here.
We work under the general dynamic model introduced in (ref). Our first goal will be to provide a reduced-form characterization of the policy gradient in dynamic thresholding designs and connect it to an RD estimand that can be interpreted as a marginal policy effect. This characterization will later serve as the foundation for developing a tractable method for estimation and inference in dynamic thresholding designs.
To meaningfully evaluate how changing the threshold $c$ affects the value function $V(\pi_c)$ defined in (ref), we seek to characterize the negative policy gradient $-\partial V(\pi_c)/\partial c$, which represents the marginal effect on total discounted welfare of infinitesimally lowering the cutoff (i.e., treating slightly more units). As in (ref) for the cross-sectional case, we show that this gradient can be expressed in terms of conditional response functions and the density of the running variable at the cutoff. However, the dynamic setting introduces an additional complexity: The relevant conditional response functions must now account for all future treatment dynamics induced by changing the treatment assignment for the current period. To formalize this, we introduce the $Q$-function (or action-value function), which plays a central role in reinforcement learning and dynamic treatment regime analysis:
In words, $Q_{c,\,t}(s_t,\,z_t,\,a_t)$ is the expected discounted sum of future rewards starting with history $s_t$, running variable $z_t$ and action $a_t$ at time $t$.
To characterize the policy gradient, all we need in addition to the basic model from Section (ref) is that the running variable have a density around the cutoff $c$ conditionally on past state, and that relevant conditional-response functions vary smoothly with the running variable. We note that both assumption will hold whenever there is non-trivial (continuously distributed and exogenous) noise in the running variable lee2008randomized,Eckles2025.
Under the above assumptions, the following result explicitly characterizes the policy gradient in our dynamic thresholding setting. Although it is conceptually similar to the standard policy-gradient theorem sutton1999policy, we note that it is not a direct corollary: The standard policy gradient quantifies the effect of changing action probabilities under overlap conditions (i.e., where treatment and control actions can both occur with positive probability in all states), whereas here we consider the effect of changing the treatment cutoff in a setting without overlap.
The above result establishes that the negative gradient of the value function $V(\pi_c)$ with respect to the threshold parameter $c$ equals a discounted sum of the $Q$-function differences at the threshold across all future time periods, weighted by the conditional density at the threshold in each period. In addition to providing a unified expression for both finite and infinite horizon settings, (ref) organically handles two core challenges that arise specifically in dynamic thresholding designs:
We also define the $Q$-function for the treatment indicators as:
In words, $Q_{c,\,t}^A(s_t,\,z_t,\,a_t)$ is the expected number of future periods in which the unit is exposed to the treatment starting with history $s_t$, running variable $z_t$ and action $a_t$ at time $t$.
An immediate application of (ref) to both the numerator and denominator of $\tau_{\,\mathrm{RD}}$ as defined in (ref) then yields the following characterization.
The expression for the causal parameter $\tau_{\,\mathrm{RD}}$ given in (ref) is explicit---but at first glance may appear unwieldy to operationalize because of its dependence on the growing-dimensional state $S_{i,t}$. Perhaps surprisingly, however, it turns out that the quantity $\tau_{\,\mathrm{RD}}$ can be estimated using a carefully designed local linear regression procedure. Although we here consider sharp thresholding designs our estimator has a ratio form typically associated with fuzzy RD designs IL2008review. The ratio form is explained by the fact that, in the dynamic design, there is some uncertainty on how moving $c$ will affect the treatment frequency (since we may not get to observe realizations of $Z_{i,t}$ for $t \geq 1$ for counterfactual thresholding policies); and, as argued in sun2021treatment, cost-benefit analyses under cost uncertainty induce statistical structure resembling that encountered in instrumental-variable analyses.
(ref) provides a unified characterization of the dynamic marginal policy effect across both finite and infinite horizons. Our estimation and inference procedures, however, will differ slightly between these two cases. The key difference lies in how one can approximate the $Q$-functions: In finite-horizon settings, we can directly use the cumulative discounted sum of rewards observed from time $t$ onward, whereas infinite-horizon problems demand separate attention to the fact that the observed trajectories in any real-world data are necessarily truncated. We first present results for the finite-horizon case.
Consider the dynamic thresholding design introduced in (ref) with $T<\infty$. Define the discounted sum of future outcomes and treatment assignments from time $t$ onward as
where $t=0,1,\dots,T$, and $i=1,2,\dots,n$. Motivated by the twice-discounted structure of the expression in (ref), where discounting appears both within the $Q$-functions and in the outer summation, we propose below a twice-discounted local linear regression procedure for estimating the causal parameter $\tau_{\,\mathrm{RD}}$. We choose a small bandwidth $h=h(n)$ (where $h(n)\to 0$ as $n\to\infty$), a weighting function $K:\mathbb{R} \to [0,\infty)$ and run the following weighted linear regression on each side of the threshold:
where $r(z):=(1,\mathbf{1}\left\{z\ge c\right\},z-c,(z-c)\mathbf{1}\left\{z\ge c\right\})^\top$ denotes the regressors, and $e_2=(0,1,0,0)^\top$. Popular choices for the weighting function $K(\cdot)$ include the window function $K(z)=\mathbf{1}\left\{|z|\le 1\right\}$ or the triangular kernel $K(z)=(1-|z|)_+$. When $K(\cdot)$ is the window function, the above can also be numerically implemented as a two-stage least squares type estimator with “instrument” $A_{i,t}$ and “treatment” $H_{i,t}$ IL2008review.
The above regression deviates from the standard local linear regression in two ways. First, we use the discounted sum of rewards $G_{i,t}$ (instead of the immediate reward $Y_{i,t}$) to approximate the $Q$-functions, thus automatically incorporating long-term downsteam effects of treatment decisions. Second, in addition to the local kernel weights, we use the temporal weighting $\gamma^t$ that mirrors our dynamic policy gradient result (cf. (ref)). Our next result illustrates that with these two simple tweaks to the standard local linear regression, we can consistently estimate the causal parameter $\tau_{\,\mathrm{RD}}$. Before stating this result, we list some regularity assumptions.
Our next goal is to show that our estimator achieves the standard nonparametric rate of $n^{-2/5}$ for local linear regression based RD estimators under appropriate smoothness conditions. Continuing the parallel between dynamic thresholding designs and classical (cross-sectional) RD designs, we impose the following regularity conditions that are natural extensions of second-order smoothness assumptions standard in classical RD literature (see, e.g., Hahn-et-al-2001).
The following result establishes the limiting distribution of the local linear regression estimator proposed in (ref) with an explicit characterization of the asymptotic variance.
(ref) reveals that our estimator achieves the same $n^{-2/5}$ rate of convergence and exhibits the same bias structure as classical local linear regression estimators in static RD designs. The leading bias term has the familiar $h^2$ form, driven by the curvature of the conditional response functions $\mu_{G,\,a}(z)$ at the threshold. The asymptotic variance in (ref) reflects the additional complexity of the dynamic setting: It aggregates uncertainty across all time periods, with contributions weighted by $\gamma^{2t}$ to account for temporal discounting. It is also interesting to note that when $T=0$, the (ref) reduces to the analogous asymptotic result for the classical RD setting (see, e.g., Hahn-et-al-2001).
Equipped with the asymptotic distribution from (ref), we can construct asymptotically valid confidence intervals following standard practice in the static RD literature. Define sandwich estimators of the asymptotic variance/covariance terms used in (ref) (namely, $V_G$, $V_H$ and $V_{GH}$, as defined in (ref)) as follows.
where $e_1=(1,0,0,0)^\top$, and
with $\wh\theta(h(n))$ and $\wh\eta(h(n))$ as defined in (ref). Define ${B}_{n,\,0}$, $M_{G,\,n,\,i,\,0}$ and $M_{H,\,n,\,i,\,0}$ analogously, with $1-A_{i,t}$ replacing $A_{i,t}$ in the above definitions. The following result provides a consistent estimator of the asymptotic variance in (ref).
The above result taken together with (ref) allows us to construct asymptotically valid confidence intervals using appropriate choices for the bandwidth $h=h(n)$, as follows.
We now consider the infinite horizon setting where our target quantity is defined via (ref) with $T=\infty$; see also (ref). This case is particularly relevant for policy evaluation in ongoing programs with no predetermined endpoint. The infinite horizon introduces an additional complication: In any real-world data, we can only observe trajectories that are long but finite, even though our target parameter involves an infinite sum. As a consequence, we cannot directly compute the discounted sum of all future rewards $G_{i,t}$ as we did in the finite horizon case. To overcome this issue, we work with discounted partial sum of rewards, as follows.
Denote by $T(n)$ the common\footnote{The methods we develop can also accomodate cases where the units have trajectories of varying lengths, provided all trajectories are sufficiently long for asymptotic approximations to hold. However, to keep the exposition simple, we use a common trajectory length $T(n)$ to state our results.} length of the observed trajectories and consider the regime where $T(n)\to\infty$ as $n\to\infty$. We fix a truncation window $\ell=\ell(n)$ and define the discounted partial sum of future outcomes and treatment assignments from time $t$ onward as:
where $t=0,1,\dots,T(n)-\ell(n)+1$, $i=1,2,\dots,n$. We carefully choose the truncation window $\ell(n)$ such that it grows with the sample size so that the truncation bias becomes asymptotically negligible, and it is small enough relative to the trajectory length $T(n)$ to maintain a sufficient effective sample size. We refer the reader to (ref) for precise conditions on the size of the truncation window.
Next, similar to the finite horizon case, we choose a small bandwidth $h=h(n)$ (where $h(n)\to 0$ as $n\to\infty$), a weighting function $K:\mathbb{R} \to [0,\infty)$ and run the following twice-discounted local linear regression on each side of the threshold:
where $r(z)=(1,\mathbf{1}\left\{z\ge c\right\},z-c,(z-c)\mathbf{1}\left\{z\ge c\right\})^\top$ denotes the regressors and $e_2=(0,1,0,0)^\top$. The following results mimics (ref) and gives us the desired consistency of the twice-discounted local-linear-regression estimator $\wh\tau_{\,\mathrm{RD}}$ as defined above.
The above result shows that, despite working with truncated trajectories and truncated reward sums, we can still consistently estimate the true infinite-horizon parameter $\tau_{\,\mathrm{RD}}$, provided that the truncation window $\ell$ grows fast enough with the sample size relative to the bandwidth $h$ (more precisely, $\ell \gg \log h/\log \gamma$) so that the discounted tail contributions (rewards beyond $\ell$-steps ahead) become asymptotically negligible relative to the statistical error induced by the bandwidth $h$.
The asymptotic distribution of our local-linear-regression estimator in the infinite-horizon case perfectly mirrors the same derived in in (ref) for the finite-horizon case, with the summations now running to infinity. The additional technical requirement $\gamma^{\min\{\ell(n), T(n)-\ell(n)\}} = o(h^3)$ ensures that both the forward truncation bias (from using $\ell$-step rewards) and the boundary effects (from stopping at $T(n) - \ell(n)$) are asymptotically negligible relative to the curvature bias.
Equipped with (ref), and following the same recipe as in (ref), we can now construct asymptotically valid confidence intervals. Define $\widehat{V}_{G,\,n}$, $\widehat{V}_{H,\,n}$, $\widehat{V}_{GH,\,n}$ as in (ref), with the quantities in (ref) modified as follows:
with $\wh\theta_n$ and $\wh\eta_n$ as defined in (ref). Define ${B}_{n,\,0}$, $M_{G,\,n,\,i,\,0}$ and $M_{H,\,n,\,i,\,0}$ analogously, with $1-A_{i,t}$ replacing $A_{i,t}$ in the above definitions. The following result provides a consistent estimator of the asymptotic variance $V_{\,\mathrm{RD}}$ in (ref).
The above result taken together with (ref) allows us to construct asymptotically valid confidence intervals using appropriate choices for the bandwidth $h=h(n)$ and the truncation window $\ell=\ell(n)$, as follows.
In this section, we validate our theoretical results using numerical experiments. Consider a simple simulation setting where we generate data from the following autoregressive process:
We consider a finite horizon $T=12$, and generate the noise $\varepsilon_t$ as i.i.d. from $\mathcal{N}(0,1)$. We use the treatment threshold $c=110$, the baseline mean $\mu_0=100$, the autocorrelation coefficient $\rho=0.9$, the treatment intensity parameter $\tau=0.1$, and consider two choices for the drift parameter $\delta$, namely $\delta=0$ (Setting 1) and $\delta=1$ (Setting 2). We initialize the autoregressive process with $Z_0\sim \mathcal{N}(0, 4/(1-\rho^2)^{1/2})$. Note that this is not a stationary initialization, even when $\delta=0$, because of the treatment-dependent drift $\tau\neq 0$.
To improve precision, we augment the local linear regression (LLR) specification in (ref) with time fixed effects. Specifically, we run the discounted LLR in (ref) with the following regressors:
where $d_t$ is a dummy variable for the time period $t$. In our simulations, we find that including these dummy variables substantially reduces the variance of our estimator (as well as that of the baseline methods we discuss below), and we plan to investigate formal properties of this modified estimator in a later draft. We contrast our proposed procedure with the following baseline methods:
The first baseline procedure is expected to recover the short-run, partial equilibrium effect (i.e., the one-period RD jump at the threshold), which in this case is given by $\tau_\text{partial eq.}=\tau (c-\mu_0)=1$. This parameter can be very different from the marginal policy effect parameter $\tau_{\,\mathrm{RD}}$. For example, in Setting 1 (where drift $\delta=0$), we find that for discount factor $\gamma=0.8$, $\tau_{\,\mathrm{RD}}\approx 3$ (obtained using $n=5{,}000{,}000$ monte-carlo simulations), which is three times the partial equilibrium effect.
The second baseline procedure tries to naively incorporate long-run dynamics by first aggregating all future outcomes and treatment decisions into the net-present quantities $G_{i,0}$ and $H_{i,0}$, and then running a single cross-sectional RD at time $t=0$. In the limit, $\widehat{\tau}_{\mathrm{RD},\,\text{naive}}(h)$ targets the marginal effect of infinitesimally lowering the threshold at the initial period on discounted outcomes per additional discounted treatment generated by that initial perturbation. However, this estimand only uses information from the first time the running variable is near the cutoff and treats the future evolution of $Z_t$ as fixed. As a result, this naive long-run LLR generally does not target the dynamic marginal policy effect and can be biased whenever the treatment policy has strong feedback effects on the future state trajectory.
The empirical performance of our proposed method, as well as the baseline procedures described above, depends crucially on the choice of the bandwidth $h$ used in the kernel function. Here we use the data-driven IK bandwidth imbens2012optimal, optimized for the standard LLR (Baeline 1), as the common bandwidth for all methods. While this bandwidth choice may not be optimal for our method, it provides a transparent comparison across all candidate approaches.
We vary the discount factor $\gamma$ in $\{0.5, 0.8, 1\}$, and use an exponentially increasing sequence of sample sizes, namely $\{1000, 2000, 4000,\dots, 128000\}$. We report in (ref) the empirical coverage and average width of $95\%$ confidence intervals constructed using our method as well as the two baseline methods, for Settings 1 (where $\delta=0$) and 2 (where $\delta=1$), respectively. The results are aggregated across $2000$ replications.
(ref) demonstrate that our proposed confidence intervals for the dynamic marginal effect parameter $\tau_{\,\mathrm{RD}}$ provide nominal coverage across all sample sizes and discount factors. In contrast, standard LLR confidence intervals severely undercover when used to make inference on $\tau_{\,\mathrm{RD}}$, with coverage rates dropping to near zero even for moderate sample sizes. This is expected since the standard LLR targets the partial equilibrium effect $\tau_{\,\text{partial eq.}}$, not the dynamic marginal policy effect $\tau_{\,\mathrm{RD}}$. Indeed, (ref) illustrate that the standard LLR confidence intervals provide nominal coverage for the parameter $\tau_{\,\text{partial eq.}}$.
We also observe that our confidence intervals are substantially wider than those from standard LLR, with widths increasing as the discount factor approaches 1. This is expected from our asymptotic theory in (ref): The asymptotic variance of our estimator $\wh\tau_{\,\mathrm{RD}}$ aggregates uncertainty across all future time periods, and the effective number of periods contributing to this variance grows as the discount factor increases. More fundamentally, the problem of estimating the long-term marginal policy effect $\tau_{\,\mathrm{RD}}$ is intrinsically more difficult than estimating the short-run partial equilibrium effect $\tau_{\,\text{partial eq.}}$, because $\tau_{\,\mathrm{RD}}$ incorporates spillovers and dynamic treatment switching across all future periods, each contributing additional uncertainty.
Our proposed method also outperforms the naive long-run LLR (Baseline 2) in terms of both coverage and width. For smaller discount factors and moderate sample sizes, this naive baseline method provides near-nominal coverage with confidence intervals that are substantially wider than our proposed ones. On the other hand, for discount factor $\gamma = 1$ and larger sample sizes, its coverage deteriorates substantially, dropping below 90%. This breakdown is expected since this naive method only exploits the first-period threshold proximity and ignores how changing the threshold affects the frequency with which units cross the threshold in future periods. In contrast, our method maintains nominal coverage across all settings by correctly pooling information from all time periods using the twice-discounted weighting scheme.