Temporal PLR DML for Time Series

Introduction

Chapters 1–4 developed DML for cross-sectional data where observations are independent. However, many economic questions involve time series data where observations exhibit temporal dependence. This chapter develops the companion repo’s current time-series implementation: scalar partially linear DML with lagged treatment controls, temporal cross-fitting, and HAC inference.

The time series challenge

Consider estimating the effect of a policy intervention $T_t$ on an outcome $Y_t$:

Two challenges arise that standard DML cannot handle:

1. Data leakage in cross-validation: random K-fold splits allow future information to “leak” into the training set, creating unrealistic advantage.
2. Invalid standard errors: influence function–based SEs assume i.i.d. observations; autocorrelation in $\varepsilon_t$ invalidates this assumption.

Chapter roadmap

• Section 5.2: time series cross-validation — purging, blocking, and embargo.
• Section 5.3: HAC standard errors — Newey–West estimation for autocorrelation.
• Section 5.4: TemporalPLRDML framework — combining the components.
• Section 5.5: Monte Carlo validation — verifying coverage and unbiasedness.
• Section 5.6: practical recommendations.
• Section 5.7: exercises.

Time series cross-validation

The data leakage problem

Standard K-fold cross-validation randomly assigns observations to folds. With time series data, this means the model predicting $Y_{100}$ might train on $\{Y_{99}, Y_{101}, Y_{102}\}$ — data that would be unavailable in a real forecasting scenario.

Fold 1: Train on {t=2,3,5,7,9}, Test on {t=1,4,6,8,10}
Fold 2: Train on {t=1,4,6,8,10}, Test on {t=2,3,5,7,9}

Fold 1: Train on {t=1,2,3,4,5}, Test on {t=6,7,8,9,10}

Time series cross-validation strategies

Three strategies are available (implemented upstream in temporalcv; the companion’s estimators select among them via cv_strategy):

Python implementation

Our TimeSeriesCrossValidator handles all three strategies:

from temporalcv import TimeSeriesCrossValidator

# Expanding window with 5 folds and 2-period gap
cv = TimeSeriesCrossValidator(
    n_splits=5,
    gap=2,          # Periods between train end and test start
    expanding=True, # True for expanding, False for sliding
    min_train_size=50
)

# Generate train/test indices
for train_idx, test_idx in cv.split(X):
    print(f"Train: {train_idx[:5]}...{train_idx[-5:]}")
    print(f"Test:  {test_idx[:5]}...{test_idx[-5:]}")

HAC standard errors

The autocorrelation problem

Even with correct point estimates, standard errors require adjustment for serial correlation. The DML influence function

has $\mathbb{E}[\psi_t] = 0$ but $\mathrm{Cov}(\psi_t, \psi_{t-k}) \neq 0$ when outcomes are autocorrelated.

Newey–West estimation

We cannot sum to infinity, so we truncate with a kernel weighting scheme:

Kernel functions

Three commonly used kernels (implemented upstream in temporalcv):

1. Bartlett (triangular):

Linear downweighting; guaranteed positive semi-definite.

2. Parzen:

Smoother than Bartlett; used for higher-order accuracy.

3. Quadratic spectral:

Optimal for minimizing MSE; non-truncating.

Comparison of HAC kernel functions — all decay to zero as lag increases, but with different rates and smoothness:

Bandwidth selection

The bandwidth $B$ controls the bias–variance tradeoff:

• $B$ too small: misses important autocovariances (bias).
• $B$ too large: includes noise from distant lags (variance).

Data-driven selection (Andrews, 1991) adapts based on estimated autocorrelation:

where $\hat{\rho}$ is the AR(1) coefficient of the influence scores.

Python implementation

from temporalcv import newey_west_se

result = newey_west_se(
    influence_scores,
    bandwidth="auto",   # floor(n^(1/3)) Newey-West rule
    kernel="bartlett",
)

# HACResult separates the three scales explicitly:
print(f"Long-run variance (Omega):      {result.long_run_variance:.4f}")
print(f"Variance of the mean (Omega/n): {result.variance:.4f}")
print(f"HAC SE = sqrt(Omega/n):         {result.se:.4f}")
print(f"Bandwidth used:                 {result.bandwidth}")

The TemporalPLRDML framework

Algorithm overview

TemporalPLRDML combines time series cross-validation with HAC inference:

Algorithm: TemporalPLRDML Estimation
Require: Outcome Y, Treatment T, Confounders X, number of lags L
Ensure: Treatment effect τ̂ with HAC standard error

 1. Create lagged features:  X̃_t = (X_t, T_{t-1}, ..., T_{t-L})
 2. Initialize: time series cross-validator with gap g
 3. for each fold k = 1, ..., K:
 4.     Train outcome model   m̂^(-k) on training data
 5.     Train treatment model ℓ̂^(-k) on training data
 6.     Predict on test fold:  Ŷ^(-k)_t, T̂^(-k)_t
 7. Exclude uncovered early rows: drop observations without valid
    temporal out-of-fold predictions
 8. Compute residuals:  Ỹ_t = Y_t - Ŷ^(-t)_t,  T̃_t = T_t - T̂^(-t)_t
 9. Estimate τ:  τ̂ = (Σ_t Ỹ_t T̃_t) / (Σ_t T̃_t²)
10. Compute influence scores:  ψ_t = (Ỹ_t - τ̂ T̃_t) T̃_t / mean(T̃²)
11. HAC variance:  Ω̂ = γ̂_0 + 2 Σ_{k=1}^B w_k γ̂_k
12. Standard error:  SE(τ̂) = sqrt(Ω̂ / n)
return τ̂, SE(τ̂)

Implementation

from dml_ts import TemporalPLRDML
import numpy as np

# Generate time series data with autocorrelation
np.random.seed(42)
n = 500
time = np.arange(n)

# Confounders with temporal structure
X = np.column_stack([
    np.random.randn(n),
    np.sin(2 * np.pi * time / 100)  # Seasonal
])

# Autocorrelated treatment
T = np.zeros(n)
T[0] = np.random.randn()
for t in range(1, n):
    T[t] = 0.3 * T[t-1] + 0.5 * X[t, 0] + np.random.randn()

# Outcome with true effect tau = 2.0
true_tau = 2.0
Y = true_tau * T + X[:, 0]**2 + np.random.randn(n)

# Fit TemporalPLRDML
model = TemporalPLRDML(
    n_lags=3,
    model_y="random_forest",
    model_t="random_forest",
    cv_strategy="time_series_split",
    hac_kernel="bartlett",
    random_state=42,
)

result = model.fit(Y, T, X, time_index=time)
print(result.summary())

Output interpretation

Temporal PLR DML Results
========================
Treatment Effect (θ):    2.1831
HAC Standard Error:      0.1060
t-statistic:             20.60
p-value:                 0.0000
95% Confidence Interval: [1.9754, 2.3907]

Sample Information:
  Observations:          500
  Used observations:     410
  Lag rows dropped:      3
  CV rows dropped:       87

Nuisance Model Diagnostics:
  Outcome R² (CV):       0.211
  Treatment R² (CV):     0.081

HAC Inference:
  Bandwidth:             7
  CV Strategy:           time_series_split

Monte Carlo validation

We validate TemporalPLRDML through simulation studies examining:

1. Unbiasedness: $\mathbb{E}[\hat{\tau}] \approx \tau_0$.
2. Coverage: the 95% CI covers $\tau_0$ in $\approx 95\%$ of simulations.
3. SE calibration: the reported SE matches the empirical standard deviation.

Data generating process

def generate_ts_dgp(n, tau=2.0, rho_T=0.3, rho_eps=0.2, seed=None):
    """
    Generate time series DGP with autocorrelated treatment and errors.

    Parameters
    ----------
    n : int
        Sample size
    tau : float
        True treatment effect
    rho_T : float
        AR(1) coefficient for treatment
    rho_eps : float
        AR(1) coefficient for outcome error

    Returns
    -------
    Y, T, X, time : arrays
    """
    rng = np.random.default_rng(seed)
    time = np.arange(n)

    # Confounders
    X = rng.standard_normal((n, 3))

    # AR(1) treatment
    T = np.zeros(n)
    T[0] = rng.standard_normal()
    for t in range(1, n):
        T[t] = rho_T * T[t-1] + 0.5 * X[t, 0] + rng.standard_normal()

    # AR(1) errors
    eps = np.zeros(n)
    eps[0] = rng.standard_normal()
    for t in range(1, n):
        eps[t] = rho_eps * eps[t-1] + rng.standard_normal()

    # Outcome
    Y = tau * T + X[:, 0]**2 + 0.5 * X[:, 1] + eps

    return Y, T, X, time

Simulation results

TemporalPLRDML Monte Carlo results ($n=500$, 200 simulations):

Practical recommendations

When to use TemporalPLRDML

1. Time series data: any setting with temporal ordering.
2. Autocorrelated outcomes: check with Durbin–Watson or Ljung–Box tests.
3. Lagged treatment controls: when prior treatments are confounders or state variables for the current scalar effect.
4. Macro/finance applications: interest rates, prices, policies over time.

Configuration guidelines

TemporalPLRDML configuration guide

Diagnostic checks

Before reporting results, verify:

1. Nuisance model fit: check cross-validated $R^2$ for both outcome and treatment models. Poor fit ($R^2 < 0.1$) suggests missing confounders or model misspecification.
2. Residual autocorrelation: plot the autocorrelation function of influence scores. Significant spikes beyond the HAC bandwidth indicate potential issues.
3. Bandwidth sensitivity: re-estimate with bandwidth $\pm 50\%$. Large SE changes suggest autocorrelation structure uncertainty.

Exercises

Exercise 5.1: time series cross-validation

Consider a dataset with $n=100$ observations indexed $t=1,\ldots,100$. You want to use 5-fold time series CV with a gap of 2 periods.

(a) What are the train and test indices for each fold using expanding window CV? (b) How many observations are in each test fold? (c) Why is the gap parameter important for forecasting applications?

Solution.

(a) With 5 folds and gap=2:

• Fold 1: Train $\{1,\ldots,16\}$, Test $\{19,\ldots,36\}$ (after gap).
• Fold 2: Train $\{1,\ldots,34\}$, Test $\{37,\ldots,54\}$.
• …continuing the pattern.

(b) Each test fold contains approximately $n/K \approx 20$ observations (exact depends on rounding).

(c) The gap prevents data leakage when there’s a delay between when information is available and when predictions are needed (e.g., quarterly economic data released with lag).

Exercise 5.2: HAC bandwidth selection

You estimate a treatment effect with $n=1000$ observations and observe AR(1) autocorrelation $\hat{\rho} = 0.4$ in the influence scores.

(a) Calculate the rule-of-thumb bandwidth $B^* = \lfloor n^{1/3} \rfloor$. (b) Calculate the Andrews data-driven bandwidth. (c) Which would you recommend and why?

Solution.

(a) $B^* = \lfloor 1000^{1/3} \rfloor = \lfloor 10 \rfloor = 10$.

(b) $B_{\text{Andrews}} = 1.1447 \left( \frac{0.4^2 \times 1000}{(1-0.4)^4} \right)^{1/3} = 1.1447 \times \left( \frac{160}{0.1296} \right)^{1/3} \approx 12$.

(c) Both suggest similar bandwidths ($\approx 12$–$13$). The Andrews method is preferred when you have a reliable AR(1) estimate, as it adapts to the actual autocorrelation structure.

Exercise 5.3: comparing standard errors

A researcher estimates $\hat{\tau} = 1.5$ using DML and reports two standard errors: an i.i.d. influence function SE of $0.10$ and a HAC (Newey–West) SE of $0.18$.

(a) Compute the 95% confidence intervals under each assumption. (b) What does the difference in SEs tell you about the data? (c) Which CI should be reported in a time series setting?

Solution.

(a) i.i.d.: $1.5 \pm 1.96 \times 0.10 = [1.30, 1.70]$. HAC: $1.5 \pm 1.96 \times 0.18 = [1.15, 1.85]$.

(b) The HAC SE being nearly twice the i.i.d. SE indicates substantial positive autocorrelation in the influence scores. The i.i.d. assumption severely understates uncertainty.

(c) The HAC confidence interval $[1.15, 1.85]$ should be reported. Using the narrower i.i.d. interval would claim false precision.

Summary

This chapter extended DML to time series settings through three key innovations:

1. Time series cross-validation: expanding/sliding windows with purging prevent data leakage while maintaining temporal ordering.
2. HAC standard errors: Newey–West estimation accounts for autocorrelation in influence scores, providing correctly sized confidence intervals.
3. TemporalPLRDML framework: integrates lagged controls, temporal CV, and HAC inference into a scalar PLR estimation procedure.

Key takeaways:

• Standard random-fold DML can create future-data leakage and misleading uncertainty for ordered data.
• Time series CV with gaps prevents future information leakage.
• HAC SEs with bandwidth $B \approx n^{1/3}$ correct for autocorrelation.
• TemporalPLRDML reports lag rows and temporal-CV rows dropped before the final residual regression.

Roadmap to Chapter 6

Chapter 6 extends our time series framework in two directions:

1. Panel DML: combining fixed effects with DML for panel data (repeated cross-sections over time).
2. Rolling Window DML: estimating time-varying treatment effects when the causal relationship is non-stationary.

These methods are essential for applications where treatment effects may differ across individuals or change over time — common in insurance pricing, macroeconomic policy, and financial markets.