"""F3-conf-03: rolling-origin ridge vs naive / expanding-mean forecasts of US IP growth (h=1). Run from workspace root: ~/wp29/venv/bin/python out/analysis.py Reads only data/; writes out/results.json, out/forecasts.csv, out/exclusions.csv. No randomness.""" import json, os import numpy as np import pandas as pd DATA = "data/ip_predictors.csv" OUT = "out" LAM = 50.0 FIRST_TRAIN = pd.Timestamp("1992-03-01") MEAN_START = pd.Timestamp("1992-04-01") ORIGINS = pd.date_range("2021-09-01", "2026-08-01", freq="MS") raw = pd.read_csv(DATA, parse_dates=["month"]) n_raw = len(raw) df = raw.sort_values("month").set_index("month") assert df.index.is_unique and (df.index == pd.date_range(df.index[0], df.index[-1], freq="MS")).all() # Feature matrix z_s (dated s; uses only values dated <= s) F = pd.DataFrame({ "ip_l0": df["ip_growth"], "ip_l1": df["ip_growth"].shift(1), "ip_l2": df["ip_growth"].shift(2), "oil_chg": df["oil_chg"], "fedfunds": df["fedfunds"], "permit_chg": df["permit_chg"], "neworder_chg": df["neworder_chg"], "tcu": df["tcu"], }) y_next = df["ip_growth"].shift(-1) # target of row s = ip_growth_{s+1} rows, excl = [], [] n_train_last = None for t in ORIGINS: avail = df.loc[:t] # information set: values dated <= t only Fa = F.loc[:t] ya = avail["ip_growth"].shift(-1) # recomputed inside the info set s_idx = Fa.index[(Fa.index >= FIRST_TRAIN) & (Fa.index <= t - pd.offsets.MonthBegin(1))] X = Fa.loc[s_idx].to_numpy(float) y = ya.loc[s_idx].to_numpy(float) assert np.isfinite(X).all() and np.isfinite(y).all() mu = X.mean(0); sd = X.std(0, ddof=1); ybar = y.mean() Z = (X - mu) / sd b = np.linalg.solve(Z.T @ Z + LAM * np.eye(Z.shape[1]), Z.T @ (y - ybar)) zt = Fa.loc[t].to_numpy(float) assert np.isfinite(zt).all() f_ridge = ybar + ((zt - mu) / sd) @ b f_naive = avail.loc[t, "ip_growth"] f_mean = avail.loc[MEAN_START:t, "ip_growth"].mean() tgt = t + pd.offsets.MonthBegin(1) actual = df.loc[tgt, "ip_growth"] # used only for scoring for m, f in (("naive", f_naive), ("mean", f_mean), ("ridge", f_ridge)): rows.append({"origin_date": t.strftime("%Y-%m-%d"), "target_date": tgt.strftime("%Y-%m-%d"), "horizon": 1, "forecast": float(f), "actual": float(actual), "model": m}) n_train_last = len(s_idx) fc = pd.DataFrame(rows) fc.to_csv(os.path.join(OUT, "forecasts.csv"), index=False, float_format="%.10f") def err(m): e = fc.loc[fc.model == m, "forecast"] - fc.loc[fc.model == m, "actual"] return float(np.sqrt((e ** 2).mean())), float(e.abs().mean()) def dm(base): """Diebold-Mariano test (squared loss, h=1, HLN small-sample correction), d = e_ridge^2 - e_base^2.""" from scipy import stats e_r = (fc.loc[fc.model == "ridge", "forecast"].to_numpy() - fc.loc[fc.model == "ridge", "actual"].to_numpy()) e_b = (fc.loc[fc.model == base, "forecast"].to_numpy() - fc.loc[fc.model == base, "actual"].to_numpy()) d = e_r ** 2 - e_b ** 2 n = len(d) stat = d.mean() / np.sqrt(d.var(ddof=0) / n) * np.sqrt((n + 1 - 2 + 0) / n) p = 2 * stats.t.sf(abs(stat), df=n - 1) return float(stat), float(p) dm_n, p_n = dm("naive"); dm_m, p_m = dm("mean") rn, an = err("naive"); rm, am = err("mean"); rr, ar = err("ridge") # Exclusion ledger (rows of the input file vs ridge training rows at the last origin) t_last = ORIGINS[-1] n1 = n_raw n2 = int((df.index >= FIRST_TRAIN).sum()) n3 = int(((df.index >= FIRST_TRAIN) & (df.index <= t_last - pd.offsets.MonthBegin(1))).sum()) pd.DataFrame([ {"step": 1, "rule": "drop months before 1992-03-01 (lags of ip_growth / neworder_chg unavailable; spec start)", "n_before": n1, "n_after": n2, "n_removed": n1 - n2, "rationale": "TASK.md training-row definition s >= 1992-03-01"}, {"step": 2, "rule": "drop months s with s+1 > last origin 2026-08-01 (target not in info set)", "n_before": n2, "n_after": n3, "n_removed": n2 - n3, "rationale": "information-availability rule at last origin"}, ]).to_csv(os.path.join(OUT, "exclusions.csv"), index=False) assert n3 == n_train_last res = { "task_id": "F3-conf-03", "n_origins": len(ORIGINS), "first_origin": ORIGINS[0].strftime("%Y-%m-%d"), "last_origin": ORIGINS[-1].strftime("%Y-%m-%d"), "lambda": LAM, "normalisation_window": "expanding", "n_features": F.shape[1], "rmse_naive": rn, "mae_naive": an, "rmse_mean": rm, "mae_mean": am, "rmse_ridge": rr, "mae_ridge": ar, "rel_improvement_vs_naive": 1 - rr / rn, "rel_improvement_vs_mean": 1 - rr / rm, "dm_stat_ridge_vs_naive": dm_n, "dm_pvalue_ridge_vs_naive": p_n, "dm_stat_ridge_vs_mean": dm_m, "dm_pvalue_ridge_vs_mean": p_m, "feasibility": "complete", "feasibility_reason": None, "data_issues": [ {"issue": "neworder_chg missing for 1992-01 and 1992-02 (series starts 1992-02)", "handling": "Not needed: first ridge training row is 1992-03 per spec; no imputation."}, {"issue": "README says the first month of every change series is empty, but ip_growth, oil_chg, permit_chg are populated for 1992-01", "handling": "Values used as given (ip_growth 1992-01 enters only as lag 2 of training row 1992-03); noted as documentation inconsistency."}, {"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", "handling": "Retained as genuine data, as the ridge specification fixes the training rows; they affect mu, sigma and coefficients."}, {"issue": "Current-vintage (revised) data; publication lags ignored", "handling": "Accepted per TASK.md; results are not real-time forecast accuracy."}, ], "n_input_rows": {"ip_predictors.csv": n_raw}, "n_analysed": n_train_last, "random_seed": None, } with open(os.path.join(OUT, "results.json"), "w") as f: json.dump(res, f, indent=2) print(json.dumps({k: res[k] for k in ["rmse_naive", "mae_naive", "rmse_mean", "mae_mean", "rmse_ridge", "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))