SHAP in Practice: Explanations Over Time

66.1 Explaining Every Real-Time Forecast (28 min)

Day 65 explained one model fitted on all data. But in the backtest, each origin's forecast came from a model refitted on data available at that time. To narrate history honestly ("in mid-2021 the model forecast X because Y"), explain each forecast with the model that made it.

Extend the harness in forecasting/evaluate.py with an optional explainer callback:

def rolling_origin(X, y, model_factory, horizons, test_start, model_name,
                   refit_every=1, min_train=60, explainer=None):
    """… as Day 61; if `explainer(model, X_row)` is given, store its output per forecast."""
    rows, explanations = [], []
    feature_ok = X.notna().all(axis=1)
    origins = X.index[X.index >= pd.Timestamp(test_start)]
    last_obs = y.index.max()

    for h in horizons:
        target = y.shift(-h)
        target_dates = X.index.shift(h, freq="MS")
        model = None
        for i, t in enumerate(origins):
            t_target = t + pd.DateOffset(months=h)
            if t_target > last_obs:
                break
            train = feature_ok & target.notna() & (target_dates <= t)
            if train.sum() < min_train or not feature_ok.loc[t]:
                continue
            if model is None or i % refit_every == 0:
                model = model_factory()
                model.fit(X[train], target[train])
            x_row = X.loc[[t]]
            pred = float(model.predict(x_row)[0])
            rows.append({"model": model_name, "h": h, "origin": t, "target_date": t_target,
                         "forecast": pred, "actual": float(y.loc[t_target])})
            if explainer is not None:
                e = explainer(model, x_row).iloc[0]
                explanations += [{"model": model_name, "h": h, "origin": t, "feature": f,
                                  "contribution": float(v)} for f, v in e.items()]

    res = pd.DataFrame(rows)
    if explainer is not None:
        res.attrs["explanations"] = pd.DataFrame(explanations)
    return res

Run it:

from forecasting.explain import xgb_contributions

res = rolling_origin(X, df["cpi_yoy"], xgb_factory, [1, 3, 6, 12], "2012-01-01",
                     "xgboost", refit_every=3, explainer=xgb_contributions)
expl = res.attrs["explanations"]
expl.to_parquet("data/forecasts/explanations.parquet", index=False)
res.to_parquet("data/forecasts/xgb_explained_forecasts.parquet", index=False)

Verify additivity per forecast: for each (h, origin), the contributions (including base_value) must sum to the stored forecast:

recon = expl.groupby(["h", "origin"])["contribution"].sum()
check = res.set_index(["h", "origin"])["forecast"]
print("max additivity gap:", (recon - check).abs().max())

66.2 Feature Families: A Forecast Decomposition Like CPI Contributions (20 min)

Twelve individual lags are too detailed for a briefing. Because SHAP values are additive, you can sum them over groups of related features, and the group totals still add up to the forecast.

Add to forecasting/explain.py:

FAMILIES = {
    "persistence": lambda f: f.startswith("yoy_lag") or f.startswith("yoy_ma"),
    "momentum":    lambda f: f.startswith("mom_"),
    "seasonality": lambda f: f == "month_of_year",
}


def family_of(feature: str) -> str:
    if feature == "base_value":
        return "base_value"
    for name, rule in FAMILIES.items():
        if rule(feature):
            return name
    return "other"


def family_decomposition(expl: pd.DataFrame) -> pd.DataFrame:
    """Wide table: one row per (h, origin), one column per family; rows sum to the forecast."""
    return (expl.assign(family=expl["feature"].map(family_of))
                .groupby(["h", "origin", "family"])["contribution"].sum()
                .unstack("family").reset_index())

Plot the h = 12 decomposition over time as stacked bars around the base value, with the forecast and the eventual outcome as lines. It's the same visual language as the notebook's CPI component chart (Module 05) and Day 19's contributions:

dec = family_decomposition(expl).query("h == 12").set_index("origin").sort_index()
fc = res.query("h == 12").set_index("origin")

fig, ax = plt.subplots(figsize=(11, 4.5))
bottom_pos = dec["base_value"].copy()
bottom_neg = dec["base_value"].copy()
for fam, colour in [("persistence", "#4A90D9"), ("momentum", "#D4924A"), ("seasonality", "#7B5EA7")]:
    v = dec.get(fam, 0)
    pos, neg = v.clip(lower=0), v.clip(upper=0)
    ax.bar(dec.index, pos, bottom=bottom_pos, width=25, color=colour, label=fam)
    ax.bar(dec.index, neg, bottom=bottom_neg, width=25, color=colour)
    bottom_pos += pos
    bottom_neg += neg
ax.plot(dec.index, dec["base_value"], color="grey", lw=1, label="base value")
ax.plot(fc.index, fc["forecast"], color="black", lw=1.5, label="forecast (12m ahead)")
ax.plot(fc.index, fc["actual"], color="#C95555", lw=1.5, ls="--", label="outcome")
ax.set_title("What drove the 12-month-ahead CPI inflation forecast?")
ax.legend(ncol=3, fontsize=8); plt.tight_layout()
plt.savefig("data/forecasts/decomposition_h12.png", dpi=120)

Questions for discussion:

• In 2021, when inflation began to rise, which family pushed the forecast up first?
• When the forecast badly missed the 2022 peak, what did the decomposition show? Was the model unable to forecast that far above its training range (Day 62.5)?

66.3 Non-Linearities and Interactions: Dependence Plots (23 min)

A dependence plot shows how a feature's SHAP value varies with the feature's value: the model's learned response curve, with interactions revealed by colour.

shap.plots.scatter(explanation[:, "yoy_lag0"], color=explanation[:, "mom_ma3"], show=False)
plt.title("Effect of current inflation on the 12-month forecast, coloured by momentum")
plt.tight_layout(); plt.savefig("data/forecasts/dependence_yoy_lag0.png", dpi=120); plt.close()

(explanation is the Day 65 object built from a full-sample model.)

How to read it:

• A straight upward line means a linear effect.
• Flattening at high values means the effect saturates (and remember that trees can't extrapolate beyond training values).
• Vertical spread with colour structure means an interaction: the effect of current inflation depends on momentum.

Try it 66.3

Make dependence plots for mom_ma3 and month_of_year. For month_of_year, are there months whose SHAP values are systematically positive or negative? Given that the target is yoy inflation (seasonally neutral), why might month-of-year still matter?

66.4 Explanation Drift: Are the Model's Drivers Stable? (20 min)

If a model's explanations change sharply from one refit to the next, its forecasts rest on unstable foundations, even if its accuracy looks fine. Monitor mean |SHAP| by family per year:

drift = (expl.assign(family=expl["feature"].map(family_of), year=expl["origin"].dt.year)
             .query("family != 'base_value' and h == 12")
             .assign(abs_c=lambda d: d["contribution"].abs())
             .groupby(["year", "family"])["abs_c"].mean().unstack("family"))
print(drift.round(2))
drift.plot(kind="line", marker="o", figsize=(9, 3.5),
           title="Average absolute contribution by family (h=12)")
plt.tight_layout(); plt.savefig("data/forecasts/explanation_drift.png", dpi=120)

Interpreting drift:

• Gradual changes can reflect genuine economic change, e.g. momentum mattering more in high-inflation periods.
• Abrupt jumps right after a refit point to instability: noisy hyperparameters (Day 64), too little data, or a data problem.
• A feature family that suddenly dominates deserves a data-quality check first (Module 5): a broken input often shows up as a strange explanation before it shows up as a large error.

This is a quality signal for models, the counterpart of Day 47's quality history for data. On Days 68–70 it becomes part of the forecasting pipeline's monitoring.

66.5 Mini-Project: A "Forecast Explained" Dashboard Page (45 min)

Make explanations a data product, visible to economists without running Python.

Step 1 — Publish the latest explained forecasts

Write forecasting/publish_explanations.py (in .venv-ml) that:

1. fits XGBoost at the latest origin for each horizon (training on all rows whose targets are known);
2. produces the forecast, an 80% interval (empirical, from the backtest errors in xgb_explained_forecasts.parquet), and SHAP contributions;
3. writes data/forecasts/latest_forecasts.parquet (one row per horizon) and data/forecasts/latest_explanations.parquet (one row per horizon × feature, with family).

Step 2 — A dashboard page

Add a function forecast_explained(...) to dashboards/build_dashboards.py (Day 60) that reads those files and draws, for the latest origin:

• a forecast path for horizons 1–12 with the 80% band, after the last 36 months of actual inflation, with a dashed line at the 2% target;
• a family decomposition bar chart per horizon (persistence, momentum, seasonality) around the base value;
• a header with the forecast origin date, the model name and its backtest relative RMSE (from the Day 63 leaderboard), and a caveat line: "Model-based statistical forecast; explanations describe the model, not causal effects."

def forecast_explained(cpi: pd.DataFrame, latest: pd.DataFrame, expl: pd.DataFrame,
                       rel_rmse: float) -> Path:
    fig = make_subplots(rows=1, cols=2, column_widths=[0.6, 0.4],
                        subplot_titles=["CPI inflation: history and forecast (%)",
                                        "What drives the forecast (pp vs base)"])
    hist = cpi["cpi_yoy"].iloc[-36:]
    fig.add_trace(go.Scatter(x=hist.index, y=hist, name="actual", line=dict(color=C["text"])),
                  row=1, col=1)
    fig.add_trace(go.Scatter(x=latest["target_date"], y=latest["hi"], line=dict(width=0),
                             showlegend=False), row=1, col=1)
    fig.add_trace(go.Scatter(x=latest["target_date"], y=latest["lo"], fill="tonexty",
                             fillcolor="rgba(74,144,217,0.25)", line=dict(width=0),
                             name="80% interval"), row=1, col=1)
    fig.add_trace(go.Scatter(x=latest["target_date"], y=latest["forecast"], name="forecast",
                             line=dict(color=C["primary"], width=3)), row=1, col=1)
    fig.add_hline(y=2, line_dash="dash", line_color=C["positive"], row=1, col=1)

    fam = (expl[expl.family != "base_value"].groupby(["h", "family"])["contribution"].sum()
               .unstack("family"))
    for f, colour in [("persistence", C["primary"]), ("momentum", C["warning"]),
                      ("seasonality", "#7B5EA7")]:
        if f in fam:
            fig.add_trace(go.Bar(x=fam.index, y=fam[f], name=f, marker_color=colour), row=1, col=2)
    fig.update_layout(**LAYOUT, barmode="relative", height=520,
                      title=dict(text=f"<b>CPI inflation forecast, explained</b><br><sup>"
                                      f"origin {latest['origin'].iloc[0]:%b %Y} · XGBoost · "
                                      f"backtest RMSE vs naive {rel_rmse:.2f} · model-based "
                                      f"statistical forecast; explanations describe the model, "
                                      f"not causal effects</sup>"))
    out = SITE / "forecast_explained.html"
    fig.write_html(out)
    return out

On Day 70, publish_explanations.py and this page become tasks in the Prefect pipeline. For now, run them by hand and view the page with python -m http.server 8090 --directory site.

Day 66 checklist

• Every backtest forecast is explained by the model that made it, with additivity verified
• I aggregate SHAP values into feature families and plot a forecast decomposition over time
• I read dependence plots for non-linearities and interactions
• I monitor explanation drift and know what abrupt changes may signal
• A dashboard page shows the latest forecast, its interval, its drivers and an explicit caveat

git add forecasting dashboards notebooks && git commit -m "Day 66: explanations over time, families, drift, forecast dashboard" && git push