1"""F3-conf-03: rolling-origin ridge vs naive / expanding mean, h=1. Run from workspace root."""2import json, numpy as np, pandas as pd3DATA = "data/ip_predictors.csv"4LAM = 50.05d = pd.read_csv(DATA, parse_dates=["month"])6n_rows = len(d)7d = d.sort_values("month").set_index("month")8assert d.index.to_series().diff().dropna().dt.days.between(28, 31).all() # contiguous monthly9issues = []10na = d.isna().sum()11issues.append({"issue": f"Missing values: { {k: int(v) for k, v in na[na>0].items()} } (neworder_chg empty for 1992-01/02). "12 "README says the first month of each change series is empty, but ip_growth/oil_chg/permit_chg "13 "for 1992-01 are populated in the file.",14 "handling": "No imputation needed: ridge training starts at s=1992-03 (lags reach 1992-01 for ip_growth only), "15 "expanding mean starts 1992-04 as specified; populated 1992-01 values used as given."})16issues.append({"issue": "Large COVID-19 outliers in 2020 (ip_growth min %.2f, max %.2f)." % (d.ip_growth.min(), d.ip_growth.max()),17 "handling": "Kept as given; no winsorisation (not in specification)."})18ip = d["ip_growth"]19F = pd.DataFrame({"ip0": ip, "ip1": ip.shift(1), "ip2": ip.shift(2), "oil": d.oil_chg, "ff": d.fedfunds,20 "perm": d.permit_chg, "no": d.neworder_chg, "tcu": d.tcu})21Y = ip.shift(-1) # y_{s+1} aligned to s22origins = pd.date_range("2021-09-01", "2026-08-01", freq="MS")23rows = []; n_train_last = None24for t in origins:25 tgt = t + pd.DateOffset(months=1)26 actual = float(ip.loc[tgt])27 rows.append((t, tgt, float(ip.loc[t]), actual, "naive"))28 rows.append((t, tgt, float(ip.loc["1992-04-01":t].mean()), actual, "mean"))29 idx = F.index[(F.index >= "1992-03-01") & (F.index + pd.DateOffset(months=1) <= t)]30 X = F.loc[idx].to_numpy(float); y = Y.loc[idx].to_numpy(float)31 assert not np.isnan(X).any() and not np.isnan(y).any()32 mu = X.mean(0); sd = X.std(0, ddof=1); ybar = y.mean()33 Z = (X - mu) / sd34 b = np.linalg.solve(Z.T @ Z + LAM * np.eye(Z.shape[1]), Z.T @ (y - ybar))35 zt = (F.loc[t].to_numpy(float) - mu) / sd36 rows.append((t, tgt, float(ybar + zt @ b), actual, "ridge"))37 n_train_last = len(idx)38fc = pd.DataFrame(rows, columns=["origin_date", "target_date", "forecast", "actual", "model"])39fc.insert(2, "horizon", 1)40fc["origin_date"] = fc.origin_date.dt.strftime("%Y-%m-%d"); fc["target_date"] = fc.target_date.dt.strftime("%Y-%m-%d")41fc.to_csv("out/forecasts.csv", index=False)42err = {}43for m in ["naive", "mean", "ridge"]:44 e = fc[fc.model == m].eval("forecast-actual").to_numpy()45 err[m] = e46res = {"task_id": "F3-conf-03", "n_origins": len(origins),47 "first_origin": origins[0].strftime("%Y-%m-%d"), "last_origin": origins[-1].strftime("%Y-%m-%d"),48 "lambda": LAM, "normalisation_window": "expanding", "n_features": F.shape[1]}49for m, e in err.items():50 res[f"rmse_{m}"] = float(np.sqrt(np.mean(e**2))); res[f"mae_{m}"] = float(np.mean(np.abs(e)))51res["rel_improvement_vs_naive"] = 1 - res["rmse_ridge"] / res["rmse_naive"]52res["rel_improvement_vs_mean"] = 1 - res["rmse_ridge"] / res["rmse_mean"]53# Supplementary: Diebold-Mariano (squared loss, h=1, HLN small-sample correction)54from scipy import stats55def dm(e1, e2):56 dd = e1**2 - e2**2; n = len(dd); s = dd.mean() / np.sqrt(dd.var(ddof=0) / n)57 s *= np.sqrt((n - 1) / n); return float(s), float(2 * stats.t.sf(abs(s), n - 1))58supp = {"dm_ridge_vs_naive": dm(err["ridge"], err["naive"]), "dm_ridge_vs_mean": dm(err["ridge"], err["mean"])}59# bootstrap CI (iid over origins) for relative improvements60rng = np.random.default_rng(20261030); B = 5000; n = len(origins)61bs = {"naive": [], "mean": []}62for _ in range(B):63 i = rng.integers(0, n, n); r = np.sqrt(np.mean(err["ridge"][i]**2))64 for m in bs: bs[m].append(1 - r / np.sqrt(np.mean(err[m][i]**2)))65supp["boot95_rel_impr_vs_naive"] = [float(np.percentile(bs["naive"], q)) for q in (2.5, 97.5)]66supp["boot95_rel_impr_vs_mean"] = [float(np.percentile(bs["mean"], q)) for q in (2.5, 97.5)]67res.update({"feasibility": "complete", "feasibility_reason": None, "data_issues": issues,68 "n_input_rows": {"ip_predictors.csv": n_rows}, "n_analysed": n_train_last, "random_seed": 20261030,69 "supplementary": supp})70json.dump(res, open("out/results.json", "w"), indent=2)71print(json.dumps(res, indent=2))72