1"""F3-dev-01: rolling-origin (expanding window) forecasts of 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. Deterministic (no randomness)."""4import json5import os6import numpy as np7import pandas as pd89DATA = "data/UNRATE.csv"10OUT = "out"11P = 2 # AR lag order on first differences (fixed a priori)12HORIZONS = [1, 3]13MODEL_NAME = "AR2_diff"14ORIGINS = pd.date_range("2023-07-01", "2026-06-01", freq="MS")1516raw = pd.read_csv(DATA, dtype={"observation_date": str})17n_raw = len(raw)18raw["date"] = pd.to_datetime(raw["observation_date"], format="%Y-%m-%d")19s = raw.set_index("date")["UNRATE"].astype(float)20assert s.index.is_monotonic_increasing and not s.index.duplicated().any()21assert (s.index == pd.date_range(s.index[0], s.index[-1], freq="MS")).all()22missing_dates = [d.strftime("%Y-%m-%d") for d in s.index[s.isna()]]232425def fit_ar_diff(y):26 """OLS of dy_t on const + dy_{t-1..t-P}, complete cases only. y: level series (may contain NaN)."""27 dy = y.diff()28 X = pd.concat([dy.shift(k) for k in range(1, P + 1)], axis=1)29 df = pd.concat([dy, X], axis=1).dropna()30 Y = df.iloc[:, 0].to_numpy()31 Xm = np.column_stack([np.ones(len(df)), df.iloc[:, 1:].to_numpy()])32 beta, *_ = np.linalg.lstsq(Xm, Y, rcond=None)33 return beta, len(df)343536def forecast_ar_diff(y, origin, hmax):37 """Forecast levels at origin+1..origin+hmax using only y dated <= origin.38 Missing levels inside the window (and the origin month if missing) are filled sequentially39 with the model's own one-step forecast from earlier data (no information after the origin)."""40 y = y.loc[:origin]41 beta, nobs = fit_ar_diff(y)42 lv = y.copy()43 # fill internal missing values (needed only if they are among the last P+1 levels used)44 for i in np.where(lv.isna().to_numpy())[0]:45 if i >= len(lv) - (P + 1) - hmax: # only recent gaps matter for the recursion46 d = lv.diff().iloc[i - P:i].to_numpy()[::-1] # dy_{i-1},...,dy_{i-P}47 lv.iloc[i] = lv.iloc[i - 1] + beta[0] + beta[1:] @ d48 assert not lv.iloc[-(P + 1):].isna().any()49 hist_d = list(lv.diff().iloc[-P:].to_numpy()[::-1]) # most recent first50 level = lv.iloc[-1]51 out = []52 for _ in range(hmax):53 d = beta[0] + beta[1:] @ np.array(hist_d[:P])54 level = level + d55 hist_d.insert(0, d)56 out.append(level)57 return out, beta, nobs585960rows = []61last_window_n = None62for t in ORIGINS:63 window = s.loc[:t]64 naive = window.dropna().iloc[-1]65 fc, beta, nobs = forecast_ar_diff(s, t, max(HORIZONS))66 last_window_n = int(window.notna().sum())67 for h in HORIZONS:68 target = t + pd.DateOffset(months=h)69 actual = s.get(target, np.nan)70 for model, f in (("naive", naive), (MODEL_NAME, fc[h - 1])):71 rows.append(dict(origin_date=t.strftime("%Y-%m-%d"), target_date=target.strftime("%Y-%m-%d"),72 horizon=h, forecast=round(float(f), 10),73 actual=(None if pd.isna(actual) else float(actual)), model=model))7475fc_df = pd.DataFrame(rows)76assert len(fc_df) == len(ORIGINS) * len(HORIZONS) * 277assert (pd.to_datetime(fc_df.target_date) > pd.to_datetime(fc_df.origin_date)).all()78fc_df.to_csv(os.path.join(OUT, "forecasts.csv"), index=False, float_format="%.10g")7980res = {"task_id": "F3-dev-01", "model_name": MODEL_NAME,81 "model_description": ("AR(2) with intercept on first differences of UNRATE, levels forecast by iterating the "82 "difference equation; estimated by OLS on complete cases over an expanding window from "83 "1948 (first usable dependent obs 1948-04) to the origin; lag order fixed a priori; a "84 "missing in-window level is replaced for forecasting by the model's own forecast from "85 "earlier data."),86 "n_origins": len(ORIGINS), "first_origin": ORIGINS[0].strftime("%Y-%m-%d"),87 "last_origin": ORIGINS[-1].strftime("%Y-%m-%d")}88for h in HORIZONS:89 sub = fc_df[(fc_df.horizon == h) & fc_df.actual.notna()]90 n_ev = sub[sub.model == "naive"].shape[0]91 res[f"n_evaluated_h{h}"] = int(n_ev)92for h in HORIZONS:93 for model, key in (("naive", "naive"), (MODEL_NAME, "model")):94 sub = fc_df[(fc_df.horizon == h) & (fc_df.model == model) & fc_df.actual.notna()]95 e = sub.actual.to_numpy() - sub.forecast.to_numpy()96 res[f"rmse_{key}_h{h}"] = float(np.sqrt(np.mean(e ** 2)))97 res[f"mae_{key}_h{h}"] = float(np.mean(np.abs(e)))98for h in HORIZONS:99 res[f"rel_improvement_rmse_h{h}"] = 1.0 - res[f"rmse_model_h{h}"] / res[f"rmse_naive_h{h}"]100101excluded_pairs = fc_df[fc_df.actual.isna() & (fc_df.model == "naive")][["origin_date", "horizon"]]102res.update({103 "feasibility": "complete", "feasibility_reason": None,104 "data_issues": [{"issue": f"Missing UNRATE value at {d} (empty field).",105 "handling": ("Kept as missing. Naive forecast at origin 2025-10-01 uses 2025-09-01. AR estimation "106 "uses complete cases only; for forecasting, the missing level is replaced by the "107 "model's own forecast from data before it. (origin, horizon) pairs whose target is "108 "this month are kept in forecasts.csv with empty actual and excluded from error "109 "statistics: " + "; ".join(f"{r.origin_date} h={r.horizon}"110 for r in excluded_pairs.itertuples()) + ".")}111 for d in missing_dates] + [112 {"issue": "Data are the final (revised) vintage, not real-time vintages.",113 "handling": "Errors are pseudo-out-of-sample on revised data; no real-time validity is claimed."}],114 "n_input_rows": {"UNRATE.csv": int(n_raw)},115 "n_analysed": last_window_n,116 "random_seed": None,117})118# diagnostic (not a required key): coefficients at last origin119res["diag_ar_coef_last_origin"] = [float(b) for b in beta]120res["diag_ar_nobs_last_origin"] = int(nobs)121122123# Diebold-Mariano test (squared-error loss), HAC long-run variance with rectangular kernel, h-1 lags,124# Harvey-Leybourne-Newbold small-sample correction; t(n-1) reference. Diagnostic only.125from scipy import stats126for h in HORIZONS:127 a = fc_df[(fc_df.horizon == h) & fc_df.actual.notna()]128 en = (a[a.model == "naive"].actual - a[a.model == "naive"].forecast).to_numpy()129 em = (a[a.model == MODEL_NAME].actual - a[a.model == MODEL_NAME].forecast).to_numpy()130 d = em ** 2 - en ** 2 # positive => model worse131 n = len(d); dbar = d.mean(); dc = d - dbar132 lrv = np.sum(dc ** 2) / n + 2 * sum(np.sum(dc[k:] * dc[:-k]) / n for k in range(1, h))133 dm = dbar / np.sqrt(lrv / n)134 hln = dm * np.sqrt((n + 1 - 2 * h + h * (h - 1) / n) / n)135 res[f"diag_dm_stat_h{h}"] = float(hln)136 res[f"diag_dm_pvalue_two_sided_h{h}"] = float(2 * stats.t.sf(abs(hln), df=n - 1))137138with open(os.path.join(OUT, "results.json"), "w") as f:139 json.dump(res, f, indent=2)140141# exclusion ledger (observation-level)142n1 = int(s.notna().sum())143n2 = int(s.loc[:ORIGINS[-1]].notna().sum())144pd.DataFrame([145 dict(step=1, rule="drop rows with missing UNRATE", n_before=n_raw, n_after=n1, n_removed=n_raw - n1,146 rationale="empty value field (" + ", ".join(missing_dates) + ")"),147 dict(step=2, rule="drop observations dated after last origin 2026-06-01 from training", n_before=n1, n_after=n2,148 n_removed=n1 - n2, rationale="information rule: training uses only data dated <= origin (used only as targets)"),149]).to_csv(os.path.join(OUT, "exclusions.csv"), index=False)150151for h in HORIZONS:152 print(f"h={h}: n={res[f'n_evaluated_h{h}']} RMSE naive={res[f'rmse_naive_h{h}']:.4f} "153 f"model={res[f'rmse_model_h{h}']:.4f} rel.impr={res[f'rel_improvement_rmse_h{h}']:.4f}")154