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 data/ only; writes out/results.json, out/forecasts.csv, out/exclusions.csv. Deterministic (no randomness)."""4import json5import os6import numpy as np7import pandas as pd8from scipy import stats910DATA = "data/ip_predictors.csv"11OUT = "out"12LAMBDA = 50.013TRAIN_START = pd.Timestamp("1992-03-01")14MEAN_START = pd.Timestamp("1992-04-01")15ORIGINS = pd.date_range("2021-09-01", "2026-08-01", freq="MS")16FEATURES = ["ip_l0", "ip_l1", "ip_l2", "oil_chg", "fedfunds", "permit_chg", "neworder_chg", "tcu"]1718raw = pd.read_csv(DATA, parse_dates=["month"])19n_raw = len(raw)20df = raw.sort_values("month").set_index("month")21assert df.index.is_unique22assert (df.index == pd.date_range(df.index[0], df.index[-1], freq="MS")).all(), "non-contiguous months"2324# Feature vector z_s (all values dated <= s) and the h=1 target y_{s+1}25feat = pd.DataFrame(index=df.index)26feat["ip_l0"] = df["ip_growth"]27feat["ip_l1"] = df["ip_growth"].shift(1)28feat["ip_l2"] = df["ip_growth"].shift(2)29for c in ["oil_chg", "fedfunds", "permit_chg", "neworder_chg", "tcu"]:30 feat[c] = df[c]31target_next = df["ip_growth"].shift(-1) # y_{s+1}, used only as a training label for s+1 <= t323334def ridge_forecast(t):35 """Forecast of ip_growth(t+1) using only data dated <= t."""36 train_idx = feat.index[(feat.index >= TRAIN_START) & (feat.index + pd.DateOffset(months=1) <= t)]37 X = feat.loc[train_idx, FEATURES].to_numpy(float)38 y = target_next.loc[train_idx].to_numpy(float)39 # sanity: every label date is <= t40 assert (train_idx + pd.DateOffset(months=1)).max() <= t41 assert np.isfinite(X).all() and np.isfinite(y).all()42 mu = X.mean(axis=0)43 sd = X.std(axis=0, ddof=1)44 Z = (X - mu) / sd45 ybar = y.mean()46 b = np.linalg.solve(Z.T @ Z + LAMBDA * np.eye(Z.shape[1]), Z.T @ (y - ybar))47 zt = (feat.loc[t, FEATURES].to_numpy(float) - mu) / sd48 return float(ybar + zt @ b), len(train_idx), b495051rows = []52n_train = {}53coefs = {}54for t in ORIGINS:55 tgt = t + pd.DateOffset(months=1)56 actual = float(df.loc[tgt, "ip_growth"])57 naive = float(df.loc[t, "ip_growth"])58 mean_fc = float(df.loc[MEAN_START:t, "ip_growth"].mean())59 ridge_fc, ntr, b = ridge_forecast(t)60 n_train[t] = ntr61 coefs[t] = b62 for m, f in [("naive", naive), ("mean", mean_fc), ("ridge", ridge_fc)]:63 rows.append(dict(origin_date=t.strftime("%Y-%m-%d"), target_date=tgt.strftime("%Y-%m-%d"),64 horizon=1, forecast=f, actual=actual, model=m))6566fc = pd.DataFrame(rows)67os.makedirs(OUT, exist_ok=True)68fc.to_csv(os.path.join(OUT, "forecasts.csv"), index=False, float_format="%.10f")6970err = {m: (g["forecast"] - g["actual"]).to_numpy() for m, g in fc.groupby("model", sort=False)}71rmse = {m: float(np.sqrt(np.mean(e ** 2))) for m, e in err.items()}72mae = {m: float(np.mean(np.abs(e))) for m, e in err.items()}737475def dm_test(e_model, e_base):76 """Diebold-Mariano (h=1, squared loss) with Harvey-Leybourne-Newbold small-sample correction.77 d = e_model^2 - e_base^2; negative mean favours the model."""78 d = e_model ** 2 - e_base ** 279 n = len(d)80 stat = d.mean() / np.sqrt(d.var(ddof=0) / n)81 stat *= np.sqrt((n + 1 - 2 + 0) / n) # HLN factor for h=1: sqrt((n+1-2h+h(h-1)/n)/n)82 p = 2 * stats.t.sf(abs(stat), df=n - 1)83 return float(stat), float(p)848586dm_naive = dm_test(err["ridge"], err["naive"])87dm_mean = dm_test(err["ridge"], err["mean"])8889last = ORIGINS[-1]90n_analysed = int(n_train[last])9192# Exclusion ledger (rows of the input file -> training rows at the last origin)93n_after_start = int((feat.index >= TRAIN_START).sum())94excl = pd.DataFrame([95 dict(step=1, rule="month < 1992-03-01 (before training start; ip lags / neworder_chg incomplete)",96 n_before=n_raw, n_after=n_after_start, n_removed=n_raw - n_after_start,97 rationale="Task specifies training rows start at s=1992-03-01; neworder_chg missing for 1992-01/02"),98 dict(step=2, rule="s + 1 > last origin 2026-08-01 (target not available at last origin)",99 n_before=n_after_start, n_after=n_analysed, n_removed=n_after_start - n_analysed,100 rationale="Information rule: label y_{s+1} must be dated <= origin t"),101])102excl.to_csv(os.path.join(OUT, "exclusions.csv"), index=False)103assert n_raw - excl["n_removed"].sum() == n_analysed104105results = {106 "task_id": "F3-conf-03",107 "n_origins": int(len(ORIGINS)),108 "first_origin": ORIGINS[0].strftime("%Y-%m-%d"),109 "last_origin": last.strftime("%Y-%m-%d"),110 "lambda": LAMBDA,111 "normalisation_window": "expanding",112 "n_features": len(FEATURES),113 "rmse_naive": rmse["naive"], "mae_naive": mae["naive"],114 "rmse_mean": rmse["mean"], "mae_mean": mae["mean"],115 "rmse_ridge": rmse["ridge"], "mae_ridge": mae["ridge"],116 "rel_improvement_vs_naive": 1 - rmse["ridge"] / rmse["naive"],117 "rel_improvement_vs_mean": 1 - rmse["ridge"] / rmse["mean"],118 "feasibility": "complete",119 "feasibility_reason": None,120 "data_issues": [121 {"issue": "neworder_chg missing for 1992-01 and 1992-02",122 "handling": "Not needed: training rows start at s=1992-03 as specified; no imputation."},123 {"issue": "README states the 1992-01 value of every change series is empty, but ip_growth, oil_chg and permit_chg are populated for 1992-01",124 "handling": "Values used as given (ip_growth 1992-01 enters only as the 2nd lag of the first training row); documentation discrepancy noted."},125 {"issue": "Extreme COVID-era observations (ip_growth 2020-04 = -14.1%, 2020-06; oil_chg 2020-03..05) lie in every training window",126 "handling": "Retained: the model is fully specified by the task; no outlier treatment. Discussed as a limitation."},127 {"issue": "Current-vintage (revised) data and ignored publication lags",128 "handling": "Results interpreted as pseudo-out-of-sample on revised data, not real-time accuracy."},129 ],130 "n_input_rows": {"ip_predictors.csv": int(n_raw)},131 "n_analysed": n_analysed,132 "random_seed": None,133 # supplementary keys134 "n_train_first_origin": int(n_train[ORIGINS[0]]),135 "n_train_last_origin": n_analysed,136 "first_target": (ORIGINS[0] + pd.DateOffset(months=1)).strftime("%Y-%m-%d"),137 "last_target": (last + pd.DateOffset(months=1)).strftime("%Y-%m-%d"),138 "rel_improvement_mae_vs_naive": 1 - mae["ridge"] / mae["naive"],139 "rel_improvement_mae_vs_mean": 1 - mae["ridge"] / mae["mean"],140 "dm_stat_ridge_vs_naive": dm_naive[0], "dm_pvalue_ridge_vs_naive": dm_naive[1],141 "dm_stat_ridge_vs_mean": dm_mean[0], "dm_pvalue_ridge_vs_mean": dm_mean[1],142 "ridge_coef_last_origin": {k: float(v) for k, v in zip(FEATURES, coefs[last])},143}144with open(os.path.join(OUT, "results.json"), "w") as f:145 json.dump(results, f, indent=2)146print(json.dumps({k: results[k] for k in ["rmse_naive", "rmse_mean", "rmse_ridge", "mae_naive", "mae_mean",147 "mae_ridge", "rel_improvement_vs_naive", "rel_improvement_vs_mean",148 "dm_pvalue_ridge_vs_naive", "dm_pvalue_ridge_vs_mean", "n_analysed"]}, indent=1))149