1"""F3-conf-03: rolling-origin ridge vs naive / expanding-mean forecasts of US IP growth (h=1).2Run from workspace root: ~/wp29/venv/bin/python out/analysis.py3Reads only data/; writes out/results.json, out/forecasts.csv, out/exclusions.csv. No randomness."""4import json, os5import numpy as np6import pandas as pd78DATA = "data/ip_predictors.csv"9OUT = "out"10LAM = 50.011FIRST, LAST = pd.Timestamp("2021-09-01"), pd.Timestamp("2026-08-01")12TRAIN_START = pd.Timestamp("1992-03-01")13MEAN_START = pd.Timestamp("1992-04-01")14FEATS = ["ip_growth", "ip_lag1", "ip_lag2", "oil_chg", "fedfunds", "permit_chg", "neworder_chg", "tcu"]1516raw = pd.read_csv(DATA, parse_dates=["month"])17n_raw = len(raw)18d = raw.sort_values("month").reset_index(drop=True)19assert d["month"].is_unique20assert (d["month"].values == pd.date_range(d.month.min(), d.month.max(), freq="MS").values).all()21d["ip_lag1"] = d["ip_growth"].shift(1)22d["ip_lag2"] = d["ip_growth"].shift(2)23d["y_next"] = d["ip_growth"].shift(-1) # target y_{s+1}, aligned to row s24d = d.set_index("month")2526origins = pd.date_range(FIRST, LAST, freq="MS")27rows, ntrain = [], {}28for t in origins:29 tgt = t + pd.DateOffset(months=1)30 actual = float(d.loc[tgt, "ip_growth"])31 # naive32 f_naive = float(d.loc[t, "ip_growth"])33 # expanding mean 1992-04..t34 f_mean = float(d.loc[MEAN_START:t, "ip_growth"].mean())35 # ridge: training rows s in [1992-03, t-1]; y_{s+1} dated <= t36 tr = d.loc[TRAIN_START:t - pd.DateOffset(months=1)]37 X = tr[FEATS].to_numpy(float); y = tr["y_next"].to_numpy(float)38 assert np.isfinite(X).all() and np.isfinite(y).all()39 assert tr.index.max() + pd.DateOffset(months=1) <= t40 mu = X.mean(0); sd = X.std(0, ddof=1); ybar = y.mean()41 Z = (X - mu) / sd42 b = np.linalg.solve(Z.T @ Z + LAM * np.eye(len(FEATS)), Z.T @ (y - ybar))43 zt = (d.loc[t, FEATS].to_numpy(float) - mu) / sd44 f_ridge = float(ybar + zt @ b)45 ntrain[t] = len(tr)46 for m, f in (("naive", f_naive), ("mean", f_mean), ("ridge", f_ridge)):47 rows.append(dict(origin_date=t.strftime("%Y-%m-%d"), target_date=tgt.strftime("%Y-%m-%d"),48 horizon=1, forecast=f, actual=actual, model=m))4950fc = pd.DataFrame(rows)51fc.to_csv(os.path.join(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(np.mean(e ** 2))), float(np.mean(np.abs(e)))5657r = {m: metrics(m) for m in ("naive", "mean", "ridge")}58n_last = ntrain[origins[-1]]5960# Diebold-Mariano test (squared loss, h=1, HLN small-sample correction, t_{n-1} reference)61from scipy import stats62def dm(base):63 e_b = (fc.loc[fc.model == base, "forecast"].values - fc.loc[fc.model == base, "actual"].values)64 e_r = (fc.loc[fc.model == "ridge", "forecast"].values - fc.loc[fc.model == "ridge", "actual"].values)65 dd = e_b ** 2 - e_r ** 2 # >0 means ridge better66 n = len(dd); var = dd.var(ddof=0) / n67 stat = dd.mean() / np.sqrt(var) * np.sqrt((n - 1) / n) # HLN factor for h=168 p = 2 * stats.t.sf(abs(stat), df=n - 1)69 return float(stat), float(p)70dm_naive, dm_mean = dm("naive"), dm("mean")7172# exclusion ledger: raw rows -> rows usable as ridge training rows at last origin73n_after_start = int((d.index >= TRAIN_START).sum())74n_after_end = n_last75excl = pd.DataFrame([76 dict(step=1, rule="drop months before 1992-03-01 (missing neworder_chg / lags unavailable; spec training start)",77 n_before=n_raw, n_after=n_after_start, n_removed=n_raw - n_after_start,78 rationale="TASK.md fixes training start s >= 1992-03-01"),79 dict(step=2, rule="drop months s with s+1 > last origin 2026-08-01 (target not available at last origin)",80 n_before=n_after_start, n_after=n_after_end, n_removed=n_after_start - n_after_end,81 rationale="information-availability rule: y_{s+1} must be dated <= origin"),82])83excl.to_csv(os.path.join(OUT, "exclusions.csv"), index=False)8485res = {86 "task_id": "F3-conf-03",87 "n_origins": int(len(origins)),88 "first_origin": origins[0].strftime("%Y-%m-%d"),89 "last_origin": origins[-1].strftime("%Y-%m-%d"),90 "lambda": LAM,91 "normalisation_window": "expanding",92 "n_features": len(FEATS),93 "rmse_naive": r["naive"][0], "mae_naive": r["naive"][1],94 "rmse_mean": r["mean"][0], "mae_mean": r["mean"][1],95 "rmse_ridge": r["ridge"][0], "mae_ridge": r["ridge"][1],96 "rel_improvement_vs_naive": 1 - r["ridge"][0] / r["naive"][0],97 "rel_improvement_vs_mean": 1 - r["ridge"][0] / r["mean"][0],98 "dm_stat_vs_naive": dm_naive[0], "dm_pvalue_vs_naive": dm_naive[1],99 "dm_stat_vs_mean": dm_mean[0], "dm_pvalue_vs_mean": dm_mean[1],100 "feasibility": "complete",101 "feasibility_reason": None,102 "data_issues": [103 {"issue": "neworder_chg missing for 1992-01 and 1992-02 (NEWORDER starts 1992-02)",104 "handling": "not used: ridge training rows start at 1992-03-01 as specified"},105 {"issue": "data/README.md says the first month (1992-01) of every change series is empty, but the file contains values for ip_growth, oil_chg and permit_chg in 1992-01",106 "handling": "values used as given; ip_growth 1992-01 enters only as the 2nd lag of the 1992-03 training row, as implied by the TASK.md specification"},107 {"issue": "extreme COVID-19 observations in 2020-03..2020-06 (ip_growth -14.1 to +6.3 percent) inside every training window",108 "handling": "kept, as the specification defines the training window without exclusions; noted as a limitation"},109 {"issue": "current-vintage (revised) data with publication lags ignored",110 "handling": "results reported as pseudo-out-of-sample on revised data, not real-time"},111 ],112 "n_input_rows": {"ip_predictors.csv": int(n_raw)},113 "n_analysed": int(n_last),114 "random_seed": None,115}116with open(os.path.join(OUT, "results.json"), "w") as f:117 json.dump(res, f, indent=2)118print(json.dumps({k: res[k] for k in res if k.startswith(("rmse", "mae", "rel", "n_"))}, indent=2))119