Day 76: DML for Macroeconomic Time Series · Applied Decision Intelligence & Data Architecture

Written lesson

Time	Block
0:00 – 0:20	76.1 What breaks with dependent data
0:20 – 0:48	76.2 Blocked cross-fitting with buffers
0:48 – 0:55	Break
0:55 – 1:20	76.3 HAC inference for overlapping horizons
1:20 – 1:45	76.4 Monte Carlo: iid vs time-series-aware DML
1:45 – 1:53	Break
1:53 – 2:10	76.5 Custom splits in `doubleml`
2:10 – 2:30	76.6 Mini-project: `fit_plr_ts`, a time-series DML wrapper

76.1 What Breaks with Dependent Data (20 min)

The DML theory of Day 75 assumes independent observations. Monthly macro data violate this in three ways:

Feature of our data	Consequence for DML
persistent controls (inflation, unemployment move slowly)	neighbouring months are near-duplicates. Random folds put a month's twin in the training set, so "out-of-fold" predictions are effectively in-sample and overfitting bias returns
overlapping outcomes (housing growth from t−1 to t+h shares h months with t+1's outcome)	residuals are serially correlated by construction (an MA(h) structure)
gradual policy (rate changes come in series)	treatment residuals are autocorrelated too

Two fixes, both standard in time-series econometrics:

blocked cross-fitting with buffers: folds are contiguous blocks of time, and months near the test block are removed from training;
HAC (Newey–West) standard errors for the final residual regression, which account for serial correlation (you used them on Day 63 for Diebold–Mariano tests).

76.2 Blocked Cross-Fitting with Buffers (28 min)

time →   [ fold 1 ][ fold 2 ][ fold 3 ][ fold 4 ][ fold 5 ]

test fold 3:  train = folds 1–5 EXCEPT fold 3 and a buffer of b months on each side
              [ fold 1 ][ fold 2 ]xx[  TEST  ]xx[ fold 4 ][ fold 5 ]
                                ↑ buffer      ↑ buffer

The buffer removes the months whose information most closely overlaps with the test block. A natural choice is b ≥ h (the outcome overlap), plus a margin for the persistence of controls (e.g. 12 months).

Add to causal/dml.py:

def blocked_splits(n: int, n_folds: int = 5, buffer: int = 12, offset: int = 0):
    """Contiguous test blocks; training excludes the block ± `buffer` observations.

    `offset` rotates block boundaries, so repeated cross-fitting can use different splits.
    """
    idx = np.arange(n)
    rolled = np.roll(idx, -offset)
    splits = []
    for block in np.array_split(rolled, n_folds):
        test = np.sort(block)
        # a rotated block may wrap around the end; handle each contiguous piece
        pieces = np.split(test, np.where(np.diff(test) != 1)[0] + 1)
        exclude = set()
        for p in pieces:
            exclude |= set(range(max(0, p[0] - buffer), min(n, p[-1] + buffer + 1)))
        train = np.array([i for i in idx if i not in exclude])
        splits.append((train, test))
    return splits

Check the splits visually: for each fold, print the test range and the gap to the nearest training observation. Every test month should be at least buffer months away from any training month.

Trade-off: larger buffers protect better against dependence, but train nuisances on less data. With ~240 monthly observations, 5 folds and a 12-month buffer, each nuisance model trains on roughly 75% of the sample. Report the buffer, and check results with a smaller and a larger one (Day 77).

76.3 HAC Inference for Overlapping Horizons (25 min)

plr_from_residuals(y_res, d_res, hac_lags=L) (Day 75) already supports Newey–West standard errors. The question is the lag length L.

Source of serial correlation	Implied minimum lag
overlapping outcomes over h months	h (residuals are MA(h) by construction)
persistence of controls and policy	more; data-driven rules help

A practical rule: L = max(h, ⌊4 · (n/100)^(2/9)⌋), using the usual Newey–West automatic bandwidth as a floor. Add it:

def hac_lags_for(n: int, h: int) -> int:
    return int(max(h, np.floor(4 * (n / 100) ** (2 / 9))))

Why not trust DoubleML's default standard errors here? DoubleML's standard errors assume independent observations. With overlapping 12-month outcomes, they can be far too small, making a noisy estimate look precise. The next block measures how much.

76.4 Monte Carlo: iid vs Time-Series-Aware DML (25 min)

Build a world with the features of 76.1. Add to causal/simulate.py:

def ts_world(T: int = 300, h: int = 12, true_effect: float = -5.0, seed: int = 0) -> pd.DataFrame:
    """Monthly world: persistent confounders, gradual policy, overlapping h-month outcomes."""
    rng = np.random.default_rng(seed)
    n = T + h + 1
    c = np.zeros(n)
    s = np.zeros(n)
    v = np.zeros(n)
    for t in range(1, n):
        c[t] = 0.95 * c[t - 1] + rng.normal(0, 0.5)          # slow business cycle
        s[t] = 0.8 * s[t - 1] + rng.normal(0, 0.5)           # second persistent state
        v[t] = 0.5 * v[t - 1] + rng.normal(0, 0.15)          # gradual policy surprises
    c_lag, s_lag = np.roll(c, 1), np.roll(s, 1)
    d = 0.4 * np.tanh(2 * c_lag) + 0.1 * s_lag ** 2 + v      # reaction to last month's state
    e = rng.normal(0, 1.0, n)
    overlap = np.array([e[t + 1:t + h + 1].sum() for t in range(n - h)])
    y = (true_effect * d[:n - h] + 6 * np.sin(c_lag[:n - h]) + 2 * s_lag[:n - h] ** 2
         + overlap)
    df = pd.DataFrame({"c_lag": c_lag[:n - h], "s_lag": s_lag[:n - h],
                       "rate_change": d[:n - h], "housing_growth_h": y})
    return df.iloc[1:T + 1].reset_index(drop=True)

Compare three versions in causal/montecarlo.py:

from causal.dml import blocked_splits, dml_plr, hac_lags_for
from causal.simulate import ts_world


def ts_coverage(n_sims=100, T=300, h=12, true_effect=-5.0, learner=RF):
    rows = []
    for s in range(n_sims):
        df = ts_world(T=T, h=h, true_effect=true_effect, seed=2000 + s)
        X, y, d = df[["c_lag", "s_lag"]], df["housing_growth_h"], df["rate_change"]
        variants = {
            "iid folds + robust SE": dml_plr(y, d, X, learner, learner, seed=s),
            "iid folds + HAC SE": dml_plr(y, d, X, learner, learner, seed=s,
                                          hac_lags=hac_lags_for(T, h)),
            "blocked folds + HAC SE": dml_plr(y, d, X, learner, learner,
                                              splits=blocked_splits(T, 5, buffer=h + 6),
                                              hac_lags=hac_lags_for(T, h)),
        }
        for name, r in variants.items():
            lo, hi = r.ci()
            rows.append({"variant": name, "theta": r.theta, "se": r.se,
                         "covered": lo <= true_effect <= hi})
    res = pd.DataFrame(rows)
    return res.groupby("variant").agg(mean_theta=("theta", "mean"), sd_theta=("theta", "std"),
                                      mean_se=("se", "mean"), coverage_95=("covered", "mean"))

Typical results, which you should verify:

Variant	Typical coverage	Why
iid folds + robust SE	well below 95%	SEs ignore serial correlation; neighbours leak across folds
iid folds + HAC SE	better, often still below	SEs fixed; some leakage across folds remains
blocked folds + HAC SE	close to 95%	both problems addressed

Also compare mean_se with sd_theta for each variant: the iid version's reported SE is typically much smaller than the actual sampling variability. That's an overconfident result, the causal-inference version of Day 63's under-covering forecast intervals.

The practical lesson for Days 78–80: with monthly macro data and multi-month outcomes, always use blocked folds with buffers and HAC standard errors, and expect wider intervals than an iid analysis would suggest. Those wider intervals are the honest ones.

76.5 Custom Splits in `doubleml` (18 min)

DoubleML accepts your own splits through set_sample_splitting, which lets you combine its machinery (repetitions, model variants, sensitivity analysis on Day 79) with blocked folds:

import doubleml as dml
from sklearn.base import clone

from causal.dml import blocked_splits
from causal.simulate import ts_world

df = ts_world(T=300, h=12, seed=7)
data = dml.DoubleMLData(df, y_col="housing_growth_h", d_cols="rate_change",
                        x_cols=["c_lag", "s_lag"])

reps = [blocked_splits(len(df), n_folds=5, buffer=18, offset=o) for o in (0, 11, 23)]
plr = dml.DoubleMLPLR(data, ml_l=clone(RF), ml_m=clone(RF), score="partialling out")
plr.set_sample_splitting(all_smpls=reps)          # 3 repetitions × 5 blocked folds
plr.fit(store_predictions=True)
print(plr.summary)

DoubleML's reported SE is still the independence version. Compute HAC inference from its stored, cross-fitted predictions:

import numpy as np

from causal.dml import aggregate_reps, hac_lags_for, plr_from_residuals

y = df["housing_growth_h"].to_numpy()
d = df["rate_change"].to_numpy()
thetas, ses = [], []
for r in range(plr.n_rep):
    y_res = y - plr.predictions["ml_l"][:, r, 0]
    d_res = d - plr.predictions["ml_m"][:, r, 0]
    t, s = plr_from_residuals(y_res, d_res, hac_lags=hac_lags_for(len(y), 12))
    thetas.append(t); ses.append(s)
print("θ̂ (median over reps), HAC SE:", aggregate_reps(thetas, ses))

Version note: some DoubleML versions require that each repetition's test folds partition the sample (ours do) and may validate the train/test structure. If your version rejects buffered training sets, use causal/dml.py (which supports any splits) as the primary estimator, and DoubleML with ordinary blocked folds (buffer=0) as a cross-check.

76.6 Mini-Project: `fit_plr_ts`, a Time-Series DML Wrapper (20 min)

Add a single entry point to causal/dml.py that applies everything from today, so Day 77–80 code stays short:

def fit_plr_ts(y, d, X, learner_l, learner_m, h: int, n_folds: int = 5,
               buffer: int | None = None, n_rep: int = 5) -> dict:
    """Time-series DML for the PLR: blocked folds with buffers, HAC SEs, repeated splits."""
    y, d, X = np.asarray(y, float), np.asarray(d, float), np.asarray(X, float)
    n = len(y)
    buffer = h + 6 if buffer is None else buffer
    lags = hac_lags_for(n, h)
    step = max(1, n // (n_folds * n_rep))
    thetas, ses, r2l, r2m = [], [], [], []
    for r in range(n_rep):
        splits = blocked_splits(n, n_folds, buffer, offset=r * step)
        y_res, d_res = cross_fit_residuals(y, d, X, learner_l, learner_m, splits)
        t, s = plr_from_residuals(y_res, d_res, hac_lags=lags)
        thetas.append(t); ses.append(s); r2l.append(r2(y, y_res)); r2m.append(r2(d, d_res))
    theta, se = aggregate_reps(thetas, ses)
    return {"theta": theta, "se": se, "ci_lo": theta - 1.96 * se, "ci_hi": theta + 1.96 * se,
            "r2_l": float(np.median(r2l)), "r2_m": float(np.median(r2m)),
            "n": n, "h": h, "buffer": buffer, "hac_lags": lags, "n_rep": n_rep}

Validate it with ts_world (mean estimate ≈ −5, coverage near 95% over 100 simulations), then add a unit test in tests/test_causal.py that runs one small simulation and checks that the true effect lies inside the estimated interval for a fixed seed. Tests protect methods as well as data.

Day 76 checklist

I can explain how persistence, overlapping outcomes and gradual policy break iid DML
Blocked cross-fitting with buffers keeps nuisance predictions genuinely out-of-sample
Final-stage HAC standard errors use a lag length tied to the horizon
My Monte Carlo shows iid DML over-confident and time-series-aware DML near nominal coverage
I can pass custom splits to doubleml and compute HAC inference from its predictions
fit_plr_ts packages the method, validated by simulation and a unit test

python -m pytest -q tests/test_causal.py
git add causal tests notebooks && git commit -m "Day 76: blocked cross-fitting, HAC inference, ts coverage, fit_plr_ts" && git push

DML for Macroeconomic Time Series

What you'll learn today

Today's plan