1"""F3-dev-01: rolling-origin forecasts of UNRATE, naive vs ARIMA(2,1,0).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.4"""5import json, os, warnings6import numpy as np7import pandas as pd8from statsmodels.tsa.statespace.sarimax import SARIMAX910warnings.filterwarnings("ignore")11ROOT = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))12OUT = os.path.join(ROOT, "out")13MODEL_NAME = "ARIMA210"14HORIZONS = [1, 3]15SEED = None # no randomness used1617raw = pd.read_csv(os.path.join(ROOT, "data", "UNRATE.csv"))18n_raw = len(raw)19raw["observation_date"] = pd.to_datetime(raw["observation_date"])20s = raw.set_index("observation_date")["UNRATE"].astype(float)21s.index = pd.DatetimeIndex(s.index, freq="MS")22assert s.index.is_monotonic_increasing and not s.index.duplicated().any()2324origins = pd.date_range("2023-07-01", "2026-06-01", freq="MS")25missing = s[s.isna()].index2627rows = []28for t in origins:29 train = s.loc[:t] # only observations dated <= t30 naive = train.dropna().iloc[-1]31 mod = SARIMAX(train, order=(2, 1, 0), trend="n")32 res = mod.fit(disp=False, maxiter=200)33 fc = res.forecast(steps=max(HORIZONS))34 last_params = res.params # forecasts for t+1..t+3 (Kalman state handles gaps)35 for h in HORIZONS:36 tgt = t + pd.DateOffset(months=h)37 act = s.get(tgt, np.nan)38 rows.append(dict(origin_date=t, target_date=tgt, horizon=h, forecast=float(naive), actual=act, model="naive"))39 rows.append(dict(origin_date=t, target_date=tgt, horizon=h, forecast=float(fc.loc[tgt]), actual=act, model=MODEL_NAME))40led = pd.DataFrame(rows)41assert (led.target_date > led.origin_date).all()42out = led.copy()43out["origin_date"] = out.origin_date.dt.strftime("%Y-%m-%d")44out["target_date"] = out.target_date.dt.strftime("%Y-%m-%d")45out["forecast"] = out.forecast.round(10)46out.to_csv(os.path.join(OUT, "forecasts.csv"), index=False, float_format="%.10g")4748def err(model, h):49 d = led[(led.model == model) & (led.horizon == h)].dropna(subset=["actual"])50 e = d.actual - d.forecast51 return len(d), float(np.sqrt(np.mean(e ** 2))), float(np.mean(np.abs(e)))5253R = {"task_id": "F3-dev-01", "model_name": MODEL_NAME,54 "model_description": ("ARIMA(2,1,0) without constant (AR(2) on first differences of the level), "55 "estimated by Gaussian maximum likelihood (statsmodels SARIMAX, Kalman filter; missing months "56 "treated as missing) on an expanding window starting 1948-01-01 and ending at each origin, "57 "re-estimated at every origin; h-step forecasts by recursive iteration."),58 "n_origins": int(len(origins)),59 "first_origin": origins[0].strftime("%Y-%m-%d"), "last_origin": origins[-1].strftime("%Y-%m-%d")}60for h in HORIZONS:61 n, rn, mn = err("naive", h)62 n2, rm, mm = err(MODEL_NAME, h)63 assert n == n264 R[f"n_evaluated_h{h}"] = n65 R[f"rmse_naive_h{h}"], R[f"mae_naive_h{h}"] = rn, mn66 R[f"rmse_model_h{h}"], R[f"mae_model_h{h}"] = rm, mm67 R[f"rel_improvement_rmse_h{h}"] = 1 - rm / rn68order = ["task_id", "model_name", "model_description", "n_origins", "first_origin", "last_origin",69 "n_evaluated_h1", "n_evaluated_h3", "rmse_naive_h1", "mae_naive_h1", "rmse_naive_h3", "mae_naive_h3",70 "rmse_model_h1", "mae_model_h1", "rmse_model_h3", "mae_model_h3",71 "rel_improvement_rmse_h1", "rel_improvement_rmse_h3"]72# Diebold-Mariano test (squared-error loss, HAC variance with h-1 Newey-West lags,73# Harvey-Leybourne-Newbold small-sample correction, t(n-1) reference). Positive stat = model worse.74from scipy import stats75for h in HORIZONS:76 a = led[(led.model == "naive") & (led.horizon == h)].reset_index(drop=True)77 b = led[(led.model == MODEL_NAME) & (led.horizon == h)].reset_index(drop=True)78 ok = a.actual.notna().values79 d = ((b.actual - b.forecast) ** 2 - (a.actual - a.forecast) ** 2).values[ok]80 n = len(d); dbar = d.mean(); dc = d - dbar81 lrv = np.sum(dc ** 2) / n82 for k in range(1, h):83 lrv += 2 * (1 - k / h) * np.sum(dc[k:] * dc[:-k]) / n84 dm = dbar / np.sqrt(lrv / n)85 hln = np.sqrt((n + 1 - 2 * h + h * (h - 1) / n) / n)86 R[f"dm_stat_h{h}"] = float(dm * hln)87 R[f"dm_pvalue_h{h}"] = float(2 * stats.t.sf(abs(dm * hln), df=n - 1))88R["ar1_coef_last_window"] = float(last_params["ar.L1"])89R["ar2_coef_last_window"] = float(last_params["ar.L2"])90order += ["ar1_coef_last_window", "ar2_coef_last_window"]91order += ["dm_stat_h1", "dm_pvalue_h1", "dm_stat_h3", "dm_pvalue_h3"]92R = {k: R[k] for k in order}93miss_pairs = led[(led.model == "naive") & led.actual.isna()]94R["feasibility"] = "complete"95R["feasibility_reason"] = None96R["data_issues"] = [{"issue": f"Missing UNRATE value(s) at {', '.join(d.strftime('%Y-%m-%d') for d in missing)}",97 "handling": ("Kept as missing (no imputation). Naive at a missing origin uses the latest earlier "98 "observation; the ARIMA Kalman filter skips the missing month in estimation and forecasting; "99 "(origin,horizon) pairs with a missing target ("100 + "; ".join(f"{o.strftime('%Y-%m-%d')} h={h}" for o, h in zip(miss_pairs.origin_date, miss_pairs.horizon))101 + ") are in the ledger with empty actual and excluded from error statistics.")},102 {"issue": "Final (revised) vintage data, not real-time vintages",103 "handling": "Errors measure forecastability of the revised series; no real-time validity claimed."}]104R["n_input_rows"] = {"UNRATE.csv": int(n_raw)}105last_train = s.loc[:origins[-1]]106R["n_analysed"] = int(last_train.notna().sum())107R["random_seed"] = SEED108with open(os.path.join(OUT, "results.json"), "w") as f:109 json.dump(R, f, indent=2)110111# exclusion ledger (observations in the last training window)112n1 = int(s.notna().sum())113excl = pd.DataFrame([114 dict(step=1, rule="drop rows with missing UNRATE value", n_before=n_raw, n_after=n1, n_removed=n_raw - n1,115 rationale="empty value field in source (2025-10-01); treated as missing, not imputed"),116 dict(step=2, rule="keep observations dated <= last origin (2026-06-01) for the last training window",117 n_before=n1, n_after=R["n_analysed"], n_removed=n1 - R["n_analysed"],118 rationale="information-availability rule; later months used only as actuals"),119])120excl.to_csv(os.path.join(OUT, "exclusions.csv"), index=False)121for k in order[6:]:122 print(f"{k}: {R[k]}")123