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.
78,482 characters · 16 sections · 55 citation commands
Unifying regression-based and design-based causal inference in time-series experiments
{5pt} {5pt} {5pt} {5pt} \hypersetup{colorlinks,breaklinks,urlcolor=blue,linkcolor=blue}
{\bf Keywords}: carryover effect, randomization inference, robust standard error, switchback experiment.
Time-series experiments, also referred to as switchback experiments or N-of-1 trials, have attracted increasing attention across a variety of domains, including clinical trials lillie2011n, online platforms bojinov2023design, and other applied settings gabler2011n,mirza2017history,hawksworth2024methodological. Unlike traditional randomized controlled trials, which assign treatments at the unit level and maintain a fixed assignment throughout the study, time-series experiments repeatedly randomize treatment assignments over time for the same, single experimental unit. Time-series experiments are closely related to longitudinal experiments robins1986new, but focus on dynamic treatment effects on the same unit rather than average effects over a population. This experiment design allows for capturing dynamic causal effects in environments where both treatments and outcomes evolve over time and where treatment effects vary across time periods. Such designs are particularly attractive in modern applications where repeated interventions are feasible, the number of experimental units is limited, and the research interest is to estimate time-varying or lagged treatment effects. bojinov2019time provided a concrete example of a time-series experiment conducted at AHL Partners LLP, a quantitative hedge fund that randomly assigned large trading orders to either human traders or computer algorithms to compare execution performance. Their example is a time-series experiment with sequential treatment assignments, where the outcome of interest is the relative performance across 10 equity-index futures markets over the course of a year.
Despite the wide applications of time-series experiments, the statistical analysis presents challenges, mainly due to the temporal dependence of outcomes and the risk of model misspecification. In time-series experiments, the treatment effects may persist beyond the period of assignment, inducing the dependence of the observed outcomes on past treatments. Consequently, the methods for analyzing traditional randomized experiments is not applicable. The design-based perspective offers a principled way to avoid strong assumptions on the outcome model. By viewing the treatment assignment as the sole source of randomness, we propose treatment effect estimates with desirable properties without strong model assumptions. This perspective builds on a rich literature on design-based analysis for various experimental settings, including regression analysis for completely randomized experiments lin2013agnostic, matched-pairs experiments fogarty2018regression, split-plot designs zhao2022reconciling,mukerjee2022causal, cluster-randomized experiments middleton2015unbiased, su2021model, and network experiments aronow2017estimating,leung2022causal, gao2023causal. Recent comprehensive reviews are available in shi2024some and ding2025randomization. However, the design-based analysis of regression-based estimators in time-series experiments remains under explored.
In this paper, we develop a unified design-based framework for analyzing regression-based estimators in time-series experiments. We show that regressing the observed outcomes on current and past treatments using an appropriately specified working model and transformed regressors yields estimates of the lagged treatment effects of interest, even when the number of regressors diverges with the length of the time series. A key feature of our working model is the centering and scaling of treatment indicators, which ensures the regression coefficients to recover the causal estimands of interest. We establish the consistency and asymptotic normality of the ordinary least squares (OLS) and show that heteroskedasticity and autocorrelation consistent (HAC) variance estimators provide conservative variance estimates. Our results do not require the regression model to be correctly specified, reflecting the robustness inherent to the design-based analysis. Our work is related to bojinov2019time and liang2025randomization, who also develop design-based analyses for time-series experiments. However, our paper differs from theirs at two levels. First, at the level of estimands, we focus on treatment effects defined independently of the realized treatment path and outcome model assumptions. Second, at the level of estimators, we focus on simpler, regression-based point and variance estimators. We will provide a detailed comparison in the next section.
The remainder of the paper is organized as follows. Section (ref) introduces the setup of the time-series experiment, defines the causal estimand of interest, and presents our regression-based estimator. Section (ref) develops the asymptotic theory, establishing consistency, asymptotic normality, and variance estimation. Section (ref) discusses the extension to the continuous treatment. Section (ref) presents empirical studies, including simulation evidence and an application to the trading experiment data previously analyzed by bojinov2019time. Section (ref) concludes with a discussion of possible directions for future research. The online supplementary materials contain extensions to the regression with interaction terms and general exposure mapping, as well as the technical details.
Notation. For a vector $a = (a_1,\ldots,a_d)^\top \in \mathbb{R}^d$, let $\lVert a \rVert_1 = \sum_{i=1}^d \lvert a_i \rvert$, and $\lVert a \rVert_2 = ( \sum_{i=1}^d a_i^2 )^{1/2}$. For a matrix $A = (a_{ij}) \in \mathbb{R}^{m \times n}$, let $\|A\|_{\mathrm{F}} = ( \sum_{i=1}^m \sum_{j=1}^n a_{ij}^2 )^{1/2}$ denote the Frobenius norm, and $\lVert A \rVert_2$ denote the spectral norm. Following the R syntax, we use $\textbf{lm}(Y \sim X)$ to denote the linear regression model with outcome $Y$ and regressors $X$.
Consider a time-series experiment with a binary treatment and a single unit repeatedly measured at $T$ time points. At each time point $t=1,\ldots,T$, we apply $Z_t=0$ for the control and $Z_t=1$ for the treatment. For each time point $t$, define the treatment history up to $t$ as $z_{t:1} = (z_t, z_{t-1}, \ldots, z_1)$. Let $\{Y_t(z_T, \ldots, z_1):z_T, \ldots, z_1 = 0 \text{ or } 1\}_{t=1}^T$ denote the potential outcomes, which can depend on the whole treatment history. We first simplify the potential outcomes by imposing the no-anticipation assumption below.
Under the no-anticipation assumption (Assumption (ref)), the observed outcome simplifies to $Y_t=Y_t(z_{t:1})$, where $z_{t:1} = (z_t,\ldots,z_1)$ represents the treatment history up to and including time $t$, following the notation in bojinov2019time. The no-anticipation assumption holds when treatments are rapidly randomized (e.g., online experiments), or treatments are administered sequentially (e.g., clinical N-of-1 trials). However, when there exists forward-looking behavior, where the experimental unit can predict future treatments, the no-anticipation assumption would be violated. We do not consider violations of the no-anticipation assumption in this paper.
The following example illustrates the notation.
Time-series experiments and network experiments are closely related in that both extend the classical potential outcomes framework to accommodate structured interference. In network experiments, a unit’s potential outcome may depend not only on its own treatment but also on the treatments of its neighbors in a network aronow2017estimating. Similarly, in time-series experiments with a single unit observed over multiple periods, outcomes may depend not only on the current treatment but also on past or future treatments, giving rise to temporal interference. The no-anticipation assumption specifies a form of temporal interference: outcomes at time $t$ may depend on the entire treatment history up to $t$, but not on future treatments bojinov2019time. In this sense, no-anticipation can be seen as a generalization of the Stable Unit Treatment Value Assumption (SUTVA) rubin1980randomization, which requires that there are no hidden versions of treatment and there is no interference between units.
We introduce the following assumption to describe the design of the time-series experiment considered in this paper.
To estimate the lagged treatment effects, which capture how treatments administered in earlier periods continue to influence outcomes at later times, we consider regressing the observed outcomes on some past treatments. We motivate the working regression model by the standard causal inference setting with only one time point. Assume there is only one binary treatment $Z$ with ${\mathrm P}(Z) = p$ known, and one pair of potential outcomes $(Y(0),Y(1))$ and observe $Y = Y(Z)$. Then the standard inverse propensity score weighted estimator for the average treatment effect ${\mathrm E}[Y(1)-Y(0)]$ is the sample mean version of ${\mathrm E}[(ZY/p) - (1-Z)Y/(1-p)] = {\mathrm E}[((Z-p)/p(1-p))Y]$. Extending the idea to time-series experiments with multiple time points, we normalize the treatment indicator $Z_t$ as
Suppose we are interested in the lagged treatment effects up to the number of lags $K$. Now we consider the OLS regression without the intercept:
using observations from $t = K+1,\ldots,T$. Let $\tilde\tau = (\tilde\tau_0,\ldots,\tilde\tau_K)^\top \in \bR^{K+1}$ denote the OLS coefficient vector from (ref).
Several key points should be noted. First, the centering of $Z_t$ in ${\tilde Z}_t$ by ${\mathrm E}[Z_t]$ ensures that the estimand is invariant to location shifts in potential outcomes. Specifically, if the potential outcome $(Y_1,\ldots,Y_T)$ is shifted to $(Y_1+c,\ldots,Y_T+c)$ for some constant $c$, the estimand remains unchanged. This transformation invariance property guarantees that the estimand captures causal effects rather than the absolute level of the outcomes. Second, the scaling of $Z_t$ in ${\tilde Z}_t$ by $\Var[Z_t]$ ensures equal weighting across all time points in the estimand. If differential weighting across lags is desired to reflect differences in treatment assignment probabilities, the scaling can be adjusted accordingly. See the discussion about general treatment effects and the connection with bojinov2019time in Section (ref).
To ensure that the OLS estimator matches the scale of the causal estimand of interest, where the difference in scale arises from dividing $Z_t$ by $\Var[Z_t]$, we apply a linear transformation to $\tilde\tau$. Intuitively, this step rescales the regression coefficients back into the natural scale of the treatment effect. Since ${\mathrm E}[{\tilde Z}_t^2] = 1/\Var[Z_t]$ for each $t = 1,\ldots,T$, we define
where $w_k$ is the harmonic mean of the variances of the lag-$k$ treatment indicators:
By Assumption (ref), the treatment probabilities are fixed by design, and thus $w_k$ is a fixed constant. Let $$\hat\tau = (\hat\tau_0,\ldots,\hat\tau_K)^\top = (w_0^{-1} \tilde\tau_0,\ldots,w_K^{-1} \tilde\tau_K)^\top \in \bR^{K+1}$$ denote the OLS estimate from (ref) after the linear transformation (ref). In this paper, we focus on studying the design-based properties of $\hat\tau$. Under the design-based perspective, we treat the potential outcomes as fixed and the randomness solely comes from the treatment assignments.
To write $\hat\tau$ explicitly, define the outcome vector $\mY_K = (Y_{K+1},\ldots,Y_T)^\top \in \bR^{T-K}$ and the regressor matrix $\mZ_K = ({\tilde Z}_{K+1:1}^\top,\ldots,{\tilde Z}_{T:T-K}^\top)^\top \in \bR^{(T-K) \times (K+1)}$ where ${\tilde Z}_{t:t'} = ({\tilde Z}_t,{\tilde Z}_{t-1},\ldots,{\tilde Z}_{t'})^\top$ for $t>t'$. Let $\mW_K \in \bR^{(K+1)\times(K+1)}$ be the diagonal weight matrix $\mW_K = {\rm diag}(w_k:k=0,\ldots,K)$. Then, the OLS estimator in (ref) after the linear transformation (ref) can be expressed as
The OLS regression is intuitive for estimating lagged treatment effects because it directly links the outcome to both current and past treatments. Each coefficient $\hat\tau_k$ captures the average impact of a treatment administered $k$ periods earlier on the current outcome. In this way, the estimator naturally decomposes the outcome into lag-specific contributions, providing a clear measure of how treatments influence present and future outcomes over time. We will show that $\hat\tau$ can recover the causal estimand in the next subsection.
We define the treatment effect of lag-$k$ treatment on the outcome at time $t$ as:
where the expectation is taken over the joint distribution of $(Z_{t:t-k+1}, Z_{t-k-1:1})$, holding $Z_{t-k}$ at 1 or 0. This estimand $\tau_{t,k}$ represents the average treatment effect on the outcome at time $t$ by setting the treatment at time $t-k$ to 1 rather than 0, averaging over all other treatment assignments. We then define the average treatment effect of lag-$k$ treatment on the outcome over time as:
Let $\tau = (\tau_0,\ldots,\tau_K)^\top \in \bR^{K+1}$ be the vector of lagged treatment effects. In this section, we focus on the interpretation of $\tau$, and we establish the asymptotic behavior of $\hat{\tau}$ as an estimator of $\tau$ in Section (ref).
The lag-$k$ treatment effect $\tau_k$ in (ref) does not rely on any specific model assumptions about the potential outcomes. However, its interpretation becomes more explicit under additional modeling assumptions. We consider a linear model for the potential outcome to facilitate the interpretation of $\tau_k$.
The linear model in Example (ref) contains several important time series models, e.g., the autoregressive model and the moving-average model. In the following two examples, we illustrate the interpretation of $\tau_k$ under these two special models.
In this section, we analyze the asymptotic properties of $\hat\tau$ as $T \to \infty$ bojinov2019time,liang2025randomization, analogous to the asymptotic analysis in traditional randomized experiments by letting the sample size grow to infinity freedman2008regression, lin2013agnostic,li2017general. By the representation of $\hat\tau$ in (ref), we decompose $\hat\tau$ as:
Therefore, the key is to analyze the asymptotic behaviors of $(T-K)^{-1}\mZ_K^\top \mZ_K$ and $(T-K)^{-1}\mZ_K^\top (\mY_K - \mZ_K \mW_K \tau)$, since $\mW_K$ is a fixed diagonal matrix. We impose the following assumption on the treatment probabilities of the time-series experiment.
Assumption (ref) assumes that the treatment probabilities are uniformly bounded away from 0 and 1 for all time points, which rules out extreme treatment probabilities. This corresponds to the classic overlap assumption for causal inference with observational studies d2021overlap, which assume the propensity score are bounded away from 0 and 1 for almost all covariates.
Then we impose weak dependence conditions on potential outcomes sequence $[Y_t = Y_t(Z_{t:1})]_{t=1}^\infty$, which rule out strong dependence of the observed outcomes on the past treatments. We consider two types of weak dependence: decaying carryover effects and $m$-dependence. Our results can be extended to more general weak dependence structures using the same proof techniques.
We first introduce the decaying carryover effects assumption, which assumes the carryover effects of the treatment decay over time.
We then introduce the $m$-dependence assumption hoeffding1948central, which restricts the influence of past treatments to some given number of time points.
Assumption (ref) requires that the carryover effects of the treatment decay sufficiently fast. Similar conditions appear in the previous time-series experiments literature ni2023design, liang2025randomization. The decaying carryover effects assumption is analogous to the approximate neighborhood interference assumption in network experiments leung2022causal, and can be viewed as a causal analogue of mixing conditions in traditional time series analysis rosenblatt1956central, bradley2005basic. Assumption (ref) requires that the observed outcome at time $t$ depends only on the most recent $m$ treatments. Assumptions (ref) and (ref) are related but do not necessarily imply one another. Assumption (ref) states that the carryover effects of the treatment decay over time, while still allowing dependence on all prior treatments. In contrast, Assumption (ref) permits arbitrary dependence on past treatments but restricts this dependence to a given number of previous time points. In both Assumptions (ref) and (ref), the first condition is weaker than the second condition, as first condition is required for consistency and the second condition is required for asymptotic normality.
In the decomposition (ref), we can write
Consider the sequence $[{\tilde Z}_{t:t-K} (Y_t - {\tilde Z}_{t:t-K}^\top \mW_K \tau)]_{t=K+1}^\infty$. Under Assumption (ref), for a fixed $K$, the carryover effects of the sequence decay with the same asymptotic order as $\theta_{T,k}$. Under Assumption (ref), the sequence is $(K \vee m)$-dependent. Therefore, we can establish the consistency and apply the central limit theorem (CLT) for the weak dependence sequence chen2004normal, chandrasekhar2023general.
To ensure well-behaved asymptotic properties of $\hat\tau$, we impose the following standard moment conditions on the potential outcomes.
In Assumption (ref), if we further assume $K$ is fixed, the uniformly boundedness of $\lVert \tau \rVert_1$ follows directly from the uniformly boundedness of the potential outcomes $Y_t$. However, when $K$ is diverging, the uniformly boundedness of $\lVert \tau \rVert_1$ implies that $\tau_k$ decays to zero as $k \to \infty$, which indicates that the treatment effect decays over time. It is possible to relax the uniformly boundedness assumption on $\tau$, but it will complicate the asymptotic analysis. We impose Assumption (ref) for the simplicity of presentation.
We first show the consistency of $\hat\tau$ in the following theorem.
Since the dimension of $\tau$ can potentially grow with $T$, we establish the CLT for a linear projection of $\hat\tau$ rather than $\hat\tau$ itself. To characterize the asymptotic variance, we introduce a covariance matrix $\mV \in \bR^{(K+1) \times (K+1)}$ below. Since $\hat\tau$ is the estimate of $\tau$ and $\tilde\tau = \mW_K \hat\tau$ by linear transformation (ref), $\tilde\tau$ is the estimate of $\mW_K \tau$. We define the oracle residual of OLS regression (ref) at time $t$ as
Then by (ref), for $k,k'=1,\ldots,K+1$, define the $k,k'$-th element of $\mV$ as:
By definition of $\mV$ in (ref), $\mV$ depends on $T$ but not on the observed treatment path and potential outcomes path. This covariance matrix captures the dependence structure of the residuals across different lags. To ensure the limit distirbution of $\hat\tau$ does not degenerate, we need the following assumption on the smallest eigenvalue of $\mV$.
We can then establish the CLT for the linear projection of $\hat\tau$ in the following theorem.
By the Cramér–Wold theorem, Theorem (ref) directly implies a multivariate CLT for $\hat\tau$ when $K$ is fixed.
In Appendix, we establish three lemmas for proving Theorem (ref) (Consistency) and Theorem (ref) (Asymptotic normality). Recall the decomposition of $\hat\tau$ in (ref). The first lemma shows that $(T-K)^{-1} \mZ_K^\top \mZ_K$ asymptotically converges to $\mW_K^{-1}$ given Assumption (ref). The second and third lemmas show the asymptotic behavior of $(T-K)^{-1} \mZ_K^\top (\mY_K - \mZ_K \mW_K \tau)$ given all other assumptions.
Recall the oracle residual $U_t$ defined in (ref). Let
be the sample residual of OLS regression (ref) at time $t$. Let $\mU_K = {\rm diag}(U_t:t=K+1,\ldots,T) \in \bR^{(T-K) \times (T-K)}$ and $\hat \mU_K = {\rm diag}(\hat U_t:t=K+1,\ldots,T) \in \bR^{(T-K) \times (T-K)}$ be the diagonal matrices of the oracle and sample residuals, respectively.
To account for autocorrelation, we introduce a kernel weighting matrix $\mQ_K \in \bR^{(T-K) \times (T-K)}$. Specifically, we use the Bartlett kernel bartlett1950periodogram, defined entrywise as \[ [\mQ_K]_{ij} =
\] for $1 \le i,j \le T-K$, where $L$ is the bandwidth parameter controlling the effective dependence range. A common choice is $L = \lfloor T^{1/4} \rfloor$ greene2003econometric. By construction, $\mQ_K$ is positive semidefinite. The HAC estimator with the Bartlett kernel corresponds to Newey–West estimator newey1987simple. Our results can be generalized to HAC estimators with other kernels, see andrews1991heteroskedasticity for different choices of kernels.
To construct the HAC variance estimator for $\hat\tau$, we first define the HAC variance estimator for $\tilde\tau$ by OLS regression (ref) as:
and then define the HAC variance estimator for $\hat\tau$ by linear transformation (ref) as:
Intuitively, the HAC variance estimator corrects the usual OLS variance formula by allowing residuals at nearby time points to be correlated. The kernel matrix $\mQ_K$ assigns larger weights to correlations between residuals that are closer in time and gradually downweights correlations as the lag length increases, with the bandwidth $L$ determining how far into the past these correlations are taken into account. This adjustment ensures that the estimated variance remains consistent under general forms of heteroskedasticity and autocorrelation in the sampling-based inference literature. We borrow their intuition and analyze $\hat\mV$ from the design-based framework.
Recall the lag-$k$ treatment effect defined in (ref). Define $\mb_{K} \in \bR^{(T-K) \times (K+1)}$ with $$[\mb_{K}]_{t-K,k+1} = \tau_{t,k} - \frac{w_k}{p_{t-K}(1-p_{t-K})}\tau_k,$$ for $t=K+1,\ldots,T$ and $k=0,\ldots,K$, which measures the heterogeneity of the treatment effect and treatment probability across time points. We then define the bias term as \[ \mB_K = \mb_K^\top \mQ_K \mb_K. \] By construction, $\mB_K$ is positive semidefinite because $\mQ_K$ is positive semidefinite under the Bartlett kernel. Intuitively, $\mB_K$ captures the bias component in the HAC variance estimator that arises when treatment probabilities vary over time or when treatment effects are heterogeneous across time points. The following theorem establishes that $\hat\mV$ is a conservative estimator of $\mV$.
Theorem (ref) generalizes the classical variance results of neyman1923application for completely randomized experiments. In the special case with constant treatment probabilities and homogeneous treatment effects, the additional term $\mB_K$ vanishes, and the HAC variance estimator consistently estimates the true variance. This mirrors neyman1923application’s result that the randomization-based variance estimator is consistent under constant effects. When treatment effects or assignment probabilities vary over time, however, the extra positive semidefinite term $\mB_K$ inflates the variance, making the HAC estimator conservative, just as the randomization-based variance estimator is conservative without the constant effects assumption. In this sense, our framework extends Neyman’s design-based analysis from randomized experiments to time-series experiments.
The Fisher randomization test (FRT) is a powerful tool for analyzing randomized experiments imbens2015causal,ding2024first. We can consider two types of hypothesis tests: the sharp null hypothesis and the weak null hypothesis. The sharp null hypothesis states that the potential outcomes are identical under all treatment assignments: \[ H_0: Y_t(z_{t:1}) = Y_t(z'_{t:1}) \text{ for any } z_{t:1},z'_{t:1} \in \{0,1\}^t \text{ and for all } t=\{1,\ldots,T\}, \] which allows for an exact FRT. The weak null hypothesis tests whether the treatment effects are zero only for a subset of lags: \[ H_0: \tau_k = 0 \text{ for all } k \in \cS, \] where $\cS \subset \{0,\ldots,K\}$ denotes an arbitrary set of lags we wish to test, which leads to an asymptotic test by using the asymptotic distribution of $\hat\tau$ (Theorem (ref)) and the HAC variance estimator $\hat\mV$ (Theorem (ref)). One can run FRTs with studentized statistics based on our regression estimator $\hat\tau$ together with the HAC variance estimator $\hat\mV$ wu2021randomization. Studentization works here because a conservative variance estimator $\hat\mV$ is available in general, which is also consistent under the sharp null hypothesis. This highlights the value of our asymptotic analysis: even when the goal is to perform a FRT, having a conservative variance estimator is crucial to apply studentization.
Importantly, the definition of the weak null hypothesis differs across frameworks. In bojinov2019time, the weak null hypothesis is formulated in terms of their path-dependent causal estimands and focuses on a single lag, reflecting the estimands defined in their design-based analysis. In our setting, the estimands are multiple lagged treatment effects, so the weak null hypothesis naturally takes a different form, where we test multiple lag-specific effects simultaneously. Our contribution is to show how such regression-based estimators can be used to construct asymptotic tests for lag-specific treatment effects, providing a complementary perspective that is closely aligned with standard regression practice.
Although the OLS procedure in (ref) involves regressing $Y_t$ on the full set of transformed treatment indicators ${\tilde Z}_t,\ldots,{\tilde Z}_{t-K}$, all of our results remain valid when regressing $Y_t$ on a subset of these transformed treatment indicators. Suppose we are interested in estimating the treatment effect at lag $k$. There are two natural approaches. One approach is the full OLS in (ref), where we regress $Y_t$ on ${\tilde Z}_t,\ldots,{\tilde Z}_{t-K}$ to obtain the estimate $\tilde\tau_k$, and then apply the linear transformation (ref) to obtain the estimate $\hat\tau_k = w_k^{-1} \tilde\tau_k$. The other approach is the marginal OLS, where we regress $Y_t$ on ${\tilde Z}_{t-k}$, i.e.,
to obtain the estimate $\tilde\tau_{{\rm marginal},k}$ and then apply the linear transformation (ref) to obtain the estimate $\hat\tau_{{\rm marginal},k} = w_k^{-1} \tilde\tau_{{\rm marginal},k}$. Both approaches yield consistent and asymptotic normal estimators ($\hat\tau_k$ and $\hat\tau_{{\rm marginal},k}$) for $\tau_k$ with conservative variance estimators by HAC variance estimation. The key differences lie in their efficiency and variance estimation. First, the two approaches have different asymptotic variances. Second, the two approaches have different variance estimators, since the variance estimator in the marginal OLS is based on noisier residuals, which inflates the standard errors. As a direct result from Theorem (ref), both approaches share the same bias term $\mB_K$ in HAC variance estimation. Thus, comparing the HAC variance estimators $\hat\mV$ reduces to comparing the asymptotic variances $\mV$, since the two approaches have the same level of conservativeness. See Appendix for more technical details on this discussion. Therefore, we focus on comparing the asymptotic variances below.
By Theorem (ref), the asymptotic variance from the full OLS (ref) is $(T-K)^{-1} \Var[\sum_{t=K+1}^T {\tilde Z}_{t-k} (Y_t - \sum_{\ell=0}^K {\tilde Z}_{t-\ell} w_\ell \tau_\ell)],$ while the asymptotic variance from the marginal OLS (ref) is $(T-K)^{-1} \Var[\sum_{t=K+1}^T {\tilde Z}_{t-k} (Y_t - {\tilde Z}_{t-k} w_k \tau_k)]$. Intuitively, the asymptotic variance depends on the residual $Y_t - \sum_{\ell=0}^K {\tilde Z}_{t-\ell} w_\ell \tau_\ell$ for full OLS and $Y_t - {\tilde Z}_{t-k} w_k \tau_k$ for marginal OLS. By adding more regressors in full OLS, the residual should have lower variance, which imples the asymptotic variance will be smaller. To compare these asymptotic variances explicitly, we consider the linear model on the potential outcomes and establish the following result.
By Proposition (ref), under the linear model with homogeneous effects and nonrandom error, full OLS has a smaller asymptotic variance than marginal OLS. The variance difference is $(T-K)^{-1} \sum_{t=K+1}^T \sum_{\ell = 0,\ell \neq k}^K \beta_\ell^2 = \sum_{\ell = 0,\ell \neq k}^K \beta_\ell^2$ which is typically large, as smaller lags often have strong effects. Furthermore, full OLS achieves lower variance unless $\beta_\ell = 0$ for all $\ell > K$. However, including too many lags may lead to high model complexity. Therefore, it is crucial for practitioners to carefully choose the number of lags in the model. This selection can be guided by standard model selection techniques from the linear regression literature.
bojinov2019time defined the general treatment effect of lag $k$ at time $t$ as:
for some non-stochastic weights satisfying $\sum_{z_{t:t-k+1}, z_{t-k-1:1}}a_{z_{t:t-k+1}, z_{t-k-1:1}} = 1$ and $a_{z_{t:t-k+1}, z_{t-k-1:1}} \ge 0$. Then the general average treatment effect is defined as the corresponding average over all time points, analogous to the definition of $\tau_k$ based on $\tau_{t,k}$.
By the design under Assumption (ref), we can write $\tau_{t,k}$ in (ref) explicitly as
where \[ a_{z_{t:t-k+1}, z_{t-k-1:1}} = \prod_{s=1, s \neq t-k}^t p_s^{z_s} (1-p_s)^{1-z_s}. \] Therefore, $\tau_{t,k}$ is a special case of the general treatment effect in bojinov2019time. However, bojinov2019time argued that the general average treatment effect cannot be estimated without strong assumptions, so they changed the estimand to the causal effect as a function of the observed treatment path, which is random (See Section 3.2 in bojinov2019time). In contrast, we focus on a special class of the general treatment effects, which does not depend on the observed treatment path, and can be estimated by OLS under time-series experiments with independent treatments.
In general, we can adjust the weights in ${\tilde Z}_t$ to estimate a broad class of average treatment effects. Consider a new normalization: ${\tilde Z}_{h,t} = h_t(p_1,\ldots,p_t)(Z_t - {\mathrm E}[Z_t])/\Var[Z_t]$ for some function $h_t:\{0,1\}^t \to \bR$ at each time $t$. Define the weight as $w_{h,k} = [(T-K)^{-1} \sum_{t=K+1}^T h_t(p_1,\ldots,p_t)^2 (p_{t-k}(1-p_{t-k}))^{-1}]^{-1} $. We then consider the OLS regression without an intercept: $\textbf{lm}(Y_t \sim {\tilde Z}_{h,t} + {\tilde Z}_{h,t-1} + \cdots + {\tilde Z}_{h,t-K})$ using observations from $t = K+1,\ldots,T$, and let $(\tilde\tau_{h,0},\ldots,\tilde\tau_{h,K})$ denote the OLS coefficients. For each $k=0,\ldots,K$, we define $ \hat\tau_{h,k} = w_{h,k}^{-1} \tilde\tau_{h,k}$. The corresponding estimands of this procedure are $(\tau_{h,0},\ldots,\tau_{h,K})$, where
liang2025randomization also studied time-series experiments and the corresponding design-based analysis. Our paper differs from theirs in two key ways. First, their analysis assumes a linear convolution model (their Section 2) and a circular convolution model (their Section 3), along with a Bernoulli randomization scheme where $p_t = 1/2$ for all $t$. In contrast, our framework does not impose any structural assumptions on the potential outcomes, allowing for arbitrary and unknown heterogeneous effects over time. Second, we consider more general experimental designs with time-varying treatment probabilities. Third, liang2025randomization proposed a method-of-moments estimator, established a central limit theorem, and proposed a variance estimator. In our work, we focus on OLS estimation and analyze its asymptotic properties, along with the performance of the HAC variance estimator derived from the same OLS procedure. Both OLS and HAC variance estimation are widely used in practice, and are implemented in standard statistical software routines. We provide theoretical support for their application in time-series experiments.
Many modern experimental settings feature treatments that take discrete or continuous values rather than binary values. For instance, in online platforms, treatment intensity may correspond to the proportion of user exposure to a new algorithm or the duration for which an experimental interface is displayed. In clinical trials, treatments often involve dosage levels or exposure durations, which are naturally measured on a continuous scale. In this section, we extend our framework to accommodate treatments with continuous values.
In settings where the experiment involves multiple discrete treatment levels in addition to a control group, one natural approach is to run separate OLS regressions for each treatment level relative to control. Our results then can be applied. Therefore, we only focus on the case where the treatment variable $Z_t \in \bR$ is continuous. We introduce the following assumption, which is the continuous treatment analogue of Assumption (ref).
As in the binary treatment case, we normalize the treatment variable by considering ${\tilde Z}_t = (Z_t - {\mathrm E}[Z_t])/\Var[Z_t]$. We then fit the following OLS regression model: $$\textbf{lm}(Y_t \sim {\tilde Z}_t + \cdots + {\tilde Z}_{t-K}),$$ where the coefficient vector is denoted by $\tilde\tau$ and its transformed version by $\hat\tau = \mW_K^{-1} \tilde\tau$ as (ref) before with the same $\mW_K$. The causal estimand is the vector of lagged treatment effects $\tau= (\tau_0,\ldots,\tau_K)^\top \in \bR^{K+1}$, where the treatment effect of lag-$k$ treatment on the outcome at time $t$ is defined as
and the average treatment effect of lag-$k$ treatment on the outcome over time is defined as
which extend the notion of the treatment effect of lag-$k$ in (ref) and (ref) to the continuous treatment setting. This estimand $\tau_k$ represents the average slope of the potential outcome with respect to the treatment assigned $k$ periods earlier, while averaging over all other past and concurrent treatment assignments.
In the literature on continuous treatments, the potential outcome under treatment intensity $z$ is referred to as the (causal) dose-response curve, denoted $Y(z)$, while its derivative $({\mathrm d}/{\mathrm d} z) Y(z)$ is known as the causal derivative effect; see zhang2025doubly and references therein. Building on this perspective, we adapt the notion of derivative effects to the setting of time-series experiments. Specifically, we define the effect of treatment at lag-$k$ on the outcome at time $t$ as
where the expectation is taken over the joint distribution of $(Z_{t:t-k+1}, Z_{t-k-1:1})$. Intuitively, $\tau_{t,k}'(z)$ measures the local sensitivity of the expected outcome at time $t$ to a marginal change in treatment intensity at lag-$k$, holding the rest of the treatment path fixed. Aggregating across time, we can equivalently represent the estimand $\tau_k$ as a weighted average of these local effects:
where the weighting function $w_{t,k}(z)$ is given by
Thus, the estimand $\tau_k$ can be interpreted as the time-averaged, weighted average derivative effect at lag $k$. This places our estimand within the broader class of weighted average derivatives newey1993efficiency, and connects naturally to discussions of how OLS and instrument variable identify such objects graham2022semiparametrically, borusyak2024negative. Importantly, both the dose–response curve and the derivative effect are nonregular parameters and cannot generally be estimated at parametric rates. In contrast, we consider the weighted average of the derivative effect, which admits parametric convergence rates.
We next impose regularity conditions on the treatment distribution and outcome to extend our asymptotic results to the continuous treatment setting. We consider the following assumption on the distribution of $Z_t$, which serves as a continuous treatment analogue of Assumption (ref).
We also consider the following assumption, which serves as a continuous treatment analogue of Assumption (ref).
With these conditions in place, we obtain the consistency and asymptotic normality of the OLS estimator. These results are the continuous treatment counterparts to Theorem (ref) and Theorem (ref).
To estimate the asymptotic variance $\mV$ in Theorem (ref), we can use exactly the same HAC variance estimator $\hat\mV$ in (ref). Define $\mb_{K} \in \bR^{(T-K) \times (K+1)}$ with $$[\mb_{K}]_{t-K,k+1} = \tau_{t,k} - \frac{w_k}{\Var[Z_{t-k}]}\tau_k,$$ for $t=K+1,\ldots,T$ and $k=0,\ldots,K$, and define the bias term as \[ \mB_K = \mb_K^\top \mQ_K \mb_K. \] The following theorem establish that $\hat\mV$ is a conservative estimator of $\mV$, which is parallel to Theorem (ref).
In this section, we conduct simulations to examine the finite-sample performance of $\hat\tau$ and the HAC variance estimator $\hat\mV$ for estimating $\mV$. We then illustrate the method on data from a trading experiment.
We consider the autoregressive model in Example (ref). Specifically, we take a first-order autoregressive model, \[ Y_t(z_{t:1}) = \mu_t(z_{t:1}) + \phi Y_{t-1}(z_{t-1:1}) + \epsilon_t(z_{t:1}), \] for any treatment history $z_{t:1} \in \{0,1\}^t$. We set $\mu_t(z_{t:1}) = \mu_t(z_t) = \mu(z_t)$ with $\mu(1)=0.5$ and $\mu(0)=0$, take $\phi = 0.5$, and let $\epsilon_t(z_{t:1}) = \epsilon_t \sim N(0,1)$ i.i.d across $t$. Under this specification, $Y_t(z_{t:1}) = \sum_{k=0}^{t-1} \phi^k \mu(z_{t-k}) + \sum_{k=0}^{t-1} \phi^k \epsilon_{t-k}$. Then the true lag-$k$ treatment effect is $\tau_k = \phi^k (\mu(1) - \mu(0))$.
We vary the time horizon $T \in \{100,200,500,1000,2000,5000,10000\}$, fix the number of lags at $K=5$, and use a constant treatment probability $p_t=p=0.5$. Adopting a design-based perspective, randomness arises only from the assignment path $Z_1,\ldots,Z_T$. For each $T$, we therefore fix a realization $[\epsilon_1,\ldots,\epsilon_T]$ and repeatedly resample $Z_1,\ldots,Z_T$. For each draw we compute the OLS estimator $\hat\tau$ and the HAC variance estimator $\hat\mV$ with the Bartlett kernel. A common choice is $L = \lfloor T^{1/4} \rfloor$ greene2003econometric. Since the treatment probabilities are constant across time and the treatment effects are homogeneous over time, the HAC variance estimator $\hat\mV$ is consistent for the asymptotic variance $\mV$ by Theorem (ref).
Figure (ref) presents histograms of $\hat\tau_k - \tau_k$ for $k=0,\ldots,5$ across replications under $T=1000$. The empirical distributions are close to the normal distribution, supporting our central limit theorem (Theorem (ref)). This provides empirical evidence that the regression-based estimator achieves approximate normality even at moderate sample sizes.
Figure (ref) displays violin plots of $\lVert \hat\tau - \tau \rVert_2$ as a function of $T$. The $L_2$ error decreases with larger sample sizes, aligning with the $T^{-1/2}$ convergence rate established in Theorem (ref).
Figure (ref) investigates the performance of the HAC variance estimator $\hat\mV$. The plot shows the normalized Frobenius error $\lVert \hat\mV - \mV \rVert_{\mathrm F}/\lVert \mV \rVert_{\mathrm F}$ for different kernel bandwidth choices $L \in \{0,1,5,\lfloor T^{1/4} \rfloor\}$, averaged across replications. When no autocorrelation adjustment is made ($L=0$), the estimator remains biased and does not improve with $T$. By contrast, increasing the bandwidth significantly reduces error, with $L = T^{1/4}$ yielding the best finite-sample performance and approaching zero error as $T$ grows. This confirms our theoretical variance results (Theorem (ref)) and demonstrates the importance of selecting an appropriate increasing bandwidth.
Table (ref) reports the empirical coverage probabilities of nominal 95% confidence intervals for $\hat\tau_k$, constructed using the HAC variance estimator $\hat\mV$ with bandwidth $L = T^{1/4}$. Across all lags $k=0,\ldots,5$ and all sample sizes, the coverage probabilities remain very close to the nominal level, even for relatively small $T=100$. As $T$ increases, the coverage stabilizes tightly around 0.95, confirming that the asymptotic normal approximation in Theorem (ref) provides accurate inference in finite samples.
We analyze the trading experiment described in bojinov2019time, which was conducted by AHL Partners, a quantitative trading firm specializing in financial futures. The firm executes orders in one of two ways: using a human trader or using an algorithmic trader. For the experiment, AHL restricted attention to orders whose sizes lay between fixed lower and upper bounds (time-invariant within each market, but varying across markets). These eligible orders were randomly assigned to either a human trader or an algorithmic trader. Therefore, treatments were assigned independently across time, following an i.i.d. Bernoulli design, ensuring that the setup satisfied our no-anticipation assumption. AHL provided 10 independent datasets from 2016, each corresponding to a different equity index futures market spanning the US, Europe, and Asia. We analyze these markets separately and regard each as a distinct experimental unit. To protect confidentiality, AHL did not reveal which of methods “A” and “B” corresponds to human versus algorithmic traders, nor did they disclose the identities of the individual markets. For consistency, we label A as the control group and B as the treatment group. The outcome of interests $Y$ is the slippage, which is a simple function of the trade prices and volumes at which an order was executed. Lower slippage indicates better execution quality and thus serves as a proxy for higher trading performance.
With $K=4$, Table (ref) reports, for each market, the sample length, mean outcomes by group (A = control, B = treatment), the joint p-value for $H_0:\tau=\mathbf{0}$, and the per-lag estimates $\hat\tau_k$ with their corresponding p-values. We find clear heterogeneity across markets: the null is rejected in Markets 4, 6, 7, and 9, indicating the presence of dynamic treatment effects, while other markets show no joint evidence of impact. In both Market 6 ($\hat\tau_0=5.65,\,p_0=0.00$) and Market 9 ($\hat\tau_0=2.43,\,p_0=0.00$), “A” (control) performed significantly better than “B” (treatment), and there is no market where “B” performed significantly better than “A”, consistent with the findings of bojinov2019time. Market 4 shows lagged effects at longer lags ($\hat\tau_2=4.89,p_2=0.05;\ \hat\tau_4=6.69,p_4=0.02$), whereas Market 7 shows a negative lag-2 effect ($\hat\tau_2=-3.77,p_2=0.00$). Some isolated single-lag signals (e.g., Market 5’s positive lag-1 estimate) do not lead to a joint rejection.
This paper develops a unified design-based framework for regression-based estimators in time-series experiments, establishing consistency, asymptotic normality, and conservative variance estimation without requiring correct model specification. We also extend our framework to the continuous treatment. In this section, we discuss the extensions we have done and some future directions.
In Appendix (ref), we extend our framework to regression with interaction terms. Beyond estimating the main effects of treatments at different lags, it is often of interest to investigate interaction effects between treatments administered at distinct time periods. Including interaction terms is particularly useful in settings where the effect of a treatment depends on its recent history or on combinations of past interventions. In online platforms, repeated exposures (e.g., showing a user the same advertisement or recommendation on consecutive days) may produce reinforcement or fatigue, leading to nonlinear cumulative effects. In clinical and behavioral studies, the efficacy of an intervention (e.g., medication, exercise, or reminders) often depends on treatment consistency, with consecutive applications generating biological accumulation or tolerance effects.
In Appendix (ref), we extend our framework to general exposure mapping. So far, we have proposed working models for the treatment effects through individual lags or interactions of consecutive lags. While this approach captures direct and pairwise temporal dependencies, it may be restrictive in settings where outcomes depend on more general features of the treatment path. To address this, we consider exposure mappings that summarize treatment assignments over multiple periods into lower-dimensional summary statistics. Specifically, we divide the past treatment assignments into several time blocks, and the treatments within each block are aggregated into a single exposure mapping. For example, in behavioral or clinical studies, treatments (e.g., reminders or interventions) often have cumulative effects that depend on how consistently the treatment has been received in the past rather than on any single recent treatment. This allows estimation of short-term, medium-term and long-term persistence in treatment effects.
We now discuss some future directions. First, while our focus has been on the single-unit, time-series setting, a natural next step is to consider panel experiments bojinov2021panel, where multiple units are observed over time and treatments are repeatedly randomized both across units and across time. Panel experiments introduce additional layers of complexity beyond single-unit time-series experiments. Interference may occur both across time within each unit and across units at a given time, leading to a rich class of potential outcomes. Researchers may be interested not only in lagged effects within units, but also in spillover effects across units and dynamic causal effects aggregated across the population. Recent work has begun to formalize and analyze these designs bojinov2021panel,ni2023design,han2024population. Our results provide a foundation for extending design-based regression analysis beyond single-unit time-series experiments. The working model we specify, together with HAC variance estimation, can be naturally generalized to panel settings. Such an extension would allow for the joint estimation of both lagged effects within units and spillover effects across units, thereby offering a unified design-based framework for dynamic causal inference in panel experiments.
Second, while our analysis assumes that treatment assignments are independent across time, many practical experiments, particularly in online platforms and sequential decision-making contexts, employ adaptive designs, where treatment assignment probabilities depend on past outcomes or treatment history. Our framework naturally extends to such adaptive experimental settings. The key observation is that one can normalize the treatment indicators conditionally on the past history instead of unconditionally. With this normalization, the same regression-based procedure and design-based analysis carry through, and asymptotic properties can be established using martingale theory.
Third, the regression-based procedure allows for incorporating lagged outcomes and additional covariates to improve efficiency. This type of regression adjustment parallels the use of covariate adjustment in classical randomized experiments. In the time-series experiments, one can augment the regression with covariates that are unaffected by the treatment of interest. Incorporating lagged outcomes is particularly relevant in time series literature, as it aligns with the structure of autoregressive models and helps account for temporal dependence in the outcome. By absorbing predictable variation, such adjustments improve efficiency without introducing bias to the treatment effect estimation.
We thank Iavor Bojinov for providing the data and Iavor Bojinov, Avi Feller and participants of the 2025 NBER-NSF Time Series Conference for helpful comments. Lin was partially supported by the Two Sigma PhD Fellowship. Ding was supported by the U.S. National Science Foundation (1945136, 2514234).