1"""F3-conf-03: rolling-origin ridge vs naive / expanding-mean forecasts of IP growth (h=1).2Run from workspace root: ~/wp29/venv/bin/python out/analysis.py"""3import json, numpy as np, pandas as pd4from pathlib import Path56DATA = Path("data/ip_predictors.csv"); OUT = Path("out"); OUT.mkdir(exist_ok=True)7LAM = 50.08d = pd.read_csv(DATA, parse_dates=["month"])9n_rows = len(d)10d = d.sort_values("month").set_index("month")11issues = []12assert not d.index.duplicated().any()13assert (d.index.to_series().diff().dropna().dt.days.between(28, 31)).all(), "gaps in months"14# README says first month of each change series is empty; check15first_filled = [c for c in ["ip_growth", "oil_chg", "permit_chg"] if pd.notna(d.iloc[0][c])]16if first_filled:17 issues.append({"issue": f"README states the first month (1992-01) of each change series is empty, but {first_filled} are populated for 1992-01.",18 "handling": "Values used as given. Only ip_growth 1992-01 enters (as the lag-2 feature of the 1992-03 training row); the task spec explicitly starts training rows at 1992-03, implying this is intended."})19na = d.isna().sum(); na = na[na > 0]20issues.append({"issue": f"Missing values: {na.to_dict()} (neworder_chg 1992-01/02, documented in README).",21 "handling": "Not used: training rows start at 1992-03 per spec, so no imputation needed."})22big = d.loc[d.ip_growth.abs() > 4, "ip_growth"]23issues.append({"issue": "Extreme ip_growth values (|x|>4; 2008-09 hurricanes/financial crisis, 2020 COVID-19) (" + ", ".join(f"{k:%Y-%m}: {v:.2f}" for k, v in big.items()) + ").",24 "handling": "Genuine observations; kept, no winsorising (spec fixes the model). They lie in every training window and influence standardisation and coefficients."})25issues.append({"issue": "Current-vintage (revised) data; publication lags ignored.",26 "handling": "Results interpreted as pseudo-out-of-sample on revised data only, not real-time."})2728y = d["ip_growth"]29F = pd.DataFrame({"ip0": y, "ip1": y.shift(1), "ip2": y.shift(2), "oil": d.oil_chg,30 "ff": d.fedfunds, "perm": d.permit_chg, "no": d.neworder_chg, "tcu": d.tcu})31origins = pd.date_range("2021-09-01", "2026-08-01", freq="MS")32months = d.index33rows = []; n_train_last = None34for t in origins:35 tgt = t + pd.DateOffset(months=1)36 actual = float(y.loc[tgt])37 rows.append((t, tgt, float(y.loc[t]), actual, "naive"))38 rows.append((t, tgt, float(y.loc["1992-04-01":t].mean()), actual, "mean"))39 # training rows s: 1992-03 <= s, s+1 <= t40 s_idx = months[(months >= "1992-03-01") & (months <= t - pd.DateOffset(months=1))]41 X = F.loc[s_idx].to_numpy()42 yt = y.reindex(s_idx + pd.DateOffset(months=1)).to_numpy()43 assert not np.isnan(X).any() and not np.isnan(yt).any()44 assert (s_idx + pd.DateOffset(months=1)).max() <= t # no future info45 mu = X.mean(0); sd = X.std(0, ddof=1); ybar = yt.mean()46 Z = (X - mu) / sd47 b = np.linalg.solve(Z.T @ Z + LAM * np.eye(Z.shape[1]), Z.T @ (yt - ybar))48 zt = (F.loc[t].to_numpy() - mu) / sd49 rows.append((t, tgt, float(ybar + zt @ b), actual, "ridge"))50 n_train_last = len(s_idx)5152fc = pd.DataFrame(rows, columns=["origin_date", "target_date", "forecast", "actual", "model"])53fc.insert(2, "horizon", 1)54fc["origin_date"] = fc.origin_date.dt.strftime("%Y-%m-%d"); fc["target_date"] = fc.target_date.dt.strftime("%Y-%m-%d")55fc.to_csv(OUT / "forecasts.csv", index=False)5657res = {"task_id": "F3-conf-03", "n_origins": len(origins),58 "first_origin": f"{origins[0]:%Y-%m-%d}", "last_origin": f"{origins[-1]:%Y-%m-%d}",59 "lambda": LAM, "normalisation_window": "expanding", "n_features": F.shape[1]}60err = {}61for m in ["naive", "mean", "ridge"]:62 e = (fc[fc.model == m].forecast - fc[fc.model == m].actual).to_numpy(); err[m] = e63 res[f"rmse_{m}"] = float(np.sqrt(np.mean(e**2))); res[f"mae_{m}"] = float(np.mean(np.abs(e)))64res["rel_improvement_vs_naive"] = 1 - res["rmse_ridge"] / res["rmse_naive"]65res["rel_improvement_vs_mean"] = 1 - res["rmse_ridge"] / res["rmse_mean"]66res.update({"feasibility": "complete", "feasibility_reason": None, "data_issues": issues,67 "n_input_rows": {"ip_predictors.csv": n_rows}, "n_analysed": n_train_last, "random_seed": 0})6869# Supplementary uncertainty: Diebold-Mariano (squared loss, h=1, HAC lag 0 -> also report Newey-West lag 4) with HLN correction70from scipy import stats71def dm(e1, e2, L):72 dd = e1**2 - e2**2; T = len(dd); dbar = dd.mean(); u = dd - dbar73 v = u @ u / T + 2 * sum((1 - k/(L+1)) * (u[k:] @ u[:-k]) / T for k in range(1, L+1))74 stat = dbar / np.sqrt(v / T); stat *= np.sqrt((T - 1) / T) # HLN for h=175 return float(stat), float(2 * stats.t.sf(abs(stat), T - 1))76supp = {}77for base in ["naive", "mean"]:78 for L in [0, 4]:79 s, p = dm(err[base], err["ridge"], L)80 supp[f"dm_{base}_vs_ridge_lag{L}"] = {"stat": s, "p_value": p}81rng = np.random.default_rng(0)82T = len(err["ridge"]); bl = 6; B = 500083for base in ["naive", "mean"]:84 vals = []85 for _ in range(B):86 st = rng.integers(0, T - bl + 1, size=int(np.ceil(T / bl)))87 idx = np.concatenate([np.arange(s, s + bl) for s in st])[:T]88 vals.append(1 - np.sqrt(np.mean(err["ridge"][idx]**2)) / np.sqrt(np.mean(err[base][idx]**2)))89 supp[f"block_bootstrap_rel_impr_vs_{base}_95ci"] = [float(np.quantile(vals, .025)), float(np.quantile(vals, .975))]90json.dump(supp, open(OUT / "supplementary_uncertainty.json", "w"), indent=2)91json.dump(res, open(OUT / "results.json", "w"), indent=2)92print(json.dumps(res, indent=2)); print(json.dumps(supp, indent=2))93