1"""F3-conf-03: rolling-origin ridge vs naive / expanding-mean forecasts of IP growth (h=1).2Run from workspace root: ~/wp29/venv/bin/python out/analysis.py"""3import json, os4import numpy as np5import pandas as pd67LAM = 50.08SEED = None # no randomness used9FIRST, LAST = pd.Timestamp("2021-09-01"), pd.Timestamp("2026-08-01")10TRAIN_START = pd.Timestamp("1992-03-01")11MEAN_START = pd.Timestamp("1992-04-01")12FEATS = ["ip_growth", "ip_l1", "ip_l2", "oil_chg", "fedfunds", "permit_chg", "neworder_chg", "tcu"]1314os.makedirs("out", exist_ok=True)15raw = pd.read_csv("data/ip_predictors.csv", parse_dates=["month"])16n_raw = len(raw)17d = raw.sort_values("month").set_index("month")18assert d.index.is_unique19d["ip_l1"] = d["ip_growth"].shift(1) # lags use only past values20d["ip_l2"] = d["ip_growth"].shift(2)21d["y_next"] = d["ip_growth"].shift(-1) # target of row s = ip_growth_{s+1}2223origins = pd.date_range(FIRST, LAST, freq="MS")24rows, n_train_last, excl = [], None, []25for t in origins:26 tgt = t + pd.DateOffset(months=1)27 actual = float(d.loc[tgt, "ip_growth"])28 avail = d.loc[:t] # information set: dated <= t29 # baselines30 naive = float(avail.loc[t, "ip_growth"])31 mean = float(avail.loc[MEAN_START:t, "ip_growth"].mean())32 # ridge training rows: TRAIN_START <= s, s+1 <= t33 tr = avail.loc[TRAIN_START:t - pd.DateOffset(months=1)]34 Xr = tr[FEATS].to_numpy(float)35 # y_{s+1} with s+1 <= t, taken from the available block36 yr = avail["ip_growth"].shift(-1).loc[tr.index].to_numpy(float)37 assert np.isfinite(Xr).all() and np.isfinite(yr).all(), t38 mu, sd = Xr.mean(0), Xr.std(0, ddof=1)39 ybar = yr.mean()40 Z = (Xr - mu) / sd41 b = np.linalg.solve(Z.T @ Z + LAM * np.eye(Z.shape[1]), Z.T @ (yr - ybar))42 zt = avail.loc[t, FEATS].to_numpy(float)43 assert np.isfinite(zt).all()44 ridge = float(ybar + ((zt - mu) / sd) @ b)45 for m, f in (("naive", naive), ("mean", mean), ("ridge", ridge)):46 rows.append(dict(origin_date=t.strftime("%Y-%m-%d"), target_date=tgt.strftime("%Y-%m-%d"),47 horizon=1, forecast=f, actual=actual, model=m))48 n_train_last = len(tr)4950fc = pd.DataFrame(rows)51fc.to_csv("out/forecasts.csv", index=False, float_format="%.10f")5253def metrics(m):54 e = fc.loc[fc.model == m, "forecast"] - fc.loc[fc.model == m, "actual"]55 return float(np.sqrt((e ** 2).mean())), float(e.abs().mean())5657res = {"task_id": "F3-conf-03", "n_origins": len(origins),58 "first_origin": FIRST.strftime("%Y-%m-%d"), "last_origin": LAST.strftime("%Y-%m-%d"),59 "lambda": LAM, "normalisation_window": "expanding", "n_features": len(FEATS)}60for m in ("naive", "mean", "ridge"):61 res[f"rmse_{m}"], res[f"mae_{m}"] = metrics(m)62res["rel_improvement_vs_naive"] = 1 - res["rmse_ridge"] / res["rmse_naive"]63res["rel_improvement_vs_mean"] = 1 - res["rmse_ridge"] / res["rmse_mean"]64res["baseline_best_rmse"] = min(res["rmse_naive"], res["rmse_mean"])6566# Uncertainty: Diebold-Mariano tests (squared-error loss, h=1, normal approx, Harvey et al. small-sample correction)67from scipy import stats68def dm(base):69 e_r = (fc.loc[fc.model == "ridge", "forecast"].to_numpy() - fc.loc[fc.model == "ridge", "actual"].to_numpy())70 e_b = (fc.loc[fc.model == base, "forecast"].to_numpy() - fc.loc[fc.model == base, "actual"].to_numpy())71 dd = e_b ** 2 - e_r ** 2 # >0 favours ridge72 n = len(dd)73 stat = dd.mean() / np.sqrt(dd.var(ddof=0) / n)74 stat *= np.sqrt((n + 1 - 2 + 0) / n) # HLN correction for h=175 p = 2 * stats.t.sf(abs(stat), df=n - 1)76 return float(stat), float(p)77res["dm_stat_ridge_vs_naive"], res["dm_pvalue_ridge_vs_naive"] = dm("naive")78res["dm_stat_ridge_vs_mean"], res["dm_pvalue_ridge_vs_mean"] = dm("mean")79res.update({"feasibility": "complete", "feasibility_reason": None,80 "data_issues": [81 {"issue": "neworder_chg missing for 1992-01 and 1992-02 (series starts 1992-02).",82 "handling": "Not used: earliest neworder_chg needed is 1992-03 (first training row)."},83 {"issue": "README states the first month (1992-01) of every change series is empty, but ip_growth, oil_chg and permit_chg are populated for 1992-01.",84 "handling": "Values used as given; only ip_growth 1992-01 enters (lag-2 feature of training row 1992-03)."},85 {"issue": "Extreme COVID-19 observations (e.g. ip_growth 2020-04 = -14.1, oil_chg 2020-03..05) lie in every training window.",86 "handling": "Kept, as the task prescribes the training window; no outlier rule applied."},87 {"issue": "Current-vintage data; publication lags ignored.",88 "handling": "Errors computed on revised data as instructed; no real-time validity claimed."}],89 "n_input_rows": {"ip_predictors.csv": n_raw},90 "n_analysed": int(n_train_last), "random_seed": SEED})91with open("out/results.json", "w") as f:92 json.dump(res, f, indent=2)9394# exclusion ledger: raw rows -> training rows at last origin95n_after_start = int((d.index >= TRAIN_START).sum())96n_final = n_train_last97pd.DataFrame([98 dict(step=1, rule="drop months before 1992-03-01 (first training row; lag/neworder availability)",99 n_before=n_raw, n_after=n_after_start, n_removed=n_raw - n_after_start,100 rationale="TASK.md ridge spec: training rows start at s=1992-03-01"),101 dict(step=2, rule="drop months s with s+1 > last origin 2026-08-01 (target not available at origin)",102 n_before=n_after_start, n_after=n_final, n_removed=n_after_start - n_final,103 rationale="information rule: training target y_{s+1} must be dated <= t"),104]).to_csv("out/exclusions.csv", index=False)105print(json.dumps({k: res[k] for k in res if k.startswith(("rmse", "mae", "rel_", "n_", "dm_"))}, indent=1))106