The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
71,240 characters
Non-parametric Causal Inference in Dynamic Thresholding Designs
\allowdisplaybreaks
\maketitle
\begin{abstract}
Consider a setting where we regularly monitor patients' fasting blood sugar, and declare them to have prediabetes (and encourage preventative care) if this number crosses a pre-specified threshold. The sharp, threshold-based treatment policy suggests that we should be able to estimate the long-term benefit of this preventative care by comparing the health trajectories of patients with blood sugar measurements right above and below the threshold. A naive regression-discontinuity analysis, however, is not applicable here, as it ignores the temporal dynamics of the problem where, e.g., a patient just below the threshold on one visit may become prediabetic (and receive treatment) following their next visit. Here, we study thresholding designs in general dynamic systems, and show that simple reduced-form characterizations remain available for a relevant causal target, namely a dynamic marginal policy effect at the treatment threshold. We develop a local-linear-regression approach for estimation and inference of this estimand, and demonstrate promise of our approach in numerical experiments.
\end{abstract}
\section{Introduction}
Dynamic threshold-based rules
\begingroup
\renewcommand{}
\footnote{\hspace{-6mm}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,
\citet{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 \citet{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 \citep{thistlethwaite1960regression}, and recent decades have seen a flurry of work on
regression discontinuity designs following this insight
\citep{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 \citep{imbens2015causal}.
In contrast, we are interested in settings where each unit is eligible for treatment multiple times
(e.g., in the setting of \citet{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 \citep{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 \citet{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
\citep{carneiro2010evaluating}. In the setting of \citet{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
\citep{SuttonBarto2018}; our analysis is in particular motivated by the policy gradient theorem
and the work of \citet{sutton1999policy} on first-order optimization of reinforcement learning models.
\subsection{Regression Discontinuities as Marginal Policy Effects}\label{sec:classical-RD-as-policy-gradient}
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 \citet{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$ \citep{Hahn-et-al-2001},
\begin{equation}\label{eqn:tau-definition}
\tau_{\,\mathrm{RD}} := \tau(c) = \lim_{z\,\downarrow\, c} \, \mathbb{E}[Y_i \,|\, Z_i = z] - \lim_{z\,\uparrow\, c} \, \mathbb{E}[Y_i \,|\, Z_i = z],
\end{equation}
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 \citet{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 \citet{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
\begin{equation*}
V(\pi_c):=\mathbb{E}[Y_i(\pi_c(Z_i))]=\mathbb{E}_{\pi_c}[Y_i].
\end{equation*}
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 \citet{dong2015identifying}.
\begin{lemma}\label{lemma:classical-RD-as-policy-gradient}
Suppose that the distribution of $\{Y_i(0), \, Y_i(1), \, Z_i\}$ is exogenous to the chosen cutoff $c$.
Suppose furthermore that the running variable $Z_i$ has density $f(\cdot)$ which is continuous and positive at $c$, and that the conditional response functions $\mu_{a}(z):=\mathbb{E}[Y(a)\mid Z = z]$ $(a=0,1)$ are continuous at $c$. Then,
\begin{equation}
\label{eq:mpe_static}
\tau_{\,\mathrm{RD}} = \frac{\partial}{\partial c} V(\pi_c) \, \bigg/\, \frac{\partial}{\partial c} \mathbb{E}_{\pi_c}[A_i], \ \ \ \
\frac{\partial}{\partial c} \mathbb{E}_{\pi_c}[A_i] = -f(c).
\end{equation}
\end{lemma}
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.
\subsection{Modeling Dynamics}
\label{sec:model}
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 \citep{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 \citep{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 \citep{SuttonBarto2018}; recent textbook
discussions are given in \citet{HernanRobins2020} and \citet{wager2024causal}. Throughout, we will assume
that conditions required for the $g$-formula to hold are satisfied.
\begin{assumption}\label{assump:data-collected-under-thresholding-policy}
We observe data collected under a dynamic thresholding policy $\pi_c$ for some $c \in \mathbb{R}$, i.e.,
actions are taken according to $A_{i,t} = \mathbf{1}\left\{Z_{i,t}\ge c\right\}$.
\end{assumption}
\begin{assumption}\label{assump:g-formula}
Under policy $\pi_c$, observation sequences for each unit $i = 1, \, \ldots, \, n$ are sampled IID
from a distribution $\mathbb{P}_{\pi_c}$ which factors according to the $g$-formula,
\begin{equation*}
\mathbb{P}_{\pi_c}\left[Z_{i,0}, \, Y_{i,0}, \, \ldots, \, Z_{i,T}, \, Y_{i,T}\right]
= \prod_{t = 0}^{T} \mathbb{P}\left[Z_{i,t} \,\big|\, S_{i,t}\right] \mathbb{P}\left[Y_{i,t} \,\big|\, S_{i,t}, \, Z_{i,t}, \, A_{i,t} = \mathbf{1}\left\{Z_{i,t}\ge c\right\}\right],
\end{equation*}
where $S_{i,t} = \{Z_{i,0}, \, A_{i,0}, \, Y_{i,0}, \, \ldots, \, Z_{i,t-1}, \, A_{i,t-1}, \, Y_{i,t-1}\}$ denotes observation history up
to time $t$ and $S_{i,0} = \emptyset$,
and we emphasize that all conditional probabilities on the right-hand side of the $g$-formula are
policy-independent.
\end{assumption}
Our main question of interest is how treatment---as determined by dynamic thresholding as in
Assumption \ref{assump:data-collected-under-thresholding-policy}---affects net-present expected welfare and treatment frequency,
\begin{equation}
\label{eq:value_dynamic}
V(\pi_c):=\mathbb{E}_{\pi_c}\left[\sum_{t = 0}^T \gamma^t\, Y_{i,t} \right], \ \ \ \
V^A(\pi_c):=\mathbb{E}_{\pi_c}\left[\sum_{t = 0}^T \gamma^t\, A_{i,t} \right],
\end{equation}
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 \eqref{eqn:tau-definition}
does not immediately enable meaningful program evaluation. However, perhaps surprisingly,
we will find that marginal policy effect characterizations of the form \eqref{eq:mpe_static} remain
useful: An RD estimand defined as
\begin{equation}
\label{eq:mpe_dynamic}
\tau_{\,\mathrm{RD}} := \frac{\partial}{\partial c} V(\pi_c) \,\bigg/\, \frac{\partial}{\partial c} V^A(\pi_c)
\end{equation}
can still be effectively estimated in a sharp RD design---and can still be used to resolve
policy-relevant cost-benefit tradeoffs.
\subsection{Related Work}
The modern literature on regression discontinuity designs goes back to \citet{Hahn-et-al-2001};
influential contributions to this literature include \citet{IL2008review},
\citet{imbens2012optimal}, \citet{CCT2014robustCI} and \citet{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 \citet{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 \citep{heckman2016dynamic}.
One notable exception is \citet{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 \citet{cellini2010value}, however, makes crucial
use of a linear parametric model whereby
\begin{equation}
\label{eq:cellini}
Y_{i,t} = \sum_{t' \leq t} \theta_{t - t'} A_{i,t'} + \varepsilon_{i,t},
\end{equation}
i.e., treatment effects are homogeneous and decay uniformly over time. And, as shown by
\citet{hsu2024dynamic}, their approach no longer recovers an ATT if we allow for treatment
heterogeneity. \citet{hsu2024dynamic} propose an alternative analysis that avoids \eqref{eq:cellini}.
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 \citet[p.~1049]{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 \citet{Hahn-et-al-2001}.
Our flexible potential-outcomes based model for dynamic causal inference goes back to \citet{robins1986new}.
This model is widely used in biostatistics in the context of, e.g., marginal structural models \citep{robins2000marginal}
and optimal treatment regimes \citep{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 \citet{hsu2024dynamic},
who pair the model of \citet{robins1986new} with their potentially restrictive CIA assumption to make progress.
Our approach is motivated by results from the reinforcement learning literature \citep{SuttonBarto2018}, and
especially the policy-gradient theorem. The policy-gradient theorem is widely used for optimizing reinforcement-learning
systems via first-order algorithms \citep{sutton1999policy}, and has recently been deployed for estimating global treatment effects in nonstationary Markovian A/B tests with temporal interference \citep{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.
\section{Characterizing the Dynamic RD Estimand}
We work under the general dynamic model introduced in \cref{sec:model}. 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 \eqref{eq:value_dynamic}, 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 \cref{lemma:classical-RD-as-policy-gradient} 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:
\begin{equation}
\label{eq:Qfn}
Q_{c,\,t}\left(s_t,\,z_t,\,a_t\right):=\mathbb{E}_{\pi_c}\left[\sum_{j=0}^{T-t} \gamma^j\,Y_{i,t+j}\,\bigg|\, S_{i,t}=s_t, \, Z_{i,t}=z_t,\,A_{i,t}=a_t\right].
\end{equation}
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{sec:model}
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
\citep{lee2008randomized,Eckles2025}.
\begin{assumption}\label{assump:condtional-density}
Conditional on the history $S_{i,t}$, the running variable $Z_{i,t}$ has a density $f_t(\,\cdot\mid S_{i,t})$ that is continuous and strictly positive at $c$ for almost every $S_{i,t}$ and for each $t\ge 0$.
\end{assumption}
\begin{assumption}\label{assump:continuous-Q}
For each $t\ge 0$, the $Q$-functions $Q_{c,\,t}(s_t,\, z_t\,,\, a_t)$ $(a_t=0, \, 1)$ are continuous at $z_t = c$ for almost every $s_t$ and for $a_t = 0, \, 1$.
\end{assumption}
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 \citep{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.
\begin{theorem}\label{thm:expression-for-policy-gradient}
Suppose that \cref{assump:data-collected-under-thresholding-policy,assump:g-formula,assump:condtional-density,assump:continuous-Q} hold, and assume furthermore that the
following intergrability conditions hold: For some $\eta>0$,
\begin{equation}
\label{eq:integrability}
\begin{split}
&\sup_{t\ge 0}\,\sup_{s_t}\,\sup_{|c'-c|\le\eta}\,\mathbb{E}_{\pi_{c'}}\left[\sum_{j=0}^{T-t}\gamma^j\,|Y_{i,t+j}|\,\bigg|\, S_{i,t}=s_t\right]<\infty, \\
&\sup_{t\ge 0}\,\sup_{s_t}\,\sup_{|z-c|\le\eta} \max\left\{\left|Q_{c,\,t}(s_t,\,z,\,1)-Q_{c,\,t}(s_t,\,z,\,0)\right|, 1\right\}f_t(z\mid s_t)<\infty.
\end{split}
\end{equation}
Then, the gradient of the total discounted reward under
the threshold-based policy $\pi_c$ with respect to the threshold parameter $c$ is given by
\begin{equation*}
-\frac{\partial}{\partial c} V(\pi_c)=\sum_{t=0}^T \gamma^t\, \mathbb{E}_{\pi_c}\left[\left(Q_{c,\,t}(S_{i,t},\,c,\,1)-Q_{c,\,t}(S_{i,t},\,c,\,0)\right)f_t(c\mid S_{i,t})\right],
\end{equation*}
provided that either $\gamma < 1$ or $T < \infty$.
\end{theorem}
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, \cref{thm:expression-for-policy-gradient} organically handles two core challenges that arise specifically in dynamic thresholding designs:
\begin{enumerate}
\item \emph{Carryover effects}: The $Q$-function differences $Q_{c,\,t}(S_{i,t}, c, 1) - Q_{c,\,t}(S_{i,t}, c, 0)$ automatically incorporate downstream effects on all future outcomes; we do not need to separately model how treatment decision at time $t$ affects the outcomes at times $t' > t$.
\item \emph{Dynamic threshold proximity}: The density factors $f_t(c\mid S_{i,t})$ combined with the temporal discounting $\gamma^t$ accounts for the evolving frequency with which units appear near the threshold over time, thereby accommodating the circumstances where individuals may repeatedly switch between treatment and control groups.
\end{enumerate}
We also define the $Q$-function for the treatment indicators as:
\begin{equation}
\label{eq:Qfn-A}
Q_{c,\,t}^A\left(s_t,\,z_t,\,a_t\right):=\mathbb{E}_{\pi_c}\left[\sum_{j=0}^{T-t} \gamma^j\,A_{i,t+j}\,\bigg|\, S_{i,t}=s_t, \, Z_{i,t}=z_t,\,A_{i,t}=a_t\right].
\end{equation}
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$.
\begin{assumption}\label{assump:continuous-QA}
For each $t\ge 0$, the $Q$-functions $Q_{c,\,t}^A(s_t,\, z_t\,,\, a_t)$ $(a_t=0, \, 1)$ are continuous at $z_t = c$ for almost every $s_t$ and for $a_t = 0, \, 1$. Furthermore, assume that $$\mathbb{E}_{\pi_c}\left[\sum_{t=0}^T\gamma^t\, (Q_{c,t}^A(S_{i,t},\,c,\,1) - Q_{c,t}^A(S_{i,t},\,c,\,0)) f_t(c\mid S_{i,t})\right]>0.$$
\end{assumption}
An immediate application of \cref{thm:expression-for-policy-gradient} to both the numerator and
denominator of $\tau_{\,\mathrm{RD}}$ as defined in \eqref{eq:mpe_dynamic} then yields the following characterization.
\begin{corollary}\label{coro:mpe_dynamic}
Under the conditions of Theorem \ref{thm:expression-for-policy-gradient} as well as Assumption \ref{assump:continuous-QA}, the dynamic marginal policy effect parameter $\tau_{\,\mathrm{RD}}$ introduced in \eqref{eq:mpe_dynamic} can be equivalently expressed as
\begin{equation}\label{eq:tauRD}
\tau_{\,\mathrm{RD}}=\frac{\sum_{t=0}^T \gamma^t\, \mathbb{E}_{\pi_c}\left[(Q_{c,\,t}(S_{i,t},\, c,\, 1)-Q_{c,\,t}(S_{i,t},\, c,\, 0))f_t(c\mid S_{i,t})\right]}{\sum_{t=0}^T \gamma^t\, \mathbb{E}_{\pi_c}\left[(Q^A_{c,\,t}(S_{i,t},\, c,\, 1)-Q^A_{c,\,t}(S_{i,t},\, c,\, 0))f_t(c\mid S_{i,t})\right]},
\end{equation}
provided that either $\gamma < 1$ or $T < \infty$.
\end{corollary}
\section{Inference via Local Linear Regression}
The expression for the causal parameter $\tau_{\,\mathrm{RD}}$ given in \cref{coro:mpe_dynamic}
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 \citep{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 \citet{sun2021treatment}, cost-benefit
analyses under cost uncertainty induce statistical structure resembling that encountered
in instrumental-variable analyses.
\cref{coro:mpe_dynamic} 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.
\subsection{Finite-horizon Case}\label{sec:finite-horizon}
Consider the dynamic thresholding design introduced in \cref{sec:model} with $T<\infty$. Define the discounted sum of future outcomes and treatment assignments from time $t$ onward as
\begin{equation}\label{eqn:define-G-finite-horizon}
G_{i,t} =\sum_{j=0}^{T-t} \gamma^j \,Y_{i,\,t+j}, \ \ \ \
H_{i,t} =\sum_{j=0}^{T-t} \gamma^j \,A_{i,\,t+j}
\end{equation}
where $t=0,1,\dots,T$, and $i=1,2,\dots,n$.
Motivated by the twice-discounted structure of the expression in \cref{coro:mpe_dynamic}, 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:
\begin{equation}
\label{pooled-LLR-finite-horizon}
\begin{split}
&\widehat{\theta}(h):=\mathop{\rm argmin}_{\theta=(\alpha,\,\tau_G,\,\beta_0,\,\beta_1)}\frac{1}{n}\sum_{i=1}^n\sum_{t=0}^{T }\gamma^t \,K \left(\frac{|Z_{i,t}-c|}{h}\right)\left(G_{i,t} -r(Z_{i,t})^\top\theta\right)^2, \\[2mm]
&\widehat{\eta}(h):=\mathop{\rm argmin}_{\theta=(\alpha,\,\tau_H,\,\beta_0,\,\beta_1)}\frac{1}{n}\sum_{i=1}^n\sum_{t=0}^{T }\gamma^t \,K \left(\frac{|Z_{i,t}-c|}{h}\right)\left(H_{i,t} -r(Z_{i,t})^\top\eta\right)^2, \\[2mm]
&\wh\tau_{\,\mathrm{RD}}(h):=e_2^\top \widehat{\theta}(h) \,\big/\, e_2^\top \widehat{\eta}(h).
\end{split}
\end{equation}
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}$ \citep{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.~\cref{thm:expression-for-policy-gradient}).
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.
\begin{assumption}\label{assump5:regularity-conditions-for-consistency}
\hspace{1mm}
\begin{enumerate}[label=(\roman*)]
\item The kernel $K:\mathbb{R} \to [0,\infty)$ is supported on $[-1,1]$, symmetric, bounded, and integrable. Moreover, $\kappa_2\kappa_0-\kappa_1^2>0$ where $\kappa_j:=\int_0^1 u^j K(u)du$, $j\ge 0$.
\item The second-moment functions $z\mapsto m_{2,\,t}(s,\,z,\,a):=\mathbb{E}_{\pi_c}[G_{i,t}^2\mid S_{i,t}=s,\, Z_{i,t}=z,\, A_{i,t}=a]$ and $z\mapsto m_{2,\,t}^A(s,\,z,\,a):=\mathbb{E}_{\pi_c}[H_{i,t}^2\mid S_{i,t}=s,\, Z_{i,t}=z,\, A_{i,t}=a]$ $(a=0,1)$ are continuous at $c$ for a.e.~$s$, for every $t\ge 0$.
\item The densities $f_t$ and second-moments $m_{2,t}$ are locally bounded: For some $\eta>0$, there exists measurable envelopes $B_{f,\,t}$ and $B_{m_2,\,t}$ such that for a.e.~$s$, and for $a=0,1$,
$$\sup_{|z-c|\le \eta}f_t(z\mid s)\vee 1\le B_{f,\,t}(s),\quad \sup_{|z-c|\le \eta} m_{2,\,t}(s,\, z,\, a)\vee 1\le B_{m_2,\,t}(s).$$
Moreover, assume that $\sum_{t= 0}^T\gamma^t\,\mathbb{E}_{\pi_c}[ B_{m_2,\,t}(S_{i,t})\, B_{f,\,t}(S_{i,t})]<\infty$.
\item $\sup_{t\ge 0}\mathbb{E}_{\pi_c}[Y_t^2]<\infty$.
\end{enumerate}
\end{assumption}
\begin{theorem}\label{thm:consistency-finite-horizon}
Suppose that
\cref{assump:data-collected-under-thresholding-policy,assump:g-formula,assump:condtional-density,assump:continuous-Q,assump:continuous-QA,assump5:regularity-conditions-for-consistency} hold true, and that we run the twice-discounted local linear regression \eqref{pooled-LLR-finite-horizon} with bandwidth $h=h(n)$ that satisfies $h(n)\to 0$ and $n h(n)\to \infty$. Then $\wh\tau_{\,\mathrm{RD}}(h(n))$ converges in probability to the causal parameter $\tau_{\,\mathrm{RD}}$ defined in \eqref{eq:mpe_dynamic}.
\end{theorem}
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., \citet{Hahn-et-al-2001}).
\begin{assumption}\label{assump6:bounded-second-derivatives}
The functions $f_t(\,\cdot\mid s)$, $Q_{c,\,t}(s,\, \cdot\,,\, a)$ and $Q_{c,\,t}^A(s,\, \cdot\,,\, a)$ $(a=0,1)$ have second derivatives in $[c-\eta,c+\eta]$, and there exist measurable envelopes $B_{f'',\,t}(s)$ and $B_{Q'',\,t}(s)$ such that
\begin{equation*}
\begin{split}
&\sup_{|z-c|\le \eta}\max\left\{\left|\frac{\partial^2}{\partial z^2} \,f_t(z\mid s)\right|, 1\right\}\le B_{f'',\,t}(s),\\
&\sup_{|z-c|\le \eta}\max\left\{ \left|\frac{\partial^2}{\partial z^2} \,Q_{c,\,t}(s,\, z,\, a)\right|, \left|\frac{\partial^2}{\partial z^2} \,Q_{c,\,t}^A(s,\, z,\, a)\right|, 1\right\}\le B_{Q'',\,t}(s),\\
& \sum_{t=0}^T\gamma^t\, \mathbb{E}_{\pi_c}\left[B_{Q'',\,t}(S_t)\,B_{f'',\,t}(S_t)\right]<\infty,\\
&\sum_{t=0}^T\sum_{t'=t+1}^{T}\gamma^{t+t'}\, \mathbb{E}_{\pi_c}\left[B_{m_2,\,t}^{1/2}(S_t)\,B_{m_2,\,t'}^{1/2}(S_{t'})\,B_{f,\,t}(S_t)\,B_{f,\,t'}(S_{t'})\right]<\infty,
\end{split}
\end{equation*}
where the envelopes $B_{m_2,\,t}$ and $B_{f,\,t}$ are as defined in \cref{assump5:regularity-conditions-for-consistency}.
\end{assumption}
The following result establishes the limiting distribution of the local linear regression estimator proposed in \eqref{pooled-LLR-finite-horizon} with an explicit characterization of the asymptotic variance.
\begin{theorem}\label{thm:clt-for-twice-discounted-llr-finite-horizon}
Suppose that \cref{assump:data-collected-under-thresholding-policy,assump:g-formula,assump:condtional-density,assump:continuous-Q,assump:continuous-QA,assump5:regularity-conditions-for-consistency,assump6:bounded-second-derivatives} hold true, and that the twice-discounted local linear regression \eqref{pooled-LLR-finite-horizon} is run with bandwidth $h=h(n)$ that satisfies $h=O(n^{-1/5})$. Then the asymptotic distribution of the twice-discounted local linear regression estimator $\wh\tau_{\,\mathrm{RD}}=\wh\tau_{\,\mathrm{RD}}(h(n))$ is given by
\begin{equation*}
\sqrt{nh(n)}\left(\wh\tau_{\,\mathrm{RD}} -\tau_{\,\mathrm{RD}}-\frac{1}{2}h^2\xi_1\frac{\Delta\mu_G''(c)-\tau_{\,\mathrm{RD}}\Delta\mu_H''(c)}{\Delta\mu_H(c)}\right)\stackrel{d}{\longrightarrow} \mathcal{N}\left(0, V_{\,\mathrm{RD}} \right),
\end{equation*}
where $\xi_1:=(\kappa_2^2-\kappa_1\kappa_3)/(\kappa_0\kappa_2-\kappa_1^2)$ with $\kappa_j:=\int_0^1 u^j K(u)du$, $\Delta\mu_G(z):=\mu_{G,\,1}(z)-\mu_{G,\,0}(z)$ and $\Delta\mu_H(z):=\mu_{H,1}(z)-\mu_{H,0}(z)$, where $$\mu_{G,\,a}(z):=\frac{\sum_{t=0}^T \gamma^t\,\mathbb{E}_{\pi_c}[Q_{c,\,t}(S_t,\, z,\, a)f_t(z\mid S_t)]}{\sum_{t=0}^T \gamma^t\,\mathbb{E}_{\pi_c}[f_t(z\mid S_t)]},\ \ \mu_{H,\,a}(z):=\frac{\sum_{t=0}^T \gamma^t\,\mathbb{E}_{\pi_c}[Q_{c,\,t}^A(S_t,\, z,\, a)f_t(z\mid S_t)]}{\sum_{t=0}^T \gamma^t\,\mathbb{E}_{\pi_c}[f_t(z\mid S_t)]},$$
and the asymptotic variance is given by \begin{equation}\label{eq:asymp-var-ratio-finite-horizon}
V_{\,\mathrm{RD}}:=\frac{V_G + \tau_{\,\mathrm{RD}}^2 V_H - 2\tau_{\,\mathrm{RD}} V_{GH}}{(\Delta\mu_H(c))^2},
\end{equation} where the quantities $V_G$, $V_H$ and $V_{GH}$ are defined as follows:
\begin{align*}
V_{G}&:=\sum_{a=0,1}\,\xi_2\,\frac{\sum_{t=0}^T \gamma^{2 t}\, \mathbb{E}_{\pi_c}\left[\mathbb{E}_{\pi_c}\left[\left(G_t-\mu_{G,\,a}(c)\right)^2\mid S_t,\,Z_t=c,\,A_t=a\right] f_t(c \mid S_t)\right]}{\left(\sum_{t= 0}^T\gamma^t \,\mathbb{E}_{\pi_c}\left[f_t(c\mid S_t)\right]\right)^2},\\
V_{H}&:=\sum_{a=0,1}\,\xi_2\,\frac{\sum_{t=0}^T \gamma^{2 t}\, \mathbb{E}_{\pi_c}\left[\mathbb{E}_{\pi_c}\left[\left(H_t-\mu_{H,\,a}(c)\right)^2\mid S_t,\,Z_t=c,\,A_t=a\right] f_t(c \mid S_t)\right]}{\left(\sum_{t= 0}^T\gamma^t \,\mathbb{E}_{\pi_c}\left[f_t(c\mid S_t)\right]\right)^2},\addtocounter{equation}{1}\tag{\@Alph\c@section.\arabic{equation}}\label{eq:VGVHCGH}\\
V_{GH}&:=\sum_{a=0,1}\,\xi_2\,\frac{\sum_{t=0}^T \gamma^{2 t}\, \mathbb{E}_{\pi_c}\left[\mathbb{E}_{\pi_c}\left[\left(G_t - \mu_{G,\,a}(c)\right)\left(H_t - \mu_{H,\,a}(c)\right)\mid S_t,\,Z_t=c,\,A_t=a\right] f_t(c \mid S_t)\right]}{\left(\sum_{t= 0}^T\gamma^t \,\mathbb{E}_{\pi_c}\left[f_t(c\mid S_t)\right]\right)^2},
\end{align*}
with $\xi_2:=(\kappa_2^2\rho_0 - 2\kappa_1\kappa_2\rho_1+\kappa_0^2\rho_2)/(\kappa_0\kappa_2-\kappa_1^2)^2$, $\rho_j:=\int_0^1 u^j\,K^2(u)\, du$, $j\ge 0$.
\end{theorem}
\cref{thm:clt-for-twice-discounted-llr-finite-horizon} 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 \cref{thm:clt-for-twice-discounted-llr-finite-horizon} 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 \cref{thm:clt-for-twice-discounted-llr-finite-horizon} reduces to the analogous asymptotic result for the classical RD setting (see, e.g., \citet[Theorem 4]{Hahn-et-al-2001}).
Equipped with the asymptotic distribution from \cref{thm:clt-for-twice-discounted-llr-finite-horizon}, 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 \eqref{eq:asymp-var-ratio-finite-horizon} (namely, $V_G$, $V_H$ and $V_{GH}$, as defined in \eqref{eq:VGVHCGH}) as follows.
\begin{equation} \label{eq:var-est-finite-horizon}
\begin{split}
\widehat{V}_{G,\,n} &:=\sum_{a=0,1} e_1^\top {B}_{n,\,a}^{-1} \left(\frac{1}{n}\sum_{i=1}^n M_{G,\,n,\,i,\,a}^{\otimes2} \right){B}_{n,\,a}^{-1}\,e_1,\\[2mm]
\widehat{V}_{H,\,n} &:=\sum_{a=0,1} e_1^\top {B}_{n,\,a}^{-1} \left(\frac{1}{n}\sum_{i=1}^n M_{H,\,n,\,i,\,a}^{\otimes2} \right){B}_{n,\,a}^{-1}\,e_1,\\[2mm]
\widehat{V}_{GH,\,n} &:=\sum_{a=0,1} e_1^\top {B}_{n,\,a}^{-1} \left(\frac{1}{n}\sum_{i=1}^n M_{G,\,n,\,i,\,a}\,M_{H,\,n,\,i,\,a}^{\top} \right){B}_{n,\,a}^{-1}\,e_1,
\end{split}
\end{equation}
where $e_1=(1,0,0,0)^\top$, and
\begin{equation} \label{eq:BM-finite-horizon}
\begin{split}
{B}_{n,\,1}&:=\frac{1}{n}\sum_{i=1}^n\sum_{t=0}^{T} \gamma^t \,K \left(\frac{|Z_{i,t}-c|}{h(n)}\right)A_{i,t}\begin{bmatrix}
1\\ Z_{i,t}-c
\end{bmatrix}^{\otimes2},\\[2mm]
M_{G,\,n,\,i,\,1} &:= \sqrt{h(n)}\,\sum_{t=0}^{T} \gamma^t\, K \left(\frac{|Z_{i,t}-c|}{h(n)}\right)A_{i,t}\,\begin{bmatrix}
1\\ Z_{i,t}-c
\end{bmatrix}\left(G_{i,t}-r(Z_{i,t})^\top\widehat{\theta}(h(n))\right),\\[2mm]
M_{H,\,n,\,i,\,1} &:= \sqrt{h(n)}\,\sum_{t=0}^{T} \gamma^t\, K \left(\frac{|Z_{i,t}-c|}{h(n)}\right)A_{i,t}\,\begin{bmatrix}
1\\ Z_{i,t}-c
\end{bmatrix}\left(H_{i,t}-r(Z_{i,t})^\top\widehat{\eta}(h(n))\right),
\end{split}
\end{equation}
with $\wh\theta(h(n))$ and $\wh\eta(h(n))$ as defined in \eqref{pooled-LLR-finite-horizon}. 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 \eqref{eq:asymp-var-ratio-finite-horizon}.
\begin{proposition}\label{propo:sandwich-estimator-finite-horizon}
Under the conditions of \cref{thm:clt-for-twice-discounted-llr-finite-horizon}, it holds that $\widehat{V}_{G,\,n}\,{\rightarrow}_p \,V_G$, $\widehat{V}_{H,\,n}\,{\rightarrow}_p \,V_H$, and $\widehat{V}_{GH,\,n}\,{\rightarrow}_p \,V_{GH}$. As a consequence,
$$
\widehat{V}_{\,\mathrm{RD},\,n} := \frac{\widehat{V}_{G,\,n} +\wh\tau_{\,\mathrm{RD}}^2\widehat{V}_{H,\,n} - 2\wh\tau_{\,\mathrm{RD}}\widehat{V}_{GH,\,n}}{(\wh\tau_{H,\,n})^2}\ {\rightarrow}_p\
V_{\,\mathrm{RD}},
$$
where $V_{\,\mathrm{RD}}$ is defined in \eqref{eq:asymp-var-ratio-finite-horizon}, and $\wh\tau_{H,\,n}=e_2^\top\eta(h(n))$ is obtained from \eqref{pooled-LLR-finite-horizon}.
\end{proposition}
The above result taken together with \cref{thm:clt-for-twice-discounted-llr-finite-horizon} allows us to construct asymptotically valid confidence intervals using appropriate choices for the bandwidth $h=h(n)$, as follows.
\begin{corollary}
Suppose that \cref{assump:data-collected-under-thresholding-policy,assump:g-formula,assump:condtional-density,assump:continuous-Q,assump:continuous-QA,assump5:regularity-conditions-for-consistency,assump6:bounded-second-derivatives} hold true, and that the twice-discounted local linear regression \eqref{pooled-LLR-finite-horizon} is run with bandwidth $h=h(n)$ satisfying $h=o(n^{-1/5})$. Then, with $\widehat{V}_{\,\mathrm{RD},\,n}$ as defined in \cref{propo:sandwich-estimator-finite-horizon}, it holds for any $\alpha\in (0,1)$ that
$$\lim_{n\to\infty} \mathbb{P}_{\pi_c}\left(\tau_{\,\mathrm{RD}}\in\left[\wh\tau_{\,\mathrm{RD}}\pm z_{1-\alpha/2}\,\widehat{V}_{\,\mathrm{RD},\,n}^{1/2} (nh)^{-1/2}\right]\right)= 1-\alpha.$$
\end{corollary}
\subsection{Infinite-horizon Case}
We now consider the infinite horizon setting where our target quantity is defined via \eqref{eq:mpe_dynamic} with $T=\infty$; see also \cref{coro:mpe_dynamic}. 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:
\begin{equation}\label{eqn:define-G}
G_{i,t}(\ell):=\sum_{j=0}^{\ell-1} \gamma^j\, Y_{i,t+j},\quad H_{i,t}(\ell):=\sum_{j=0}^{\ell-1} \gamma^j\, A_{i,t+j},
\end{equation}
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 \cref{thm:clt-for-twice-discounted-llr,thm:consistency} 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:
\begin{equation} \label{pooled-LLR}
\begin{split}
\widehat{\theta}_n&:=\mathop{\rm argmin}_{\theta=(\alpha,\,\tau_G,\,\beta_0,\,\beta_1)}\frac{1}{n}\sum_{i=1}^n\sum_{t=0}^{T(n)-\ell(n)}\gamma^t \,K \left(\frac{|Z_{i,t}-c|}{h(n)}\right)\left(G_{i,t}(\ell(n))-r(Z_{i,t})^\top\theta\right)^2,\\[2mm]
\widehat{\eta}_n&:=\mathop{\rm argmin}_{\eta=(\alpha,\,\tau_H,\,\beta_0,\,\beta_1)}\frac{1}{n}\sum_{i=1}^n\sum_{t=0}^{T(n)-\ell(n)}\gamma^t \,K \left(\frac{|Z_{i,t}-c|}{h(n)}\right)\left(H_{i,t}(\ell(n))-r(Z_{i,t})^\top\eta\right)^2,\\[2mm]
\wh\tau_{\,\mathrm{RD}} &=\wh\tau_{\,\mathrm{RD}}(h(n),\,T(n),\,\ell(n)):=e_2^\top \widehat{\theta}_n/e_2^\top \widehat{\eta}_n,
\end{split}
\end{equation}
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 \cref{thm:consistency-finite-horizon} and gives us the desired consistency of the twice-discounted local-linear-regression estimator $\wh\tau_{\,\mathrm{RD}}$ as defined above.
\begin{theorem}\label{thm:consistency}
Suppose that
\cref{assump:data-collected-under-thresholding-policy,assump:g-formula,assump:condtional-density,assump:continuous-Q,assump:continuous-QA,assump5:regularity-conditions-for-consistency} hold for $T=\infty$. Assume that we run the twice-discounted local linear regression \eqref{pooled-LLR} with bandwidth $h=h(n)$ and truncation window $\ell=\ell(n)$ satisfying $h(n)\to 0$, $\ell(n)\to \infty$, $T(n)-\ell(n)\to\infty$, $n h(n)\to \infty$ and $\gamma^{\ell(n)}=o(h(n))$. Then the estimator $\wh\tau_{\,\mathrm{RD}}=\wh\tau_{\,\mathrm{RD}}(h(n), T(n), \ell(n))$ defined in \eqref{pooled-LLR} converges in probability to the causal parameter $\tau_{\,\mathrm{RD}}$ defined in \eqref{eq:mpe_dynamic} with $T=\infty$.
\end{theorem}
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$.
\begin{theorem}\label{thm:clt-for-twice-discounted-llr}
Suppose that \cref{assump:data-collected-under-thresholding-policy,assump:g-formula,assump:condtional-density,assump:continuous-Q,assump:continuous-QA,assump5:regularity-conditions-for-consistency,assump6:bounded-second-derivatives} hold for $T=\infty$, and that the twice-discounted local linear regression \eqref{pooled-LLR} is run with bandwidth $h=h(n)$ and truncation window $\ell=\ell(n)$ satisfying $\ell(n)\to \infty$, $T(n)-\ell(n)\to\infty$, $\gamma^{\min\{\ell(n),\,T(n)-\ell(n)\}}=o(h^3)$ and $h=O(n^{-1/5})$. Then the twice-discounted local linear regression estimator $\wh\tau_{\,\mathrm{RD}} =\wh\tau_{\,\mathrm{RD}}(h(n),\,T(n),\,\ell(n))$ defined in \eqref{pooled-LLR} satisfies
\begin{equation*}
\sqrt{nh(n)}\left(\wh\tau_{\,\mathrm{RD}} -\tau_{\,\mathrm{RD}}-\frac{1}{2}h^2\xi_1\frac{\Delta\mu_G''(c)-\tau_{\,\mathrm{RD}}\Delta\mu_H''(c)}{\Delta\mu_H(c)}\right)\stackrel{d}{\longrightarrow} \mathcal{N}\left(0, V_{\,\mathrm{RD}} \right),
\end{equation*}
where $\xi_1:=(\kappa_2^2-\kappa_1\kappa_3)/(\kappa_0\kappa_2-\kappa_1^2)$ with $\kappa_j:=\int_0^1 u^j K(u)du$, and $\Delta\mu_G$, $\Delta\mu_H$ and $V_{\,\mathrm{RD}}$ are as defined in \cref{thm:clt-for-twice-discounted-llr-finite-horizon} except with $T=\infty$.
\end{theorem}
The asymptotic distribution of our local-linear-regression estimator in the infinite-horizon case perfectly mirrors the same derived in in \cref{thm:clt-for-twice-discounted-llr-finite-horizon} 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.
\begin{remark}[Dependence of the asymptotic variance on the discount factor]
Under additional regularity conditions, the asymptotic variance $V_{\mathrm{RD}}$ in \cref{thm:clt-for-twice-discounted-llr} satisfies $V_{\mathrm{RD}} = O((1-\gamma)^{-1})$ as $\gamma \to 1$; see \cref{lemma:VRD-on-gamma} in the Appendix for a precise statement.
\end{remark}
Equipped with \cref{thm:clt-for-twice-discounted-llr}, and following the same recipe as in \cref{sec:finite-horizon}, we can now construct asymptotically valid confidence intervals. Define $\widehat{V}_{G,\,n}$, $\widehat{V}_{H,\,n}$, $\widehat{V}_{GH,\,n}$ as in \eqref{eq:var-est-finite-horizon}, with the quantities in \eqref{eq:BM-finite-horizon} modified as follows:
\begin{align*}
{B}_{n,\,1}&:=\frac{1}{n}\sum_{i=1}^n\sum_{t=0}^{T(n)-\ell(n)} \gamma^t \,K \left(\frac{|Z_{i,t}-c|}{h(n)}\right)A_{i,t}\begin{bmatrix}
1\\ Z_{i,t}-c
\end{bmatrix}^{\otimes2},\\
M_{G,\,n,\,i,\,1} &:= \sqrt{h(n)}\,\sum_{t=0}^{T(n)-\ell(n)} \gamma^t\, K \left(\frac{|Z_{i,t}-c|}{h(n)}\right)A_{i,t}\,\begin{bmatrix}
1\\ Z_{i,t}-c
\end{bmatrix}\left(G_{i,t}(\ell(n))-r(Z_{i,t})^\top\widehat{\theta}_n\right),\\
M_{H,\,n,\,i,\,1} &:= \sqrt{h(n)}\,\sum_{t=0}^{T(n)-\ell(n)} \gamma^t\, K \left(\frac{|Z_{i,t}-c|}{h(n)}\right)A_{i,t}\,\begin{bmatrix}
1\\ Z_{i,t}-c
\end{bmatrix}\left(H_{i,t}(\ell(n))-r(Z_{i,t})^\top\widehat{\eta}_n\right),
\end{align*}
with $\wh\theta_n$ and $\wh\eta_n$ as defined in \eqref{pooled-LLR}. 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 \cref{thm:clt-for-twice-discounted-llr}.
\begin{proposition}\label{propo:sandwich-estimator}
Under the conditions of \cref{thm:clt-for-twice-discounted-llr}, it holds that $\widehat{V}_{G,\,n}\,{\rightarrow}_p \,V_G$, $\widehat{V}_{H,\,n}\,{\rightarrow}_p \,V_H$, and $\widehat{V}_{GH,\,n}\,{\rightarrow}_p \,V_{GH}$, where $V_G$, $V_H$ and $V_{GH}$ are defined in \eqref{eq:VGVHCGH} with $T=\infty$. As a consequence,
\begin{equation*}
\widehat{V}_{\,\mathrm{RD},\,n} := \frac{\widehat{V}_{G,\,n} +\wh\tau_{\,\mathrm{RD}}^2\widehat{V}_{H,\,n} - 2\wh\tau_{\,\mathrm{RD}}\widehat{V}_{GH,\,n}}{(\wh\tau_{H,\,n})^2}\ {\rightarrow}_p\
V_{\,\mathrm{RD}},
\end{equation*}
where $V_{\,\mathrm{RD}}$ is the asymptotic variance in \cref{thm:clt-for-twice-discounted-llr}, and $\wh\tau_{H,\,n}=e_2^\top\wh\eta_n$ is obtained from \eqref{pooled-LLR}.
\end{proposition}
The above result taken together with \cref{thm:clt-for-twice-discounted-llr} 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.
\begin{corollary}
Suppose that \cref{assump:data-collected-under-thresholding-policy,assump:g-formula,assump:condtional-density,assump:continuous-Q,assump:continuous-QA,assump5:regularity-conditions-for-consistency,assump6:bounded-second-derivatives} hold with $T=\infty$, and that the twice-discounted local linear regression \eqref{pooled-LLR} is run with bandwidth $h=h(n)$ and truncation window $\ell=\ell(n)$ satisfying $\ell(n)\to \infty$, $T(n)-\ell(n)\to\infty$, $\gamma^{\min\{\ell(n),\,T(n)-\ell(n)\}}=o(h^3)$ and $h=o(n^{-1/5})$. Then, with $\widehat{V}_{\,\mathrm{RD},\,n}$ as defined in \cref{propo:sandwich-estimator}, it holds for any $\alpha\in (0,1)$ that
$$\lim_{n\to\infty} \mathbb{P}_{\pi_c}\left(\tau_{\,\mathrm{RD}}\in\left[\wh\tau_{\,\mathrm{RD}}\pm z_{1-\alpha/2}\,\widehat{V}_{\,\mathrm{RD},\,n}^{1/2} (nh)^{-1/2}\right]\right)= 1-\alpha.$$
\end{corollary}
\section{Numerical Experiments}
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:
\begin{equation}
\begin{split}
Z_{t+1} &= \delta + Z_t- (1-\rho)(Z_t -\mu_0) - \tau\, A_t\, (Z_t-\mu_0)_+ +4\varepsilon_t,\\ A_t &= \mathbf{1}\left\{Z_t\ge c\right\},\\ Y_t &=-Z_{t+1}, \quad t=0,1,\dots,T-1.
\end{split}
\end{equation}
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 \eqref{pooled-LLR-finite-horizon} with time fixed effects. Specifically, we run the discounted LLR in \eqref{pooled-LLR-finite-horizon} with the following regressors:
\begin{equation}
\label{time-fe-regressors}
r(z)= (1,\mathbf{1}\left\{z\ge c\right\},z-c,(z-c)\mathbf{1}\left\{z\ge c\right\}, d_1, \dots, d_{T-1})^\top,
\end{equation}
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:
\begin{enumerate}
\item Baseline 1: Run the standard LLR of $Y_{i,t}$ on $Z_{i,t}$, i.e.,
\begin{equation}
\widehat{\tau}_{\mathrm{LLR}}:=e_2^\top \mathop{\rm argmin}_{\theta=(\alpha,\,\tau_G,\,\beta_0,\,\beta_1)}\frac{1}{n}\sum_{i=1}^n\sum_{t=0}^{T } K \left(\frac{|Z_{i,t}-c|}{h}\right)\left(Y_{i,t} -r(Z_{i,t})^\top\theta\right)^2,
\end{equation}
where the regressors $r(z)$ are as defined in \eqref{time-fe-regressors}, and $e_2=(0,1,0,0)^\top$.
\item Baseline 2: Consider a naive long-run LLR approach where
we collapse each trajectory to the first-period net-present outcome $G_{i,0}$ and treatment exposure $H_{i,0}$ (as defined in~\eqref{eqn:define-G-finite-horizon}), and run separate standard LLRs of $G_{i,0}$ and $H_{i,0}$ on the first-period running variable $Z_{i,0}$. Precisely,
\begin{equation}
\label{pooled-LLR-first-period}
\begin{split}
&\widehat{\tau}_{G,\,\text{naive}}(h):=e_2^\top \mathop{\rm argmin}_{\theta=(\alpha,\,\tau_G,\,\beta_0,\,\beta_1)}\frac{1}{n}\sum_{i=1}^n K \left(\frac{|Z_{i,0}-c|}{h}\right)\left(G_{i,0} -r(Z_{i,0})^\top\theta\right)^2, \\[2mm]
&\widehat{\tau}_{H,\,\text{naive}}(h):=e_2^\top\mathop{\rm argmin}_{\theta=(\alpha,\,\tau_H,\,\beta_0,\,\beta_1)}\frac{1}{n}\sum_{i=1}^n K \left(\frac{|Z_{i,0}-c|}{h}\right)\left(H_{i,0} -r(Z_{i,0})^\top\eta\right)^2, \\[2mm]
&\widehat{\tau}_{\mathrm{RD},\,\text{naive}}(h):=\widehat{\tau}_{G,\,\text{naive}}(h) \,\big/\, \widehat{\tau}_{H,\,\text{naive}}(h).
\end{split}
\end{equation}
where the regressors $r(z)$ are as defined in \eqref{time-fe-regressors}, and $e_2=(0,1,0,0)^\top$.
\end{enumerate}
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.
\begin{table}[p]
\centering
\setlength{\tabcolsep}{4pt}
\begin{subtable}{\linewidth}
\caption{$\gamma=0.5$ ($\tau_{\,\mathrm{RD}}\approx 1.745$, $\tau_{\,\text{partial eq.}}=1$)}
\centering
\begin{tabular}{@{} c c c c c c c c c @{}}
\toprule
& \multicolumn{2}{c}{\textbf{standard LLR}}
& \multicolumn{2}{c}{\textbf{standard LLR}}
& \multicolumn{2}{c}{\textbf{naive method}}
& \multicolumn{2}{c}{\textbf{proposed method}}\\
& \multicolumn{2}{c}{\small (targeting $\tau_{\text{partial eq.}}$)}
& \multicolumn{2}{c}{\small (targeting $\tau_{\mathrm{RD}}$)}
& \multicolumn{2}{c}{\small (targeting $\tau_{\mathrm{RD}}$)}
& \multicolumn{2}{c}{\small (targeting $\tau_{\mathrm{RD}}$)}\\
\cmidrule(lr){2-3}\cmidrule(lr){4-5}\cmidrule(lr){6-7}\cmidrule(lr){8-9}
\textbf{$n$}
& \textbf{coverage} & \textbf{width}
& \textbf{coverage} & \textbf{width}
& \textbf{coverage} & \textbf{width}
& \textbf{coverage} & \textbf{width}\\
\midrule
1000 & 95.3\% & 1.380 & 42.9\% & 1.380 & 93.5\% & 7.096 & 94.8\% & 5.930 \\
2000 & 95.1\% & 1.021 & 17.6\% & 1.021 & 94.5\% & 5.275 & 94.2\% & 4.395 \\
4000 & 94.5\% & 0.758 & 3.5\% & 0.758 & 94.0\% & 3.906 & 93.9\% & 3.257 \\
8000 & 94.3\% & 0.563 & 0.1\% & 0.563 & 95.0\% & 2.891 & 94.5\% & 2.417 \\
16000 & 95.5\% & 0.417 & 0.0\% & 0.417 & 94.5\% & 2.141 & 95.0\% & 1.795 \\
32000 & 94.2\% & 0.310 & 0.0\% & 0.310 & 94.5\% & 1.589 & 95.9\% & 1.337 \\
64000 & 95.1\% & 0.230 & 0.0\% & 0.230 & 93.4\% & 1.180 & 94.5\% & 0.992 \\
128000 & 94.5\% & 0.171 & 0.0\% & 0.171 & 93.1\% & 0.877 & 95.0\% & 0.736 \\
\bottomrule
\end{tabular}\\\vspace{1mm}
\end{subtable}
\begin{subtable}{\linewidth}
\caption{$\gamma=0.8$ ($\tau_{\,\mathrm{RD}}\approx 3.014$, $\tau_{\,\text{partial eq.}}=1$)}
\centering
\begin{tabular}{@{} c c c c c c c c c @{}}
\toprule
& \multicolumn{2}{c}{\textbf{standard LLR}}
& \multicolumn{2}{c}{\textbf{standard LLR}}
& \multicolumn{2}{c}{\textbf{naive method}}
& \multicolumn{2}{c}{\textbf{proposed method}}\\
& \multicolumn{2}{c}{\small (targeting $\tau_{\text{partial eq.}}$)}
& \multicolumn{2}{c}{\small (targeting $\tau_{\mathrm{RD}}$)}
& \multicolumn{2}{c}{\small (targeting $\tau_{\mathrm{RD}}$)}
& \multicolumn{2}{c}{\small (targeting $\tau_{\mathrm{RD}}$)}\\
\cmidrule(lr){2-3}\cmidrule(lr){4-5}\cmidrule(lr){6-7}\cmidrule(lr){8-9}
\textbf{$n$}
& \textbf{coverage} & \textbf{width}
& \textbf{coverage} & \textbf{width}
& \textbf{coverage} & \textbf{width}
& \textbf{coverage} & \textbf{width}\\
\midrule
1000 & 95.3\% & 1.380 & 0.1\% & 1.380 & 93.8\% & 16.144 & 95.3\% & 10.408 \\
2000 & 95.1\% & 1.021 & 0.0\% & 1.021 & 93.5\% & 11.824 & 95.0\% & 7.602 \\
4000 & 94.5\% & 0.758 & 0.0\% & 0.758 & 95.2\% & 8.694 & 94.2\% & 5.598 \\
8000 & 94.3\% & 0.563 & 0.0\% & 0.563 & 94.5\% & 6.426 & 95.0\% & 4.134 \\
16000 & 95.5\% & 0.417 & 0.0\% & 0.417 & 94.1\% & 4.751 & 95.0\% & 3.064 \\
32000 & 94.2\% & 0.310 & 0.0\% & 0.310 & 94.2\% & 3.516 & 94.8\% & 2.284 \\
64000 & 95.1\% & 0.230 & 0.0\% & 0.230 & 93.2\% & 2.607 & 94.6\% & 1.693 \\
128000 & 94.5\% & 0.171 & 0.0\% & 0.171 & 94.0\% & 1.930 & 94.9\% & 1.256 \\
\bottomrule
\end{tabular}
\end{subtable}\\\vspace{1mm}
\begin{subtable}{\linewidth}
\caption{$\gamma=1.0$ ($\tau_{\,\mathrm{RD}}\approx 4.325$, $\tau_{\,\text{partial eq.}}=1$)}
\centering
\begin{tabular}{@{} c c c c c c c c c @{}}
\toprule
& \multicolumn{2}{c}{\textbf{standard LLR}}
& \multicolumn{2}{c}{\textbf{standard LLR}}
& \multicolumn{2}{c}{\textbf{naive method}}
& \multicolumn{2}{c}{\textbf{proposed method}}\\
& \multicolumn{2}{c}{\small (targeting $\tau_{\text{partial eq.}}$)}
& \multicolumn{2}{c}{\small (targeting $\tau_{\mathrm{RD}}$)}
& \multicolumn{2}{c}{\small (targeting $\tau_{\mathrm{RD}}$)}
& \multicolumn{2}{c}{\small (targeting $\tau_{\mathrm{RD}}$)}\\
\cmidrule(lr){2-3}\cmidrule(lr){4-5}\cmidrule(lr){6-7}\cmidrule(lr){8-9}
\textbf{$n$}
& \textbf{coverage} & \textbf{width}
& \textbf{coverage} & \textbf{width}
& \textbf{coverage} & \textbf{width}
& \textbf{coverage} & \textbf{width}\\
\midrule
1000 & 95.3\% & 1.380 & 0.0\% & 1.380 & 93.5\% & 45.699 & 94.9\% & 18.814 \\
2000 & 95.1\% & 1.021 & 0.0\% & 1.021 & 94.3\% & 32.761 & 95.3\% & 13.306 \\
4000 & 94.5\% & 0.758 & 0.0\% & 0.758 & 93.5\% & 23.631 & 95.4\% & 9.704 \\
8000 & 94.3\% & 0.563 & 0.0\% & 0.563 & 93.2\% & 17.086 & 95.0\% & 7.099 \\
16000 & 95.5\% & 0.417 & 0.0\% & 0.417 & 91.8\% & 12.475 & 94.1\% & 5.232 \\
32000 & 94.2\% & 0.310 & 0.0\% & 0.310 & 92.0\% & 9.158 & 94.8\% & 3.899 \\
64000 & 95.1\% & 0.230 & 0.0\% & 0.230 & 89.2\% & 6.757 & 95.1\% & 2.889 \\
128000 & 94.5\% & 0.171 & 0.0\% & 0.171 & 87.5\% & 4.999 & 94.8\% & 2.139 \\
\bottomrule
\end{tabular}
\end{subtable}
\caption{Empirical coverage and average width (across $2{,}000$ replications) of 95\% CIs for Setting 1 (with $\delta=0$) obtained using our proposed method and the two baseline methods: Standard LLR and a naive long-run LLR. We also report results for the standard LLR estimator targeting the partial equilibrium effect $\tau_{\text{partial eq.}}$. We vary the discount factor $\gamma\in\{0.5,0.8,1.0\}$ and use the uniform kernel with the \citet{imbens2012optimal} bandwidth.
The oracle estimates of $\tau_{\,\mathrm{RD}}$ are obtained using $5{,}000{,}000$ replications.}
\label{tab:coverage-llr}
\end{table}
\begin{table}[p]
\centering
\setlength{\tabcolsep}{4pt}
\begin{subtable}{\linewidth}
\caption{$\gamma=0.5$ ($\tau_{\,\mathrm{RD}}\approx 1.706$, $\tau_{\,\text{partial eq.}}=1$)}
\centering
\begin{tabular}{@{} c c c c c c c c c @{}}
\toprule
& \multicolumn{2}{c}{\textbf{standard LLR}}
& \multicolumn{2}{c}{\textbf{standard LLR}}
& \multicolumn{2}{c}{\textbf{naive method}}
& \multicolumn{2}{c}{\textbf{proposed method}}\\
& \multicolumn{2}{c}{\small (targeting $\tau_{\text{partial eq.}}$)}
& \multicolumn{2}{c}{\small (targeting $\tau_{\mathrm{RD}}$)}
& \multicolumn{2}{c}{\small (targeting $\tau_{\mathrm{RD}}$)}
& \multicolumn{2}{c}{\small (targeting $\tau_{\mathrm{RD}}$)}\\
\cmidrule(lr){2-3}\cmidrule(lr){4-5}\cmidrule(lr){6-7}\cmidrule(lr){8-9}
\textbf{$n$}
& \textbf{coverage} & \textbf{width}
& \textbf{coverage} & \textbf{width}
& \textbf{coverage} & \textbf{width}
& \textbf{coverage} & \textbf{width}\\
\midrule
1000 & 94.4\% & 1.336 & 45.3\% & 1.336 & 94.0\% & 7.983 & 94.3\% & 4.989 \\
2000 & 95.0\% & 0.994 & 21.1\% & 0.994 & 93.0\% & 5.885 & 94.0\% & 3.701 \\
4000 & 95.0\% & 0.739 & 3.5\% & 0.739 & 95.5\% & 4.368 & 94.3\% & 2.754 \\
8000 & 94.8\% & 0.550 & 0.2\% & 0.550 & 94.2\% & 3.235 & 94.5\% & 2.047 \\
16000 & 95.0\% & 0.409 & 0.0\% & 0.409 & 95.4\% & 2.399 & 95.6\% & 1.523 \\
32000 & 94.7\% & 0.304 & 0.0\% & 0.304 & 96.0\% & 1.788 & 95.2\% & 1.132 \\
64000 & 94.3\% & 0.226 & 0.0\% & 0.226 & 94.8\% & 1.328 & 93.9\% & 0.842 \\
128000 & 93.8\% & 0.168 & 0.0\% & 0.168 & 95.2\% & 0.985 & 94.2\% & 0.625 \\
\bottomrule
\end{tabular}\\\vspace{1mm}
\end{subtable}
\begin{subtable}{\linewidth}
\caption{$\gamma=0.8$ ($\tau_{\,\mathrm{RD}}\approx 2.819$, $\tau_{\,\text{partial eq.}}=1$)}
\centering
\begin{tabular}{@{} c c c c c c c c c @{}}
\toprule
& \multicolumn{2}{c}{\textbf{standard LLR}}
& \multicolumn{2}{c}{\textbf{standard LLR}}
& \multicolumn{2}{c}{\textbf{naive method}}
& \multicolumn{2}{c}{\textbf{proposed method}}\\
& \multicolumn{2}{c}{\small (targeting $\tau_{\text{partial eq.}}$)}
& \multicolumn{2}{c}{\small (targeting $\tau_{\mathrm{RD}}$)}
& \multicolumn{2}{c}{\small (targeting $\tau_{\mathrm{RD}}$)}
& \multicolumn{2}{c}{\small (targeting $\tau_{\mathrm{RD}}$)}\\
\cmidrule(lr){2-3}\cmidrule(lr){4-5}\cmidrule(lr){6-7}\cmidrule(lr){8-9}
\textbf{$n$}
& \textbf{coverage} & \textbf{width}
& \textbf{coverage} & \textbf{width}
& \textbf{coverage} & \textbf{width}
& \textbf{coverage} & \textbf{width}\\
\midrule
1000 & 94.4\% & 1.336 & 0.1\% & 1.336 & 94.5\% & 21.724 & 95.0\% & 7.642 \\
2000 & 95.0\% & 0.994 & 0.0\% & 0.994 & 94.8\% & 14.701 & 95.0\% & 5.615 \\
4000 & 95.0\% & 0.739 & 0.0\% & 0.739 & 95.5\% & 10.553 & 95.2\% & 4.176 \\
8000 & 94.8\% & 0.550 & 0.0\% & 0.550 & 95.9\% & 7.718 & 94.3\% & 3.107 \\
16000 & 95.0\% & 0.409 & 0.0\% & 0.409 & 95.5\% & 5.699 & 94.9\% & 2.312 \\
32000 & 94.7\% & 0.304 & 0.0\% & 0.304 & 94.5\% & 4.255 & 95.2\% & 1.715 \\
64000 & 94.3\% & 0.226 & 0.0\% & 0.226 & 93.8\% & 3.153 & 94.5\% & 1.276 \\
128000 & 93.8\% & 0.168 & 0.0\% & 0.168 & 93.7\% & 2.344 & 94.4\% & 0.947 \\
\bottomrule
\end{tabular}
\end{subtable}\\\vspace{1mm}
\begin{subtable}{\linewidth}
\caption{$\gamma=1.0$ ($\tau_{\,\mathrm{RD}}\approx 4.308$, $\tau_{\,\text{partial eq.}}=1$)}
\centering
\begin{tabular}{@{} c c c c c c c c c @{}}
\toprule
& \multicolumn{2}{c}{\textbf{standard LLR}}
& \multicolumn{2}{c}{\textbf{standard LLR}}
& \multicolumn{2}{c}{\textbf{naive method}}
& \multicolumn{2}{c}{\textbf{proposed method}}\\
& \multicolumn{2}{c}{\small (targeting $\tau_{\text{partial eq.}}$)}
& \multicolumn{2}{c}{\small (targeting $\tau_{\mathrm{RD}}$)}
& \multicolumn{2}{c}{\small (targeting $\tau_{\mathrm{RD}}$)}
& \multicolumn{2}{c}{\small (targeting $\tau_{\mathrm{RD}}$)}\\
\cmidrule(lr){2-3}\cmidrule(lr){4-5}\cmidrule(lr){6-7}\cmidrule(lr){8-9}
\textbf{$n$}
& \textbf{coverage} & \textbf{width}
& \textbf{coverage} & \textbf{width}
& \textbf{coverage} & \textbf{width}
& \textbf{coverage} & \textbf{width}\\
\midrule
1000 & 94.4\% & 1.336 & 0.0\% & 1.336 & 93.7\% & 249.885 & 95.2\% & 14.542 \\
2000 & 95.0\% & 0.994 & 0.0\% & 0.994 & 94.9\% & 89.565 & 95.3\% & 10.652 \\
4000 & 95.0\% & 0.739 & 0.0\% & 0.739 & 96.8\% & 32.740 & 94.5\% & 7.921 \\
8000 & 94.8\% & 0.550 & 0.0\% & 0.550 & 97.2\% & 21.575 & 94.7\% & 5.890 \\
16000 & 95.0\% & 0.409 & 0.0\% & 0.409 & 97.0\% & 15.445 & 94.5\% & 4.379 \\
32000 & 94.7\% & 0.304 & 0.0\% & 0.304 & 94.2\% & 11.459 & 95.3\% & 3.247 \\
64000 & 94.3\% & 0.226 & 0.0\% & 0.226 & 92.1\% & 8.505 & 95.1\% & 2.416 \\
128000 & 93.8\% & 0.168 & 0.0\% & 0.168 & 89.0\% & 6.347 & 93.9\% & 1.792 \\
\bottomrule
\end{tabular}
\end{subtable}
\caption{Empirical coverage and average width (across $2{,}000$ replications) of 95\% CIs for Setting 2 (with $\delta=1$) obtained using our proposed method and the two baseline methods: Standard LLR and a naive long-run LLR. We also report results for the standard LLR estimator targeting the partial equilibrium effect $\tau_{\text{partial eq.}}$. We vary the discount factor $\gamma\in\{0.5,0.8,1.0\}$ and use the uniform kernel with the \citet{imbens2012optimal} bandwidth.
The oracle estimates of $\tau_{\,\mathrm{RD}}$ are obtained using $5{,}000{,}000$ replications.}
\label{tab:coverage-llr-nonstationary}
\end{table}
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 \citep{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 \cref{tab:coverage-llr,tab:coverage-llr-nonstationary} 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.
\cref{tab:coverage-llr,tab:coverage-llr-nonstationary} 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, \cref{tab:coverage-llr,tab:coverage-llr-nonstationary} 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 \cref{thm:clt-for-twice-discounted-llr,thm:clt-for-twice-discounted-llr-finite-horizon}: 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.
\bibliographystyle{chicago}
\bibliography{dynamic_RD_refs}
\newpage