1"""F3-dev-01: rolling-origin forecasts of US 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 scipy import stats9from statsmodels.tsa.statespace.sarimax import SARIMAX1011warnings.filterwarnings("ignore")12OUT = "out"13os.makedirs(OUT, exist_ok=True)14MODEL = "ARIMA210_nc"15HORIZONS = [1, 3]1617raw = pd.read_csv("data/UNRATE.csv", dtype={"observation_date": str})18n_raw = len(raw)19s = pd.Series(pd.to_numeric(raw["UNRATE"], errors="coerce").values,20 index=pd.DatetimeIndex(pd.to_datetime(raw["observation_date"]), freq="MS"))21assert s.index.is_monotonic_increasing and not s.index.duplicated().any()22missing_dates = [d.strftime("%Y-%m-%d") for d in s.index[s.isna()]]2324origins = pd.date_range("2023-07-01", "2026-06-01", freq="MS")25rows = []26for t in origins:27 train = s.loc[:t] # only observations dated <= t28 # naive: latest available observation dated <= t29 naive = float(train.dropna().iloc[-1])30 # model: ARIMA(2,1,0), no constant, Gaussian MLE via Kalman filter (handles NaN)31 mod = SARIMAX(train, order=(2, 1, 0), trend="n")32 res = mod.fit(disp=False, method="lbfgs", maxiter=500)33 fc = res.forecast(steps=max(HORIZONS))34 for h in HORIZONS:35 tgt = t + pd.DateOffset(months=h)36 act = s.get(tgt, np.nan)37 for name, f in (("naive", naive), (MODEL, float(fc.iloc[h - 1]))):38 rows.append({"origin_date": t.strftime("%Y-%m-%d"), "target_date": tgt.strftime("%Y-%m-%d"),39 "horizon": h, "forecast": round(f, 10),40 "actual": (None if pd.isna(act) else float(act)), "model": name})41fc_df = pd.DataFrame(rows)42fc_df.to_csv(f"{OUT}/forecasts.csv", index=False, na_rep="")4344def err(model, h):45 d = fc_df[(fc_df.model == model) & (fc_df.horizon == h) & fc_df.actual.notna()].sort_values("origin_date")46 return (d.actual - d.forecast).to_numpy(dtype=float)4748def dm_test(e1, e2, h):49 """Diebold-Mariano (squared loss), HAC variance with h-1 rectangular lags, HLN small-sample correction."""50 d = e1 ** 2 - e2 ** 251 n = len(d); dbar = d.mean(); dc = d - dbar52 g = [np.sum(dc[k:] * dc[:n - k]) / n for k in range(h)]53 v = (g[0] + 2 * sum(g[1:])) / n54 if v <= 0:55 return None, None56 dm = dbar / np.sqrt(v)57 hln = np.sqrt((n + 1 - 2 * h + h * (h - 1) / n) / n)58 stat = dm * hln59 p = 2 * stats.t.sf(abs(stat), df=n - 1)60 return float(stat), float(p)6162R = {"task_id": "F3-dev-01", "model_name": MODEL,63 "model_description": ("ARIMA(2,1,0) without constant, i.e. AR(2) on first differences of UNRATE, "64 "estimated by Gaussian maximum likelihood (Kalman filter, statsmodels SARIMAX, missing months skipped) "65 "on an expanding window starting 1948-01-01 and ending at the origin, re-estimated at each origin; "66 "multi-step forecasts are iterated."),67 "n_origins": int(len(origins)), "first_origin": origins[0].strftime("%Y-%m-%d"),68 "last_origin": origins[-1].strftime("%Y-%m-%d")}69for h in HORIZONS:70 en, em = err("naive", h), err(MODEL, h)71 assert len(en) == len(em)72 R[f"n_evaluated_h{h}"] = int(len(en))73for h in HORIZONS:74 en = err("naive", h)75 R[f"rmse_naive_h{h}"] = float(np.sqrt(np.mean(en ** 2)))76 R[f"mae_naive_h{h}"] = float(np.mean(np.abs(en)))77for h in HORIZONS:78 em = err(MODEL, h)79 R[f"rmse_model_h{h}"] = float(np.sqrt(np.mean(em ** 2)))80 R[f"mae_model_h{h}"] = float(np.mean(np.abs(em)))81for h in HORIZONS:82 R[f"rel_improvement_rmse_h{h}"] = 1.0 - R[f"rmse_model_h{h}"] / R[f"rmse_naive_h{h}"]83for h in HORIZONS:84 st, p = dm_test(err("naive", h), err(MODEL, h), h)85 R[f"dm_stat_h{h}"] = st # >0 means model has lower squared error86 R[f"dm_pvalue_h{h}"] = p87n_missing_targets = {h: int(((fc_df.model == "naive") & (fc_df.horizon == h) & fc_df.actual.isna()).sum()) for h in HORIZONS}88R["feasibility"] = "complete"89R["feasibility_reason"] = None90R["data_issues"] = [91 {"issue": f"Missing observation(s) in UNRATE.csv: {missing_dates}",92 "handling": ("Naive at origin 2025-10-01 uses the 2025-09-01 value; the ARIMA likelihood/filter skips the missing month; "93 f"(origin,horizon) pairs targeting it are kept in forecasts.csv with empty actual and excluded from error statistics "94 f"(h1: {n_missing_targets[1]}, h3: {n_missing_targets[3]}).")},95 {"issue": "Data are the final (revised) vintage retrieved 2026-10-30, not real-time vintages.",96 "handling": "Evaluation is pseudo-out-of-sample on revised data; no real-time validity is claimed."},97 {"issue": "COVID-19 spike (2020-04 = 14.8%) lies inside every training window.",98 "handling": "Retained (genuine observation); noted as a limitation for Gaussian MLE of the AR coefficients."},99]100R["n_input_rows"] = {"UNRATE.csv": int(n_raw)}101last_train = s.loc[:origins[-1]]102R["n_analysed"] = int(last_train.notna().sum())103R["random_seed"] = None104with open(f"{OUT}/results.json", "w") as f:105 json.dump(R, f, indent=2)106107# exclusion ledger: raw rows -> observations in the last training window108n1 = int(s.notna().sum())109excl = pd.DataFrame([110 {"step": 1, "rule": "drop rows with missing UNRATE value", "n_before": n_raw, "n_after": n1,111 "n_removed": n_raw - n1, "rationale": "empty value field (2025-10-01); no data to use"},112 {"step": 2, "rule": "drop observations dated after last origin (2026-06-01) from training",113 "n_before": n1, "n_after": R["n_analysed"], "n_removed": n1 - R["n_analysed"],114 "rationale": "information rule: training uses only obs <= origin; these serve only as evaluation actuals"},115])116excl.to_csv(f"{OUT}/exclusions.csv", index=False)117print(json.dumps({k: R[k] for k in R if k.startswith(("rmse", "mae", "rel", "n_eval", "dm"))}, indent=1))118