"""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. Deterministic (no randomness).""" import json, os import numpy as np import pandas as pd DATA = "data/ip_predictors.csv" OUT = "out" LAMBDA = 50.0 TRAIN_START = pd.Timestamp("1992-03-01") MEAN_START = pd.Timestamp("1992-04-01") FIRST_ORIGIN = pd.Timestamp("2021-09-01") LAST_ORIGIN = pd.Timestamp("2026-08-01") FEATS = ["ip_growth", "ip_lag1", "ip_lag2", "oil_chg", "fedfunds", "permit_chg", "neworder_chg", "tcu"] raw = pd.read_csv(DATA, parse_dates=["month"]) n_raw = len(raw) d = raw.sort_values("month").set_index("month") assert d.index.is_unique # lags use only past values (shift forward in time) d["ip_lag1"] = d["ip_growth"].shift(1) d["ip_lag2"] = d["ip_growth"].shift(2) d["y_next"] = d["ip_growth"].shift(-1) # target of row s = ip_growth at s+1 origins = pd.date_range(FIRST_ORIGIN, LAST_ORIGIN, freq="MS") rows = [] last_ntrain = None for t in origins: tgt = t + pd.DateOffset(months=1) actual = float(d.loc[tgt, "ip_growth"]) avail = d.loc[:t] # information set: dated <= t # baselines f_naive = float(avail.loc[t, "ip_growth"]) f_mean = float(avail.loc[MEAN_START:t, "ip_growth"].mean()) # ridge training rows: TRAIN_START <= s, s+1 <= t -> s <= t - 1 month s_end = t - pd.DateOffset(months=1) tr = avail.loc[TRAIN_START:s_end] X = tr[FEATS].to_numpy(float) # target y_{s+1}: all s+1 <= t, so values are inside the information set y = avail["ip_growth"].shift(-1).loc[TRAIN_START:s_end].to_numpy(float) assert np.isfinite(X).all() and np.isfinite(y).all() mu = X.mean(axis=0) sd = X.std(axis=0, ddof=1) Z = (X - mu) / sd ybar = y.mean() b = np.linalg.solve(Z.T @ Z + LAMBDA * np.eye(Z.shape[1]), Z.T @ (y - ybar)) zt = (avail.loc[t, FEATS].to_numpy(float) - mu) / sd f_ridge = float(ybar + zt @ b) last_ntrain = len(tr) for m, f in (("naive", f_naive), ("mean", f_mean), ("ridge", f_ridge)): rows.append(dict(origin_date=t.strftime("%Y-%m-%d"), target_date=tgt.strftime("%Y-%m-%d"), horizon=1, forecast=f, actual=actual, model=m)) fc = pd.DataFrame(rows) os.makedirs(OUT, exist_ok=True) fc.to_csv(os.path.join(OUT, "forecasts.csv"), index=False, float_format="%.10f") def metrics(m): e = fc.loc[fc.model == m, "actual"].to_numpy() - fc.loc[fc.model == m, "forecast"].to_numpy() return float(np.sqrt(np.mean(e ** 2))), float(np.mean(np.abs(e))) rmse_n, mae_n = metrics("naive") rmse_m, mae_m = metrics("mean") rmse_r, mae_r = metrics("ridge") # exclusion ledger: raw rows -> training rows at the last origin s_end_last = LAST_ORIGIN - pd.DateOffset(months=1) n1 = int((d.index >= TRAIN_START).sum()) n2 = int(((d.index >= TRAIN_START) & (d.index <= s_end_last)).sum()) excl = pd.DataFrame([ dict(step=1, rule="month >= 1992-03-01 (lag-2 of ip_growth and neworder_chg required)", n_before=n_raw, n_after=n1, n_removed=n_raw - n1, rationale="1992-01/02 lack neworder_chg and/or two ip_growth lags; task fixes training start at 1992-03"), dict(step=2, rule="s + 1 <= last origin 2026-08-01 (s <= 2026-07-01)", n_before=n1, n_after=n2, n_removed=n1 - n2, rationale="information rule: targets dated after the origin may not be used in training (2026-08 row's target 2026-09 is the evaluation target; 2026-09 has no target)"), ]) excl.to_csv(os.path.join(OUT, "exclusions.csv"), index=False) assert n2 == last_ntrain # Uncertainty: Diebold-Mariano test (squared-error loss, h=1, Harvey-Leybourne-Newbold small-sample correction) from scipy import stats def dm(m1, m2): e1 = fc.loc[fc.model == m1, "actual"].to_numpy() - fc.loc[fc.model == m1, "forecast"].to_numpy() e2 = fc.loc[fc.model == m2, "actual"].to_numpy() - fc.loc[fc.model == m2, "forecast"].to_numpy() dl = e1 ** 2 - e2 ** 2 T = len(dl) stat = dl.mean() / np.sqrt(dl.var(ddof=0) / T) # h=1: no autocovariance terms stat *= np.sqrt((T + 1 - 2 + 0) / T) # HLN factor for h=1: sqrt((T+1-2h+h(h-1)/T)/T) p = 2 * stats.t.sf(abs(stat), df=T - 1) return float(stat), float(p) dm_rn, p_rn = dm("ridge", "naive") dm_rm, p_rm = dm("ridge", "mean") res = { "task_id": "F3-conf-03", "n_origins": int(len(origins)), "first_origin": origins[0].strftime("%Y-%m-%d"), "last_origin": origins[-1].strftime("%Y-%m-%d"), "lambda": LAMBDA, "normalisation_window": "expanding", "n_features": len(FEATS), "rmse_naive": rmse_n, "mae_naive": mae_n, "rmse_mean": rmse_m, "mae_mean": mae_m, "rmse_ridge": rmse_r, "mae_ridge": mae_r, "rel_improvement_vs_naive": 1 - rmse_r / rmse_n, "rel_improvement_vs_mean": 1 - rmse_r / rmse_m, "feasibility": "complete", "feasibility_reason": None, "data_issues": [ {"issue": "neworder_chg missing for 1992-01 and 1992-02", "handling": "these months are never ridge feature rows (training starts 1992-03); no imputation needed"}, {"issue": "README states 1992-01 is empty for every change series, but ip_growth, oil_chg and permit_chg have values in 1992-01", "handling": "values used as given; ip_growth 1992-01 serves as lag-2 feature of the 1992-03 training row"}, {"issue": "extreme COVID-2020 and 2008 observations (e.g. ip_growth 2020-04 = -14.1, oil_chg 2020-03..05 about +/-55)", "handling": "kept unchanged as genuine data; they enter training windows only (not the evaluation period)"}, {"issue": "current-vintage data with publication lags ignored", "handling": "results describe revised-data predictability under stylised timing, not real-time performance"}, ], "n_input_rows": {"ip_predictors.csv": int(n_raw)}, "n_analysed": int(last_ntrain), "random_seed": None, "dm_stat_ridge_vs_naive": dm_rn, "dm_pvalue_ridge_vs_naive": p_rn, "dm_stat_ridge_vs_mean": dm_rm, "dm_pvalue_ridge_vs_mean": p_rm, } 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))