1"""F3-dev-01: rolling-origin (expanding window) forecasts of UNRATE.2Model: AR(p) on monthly first differences, intercept, OLS, p chosen by BIC in 0..63at every origin, estimation sample from 1948, with impulse dummies for the42020-04 and 2020-05 COVID changes. Deterministic; no randomness.5Run from workspace root: ~/wp29/venv/bin/python out/analysis.py6"""7import json, os8import numpy as np, pandas as pd9from scipy import stats1011DATA = "data/UNRATE.csv"; OUT = "out"12MODEL = "AR_BIC_diff"13PMAX = 614DUMMY_DATES = [pd.Timestamp("2020-04-01"), pd.Timestamp("2020-05-01")]15HORIZONS = [1, 3]1617raw = pd.read_csv(DATA)18n_rows = len(raw)19df = raw.copy()20df["observation_date"] = pd.to_datetime(df["observation_date"])21df["UNRATE"] = pd.to_numeric(df["UNRATE"], errors="coerce")22df = df.set_index("observation_date").asfreq("MS") # full monthly grid23y_all = df["UNRATE"]2425origins = pd.date_range("2023-07-01", "2026-06-01", freq="MS")262728def fit_predict(y, origin, hmax, use_dummies=True):29 """y: level series on monthly grid, already truncated to dates <= origin."""30 dy = y.diff() # NaN where either month missing -> those rows dropped in estimation31 best = None32 for p in range(PMAX + 1):33 X = [np.ones(len(dy))]34 for k in range(1, p + 1):35 X.append(dy.shift(k).values)36 cols = len(X)37 if use_dummies:38 for d in DUMMY_DATES:39 if d <= origin:40 X.append((dy.index == d).astype(float))41 X = np.column_stack(X); Y = dy.values42 ok = ~np.isnan(Y) & ~np.isnan(X).any(1)43 # common estimation sample across p for comparable BIC: start after PMAX lags44 ok &= np.arange(len(Y)) > PMAX45 Xo, Yo = X[ok], Y[ok]46 b, *_ = np.linalg.lstsq(Xo, Yo, rcond=None)47 e = Yo - Xo @ b; n = len(Yo); kpar = Xo.shape[1]48 bic = n * np.log(e @ e / n) + kpar * np.log(n)49 if best is None or bic < best[0]:50 best = (bic, p, b[:cols], n)51 _, p, b, nreg = best52 # state for forecasting: interior missing levels linearly interpolated (uses only data <= origin);53 yy = y.copy()54 last_obs = yy.last_valid_index()55 yy.loc[:last_obs] = yy.loc[:last_obs].interpolate()56 dyy = list(yy.loc[:last_obs].diff().values)57 level = yy.loc[last_obs]58 steps = (origin.year - last_obs.year) * 12 + origin.month - last_obs.month + hmax59 fc = []60 for _ in range(steps):61 d = b[0] + sum(b[k] * dyy[-k] for k in range(1, p + 1))62 dyy.append(d); level = level + d; fc.append(level)63 # fc[i] is forecast for last_obs + i+1 months; return by horizon relative to origin64 off = steps - hmax65 return {h: fc[off + h - 1] for h in range(1, hmax + 1)}, p, nreg666768def run(use_dummies=True):69 rows, ps = [], []70 for t in origins:71 y = y_all.loc[:t]72 naive = y.loc[y.last_valid_index()]73 fc, p, nreg = fit_predict(y, t, max(HORIZONS), use_dummies)74 ps.append(p)75 for h in HORIZONS:76 tgt = t + pd.DateOffset(months=h)77 act = y_all.get(tgt, np.nan)78 for m, f in (("naive", naive), (MODEL, fc[h])):79 rows.append(dict(origin_date=t.strftime("%Y-%m-%d"), target_date=tgt.strftime("%Y-%m-%d"),80 horizon=h, forecast=round(float(f), 6),81 actual=None if pd.isna(act) else float(act), model=m))82 return pd.DataFrame(rows), ps, nreg838485led, ps, nreg_last = run(True)86led.to_csv(os.path.join(OUT, "forecasts.csv"), index=False, float_format="%.6f")87led_nd, _, _ = run(False) # sensitivity only (no COVID dummies)888990def errs(L, m, h):91 s = L[(L.model == m) & (L.horizon == h) & L.actual.notna()]92 e = s.actual - s.forecast93 return e.values, s.origin_date.values949596res = {"task_id": "F3-dev-01", "model_name": MODEL,97 "model_description": ("AR(p) on monthly first differences of UNRATE with intercept, p in 0..6 re-selected by BIC at "98 "each origin, OLS on an expanding window starting 1948 (first usable change 1948-02, common "99 "sample after 6 lags) with impulse dummies for the 2020-04 and 2020-05 changes; iterated "100 "multi-step forecasts cumulated to levels."),101 "n_origins": len(origins), "first_origin": origins[0].strftime("%Y-%m-%d"),102 "last_origin": origins[-1].strftime("%Y-%m-%d")}103extra = {}104for h in HORIZONS:105 en, _ = errs(led, "naive", h); em, _ = errs(led, MODEL, h)106 res[f"n_evaluated_h{h}"] = int(len(en))107 for nm, e in (("naive", en), ("model", em)):108 res[f"rmse_{nm}_h{h}"] = float(np.sqrt(np.mean(e ** 2)))109 res[f"mae_{nm}_h{h}"] = float(np.mean(np.abs(e)))110for h in HORIZONS:111 res[f"rel_improvement_rmse_h{h}"] = 1 - res[f"rmse_model_h{h}"] / res[f"rmse_naive_h{h}"]112113# Uncertainty: Diebold-Mariano (squared loss), HAC variance with h-1 Bartlett-free rectangular lags114# (Harvey-Leybourne-Newbold small-sample correction), t(n-1) reference. Note the evaluation set has one gap.115for h in HORIZONS:116 en, _ = errs(led, "naive", h); em, _ = errs(led, MODEL, h)117 d = em ** 2 - en ** 2; n = len(d); dbar = d.mean(); dc = d - dbar118 g = [dc @ dc / n] + [dc[k:] @ dc[:-k] / n for k in range(1, h)]119 var = (g[0] + 2 * sum(g[1:])) / n120 dm = dbar / np.sqrt(var) if var > 0 else np.nan121 hln = dm * np.sqrt((n + 1 - 2 * h + h * (h - 1) / n) / n)122 extra[f"dm_stat_h{h}"] = float(hln)123 extra[f"dm_pvalue_two_sided_h{h}"] = float(2 * stats.t.sf(abs(hln), n - 1))124 e2, _ = errs(led_nd, MODEL, h)125 extra[f"sensitivity_no_covid_dummies_rmse_model_h{h}"] = float(np.sqrt(np.mean(e2 ** 2)))126 extra[f"sensitivity_no_covid_dummies_rel_improvement_h{h}"] = float(1 - np.sqrt(np.mean(e2 ** 2)) / res[f"rmse_naive_h{h}"])127extra["selected_p_counts"] = {str(k): int(v) for k, v in zip(*np.unique(ps, return_counts=True))}128extra["n_regression_rows_last_window"] = int(nreg_last)129130n_missing = int(y_all.isna().sum())131res.update({132 "feasibility": "complete", "feasibility_reason": None,133 "data_issues": [134 {"issue": f"{n_missing} missing observation(s): " + ", ".join(d.strftime('%Y-%m-%d') for d in y_all.index[y_all.isna()]) + " (empty value field)",135 "handling": "Target pairs with missing actual (origin 2025-09 h1, origin 2025-07 h3) kept in ledger with empty actual and excluded from error statistics; differences touching the missing month dropped from estimation; for forecasting at later origins the interior gap is linearly interpolated from neighbours dated <= origin; at origin 2025-10 the naive forecast uses 2025-09 and the model forecasts 2025-10 as an intermediate step."},136 {"issue": "Extreme COVID-19 jump (2020-04 +10.4 pp, 2020-05 -1.6 pp) inside every training window",137 "handling": "Impulse dummies on the 2020-04 and 2020-05 monthly changes (dates known before all origins); sensitivity without dummies reported in report.md."},138 {"issue": "Final-vintage data (revised seasonal factors), not real-time vintages",139 "handling": "Errors computed on data as given; no real-time validity claimed."}],140 "n_input_rows": {"UNRATE.csv": int(n_rows)},141 "n_analysed": int(y_all.loc[:origins[-1]].notna().sum()),142 "random_seed": None,143})144with open(os.path.join(OUT, "results.json"), "w") as f:145 json.dump(res, f, indent=2)146with open(os.path.join(OUT, "supplementary.json"), "w") as f:147 json.dump(extra, f, indent=2)148print(json.dumps(res, indent=2)); print(json.dumps(extra, indent=2))149