1"""Rolling-origin (expanding window) forecasts of US UNRATE: naive vs AR(2) on first differences.2Run from workspace root: ~/wp29/venv/bin/python out/analysis.py3Reads only data/UNRATE.csv; writes out/results.json, out/forecasts.csv, out/exclusions.csv. No randomness."""4import json, os5import numpy as np, pandas as pd6from scipy import stats78DATA = "data/UNRATE.csv"; OUT = "out"9P = 2 # AR order on first differences (fixed a priori)10HORIZONS = [1, 3]11MODEL = "AR2_diff"12origins = pd.date_range("2023-07-01", "2026-06-01", freq="MS")1314raw = pd.read_csv(DATA, dtype={"observation_date": str})15n_raw = len(raw)16raw["date"] = pd.to_datetime(raw["observation_date"], format="%Y-%m-%d")17y = raw.set_index("date")["UNRATE"].astype(float)18assert y.index.is_monotonic_increasing and not y.index.duplicated().any()19assert (y.index == pd.date_range(y.index[0], y.index[-1], freq="MS")).all()20missing_dates = [d.strftime("%Y-%m-%d") for d in y.index[y.isna()]]2122def ar_forecast(s, steps):23 """s: series with obs dated <= origin (may contain NaN). Returns forecasts 1..steps months24 ahead of the last *observed* month."""25 s = s.loc[:s.last_valid_index()]26 d = s.diff()27 X = pd.concat({f"l{i}": d.shift(i) for i in range(1, P + 1)}, axis=1)28 reg = pd.concat([d.rename("y"), X], axis=1).dropna() # drop rows touching missing values29 A = np.column_stack([np.ones(len(reg)), reg[[f"l{i}" for i in range(1, P + 1)]].values])30 beta, *_ = np.linalg.lstsq(A, reg["y"].values, rcond=None)31 # forecasting state: interior missing values linearly interpolated using only obs <= origin32 si = s.interpolate(method="linear", limit_area="inside")33 lags = list(si.diff().values[-P:][::-1]) # [d_T, d_{T-1}]34 level, out = si.values[-1], []35 for _ in range(steps):36 dhat = beta[0] + np.dot(beta[1:], lags[:P])37 level += dhat; out.append(level); lags = [dhat] + lags[:P - 1]38 return np.array(out), len(reg)3940rows = []41for t in origins:42 s = y.loc[:t] # information rule: only obs dated <= t43 last_obs = s.last_valid_index()44 gap = (t.year - last_obs.year) * 12 + (t.month - last_obs.month) # months origin is ahead of last obs45 fc, nreg = ar_forecast(s, gap + max(HORIZONS))46 for h in HORIZONS:47 tgt = t + pd.DateOffset(months=h)48 act = y.get(tgt, np.nan)49 rows.append(dict(origin_date=t.strftime("%Y-%m-%d"), target_date=tgt.strftime("%Y-%m-%d"),50 horizon=h, forecast=float(s.loc[last_obs]), actual=act, model="naive"))51 rows.append(dict(origin_date=t.strftime("%Y-%m-%d"), target_date=tgt.strftime("%Y-%m-%d"),52 horizon=h, forecast=float(fc[gap + h - 1]), actual=act, model=MODEL))53fc_df = pd.DataFrame(rows)54fc_df.to_csv(f"{OUT}/forecasts.csv", index=False, float_format="%.6f")5556res = {"task_id": "F3-dev-01", "model_name": MODEL,57 "model_description": ("AR(2) with intercept on monthly first differences of UNRATE, estimated by OLS on an "58 "expanding window starting 1948-01 (regression rows touching missing values dropped), "59 "re-estimated at every origin; level forecasts by iterating the difference recursion."),60 "n_origins": len(origins), "first_origin": origins[0].strftime("%Y-%m-%d"),61 "last_origin": origins[-1].strftime("%Y-%m-%d")}62dm = {}63for h in HORIZONS:64 sub = fc_df[fc_df.horizon == h]65 ev = sub.dropna(subset=["actual"])66 e = {m: (ev[ev.model == m].set_index("origin_date")["actual"] - ev[ev.model == m].set_index("origin_date")["forecast"])67 for m in ["naive", MODEL]}68 n = len(e["naive"]); assert n == len(e[MODEL])69 res[f"n_evaluated_h{h}"] = int(n)70 for m, tag in [("naive", "naive"), (MODEL, "model")]:71 res[f"rmse_{tag}_h{h}"] = float(np.sqrt(np.mean(e[m] ** 2)))72 res[f"mae_{tag}_h{h}"] = float(np.mean(np.abs(e[m])))73 # Diebold-Mariano (squared loss), HAC variance with h-1 Bartlett-free (rectangular) lags, HLN small-sample correction74 dl = (e[MODEL] ** 2 - e["naive"] ** 2).values75 dbar = dl.mean(); dc = dl - dbar76 gam = [np.sum(dc[k:] * dc[:n - k]) / n for k in range(h)]77 var = (gam[0] + 2 * sum(gam[1:])) / n78 if var <= 0: var = gam[0] / n79 hln = np.sqrt((n + 1 - 2 * h + h * (h - 1) / n) / n)80 stat = hln * dbar / np.sqrt(var)81 dm[h] = (float(stat), float(2 * stats.t.sf(abs(stat), df=n - 1)), float(dbar))82for h in HORIZONS:83 res[f"rel_improvement_rmse_h{h}"] = 1 - res[f"rmse_model_h{h}"] / res[f"rmse_naive_h{h}"]84for h in HORIZONS:85 res[f"dm_stat_h{h}"], res[f"dm_pvalue_h{h}"], res[f"mean_loss_diff_h{h}"] = dm[h]86res.update({87 "feasibility": "complete", "feasibility_reason": None,88 "data_issues": [89 {"issue": f"Missing UNRATE value(s) at {', '.join(missing_dates)} (empty field).",90 "handling": "Origin 2025-10-01: naive uses 2025-09 value; AR forecasts iterated from last observed month (2025-09). "91 "Pairs with missing target (origin 2025-09 h1, origin 2025-07 h3) kept in ledger with empty actual and "92 "excluded from error statistics. In later training windows, regression rows touching the missing month are "93 "dropped; the forecasting state uses linear interpolation between observations dated <= origin."},94 {"issue": "Final (revised) vintage, not real-time data.",95 "handling": "Errors measure forecastability of the revised series; no real-time validity claimed."},96 {"issue": "COVID-19 spike 2020-03..2020-05 (4.4 -> 14.8 -> 13.2) inside every training window.",97 "handling": "Kept as genuine data; no outlier treatment (specification fixed a priori, not tuned on test errors)."}],98 "n_input_rows": {"UNRATE.csv": int(n_raw)},99})100last_win = y.loc[:origins[-1]]101res["n_analysed"] = int(last_win.notna().sum())102res["n_regression_rows_last_window"] = int(ar_forecast(last_win, 1)[1])103res["random_seed"] = None104with open(f"{OUT}/results.json", "w") as f:105 json.dump(res, f, indent=2)106107# exclusion ledger (training data for the last window)108n1 = int(y.notna().sum()); n2 = res["n_analysed"]109pd.DataFrame([110 dict(step=1, rule="drop rows with missing UNRATE", n_before=n_raw, n_after=n1, n_removed=n_raw - n1,111 rationale="value not available (2025-10-01)"),112 dict(step=2, rule="drop observations dated after last origin 2026-06-01", n_before=n1, n_after=n2, n_removed=n1 - n2,113 rationale="information rule: training only on obs <= origin (these are used only as evaluation actuals)"),114]).to_csv(f"{OUT}/exclusions.csv", index=False)115print(json.dumps({k: v for k, v in res.items() if k not in ("data_issues", "model_description")}, indent=1))116