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, 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")1617raw = pd.read_csv(DATA, parse_dates=["month"])18n_raw = len(raw)19d = raw.sort_values("month").set_index("month")20assert not d.index.duplicated().any()21assert (d.index == pd.date_range(d.index[0], d.index[-1], freq="MS")).all(), "non-contiguous months"2223# feature matrix z_s24F = pd.DataFrame({25 "ip0": d["ip_growth"], "ip1": d["ip_growth"].shift(1), "ip2": d["ip_growth"].shift(2),26 "oil": d["oil_chg"], "ff": d["fedfunds"], "permit": d["permit_chg"],27 "neworder": d["neworder_chg"], "tcu": d["tcu"]})28y_next = d["ip_growth"].shift(-1) # target for row s is ip_{s+1}29months = d.index3031origins = pd.date_range(FIRST_ORIGIN, LAST_ORIGIN, freq="MS")32rows = []33n_train_last = None34for t in origins:35 tgt = t + pd.offsets.MonthBegin(1)36 actual = float(d.loc[tgt, "ip_growth"])37 # naive38 rows.append((t, tgt, d.loc[t, "ip_growth"], actual, "naive"))39 # expanding mean 1992-04..t40 rows.append((t, tgt, d.loc[MEAN_START:t, "ip_growth"].mean(), actual, "mean"))41 # ridge: training rows s in [1992-03, t-1] (s+1 <= t)42 s_last = t - pd.offsets.MonthBegin(1)43 Xtr = F.loc[TRAIN_START:s_last].to_numpy(float)44 ytr = y_next.loc[TRAIN_START:s_last].to_numpy(float)45 assert np.isfinite(Xtr).all() and np.isfinite(ytr).all()46 assert y_next.loc[TRAIN_START:s_last].index.max() + pd.offsets.MonthBegin(1) <= t47 mu = Xtr.mean(0); sd = Xtr.std(0, ddof=1); ybar = ytr.mean()48 Z = (Xtr - mu) / sd49 b = np.linalg.solve(Z.T @ Z + LAM * np.eye(Z.shape[1]), Z.T @ (ytr - ybar))50 zt = (F.loc[t].to_numpy(float) - mu) / sd51 rows.append((t, tgt, float(ybar + zt @ b), actual, "ridge"))52 n_train_last = len(ytr)5354fc = pd.DataFrame(rows, columns=["origin_date", "target_date", "forecast", "actual", "model"])55fc.insert(2, "horizon", 1)56fc["origin_date"] = fc["origin_date"].dt.strftime("%Y-%m-%d")57fc["target_date"] = fc["target_date"].dt.strftime("%Y-%m-%d")58fc.to_csv(OUT / "forecasts.csv", index=False, float_format="%.8f")5960res = {"task_id": "F3-conf-03", "n_origins": len(origins),61 "first_origin": origins[0].strftime("%Y-%m-%d"), "last_origin": origins[-1].strftime("%Y-%m-%d"),62 "lambda": LAM, "normalisation_window": "expanding", "n_features": F.shape[1]}63err = {}64for m in ["naive", "mean", "ridge"]:65 e = (fc.loc[fc.model == m, "actual"] - fc.loc[fc.model == m, "forecast"]).to_numpy()66 err[m] = e67 res[f"rmse_{m}"] = float(np.sqrt(np.mean(e ** 2)))68 res[f"mae_{m}"] = float(np.mean(np.abs(e)))69res["rel_improvement_vs_naive"] = 1 - res["rmse_ridge"] / res["rmse_naive"]70res["rel_improvement_vs_mean"] = 1 - res["rmse_ridge"] / res["rmse_mean"]71res["feasibility"] = "complete"72res["feasibility_reason"] = None73res["data_issues"] = [74 {"issue": "neworder_chg missing for 1992-01 and 1992-02 (NEWORDER starts 1992-02).",75 "handling": "Outside every training window (training starts 1992-03 as specified); no imputation needed."},76 {"issue": "data/README.md says the first month (1992-01) of each change series is empty, but ip_growth, oil_chg and permit_chg have values in 1992-01.",77 "handling": "Values used as given; only ip_growth 1992-01 enters (as lag ip_{s-2} for s=1992-03), as implied by the task's feature definition."},78 {"issue": "Current-vintage (revised) data with publication lags ignored; training windows include the 2020 COVID swing (ip_growth minimum -14.14 in 2020-04) while evaluation-period actuals lie in [-1.38, 1.33].",79 "handling": "No adjustment or outlier treatment (not part of the specified model); results are not real-time and are interpreted accordingly."},80]81res["n_input_rows"] = {"ip_predictors.csv": n_raw}82res["n_analysed"] = n_train_last83res["random_seed"] = None84# supplementary (not required): Diebold-Mariano-type tests on squared-error loss, h=1, plain variance85from scipy import stats86supp = {}87for base in ["naive", "mean"]:88 dl = err[base] ** 2 - err["ridge"] ** 289 dm = dl.mean() / np.sqrt(dl.var(ddof=1) / len(dl))90 supp[f"dm_stat_ridge_vs_{base}"] = float(dm)91 supp[f"dm_pvalue_ridge_vs_{base}"] = float(2 * (1 - stats.t.cdf(abs(dm), df=len(dl) - 1)))92res["supplementary"] = supp93(OUT / "results.json").write_text(json.dumps(res, indent=2))9495# exclusion ledger: raw rows -> rows used as ridge training rows at last origin96last = origins[-1]97steps = [("raw", "all rows read", n_raw, n_raw)]98n1 = int((months >= TRAIN_START).sum())99steps.append(("1", "drop months before 1992-03 (no lags / neworder_chg missing)", n_raw, n1))100n2 = int(((months >= TRAIN_START) & (months <= last - pd.offsets.MonthBegin(1))).sum())101steps.append(("2", "drop months s with s+1 > last origin 2026-08 (target not available at last origin)", n1, n2))102ex = pd.DataFrame([(s, r, a, b, a - b) for s, r, a, b in steps],103 columns=["step", "rule", "n_before", "n_after", "n_removed"])104ex["rationale"] = ["input", "task-specified training start; features need ip lags and neworder_chg",105 "information rule: s+1 <= t at the last origin"]106ex.to_csv(OUT / "exclusions.csv", index=False)107assert n2 == n_train_last108print(json.dumps({k: v for k, v in res.items() if k.startswith(("rmse", "mae", "rel"))}, indent=2))109