Serial Correlation and Inference in Comparative Interrupted Series

Why Serial Dependence Matters

In a comparative interrupted time series (CITS), observations close together in time often share influences that are not fully represented by the segmented regression. A period with unusually high residuals may be followed by another high residual, for example. This pattern is serial dependence. It matters because conventional regression standard errors generally treat errors from different time points as uncorrelated. When that assumption fails, the estimated treatment contrast may remain the same while its reported uncertainty is too small or too large.

Here, let \(D_t\) denote the treated-series outcome minus the comparison-series outcome at time \(t\). A simple segmented model for the difference is

where \(P_t\) is 0 before the intervention and 1 from the intervention onward, and \(R_t\) is 0 before the intervention and counts elapsed post-intervention time afterward. In this parameterization, \(\tau_0\) represents an immediate change in the treated-versus-comparison difference, and \(\tau_1\) represents a change in its slope. The residual \(u_t\) is the part of the observed difference not explained by the model. If the target is the relative effect \(h\) time units after intervention, its estimated value is \(\hat{\tau}_0+h\hat{\tau}_1\); its uncertainty must account for the covariance between the two estimates as well as the uncertainty in each.

A residual time plot is a useful first check: look for long runs above or below zero, gradual waves, or repeated oscillation. The autocorrelation at lag \(k\) summarizes the association between residuals \(k\) periods apart. A residual autocorrelation plot can show whether dependence is concentrated at lag 1 or persists across several lags. These are diagnostics, not proof that a particular error model is correct. In a CITS, an omitted change in the mean pattern can also create apparent residual structure, so first check that the segmented mean model is plausible.

Worked Example: Reading the Durbin–Watson Statistic

For many settings, \(DW\) is roughly related to the lag-one residual correlation by \(DW\approx 2(1-\hat{\rho}_1)\). Thus, a value near 2 suggests little first-order correlation, a value below 2 suggests positive correlation, and a value above 2 suggests negative correlation. This is a useful orientation, not a universal testing rule. The usual Durbin–Watson reference procedures have limitations, especially when the regression includes lagged dependent variables; a single statistic also cannot describe all patterns of dependence.

Why the Conventional Standard Error Can Mislead

In matrix notation, write the fitted regression as \(D=X\beta+u\). The ordinary least squares (OLS) coefficient estimator is \(\hat{\beta}=(X'X)^{-1}X'D\). Conditional on the design matrix \(X\), its variance is

where \(\Omega=Var(u\mid X)\) is the residual covariance matrix. The conventional OLS standard-error formula uses the special case \(\Omega=\sigma^2 I\), which assumes equal error variances and no correlation across time. If errors are serially correlated, the off-diagonal entries of \(\Omega\) are not zero, and that formula generally does not estimate the actual sampling variance.

The direction of the error is not guaranteed in every design, but persistent positive dependence commonly makes conventional standard errors too small. That can produce confidence intervals that are too narrow and tests that reject \(H_0\) too often. Changing the standard error addresses uncertainty under a model; it does not repair a misspecified mean pattern, eliminate confounding, or establish that the CITS comparison identifies a causal effect.

Three Ways to Account for Dependence

A heteroskedasticity-robust standard error by itself is not generally serial-correlation robust. A heteroskedasticity-and-autocorrelation-consistent (HAC) estimator, often implemented with Newey–West weights, estimates the coefficient variance using products of residuals at nearby time points as well as same-time residuals. A schematic “meat” component is

where \(x_t\) is the model’s regressor vector at time \(t\), \(L\) is the maximum lag included, and \(w_k\) are lag weights. One common choice is \(w_k=1-k/(L+1)\). The complete sandwich estimate multiplies \(\widehat{S}\) on the left and right by \((X'X)^{-1}\). The lag choice should reflect the sampling interval and plausible persistence, while avoiding an unnecessarily large number of noisy lag terms. HAC inference is usually justified by large-sample conditions for a weakly dependent process; it is not an automatic small-sample correction.

An alternative is to model the covariance explicitly. Under a stationary AR(1) error model,

where the innovations \(\eta_t\) are uncorrelated over time and have constant variance. This implies \(Corr(u_t,u_{t-k})=\phi^k\). Generalized least squares (GLS) uses the covariance matrix implied by the fitted \(\phi\) to estimate the regression and its uncertainty. An AR(1) model is appealing when that pattern is credible, but it should not be selected merely because it yields a preferred p-value. Residual diagnostics should be repeated after fitting, and additional dependence patterns may require a different model.

Worked Example: Comparing Conventional and HAC Inference

Block-Based Resampling

Block methods offer another way to reflect temporal dependence. In a residual moving-block bootstrap, the analyst fits the segmented model, centers its residuals as appropriate, and samples consecutive residual blocks rather than individual residuals. The sampled residual sequence is combined with the fitted mean structure to form a bootstrap outcome series; the model is refitted, and the coefficient or treatment-effect contrast is recorded. Repeating this process yields an empirical distribution for the statistic. Percentile intervals can be formed from the corresponding quantiles, while the standard deviation of bootstrap estimates can estimate a standard error.

The block length is a substantive and statistical choice. Very short blocks may break important dependence; very long blocks leave few distinct blocks to resample. The residual process should be sufficiently stable for the resampling scheme. If residual behavior differs across time or across the intervention boundary, a scheme that treats those periods as exchangeable may be inappropriate. Resampling should preserve the relevant time ordering and the design’s intervention structure, not shuffle dates or erase the interruption.

Worked Example: Calculating a Block-Bootstrap Standard Error

Common Mistakes and Practical Interpretation

• Reading the Durbin–Watson direction backward. A value below 2 is associated with positive first-order dependence, not negative dependence; a value above 2 is associated with negative first-order dependence. Check the residual sequence and the statistic’s definition.
• Using a heteroskedasticity-only correction for autocorrelation. Robustness to unequal variances does not by itself account for correlation between different months. Use a method whose covariance calculation addresses serial dependence.
• Assuming every correction is interchangeable. HAC inference, AR(1) GLS, and block resampling make different assumptions. Report the method and relevant choices, and assess whether they match the residual process and series length.
• Treating a larger standard error as a causal fix. A correlation-aware uncertainty estimate does not remove bias from an inappropriate comparison, an incorrect mean model, or other identification failures.

The aim is not to choose whichever method produces the narrowest interval. Inspect the fitted residuals, state the dependence structure or resampling scheme, and compare reasonable uncertainty analyses. If substantive conclusions change across plausible approaches, that sensitivity is important evidence about how much the data support the estimated effect. The estimate and its uncertainty remain conditional on the CITS design and its identifying assumptions.