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.07d = pd.read_csv("data/ip_predictors.csv", parse_dates=["month"])8n_rows = len(d)9d = d.sort_values("month").set_index("month")10assert d.index.is_unique11months = d.index12issues = []13# data checks14gaps = pd.Series(months).diff().dt.days.dropna()15if not gaps.between(28, 31).all():16 issues.append({"issue": "non-consecutive months", "handling": "none needed check failed"})17na = d.isna().sum()18issues.append({"issue": f"missing values: neworder_chg {int(na['neworder_chg'])} (1992-01/02); other columns {int(na.drop('neworder_chg').sum())}",19 "handling": "outside every training window (training starts 1992-03, lags reach back to 1992-01 only for ip_growth); no imputation"})20if not np.isnan(d["ip_growth"].iloc[0]):21 issues.append({"issue": "README says the first month of each change series is empty, but ip_growth (and other change columns) for 1992-01 are populated",22 "handling": "values kept as given; ip_growth 1992-01 is used only as the second lag of the 1992-03 training row, as the task specification implies"})23covid = d.loc["2020-03-01":"2020-07-01", "ip_growth"]24issues.append({"issue": f"extreme COVID-19 observations in 2020 (ip_growth range {covid.min():.2f} to {covid.max():.2f})",25 "handling": "kept in all training windows as specified (no outlier treatment); they influence standardisation statistics and coefficients"})2627ip = d["ip_growth"]28feat = pd.DataFrame({29 "ip0": ip, "ip1": ip.shift(1), "ip2": ip.shift(2),30 "oil": d["oil_chg"], "ff": d["fedfunds"], "perm": d["permit_chg"],31 "nord": d["neworder_chg"], "tcu": d["tcu"]}, index=months)32target_next = ip.shift(-1) # y_{s+1} indexed by s3334origins = pd.date_range("2021-09-01", "2026-08-01", freq="MS")35rows = []36n_train_last = None37for t in origins:38 tgt = t + pd.offsets.MonthBegin(1)39 actual = float(ip.loc[tgt])40 # naive41 f_naive = float(ip.loc[t])42 # expanding mean 1992-04..t43 f_mean = float(ip.loc["1992-04-01":t].mean())44 # ridge: training s in [1992-03, t-1]45 s_idx = months[(months >= "1992-03-01") & (months <= t - pd.offsets.MonthBegin(1))]46 X = feat.loc[s_idx].to_numpy(float)47 y = target_next.loc[s_idx].to_numpy(float)48 assert np.isfinite(X).all() and np.isfinite(y).all()49 assert (s_idx + pd.offsets.MonthBegin(1)).max() <= t # no info after t50 mu = X.mean(0); sd = X.std(0, ddof=1); ybar = y.mean()51 Z = (X - mu) / sd52 b = np.linalg.solve(Z.T @ Z + LAM * np.eye(Z.shape[1]), Z.T @ (y - ybar))53 zt = (feat.loc[t].to_numpy(float) - mu) / sd54 f_ridge = float(ybar + zt @ b)55 n_train_last = len(s_idx)56 for m, f in (("naive", f_naive), ("mean", f_mean), ("ridge", f_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})5960fc = pd.DataFrame(rows)61os.makedirs("out", exist_ok=True)62fc.to_csv("out/forecasts.csv", index=False)6364def err(m):65 e = fc.loc[fc.model == m, "forecast"] - fc.loc[fc.model == m, "actual"]66 return float(np.sqrt((e ** 2).mean())), float(e.abs().mean())67r = {k: err(k) for k in ("naive", "mean", "ridge")}68res = {69 "task_id": "F3-conf-03", "n_origins": len(origins),70 "first_origin": origins[0].strftime("%Y-%m-%d"), "last_origin": origins[-1].strftime("%Y-%m-%d"),71 "lambda": LAM, "normalisation_window": "expanding", "n_features": 8,72 "rmse_naive": r["naive"][0], "mae_naive": r["naive"][1],73 "rmse_mean": r["mean"][0], "mae_mean": r["mean"][1],74 "rmse_ridge": r["ridge"][0], "mae_ridge": r["ridge"][1],75 "rel_improvement_vs_naive": 1 - r["ridge"][0] / r["naive"][0],76 "rel_improvement_vs_mean": 1 - r["ridge"][0] / r["mean"][0],77 "feasibility": "complete", "feasibility_reason": None,78 "data_issues": issues, "n_input_rows": {"ip_predictors.csv": n_rows},79 "n_analysed": n_train_last, "random_seed": 12345, # used only by the bootstrap in out/uncertainty.json80}81with open("out/results.json", "w") as f:82 json.dump(res, f, indent=2)83print(json.dumps(res, indent=2))8485# ---- Uncertainty diagnostics (written to out/uncertainty.json; not part of results.json) ----86from scipy import stats87def dm(base, model="ridge"):88 e1 = (fc.loc[fc.model == model, "forecast"].values - fc.loc[fc.model == model, "actual"].values)89 e0 = (fc.loc[fc.model == base, "forecast"].values - fc.loc[fc.model == base, "actual"].values)90 dl = e1 ** 2 - e0 ** 291 n = len(dl)92 stat = dl.mean() / np.sqrt(dl.var(ddof=0) / n)93 stat_hln = stat * np.sqrt((n + 1 - 2 + 0) / n) # Harvey-Leybourne-Newbold correction, h=194 p = 2 * stats.t.sf(abs(stat_hln), df=n - 1)95 return {"mean_loss_diff_sq": float(dl.mean()), "dm_stat_hln": float(stat_hln), "p_two_sided": float(p)}96rng = np.random.default_rng(12345)97def boot_rel(base, B=5000, block=6):98 e1 = (fc.loc[fc.model == "ridge", "forecast"].values - fc.loc[fc.model == "ridge", "actual"].values)99 e0 = (fc.loc[fc.model == base, "forecast"].values - fc.loc[fc.model == base, "actual"].values)100 n = len(e1); out = []101 for _ in range(B):102 starts = rng.integers(0, n - block + 1, size=int(np.ceil(n / block)))103 idx = np.concatenate([np.arange(s, s + block) for s in starts])[:n]104 out.append(1 - np.sqrt((e1[idx] ** 2).mean()) / np.sqrt((e0[idx] ** 2).mean()))105 return [float(np.quantile(out, 0.025)), float(np.quantile(out, 0.975))]106unc = {"dm_ridge_vs_naive": dm("naive"), "dm_ridge_vs_mean": dm("mean"),107 "boot95_rel_impr_vs_naive": boot_rel("naive"), "boot95_rel_impr_vs_mean": boot_rel("mean"),108 "bootstrap": "moving block, block=6, B=5000, seed=12345"}109with open("out/uncertainty.json", "w") as f:110 json.dump(unc, f, indent=2)111print(json.dumps(unc, indent=2))112