1"""F3-conf-03: rolling-origin ridge vs naive / expanding-mean forecasts of US IP growth (h=1).2Reads data/ip_predictors.csv only; writes out/results.json, out/forecasts.csv, out/exclusions.csv.3Deterministic; no randomness used (random_seed = null)."""4import json, os5import numpy as np6import pandas as pd78ROOT = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))9DATA = os.path.join(ROOT, "data", "ip_predictors.csv")10OUT = os.path.join(ROOT, "out")11LAM = 50.012TRAIN_START = pd.Timestamp("1992-03-01")13MEAN_START = pd.Timestamp("1992-04-01")14ORIGINS = pd.date_range("2021-09-01", "2026-08-01", freq="MS")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").reset_index(drop=True)20assert d["month"].is_unique21assert (d["month"].diff().dropna().dt.days.between(28, 31)).all(), "gap in monthly index"22d["ip_lag1"] = d["ip_growth"].shift(1)23d["ip_lag2"] = d["ip_growth"].shift(2)24d["y_next"] = d["ip_growth"].shift(-1)25d["month_next"] = d["month"].shift(-1)26d = d.set_index("month", drop=False)2728rows, n_train_last = [], None29for t in ORIGINS:30 tgt = t + pd.DateOffset(months=1)31 actual = float(d.loc[tgt, "ip_growth"])32 # naive33 f_naive = float(d.loc[t, "ip_growth"])34 # expanding mean 1992-04..t35 f_mean = float(d.loc[(d.month >= MEAN_START) & (d.month <= t), "ip_growth"].mean())36 # ridge: training rows s with TRAIN_START <= s and s+1 <= t37 tr = d[(d.month >= TRAIN_START) & (d.month_next <= t)]38 assert tr["month_next"].max() <= t39 X = tr[FEATS].to_numpy(float); y = tr["y_next"].to_numpy(float)40 assert np.isfinite(X).all() and np.isfinite(y).all()41 mu = X.mean(axis=0); sd = X.std(axis=0, ddof=1); ybar = y.mean()42 Z = (X - mu) / sd43 b = np.linalg.solve(Z.T @ Z + LAM * np.eye(Z.shape[1]), Z.T @ (y - ybar))44 zt = (d.loc[t, FEATS].to_numpy(float) - mu) / sd45 f_ridge = float(ybar + zt @ b)46 n_train_last = len(tr)47 for m, f in (("naive", f_naive), ("mean", f_mean), ("ridge", f_ridge)):48 rows.append(dict(origin_date=t.strftime("%Y-%m-%d"), target_date=tgt.strftime("%Y-%m-%d"),49 horizon=1, forecast=f, actual=actual, model=m))5051fc = pd.DataFrame(rows)52fc.to_csv(os.path.join(OUT, "forecasts.csv"), index=False, float_format="%.10f")5354def err(m):55 e = fc[fc.model == m]; r = e["forecast"] - e["actual"]56 return float(np.sqrt(np.mean(r ** 2))), float(np.mean(np.abs(r)))57rn, an = err("naive"); rm, am = err("mean"); rr, ar = err("ridge")5859# Diebold-Mariano (squared-error loss, h=1, Harvey-Leybourne-Newbold small-sample correction, t_{n-1})60from scipy import stats61def dm(m1, m2):62 e1 = fc[fc.model == m1].set_index("origin_date"); e2 = fc[fc.model == m2].set_index("origin_date")63 dl = (e1.forecast - e1.actual) ** 2 - (e2.forecast - e2.actual) ** 264 n = len(dl); stat = dl.mean() / np.sqrt(dl.var(ddof=0) / n)65 stat *= np.sqrt((n - 1) / n) # HLN correction factor for h=166 return float(stat), float(2 * stats.t.sf(abs(stat), n - 1))67dm_rn = dm("ridge", "naive"); dm_rm = dm("ridge", "mean")6869# exclusion ledger (training sample at last origin)70n1 = int((d.month >= TRAIN_START).sum())71excl = pd.DataFrame([72 dict(step=1, rule="drop months s < 1992-03-01 (lags ip_growth_{s-2} and neworder_chg undefined before)",73 n_before=n_raw, n_after=n1, n_removed=n_raw - n1, rationale="task-specified training start"),74 dict(step=2, rule="drop months s with s+1 > last origin 2026-08-01 (target not available at origin)",75 n_before=n1, n_after=n_train_last, n_removed=n1 - n_train_last, rationale="information-availability rule"),76])77excl.to_csv(os.path.join(OUT, "exclusions.csv"), index=False)7879res = dict(80 task_id="F3-conf-03", n_origins=int(len(ORIGINS)),81 first_origin=ORIGINS[0].strftime("%Y-%m-%d"), last_origin=ORIGINS[-1].strftime("%Y-%m-%d"),82 **{"lambda": LAM}, normalisation_window="expanding", n_features=len(FEATS),83 rmse_naive=rn, mae_naive=an, rmse_mean=rm, mae_mean=am, rmse_ridge=rr, mae_ridge=ar,84 rel_improvement_vs_naive=1 - rr / rn, rel_improvement_vs_mean=1 - rr / rm,85 feasibility="complete", feasibility_reason=None,86 data_issues=[87 {"issue": "neworder_chg missing for 1992-01 and 1992-02 (NEWORDER starts 1992-02)",88 "handling": "no effect: training rows start at 1992-03 as specified, so these cells are never used"},89 {"issue": "data/README.md states the first month of every change series is empty, but ip_growth, oil_chg and permit_chg are populated for 1992-01",90 "handling": "values used as given (only ip_growth 1992-01/02 enters, as lags of the 1992-03 training row, as the task specification implies)"},91 {"issue": "extreme COVID-19 observations (e.g. ip_growth 2020-04 = -14.1, 2020-06 = +6.3) inside every training window",92 "handling": "retained; the ridge specification is fixed by the task and no outlier treatment is allowed for"},93 {"issue": "current-vintage (revised) data with publication lags ignored",94 "handling": "results describe pseudo-out-of-sample accuracy on revised data, not real-time performance"},95 ],96 n_input_rows={"ip_predictors.csv": int(n_raw)}, n_analysed=int(n_train_last), random_seed=None,97 dm_stat_ridge_vs_naive=dm_rn[0], dm_pvalue_ridge_vs_naive=dm_rn[1],98 dm_stat_ridge_vs_mean=dm_rm[0], dm_pvalue_ridge_vs_mean=dm_rm[1],99)100with open(os.path.join(OUT, "results.json"), "w") as f:101 json.dump(res, f, indent=2)102print(json.dumps({k: res[k] for k in ["rmse_naive", "mae_naive", "rmse_mean", "mae_mean", "rmse_ridge",103 "mae_ridge", "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))104