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 np, pandas as pd56LAM = 50.07FEATS = ["ip0", "ip1", "ip2", "oil_chg", "fedfunds", "permit_chg", "neworder_chg", "tcu"]8os.makedirs("out", exist_ok=True)9raw = pd.read_csv("data/ip_predictors.csv")10n_rows = len(raw)11d = raw.copy()12d["month"] = pd.to_datetime(d["month"])13d = d.set_index("month").sort_index()1415issues = []16# checks17assert d.index.is_unique and d.index.is_monotonic_increasing18full = pd.date_range(d.index[0], d.index[-1], freq="MS")19assert len(full) == len(d), "gaps in monthly index"20na = d.isna().sum()21issues.append({"issue": "Missing values: " + ", ".join(f"{k}={int(v)}" for k, v in na[na>0].items()) + " (neworder_chg for 1992-01/02).",22 "handling": "Not used: earliest ridge training row is 1992-03, whose features (incl. lags to 1992-01 of ip_growth) are complete."})23if pd.notna(d["ip_growth"].iloc[0]) or pd.notna(d["oil_chg"].iloc[0]):24 issues.append({"issue": "README says the first month (1992-01) of each change series is empty, but ip_growth, oil_chg and permit_chg have values for 1992-01 (presumably computed from 1991-12 levels).",25 "handling": "Values used as given; ip_growth 1992-01 enters only as the second lag of the 1992-03 training row; the expanding-mean baseline starts at 1992-04 as specified."})26mn = d["ip_growth"].idxmin()27issues.append({"issue": f"Extreme COVID-19 observations in ip_growth (min {d['ip_growth'].min():.3f} at {mn:%Y-%m}, max {d['ip_growth'].max():.3f} at {d['ip_growth'].idxmax():%Y-%m}) lie in every training window.",28 "handling": "Kept as given (task specifies the estimator exactly); noted as a limitation."})2930X = pd.DataFrame(index=d.index)31X["ip0"] = d["ip_growth"]; X["ip1"] = d["ip_growth"].shift(1); X["ip2"] = d["ip_growth"].shift(2)32for c in ["oil_chg", "fedfunds", "permit_chg", "neworder_chg", "tcu"]:33 X[c] = d[c]34y_next = d["ip_growth"].shift(-1) # y_{s+1} aligned to s3536origins = pd.date_range("2021-09-01", "2026-08-01", freq="MS")37rows = []38n_train_last = None39for t in origins:40 tgt = t + pd.offsets.MonthBegin(1)41 actual = float(d.loc[tgt, "ip_growth"])42 # baselines43 naive = float(d.loc[t, "ip_growth"])44 mean_ = float(d.loc[pd.Timestamp("1992-04-01"):t, "ip_growth"].mean())45 # ridge training rows: 1992-03 <= s <= t-1 (so s+1 <= t)46 s_end = t - pd.offsets.MonthBegin(1)47 Ztr = X.loc[pd.Timestamp("1992-03-01"):s_end, FEATS].to_numpy(float)48 ytr = y_next.loc[pd.Timestamp("1992-03-01"):s_end].to_numpy(float)49 assert not np.isnan(Ztr).any() and not np.isnan(ytr).any()50 mu = Ztr.mean(0); sd = Ztr.std(0, ddof=1); ybar = ytr.mean()51 Z = (Ztr - mu) / sd52 b = np.linalg.solve(Z.T @ Z + LAM * np.eye(Z.shape[1]), Z.T @ (ytr - ybar))53 zt = X.loc[t, FEATS].to_numpy(float)54 ridge = float(ybar + ((zt - mu) / sd) @ b)55 n_train_last = len(ytr)56 for m, f in [("naive", naive), ("mean", mean_), ("ridge", ridge)]:57 rows.append({"origin_date": t.strftime("%Y-%m-%d"), "target_date": tgt.strftime("%Y-%m-%d"),58 "horizon": 1, "forecast": f, "actual": actual, "model": m})59fc = pd.DataFrame(rows)60fc.to_csv("out/forecasts.csv", index=False)6162met = {}63for m in ["naive", "mean", "ridge"]:64 e = fc[fc.model == m]; err = e.forecast - e.actual65 met[f"rmse_{m}"] = float(np.sqrt((err**2).mean())); met[f"mae_{m}"] = float(err.abs().mean())6667res = {"task_id": "F3-conf-03", "n_origins": len(origins),68 "first_origin": origins[0].strftime("%Y-%m-%d"), "last_origin": origins[-1].strftime("%Y-%m-%d"),69 "lambda": LAM, "normalisation_window": "expanding", "n_features": len(FEATS), **met,70 "rel_improvement_vs_naive": 1 - met["rmse_ridge"] / met["rmse_naive"],71 "rel_improvement_vs_mean": 1 - met["rmse_ridge"] / met["rmse_mean"],72 "feasibility": "complete", "feasibility_reason": None, "data_issues": issues,73 "n_input_rows": {"ip_predictors.csv": n_rows}, "n_analysed": n_train_last, "random_seed": 20261030} # used only for the supplementary bootstrap (out/uncertainty.json)74with open("out/results.json", "w") as f:75 json.dump(res, f, indent=2)76print(json.dumps(res, indent=2))7778# Supplementary uncertainty: Diebold-Mariano (squared loss, h=1, Newey-West lag 0..3 shown with79# Harvey-Leybourne-Newbold small-sample correction) and a moving-block bootstrap of RMSE ratios.80from scipy import stats81unc = {}82e = {m: (fc[fc.model == m].forecast - fc[fc.model == m].actual).to_numpy() for m in ["naive", "mean", "ridge"]}83T = len(origins)84for base in ["naive", "mean"]:85 dl = e["ridge"]**2 - e[base]**286 dbar = dl.mean()87 out = {}88 for L in [0, 3]:89 g = [np.mean((dl[k:] - dbar) * (dl[:T-k] - dbar)) for k in range(L + 1)]90 v = g[0] + 2 * sum((1 - k / (L + 1)) * g[k] for k in range(1, L + 1))91 dm = dbar / np.sqrt(v / T) * np.sqrt((T + 1 - 2 + 0) / T) # HLN factor for h=192 out[f"nw_lag{L}"] = {"dm_stat": float(dm), "p_two_sided": float(2 * stats.t.sf(abs(dm), T - 1))}93 rng = np.random.default_rng(20261030)94 B, bl = 5000, 695 ratios = []96 for _ in range(B):97 starts = rng.integers(0, T - bl + 1, size=int(np.ceil(T / bl)))98 idx = np.concatenate([np.arange(s0, s0 + bl) for s0 in starts])[:T]99 ratios.append(1 - np.sqrt(np.mean(e["ridge"][idx]**2)) / np.sqrt(np.mean(e[base][idx]**2)))100 out["mean_loss_diff_ridge_minus_base"] = float(dbar)101 out["block_bootstrap_rel_improvement_90ci"] = [float(np.percentile(ratios, 5)), float(np.percentile(ratios, 95))]102 out["bootstrap"] = {"B": B, "block_length": bl, "seed": 20261030}103 unc[f"ridge_vs_{base}"] = out104with open("out/uncertainty.json", "w") as f:105 json.dump(unc, f, indent=2)106print(json.dumps(unc, indent=2))107