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. Deterministic (no randomness)."""4import json, os5import numpy as np6import pandas as pd78DATA = "data/ip_predictors.csv"9OUT = "out"10LAMBDA = 50.011TRAIN_START = pd.Timestamp("1992-03-01")12MEAN_START = pd.Timestamp("1992-04-01")13FIRST_ORIGIN = pd.Timestamp("2021-09-01")14LAST_ORIGIN = pd.Timestamp("2026-08-01")15FEATS = ["ip_growth", "ip_lag1", "ip_lag2", "oil_chg", "fedfunds", "permit_chg", "neworder_chg", "tcu"]1617raw = pd.read_csv(DATA, parse_dates=["month"])18n_raw = len(raw)19d = raw.sort_values("month").set_index("month")20assert d.index.is_unique21# lags use only past values (shift forward in time)22d["ip_lag1"] = d["ip_growth"].shift(1)23d["ip_lag2"] = d["ip_growth"].shift(2)24d["y_next"] = d["ip_growth"].shift(-1) # target of row s = ip_growth at s+12526origins = pd.date_range(FIRST_ORIGIN, LAST_ORIGIN, freq="MS")27rows = []28last_ntrain = None29for t in origins:30 tgt = t + pd.DateOffset(months=1)31 actual = float(d.loc[tgt, "ip_growth"])32 avail = d.loc[:t] # information set: dated <= t33 # baselines34 f_naive = float(avail.loc[t, "ip_growth"])35 f_mean = float(avail.loc[MEAN_START:t, "ip_growth"].mean())36 # ridge training rows: TRAIN_START <= s, s+1 <= t -> s <= t - 1 month37 s_end = t - pd.DateOffset(months=1)38 tr = avail.loc[TRAIN_START:s_end]39 X = tr[FEATS].to_numpy(float)40 # target y_{s+1}: all s+1 <= t, so values are inside the information set41 y = avail["ip_growth"].shift(-1).loc[TRAIN_START:s_end].to_numpy(float)42 assert np.isfinite(X).all() and np.isfinite(y).all()43 mu = X.mean(axis=0)44 sd = X.std(axis=0, ddof=1)45 Z = (X - mu) / sd46 ybar = y.mean()47 b = np.linalg.solve(Z.T @ Z + LAMBDA * np.eye(Z.shape[1]), Z.T @ (y - ybar))48 zt = (avail.loc[t, FEATS].to_numpy(float) - mu) / sd49 f_ridge = float(ybar + zt @ b)50 last_ntrain = len(tr)51 for m, f in (("naive", f_naive), ("mean", f_mean), ("ridge", f_ridge)):52 rows.append(dict(origin_date=t.strftime("%Y-%m-%d"), target_date=tgt.strftime("%Y-%m-%d"),53 horizon=1, forecast=f, actual=actual, model=m))5455fc = pd.DataFrame(rows)56os.makedirs(OUT, exist_ok=True)57fc.to_csv(os.path.join(OUT, "forecasts.csv"), index=False, float_format="%.10f")5859def metrics(m):60 e = fc.loc[fc.model == m, "actual"].to_numpy() - fc.loc[fc.model == m, "forecast"].to_numpy()61 return float(np.sqrt(np.mean(e ** 2))), float(np.mean(np.abs(e)))6263rmse_n, mae_n = metrics("naive")64rmse_m, mae_m = metrics("mean")65rmse_r, mae_r = metrics("ridge")6667# exclusion ledger: raw rows -> training rows at the last origin68s_end_last = LAST_ORIGIN - pd.DateOffset(months=1)69n1 = int((d.index >= TRAIN_START).sum())70n2 = int(((d.index >= TRAIN_START) & (d.index <= s_end_last)).sum())71excl = pd.DataFrame([72 dict(step=1, rule="month >= 1992-03-01 (lag-2 of ip_growth and neworder_chg required)", n_before=n_raw, n_after=n1,73 n_removed=n_raw - n1, rationale="1992-01/02 lack neworder_chg and/or two ip_growth lags; task fixes training start at 1992-03"),74 dict(step=2, rule="s + 1 <= last origin 2026-08-01 (s <= 2026-07-01)", n_before=n1, n_after=n2,75 n_removed=n1 - n2, rationale="information rule: targets dated after the origin may not be used in training (2026-08 row's target 2026-09 is the evaluation target; 2026-09 has no target)"),76])77excl.to_csv(os.path.join(OUT, "exclusions.csv"), index=False)78assert n2 == last_ntrain7980# Uncertainty: Diebold-Mariano test (squared-error loss, h=1, Harvey-Leybourne-Newbold small-sample correction)81from scipy import stats82def dm(m1, m2):83 e1 = fc.loc[fc.model == m1, "actual"].to_numpy() - fc.loc[fc.model == m1, "forecast"].to_numpy()84 e2 = fc.loc[fc.model == m2, "actual"].to_numpy() - fc.loc[fc.model == m2, "forecast"].to_numpy()85 dl = e1 ** 2 - e2 ** 286 T = len(dl)87 stat = dl.mean() / np.sqrt(dl.var(ddof=0) / T) # h=1: no autocovariance terms88 stat *= np.sqrt((T + 1 - 2 + 0) / T) # HLN factor for h=1: sqrt((T+1-2h+h(h-1)/T)/T)89 p = 2 * stats.t.sf(abs(stat), df=T - 1)90 return float(stat), float(p)91dm_rn, p_rn = dm("ridge", "naive")92dm_rm, p_rm = dm("ridge", "mean")9394res = {95 "task_id": "F3-conf-03",96 "n_origins": int(len(origins)),97 "first_origin": origins[0].strftime("%Y-%m-%d"),98 "last_origin": origins[-1].strftime("%Y-%m-%d"),99 "lambda": LAMBDA,100 "normalisation_window": "expanding",101 "n_features": len(FEATS),102 "rmse_naive": rmse_n, "mae_naive": mae_n,103 "rmse_mean": rmse_m, "mae_mean": mae_m,104 "rmse_ridge": rmse_r, "mae_ridge": mae_r,105 "rel_improvement_vs_naive": 1 - rmse_r / rmse_n,106 "rel_improvement_vs_mean": 1 - rmse_r / rmse_m,107 "feasibility": "complete",108 "feasibility_reason": None,109 "data_issues": [110 {"issue": "neworder_chg missing for 1992-01 and 1992-02",111 "handling": "these months are never ridge feature rows (training starts 1992-03); no imputation needed"},112 {"issue": "README states 1992-01 is empty for every change series, but ip_growth, oil_chg and permit_chg have values in 1992-01",113 "handling": "values used as given; ip_growth 1992-01 serves as lag-2 feature of the 1992-03 training row"},114 {"issue": "extreme COVID-2020 and 2008 observations (e.g. ip_growth 2020-04 = -14.1, oil_chg 2020-03..05 about +/-55)",115 "handling": "kept unchanged as genuine data; they enter training windows only (not the evaluation period)"},116 {"issue": "current-vintage data with publication lags ignored",117 "handling": "results describe revised-data predictability under stylised timing, not real-time performance"},118 ],119 "n_input_rows": {"ip_predictors.csv": int(n_raw)},120 "n_analysed": int(last_ntrain),121 "random_seed": None,122 "dm_stat_ridge_vs_naive": dm_rn, "dm_pvalue_ridge_vs_naive": p_rn,123 "dm_stat_ridge_vs_mean": dm_rm, "dm_pvalue_ridge_vs_mean": p_rm,124}125with open(os.path.join(OUT, "results.json"), "w") as f:126 json.dump(res, f, indent=2)127print(json.dumps({k: res[k] for k in ["rmse_naive", "mae_naive", "rmse_mean", "mae_mean", "rmse_ridge", "mae_ridge",128 "rel_improvement_vs_naive", "rel_improvement_vs_mean", "n_analysed", "dm_stat_ridge_vs_naive", "dm_pvalue_ridge_vs_naive", "dm_stat_ridge_vs_mean", "dm_pvalue_ridge_vs_mean"]}, indent=2))129