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, os5import numpy as np6import pandas as pd78DATA = "data/ip_predictors.csv"9OUT = "out"10LAM = 50.011FIRST_TRAIN = pd.Timestamp("1992-03-01")12MEAN_START = pd.Timestamp("1992-04-01")13ORIGINS = pd.date_range("2021-09-01", "2026-08-01", freq="MS")1415raw = pd.read_csv(DATA, parse_dates=["month"])16n_raw = len(raw)17df = raw.sort_values("month").set_index("month")18assert df.index.is_unique and (df.index == pd.date_range(df.index[0], df.index[-1], freq="MS")).all()1920# Feature matrix z_s (dated s; uses only values dated <= s)21F = pd.DataFrame({22 "ip_l0": df["ip_growth"], "ip_l1": df["ip_growth"].shift(1), "ip_l2": df["ip_growth"].shift(2),23 "oil_chg": df["oil_chg"], "fedfunds": df["fedfunds"], "permit_chg": df["permit_chg"],24 "neworder_chg": df["neworder_chg"], "tcu": df["tcu"],25})26y_next = df["ip_growth"].shift(-1) # target of row s = ip_growth_{s+1}2728rows, excl = [], []29n_train_last = None30for t in ORIGINS:31 avail = df.loc[:t] # information set: values dated <= t only32 Fa = F.loc[:t]33 ya = avail["ip_growth"].shift(-1) # recomputed inside the info set34 s_idx = Fa.index[(Fa.index >= FIRST_TRAIN) & (Fa.index <= t - pd.offsets.MonthBegin(1))]35 X = Fa.loc[s_idx].to_numpy(float)36 y = ya.loc[s_idx].to_numpy(float)37 assert np.isfinite(X).all() and np.isfinite(y).all()38 mu = X.mean(0); sd = X.std(0, ddof=1); ybar = y.mean()39 Z = (X - mu) / sd40 b = np.linalg.solve(Z.T @ Z + LAM * np.eye(Z.shape[1]), Z.T @ (y - ybar))41 zt = Fa.loc[t].to_numpy(float)42 assert np.isfinite(zt).all()43 f_ridge = ybar + ((zt - mu) / sd) @ b44 f_naive = avail.loc[t, "ip_growth"]45 f_mean = avail.loc[MEAN_START:t, "ip_growth"].mean()46 tgt = t + pd.offsets.MonthBegin(1)47 actual = df.loc[tgt, "ip_growth"] # used only for scoring48 for m, f in (("naive", f_naive), ("mean", f_mean), ("ridge", f_ridge)):49 rows.append({"origin_date": t.strftime("%Y-%m-%d"), "target_date": tgt.strftime("%Y-%m-%d"),50 "horizon": 1, "forecast": float(f), "actual": float(actual), "model": m})51 n_train_last = len(s_idx)5253fc = pd.DataFrame(rows)54fc.to_csv(os.path.join(OUT, "forecasts.csv"), index=False, float_format="%.10f")5556def err(m):57 e = fc.loc[fc.model == m, "forecast"] - fc.loc[fc.model == m, "actual"]58 return float(np.sqrt((e ** 2).mean())), float(e.abs().mean())5960def dm(base):61 """Diebold-Mariano test (squared loss, h=1, HLN small-sample correction), d = e_ridge^2 - e_base^2."""62 from scipy import stats63 e_r = (fc.loc[fc.model == "ridge", "forecast"].to_numpy() - fc.loc[fc.model == "ridge", "actual"].to_numpy())64 e_b = (fc.loc[fc.model == base, "forecast"].to_numpy() - fc.loc[fc.model == base, "actual"].to_numpy())65 d = e_r ** 2 - e_b ** 266 n = len(d)67 stat = d.mean() / np.sqrt(d.var(ddof=0) / n) * np.sqrt((n + 1 - 2 + 0) / n)68 p = 2 * stats.t.sf(abs(stat), df=n - 1)69 return float(stat), float(p)7071dm_n, p_n = dm("naive"); dm_m, p_m = dm("mean")72rn, an = err("naive"); rm, am = err("mean"); rr, ar = err("ridge")7374# Exclusion ledger (rows of the input file vs ridge training rows at the last origin)75t_last = ORIGINS[-1]76n1 = n_raw77n2 = int((df.index >= FIRST_TRAIN).sum())78n3 = int(((df.index >= FIRST_TRAIN) & (df.index <= t_last - pd.offsets.MonthBegin(1))).sum())79pd.DataFrame([80 {"step": 1, "rule": "drop months before 1992-03-01 (lags of ip_growth / neworder_chg unavailable; spec start)",81 "n_before": n1, "n_after": n2, "n_removed": n1 - n2, "rationale": "TASK.md training-row definition s >= 1992-03-01"},82 {"step": 2, "rule": "drop months s with s+1 > last origin 2026-08-01 (target not in info set)",83 "n_before": n2, "n_after": n3, "n_removed": n2 - n3, "rationale": "information-availability rule at last origin"},84]).to_csv(os.path.join(OUT, "exclusions.csv"), index=False)85assert n3 == n_train_last8687res = {88 "task_id": "F3-conf-03",89 "n_origins": len(ORIGINS),90 "first_origin": ORIGINS[0].strftime("%Y-%m-%d"),91 "last_origin": ORIGINS[-1].strftime("%Y-%m-%d"),92 "lambda": LAM,93 "normalisation_window": "expanding",94 "n_features": F.shape[1],95 "rmse_naive": rn, "mae_naive": an,96 "rmse_mean": rm, "mae_mean": am,97 "rmse_ridge": rr, "mae_ridge": ar,98 "rel_improvement_vs_naive": 1 - rr / rn,99 "rel_improvement_vs_mean": 1 - rr / rm,100 "dm_stat_ridge_vs_naive": dm_n, "dm_pvalue_ridge_vs_naive": p_n,101 "dm_stat_ridge_vs_mean": dm_m, "dm_pvalue_ridge_vs_mean": p_m,102 "feasibility": "complete",103 "feasibility_reason": None,104 "data_issues": [105 {"issue": "neworder_chg missing for 1992-01 and 1992-02 (series starts 1992-02)",106 "handling": "Not needed: first ridge training row is 1992-03 per spec; no imputation."},107 {"issue": "README says the first month of every change series is empty, but ip_growth, oil_chg, permit_chg are populated for 1992-01",108 "handling": "Values used as given (ip_growth 1992-01 enters only as lag 2 of training row 1992-03); noted as documentation inconsistency."},109 {"issue": "Extreme COVID-19 observations (ip_growth 2020-04 = -14.1, 2020-06 = +6.3; oil_chg 2020-03..05 ~ +/-55) inside every training window",110 "handling": "Retained as genuine data, as the ridge specification fixes the training rows; they affect mu, sigma and coefficients."},111 {"issue": "Current-vintage (revised) data; publication lags ignored",112 "handling": "Accepted per TASK.md; results are not real-time forecast accuracy."},113 ],114 "n_input_rows": {"ip_predictors.csv": n_raw},115 "n_analysed": n_train_last,116 "random_seed": None,117}118with open(os.path.join(OUT, "results.json"), "w") as f:119 json.dump(res, f, indent=2)120print(json.dumps({k: res[k] for k in ["rmse_naive", "mae_naive", "rmse_mean", "mae_mean", "rmse_ridge", "mae_ridge",121 "rel_improvement_vs_naive", "rel_improvement_vs_mean", "n_analysed",122 "dm_stat_ridge_vs_naive", "dm_pvalue_ridge_vs_naive", "dm_stat_ridge_vs_mean", "dm_pvalue_ridge_vs_mean"]}, indent=2))123