The Partially Linear Model in Practice: A Dress Rehearsal

77.1 From design_spec.yml to an Estimation Dataset (25 min)

The Day 73 spec defines treatment, outcome, controls (with lags and transforms) and exclusions. Build the dataset from the spec, reusing Day 68's availability logic, so the estimated model is exactly the reviewed design.

Create causal/design.py:

"""Build the estimation dataset defined in causal/design_spec.yml."""
from pathlib import Path

import numpy as np
import pandas as pd
import yaml

from forecasting.availability import load_monthly_panel, transform

SPEC_PATH = Path(__file__).resolve().parent / "design_spec.yml"


def load_design(path=SPEC_PATH) -> dict:
    with open(path) as f:
        return yaml.safe_load(f)


def transform_ext(x: pd.Series, how: str | None) -> pd.Series:
    """Day 68 transforms plus log differences."""
    if how is None or how == "level":
        return x
    if how.startswith("logdiff"):
        k = int(how[7:])
        return 100 * (np.log(x) - np.log(x.shift(k)))
    return transform(x, how)


def build_controls(panel: pd.DataFrame, spec: dict) -> pd.DataFrame:
    cols = {}
    for c in spec["controls"]:
        if c["column"] not in panel:
            continue                                    # e.g. housing until Day 78
        series = panel[c["column"]].shift(c["lag"])
        cols[c["name"]] = transform_ext(series, c.get("transform"))
    return pd.DataFrame(cols, index=panel.index)


def treatment(panel: pd.DataFrame, months: int = 1) -> pd.Series:
    """Change in the month-end overnight rate over `months` months ending in t (pp)."""
    return panel["overnight_rate"].diff(months).rename("d")


def outcome_change(series: pd.Series, h: int, log: bool) -> pd.Series:
    """From t−1 to t+h: log-growth (%) for levels like housing; pp change for rates."""
    if log:
        return (100 * (np.log(series.shift(-h)) - np.log(series.shift(1)))).rename("y")
    return (series.shift(-h) - series.shift(1)).rename("y")


def estimation_frame(panel, spec, outcome_series, h, log_outcome, treat_months=1,
                     start=None, exclude=None) -> pd.DataFrame:
    df = (build_controls(panel, spec)
          .join(treatment(panel, treat_months))
          .join(outcome_change(outcome_series, h, log_outcome)))
    if start:
        df = df.loc[start:]
    if exclude:
        df = df.drop(df.loc[exclude[0]:exclude[1]].index, errors="ignore")
    return df.dropna()

Every choice (lags, transforms, sample, exclusions) is driven by the spec or by explicit arguments, so the estimate can be traced back to the reviewed design. That's the same discipline as dbt's ref() and GX's specs.

77.2 A Rehearsal Outcome: Rate Changes and the Unemployment Rate (25 min)

The housing-starts series arrives on Day 78. A professional team doesn't wait: it rehearses the whole analysis on an outcome already in the platform, to debug code and diagnostics before the real data arrive. Our rehearsal question:

Effect of a 1pp overnight-rate increase on the change in the unemployment rate over the following h months.

Standard macro reasoning predicts that tightening raises unemployment with a lag, so θ should be positive at longer horizons. It's a useful sanity check on the machinery.

from sklearn.ensemble import RandomForestRegressor

from causal.design import estimation_frame, load_design
from causal.dml import fit_plr_ts
from forecasting.availability import load_monthly_panel

spec = load_design()
panel = load_monthly_panel()
rf = RandomForestRegressor(300, min_samples_leaf=5, max_features=0.5, random_state=0, n_jobs=2)

df = estimation_frame(panel, spec, panel["unemployment_rate"], h=12, log_outcome=False)
X = df.drop(columns=["d", "y"])
res = fit_plr_ts(df["y"], df["d"], X, rf, rf, h=12)
print(f"n={res['n']}  θ̂={res['theta']:+.2f} pp per 1pp  "
      f"95% CI [{res['ci_lo']:+.2f}, {res['ci_hi']:+.2f}]  "
      f"R²(E[D|X])={res['r2_m']:.2f}  R²(E[Y|X])={res['r2_l']:.2f}")

(For the rehearsal, the spec's lagged-outcome control is housing, which doesn't exist yet. Add a lagged unemployment change control for this outcome, e.g. unemployment_rate with transform: diff6, lag: 1. The spec already includes unemployment controls.)

Interpreting the rehearsal

• Sign and size: does a tightening raise unemployment a year later? By roughly how many percentage points per 1pp of tightening?
• Precision: how wide is the HAC interval? With ~19 years of monthly data and few tightening cycles, wide intervals are expected (Day 72.5).
• R² of E[D | X]: how much of rate setting your controls explain. Is there enough residual policy variation left?

Try it 77.2: the price puzzle

Repeat with CPI inflation as the outcome (panel["headline_yoy"], pp change). Early empirical work on monetary policy often found prices rising after a tightening, the famous "price puzzle" (Sims, 1992), widely interpreted as a sign that the central bank reacted to information about future inflation that the econometrician didn't control for. What do you find? If θ̂ is positive at short horizons, which unobserved confounder from your DAG (Day 73) is the likely culprit, and how might you proxy it (e.g. oil prices, core inflation momentum, market-based expectations)?

77.3 Learners, Diagnostics and the Reaction-Function Fit (20 min)

Run the rehearsal with several learners and report diagnostics alongside estimates:

from sklearn.linear_model import LassoCV, LinearRegression
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler
from xgboost import XGBRegressor

LEARNERS = {
    "ols (classic regression)": LinearRegression(),
    "lasso": make_pipeline(StandardScaler(), LassoCV(cv=5, max_iter=5000)),
    "random_forest": rf,
    "xgboost": XGBRegressor(n_estimators=300, learning_rate=0.03, max_depth=3, subsample=0.8,
                            colsample_bytree=0.8, random_state=0, n_jobs=2),
}
rows = []
for name, L in LEARNERS.items():
    r = fit_plr_ts(df["y"], df["d"], X, L, L, h=12)
    rows.append({"learner": name, **{k: r[k] for k in ("theta", "se", "ci_lo", "ci_hi", "r2_l", "r2_m")}})
table = pd.DataFrame(rows)
print(table.round(3).to_string(index=False))

77.4 Dynamic Responses Across Horizons (20 min)

Estimating θ for several horizons traces the dynamic effect of a rate change, in the spirit of local projections (Jordà, 2005):

import matplotlib.pyplot as plt

H = [3, 6, 12, 18, 24]
dyn = []
for h in H:
    dfh = estimation_frame(panel, spec, panel["unemployment_rate"], h=h, log_outcome=False)
    r = fit_plr_ts(dfh["y"], dfh["d"], dfh.drop(columns=["d", "y"]), rf, rf, h=h)
    dyn.append({"h": h, **r})
dyn = pd.DataFrame(dyn)
dyn.to_parquet("data/causal/rehearsal_unemployment_dynamic.parquet", index=False)

fig, ax = plt.subplots(figsize=(8, 4))
ax.plot(dyn["h"], dyn["theta"], marker="o", label="θ̂ (per 1pp rate increase)")
ax.fill_between(dyn["h"], dyn["ci_lo"], dyn["ci_hi"], alpha=0.25, label="95% CI (HAC)")
ax.axhline(0, color="grey", lw=1)
ax.set_xlabel("months after the rate change"); ax.set_ylabel("change in unemployment rate (pp)")
ax.set_title("Rehearsal: dynamic effect of a 1pp tightening on unemployment (DML, PLR)")
ax.legend(); plt.tight_layout(); plt.savefig("data/causal/rehearsal_dynamic.png", dpi=120)

Note that each horizon is a separate estimate with its own sample (longer horizons lose the most recent months). Intervals at different horizons are not independent: the outcomes overlap. Read the plot as a profile, not as a set of separate tests.

77.5 Pre-Registered Robustness Checks (23 min)

Day 73 listed robustness checks before any estimates. Run them now on the rehearsal outcome, exactly as written:

def placebo_outcome(series, lead_start=12, lead_end=6):
    """Change from t−lead_start to t−lead_end: entirely before the decision at t."""
    return (series.shift(lead_end) - series.shift(lead_start)).rename("y")


checks = {}
base = estimation_frame(panel, spec, panel["unemployment_rate"], h=12, log_outcome=False)
Xb = base.drop(columns=["d", "y"])

pl = Xb.join(base["d"]).join(placebo_outcome(panel["unemployment_rate"])).dropna()
checks["placebo outcome"] = fit_plr_ts(pl["y"], pl["d"], pl.drop(columns=["d", "y"]), rf, rf, h=6)

future_d = panel["overnight_rate"].diff().shift(-12).rename("d")
nc = Xb.join(future_d).join(base["y"]).dropna()
checks["negative-control treatment"] = fit_plr_ts(nc["y"], nc["d"], nc.drop(columns=["d", "y"]),
                                                  rf, rf, h=12)

alt = estimation_frame(panel, spec, panel["unemployment_rate"], h=12, log_outcome=False,
                       treat_months=3)
checks["3-month treatment"] = fit_plr_ts(alt["y"], alt["d"], alt.drop(columns=["d", "y"]),
                                         rf, rf, h=12)

ex = estimation_frame(panel, spec, panel["unemployment_rate"], h=12, log_outcome=False,
                      exclude=("2020-03-01", "2021-06-01"))
checks["excl. pandemic"] = fit_plr_ts(ex["y"], ex["d"], ex.drop(columns=["d", "y"]), rf, rf, h=12)

print(pd.DataFrame(checks).T[["theta", "ci_lo", "ci_hi", "n"]].round(3))

How to react to a failed placebo: a clearly non-zero placebo effect means the design has a problem (remaining confounding or anticipation), and the main estimate should be reported with that caveat, not quietly dropped. Pre-registration is what makes this honest: you can't choose which checks to show after seeing them.

77.6 A Specification Curve and run_design.py (23 min)

Many defensible choices exist: learners, buffer sizes, treatment definitions, sample periods. Rather than report one favourite, show all reasonable combinations sorted by estimate, a specification curve (Simonsohn, Simmons & Nelson, 2020). If most specifications agree in sign and rough size, the conclusion is robust; if they scatter across zero, it isn't.

import itertools

grid = {
    "learner": ["lasso", "random_forest", "xgboost"],
    "treat_months": [1, 3],
    "exclude_pandemic": [False, True],
    "buffer_extra": [6, 18],
}
rows = []
for learner, tm, exc, bx in itertools.product(*grid.values()):
    dfi = estimation_frame(panel, spec, panel["unemployment_rate"], h=12, log_outcome=False,
                           treat_months=tm,
                           exclude=("2020-03-01", "2021-06-01") if exc else None)
    L = LEARNERS[learner]
    r = fit_plr_ts(dfi["y"], dfi["d"], dfi.drop(columns=["d", "y"]), L, L, h=12, buffer=12 + bx)
    rows.append({"learner": learner, "treat_months": tm, "exclude_pandemic": exc,
                 "buffer": 12 + bx, **{k: r[k] for k in ("theta", "ci_lo", "ci_hi")}})
curve = pd.DataFrame(rows).sort_values("theta").reset_index(drop=True)

fig, ax = plt.subplots(figsize=(9, 3.5))
ax.errorbar(curve.index, curve["theta"],
            yerr=[curve["theta"] - curve["ci_lo"], curve["ci_hi"] - curve["theta"]], fmt="o")
ax.axhline(0, color="grey")
ax.set_xlabel("specification (sorted)"); ax.set_ylabel("θ̂ at h=12")
ax.set_title("Specification curve (rehearsal outcome: unemployment)")
plt.tight_layout(); plt.savefig("data/causal/rehearsal_spec_curve.png", dpi=120)
print(f"{(curve.theta > 0).mean():.0%} of specifications positive; "
      f"{((curve.ci_lo > 0) | (curve.ci_hi < 0)).mean():.0%} exclude zero")

causal/run_design.py

Package everything into one script that reads the spec, takes the outcome as an argument, and writes results plus a Markdown report: main estimates by horizon, learner table, robustness checks, specification-curve summary, diagnostics and package versions. On Day 78 the only change will be the outcome argument:

python -m causal.run_design --outcome unemployment_rate --outcome-type change   # rehearsal (today)
python -m causal.run_design --outcome housing_starts --outcome-type log          # Days 78–80

"""Run the causal design end to end and write results + a report."""
import argparse
import itertools
from datetime import datetime, timezone
from pathlib import Path

import doubleml
import pandas as pd
import sklearn

from causal.design import estimation_frame, load_design
from causal.dml import fit_plr_ts
from forecasting.availability import load_monthly_panel

OUT = Path("data/causal")


def main():
    ap = argparse.ArgumentParser()
    ap.add_argument("--outcome", required=True)
    ap.add_argument("--outcome-type", choices=["log", "change"], required=True)
    ap.add_argument("--horizons", default="3,6,12,24")
    args = ap.parse_args()

    spec, panel = load_design(), load_monthly_panel()
    series = panel[args.outcome]
    log = args.outcome_type == "log"
    H = [int(h) for h in args.horizons.split(",")]
    learners = build_learners()                      # as in 77.3

    main_rows = []
    for h in H:
        df = estimation_frame(panel, spec, series, h=h, log_outcome=log)
        r = fit_plr_ts(df["y"], df["d"], df.drop(columns=["d", "y"]),
                       learners["random_forest"], learners["random_forest"], h=h)
        main_rows.append({"h": h, **r})
    main = pd.DataFrame(main_rows)
    checks = run_robustness(panel, spec, series, log, learners)   # as in 77.5
    curve = run_spec_curve(panel, spec, series, log, learners)    # as in 77.6

    tag = f"{args.outcome}_{datetime.now(timezone.utc):%Y%m%dT%H%M%S}"
    OUT.mkdir(parents=True, exist_ok=True)
    main.to_parquet(OUT / f"main_{tag}.parquet"); checks.to_parquet(OUT / f"checks_{tag}.parquet")
    curve.to_parquet(OUT / f"curve_{tag}.parquet")
    (OUT / f"report_{tag}.md").write_text("\n\n".join([
        f"# Causal design results — outcome: {args.outcome}",
        f"doubleml {doubleml.__version__} · scikit-learn {sklearn.__version__} · "
        "estimator: PLR, blocked cross-fitting, HAC SEs, median over repeated splits",
        "## Main estimates", main.round(3).to_markdown(index=False),
        "## Robustness (pre-registered)", checks.round(3).to_markdown(),
        "## Specification curve", f"{(curve.theta > 0).mean():.0%} positive",
    ]))
    print(f"✓ results written with tag {tag}")


if __name__ == "__main__":
    main()

Recap: Days 74–77

Looking ahead: Days 78–80

• Day 78: add housing starts to the platform through the full path you've used all course: source (dlt or WAP), staging model, tests, mart, GX expectations. Then rerun run_design.py with the real outcome.
• Day 79: interpret the results, run sensitivity analysis for unobserved confounders (how strong would a missing variable need to be to overturn the result?), and explore heterogeneity (e.g. high vs low inflation regimes).
• Day 80: integrate the causal analysis into the orchestrated platform and dashboards, and complete the final capstone.

Day 77 checklist

• My estimation dataset is built from design_spec.yml, with every choice traceable
• I rehearsed the full analysis on a real outcome, and can interpret sign, size, precision and reaction-function fit
• I compared learners and reported diagnostics alongside estimates
• I traced dynamic effects across horizons with HAC intervals
• I ran every pre-registered robustness check and a specification curve
• run_design.py produces results and a report, and will only need a new outcome on Day 78

git add causal notebooks && git commit -m "Day 77: design-driven estimation, rehearsal, robustness, spec curve, run_design" && git push