1"""F3-conf-03: rolling-origin ridge vs naive / expanding-mean forecasts of US IP growth (h=1).2Reads only data/ip_predictors.csv; writes out/results.json, out/forecasts.csv, out/exclusions.csv.3Deterministic: no randomness used (random_seed = null)."""4import json, hashlib5from pathlib import Path6import numpy as np7import pandas as pd89ROOT = Path(__file__).resolve().parent.parent10DATA = ROOT / "data" / "ip_predictors.csv"11OUT = ROOT / "out"12LAM = 50.013FIRST_ORIGIN, LAST_ORIGIN = pd.Timestamp("2021-09-01"), pd.Timestamp("2026-08-01")14TRAIN_START = pd.Timestamp("1992-03-01")15MEAN_START = pd.Timestamp("1992-04-01")16FEATS = ["ip_growth", "ip_lag1", "ip_lag2", "oil_chg", "fedfunds", "permit_chg", "neworder_chg", "tcu"]1718raw = pd.read_csv(DATA, parse_dates=["month"])19n_raw = len(raw)20d = raw.sort_values("month").reset_index(drop=True)21assert d["month"].is_unique22assert (d["month"].diff().dropna().dt.days.between(28, 31)).all(), "non-contiguous months"23d = d.set_index("month")24# Lags/leads by row position (months verified contiguous). Lags use only past values.25d["ip_lag1"] = d["ip_growth"].shift(1)26d["ip_lag2"] = d["ip_growth"].shift(2)27d["y_next"] = d["ip_growth"].shift(-1) # target of feature row s is ip_growth_{s+1}2829origins = d.index[(d.index >= FIRST_ORIGIN) & (d.index <= LAST_ORIGIN)]30assert len(origins) == 603132rows = []33n_train_last = None34for t in origins:35 t_next = d.index[d.index.get_loc(t) + 1]36 actual = float(d.loc[t_next, "ip_growth"])37 # Information set: rows dated <= t only.38 info = d.loc[d.index <= t]39 # naive40 f_naive = float(info.loc[t, "ip_growth"])41 # expanding mean 1992-04..t42 f_mean = float(info.loc[info.index >= MEAN_START, "ip_growth"].mean())43 # ridge training rows: TRAIN_START <= s, s+1 <= t (so target y_{s+1} dated <= t)44 s_idx = info.index[(info.index >= TRAIN_START) & (info.index < t)]45 tr = info.loc[s_idx, FEATS].copy()46 y = info["ip_growth"].shift(-1).loc[s_idx].to_numpy() # ip_growth_{s+1}, all dated <= t47 X = tr.to_numpy(dtype=float)48 assert np.isfinite(X).all() and np.isfinite(y).all()49 mu = X.mean(axis=0)50 sd = X.std(axis=0, ddof=1)51 Z = (X - mu) / sd52 ybar = y.mean()53 b = np.linalg.solve(Z.T @ Z + LAM * np.eye(Z.shape[1]), Z.T @ (y - ybar))54 zt = (info.loc[t, FEATS].to_numpy(dtype=float) - mu) / sd55 f_ridge = float(ybar + zt @ b)56 n_train_last = len(s_idx)57 for m, f in (("naive", f_naive), ("mean", f_mean), ("ridge", f_ridge)):58 rows.append(dict(origin_date=t.strftime("%Y-%m-%d"), target_date=t_next.strftime("%Y-%m-%d"),59 horizon=1, forecast=f, actual=actual, model=m))6061fc = pd.DataFrame(rows)62fc.to_csv(OUT / "forecasts.csv", index=False, float_format="%.10f")6364def metrics(m):65 e = fc.loc[fc.model == m, "forecast"] - fc.loc[fc.model == m, "actual"]66 return float(np.sqrt(np.mean(e ** 2))), float(np.mean(np.abs(e)))6768def dm_test(m1, m2):69 """Diebold-Mariano test on squared-error loss differential d = e1^2 - e2^2 (h=1:70 plain sample variance, i.e. Newey-West with 0 lags), with the Harvey-Leybourne-Newbold71 small-sample correction and t_{n-1} reference distribution. Positive stat => m1 worse."""72 from scipy import stats73 e1 = (fc.loc[fc.model == m1, "forecast"] - fc.loc[fc.model == m1, "actual"]).to_numpy()74 e2 = (fc.loc[fc.model == m2, "forecast"] - fc.loc[fc.model == m2, "actual"]).to_numpy()75 dd = e1 ** 2 - e2 ** 276 n = len(dd)77 stat = dd.mean() / np.sqrt(dd.var(ddof=0) / n)78 stat *= np.sqrt((n + 1 - 2 * 1 + 1 * 0) / n) # HLN correction for h=1: sqrt((n+1-2h+h(h-1)/n)/n)79 p = 2 * stats.t.sf(abs(stat), df=n - 1)80 return float(stat), float(p)8182rn, an = metrics("naive"); rm, am = metrics("mean"); rr, ar = metrics("ridge")8384# Exclusion ledger: raw rows -> training rows at the last origin85last = origins[-1]86n1 = int((d.index >= TRAIN_START).sum())87n2 = int(((d.index >= TRAIN_START) & (d.index < last)).sum())88assert n2 == n_train_last89excl = pd.DataFrame([90 dict(step=1, rule="drop months before 1992-03-01 (lags ip_growth_{s-2} and neworder_chg undefined/not in spec)",91 n_before=n_raw, n_after=n1, n_removed=n_raw - n1, rationale="TASK.md training window starts at s=1992-03-01"),92 dict(step=2, rule="keep feature months s with s+1 <= last origin 2026-08-01 (s <= 2026-07-01)",93 n_before=n1, n_after=n2, n_removed=n1 - n2, rationale="target y_{s+1} must be dated <= origin; later months unavailable"),94])95excl.to_csv(OUT / "exclusions.csv", index=False)9697res = {98 "task_id": "F3-conf-03",99 "n_origins": int(len(origins)),100 "first_origin": origins[0].strftime("%Y-%m-%d"),101 "last_origin": origins[-1].strftime("%Y-%m-%d"),102 "lambda": LAM,103 "normalisation_window": "expanding",104 "n_features": len(FEATS),105 "rmse_naive": rn, "mae_naive": an,106 "rmse_mean": rm, "mae_mean": am,107 "rmse_ridge": rr, "mae_ridge": ar,108 "rel_improvement_vs_naive": 1 - rr / rn,109 "rel_improvement_vs_mean": 1 - rr / rm,110 "feasibility": "complete",111 "feasibility_reason": None,112 "data_issues": [113 {"issue": "neworder_chg missing for 1992-01 and 1992-02 (NEWORDER starts 1992-02)",114 "handling": "Not used: training rows start at s=1992-03 per TASK.md; no imputation needed."},115 {"issue": "data/README.md states the first month (1992-01) of each change series is empty, but ip_growth, oil_chg and permit_chg are populated for 1992-01",116 "handling": "Values used as given; only ip_growth_1992-01 enters (as lag ip_growth_{s-2} for s=1992-03). oil_chg/permit_chg 1992-01 are not used."},117 {"issue": "Extreme ip_growth values in training windows (e.g. 2020-04 = -14.1, 2020-06 = +6.3, 2008-09 = -4.5)",118 "handling": "Kept: the ridge specification is fixed by TASK.md; no outlier treatment. Noted as a limitation."},119 {"issue": "Current-vintage (revised) data with publication lags ignored",120 "handling": "Errors computed on data as given; no real-time validity claimed."},121 ],122 "n_input_rows": {"ip_predictors.csv": int(n_raw)},123 "n_analysed": int(n_train_last),124 "random_seed": None,125 "dm_stat_naive_vs_ridge": dm_test("naive", "ridge")[0],126 "dm_pvalue_naive_vs_ridge": dm_test("naive", "ridge")[1],127 "dm_stat_mean_vs_ridge": dm_test("mean", "ridge")[0],128 "dm_pvalue_mean_vs_ridge": dm_test("mean", "ridge")[1],129 "sha256_input": hashlib.sha256(DATA.read_bytes()).hexdigest(),130 "baseline_best_rmse_model": "naive" if rn < rm else "mean",131}132with open(OUT / "results.json", "w") as f:133 json.dump(res, f, indent=2)134print(json.dumps({k: res[k] for k in ["rmse_naive", "rmse_mean", "rmse_ridge", "mae_naive", "mae_mean", "mae_ridge",135 "rel_improvement_vs_naive", "rel_improvement_vs_mean", "n_analysed"]}, indent=2))136