1"""Rolling-origin 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 and out/forecasts.csv. Deterministic (no randomness)."""4import json, os5import numpy as np, pandas as pd67ROOT = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))8OUT = os.path.join(ROOT, "out")9raw = pd.read_csv(os.path.join(ROOT, "data", "UNRATE.csv"), dtype={"observation_date": str})10n_input = len(raw)11s = pd.Series(raw["UNRATE"].values, index=pd.PeriodIndex(raw["observation_date"], freq="M"))12missing = [str(p.to_timestamp().date()) for p in s.index[s.isna()]]1314P = 2 # AR order on first differences, fixed a priori15MODEL = "AR2_diff"16origins = pd.period_range("2023-07", "2026-06", freq="M")17H = [1, 3]1819def fit_ar_diff(y):20 """OLS of dy_t on const, dy_{t-1..t-P}. y: pandas Series (no NaN, contiguous)."""21 dy = np.diff(y.values)22 X = np.column_stack([np.ones(len(dy) - P)] + [dy[P - k:len(dy) - k] for k in range(1, P + 1)])23 beta, *_ = np.linalg.lstsq(X, dy[P:], rcond=None)24 return beta, dy2526def forecast_ar_diff(beta, y_last, dy_hist, steps):27 hist = list(dy_hist[-P:]); lvl = y_last; path = []28 for _ in range(steps):29 d = beta[0] + sum(beta[k] * hist[-k] for k in range(1, P + 1))30 hist.append(d); lvl = lvl + d; path.append(lvl)31 return path3233rows = []; last_n = None34for t in origins:35 w = s[s.index <= t] # information set: dated <= t36 obs = w.dropna()37 last_obs_p = obs.index[-1] # latest observed month <= t38 naive = float(obs.iloc[-1])39 # model input: trim trailing missing, interpolate interior gaps using only data <= t40 wm = w.loc[:last_obs_p].interpolate(method="linear")41 beta, dy = fit_ar_diff(wm)42 maxsteps = max(H) + (t - last_obs_p).n43 path = forecast_ar_diff(beta, float(wm.iloc[-1]), dy, maxsteps)44 last_n = int(obs.shape[0])45 for h in H:46 tgt = t + h47 act = s.get(tgt, np.nan)48 steps = (tgt - last_obs_p).n49 for m, f in (("naive", naive), (MODEL, path[steps - 1])):50 rows.append({"origin_date": str(t.to_timestamp().date()), "target_date": str(tgt.to_timestamp().date()),51 "horizon": h, "forecast": round(float(f), 10),52 "actual": (None if pd.isna(act) else float(act)), "model": m})5354fc = pd.DataFrame(rows)55fc.to_csv(os.path.join(OUT, "forecasts.csv"), index=False, float_format="%.10g")5657res = {"task_id": "F3-dev-01", "model_name": MODEL,58 "model_description": ("AR(2) on first differences of UNRATE with intercept, estimated by OLS on the expanding window "59 "from 1948-01 (first usable difference regression row 1948-04) to the origin; levels forecast by "60 "iterating the difference recursion and cumulating; lag order fixed a priori, no selection."),61 "n_origins": len(origins), "first_origin": str(origins[0].to_timestamp().date()),62 "last_origin": str(origins[-1].to_timestamp().date())}63for h in H:64 ev = fc[(fc.horizon == h) & fc.actual.notna()]65 res[f"n_evaluated_h{h}"] = int((ev.model == "naive").sum())66 for m, key in (("naive", "naive"), (MODEL, "model")):67 e = ev[ev.model == m]; err = e.actual - e.forecast68 res[f"rmse_{key}_h{h}"] = float(np.sqrt(np.mean(err ** 2)))69 res[f"mae_{key}_h{h}"] = float(np.mean(np.abs(err)))70from scipy import stats71for h in H:72 ev = fc[(fc.horizon == h) & fc.actual.notna()]73 en = (ev[ev.model == "naive"].set_index("origin_date").eval("actual - forecast"))74 em = (ev[ev.model == MODEL].set_index("origin_date").eval("actual - forecast")).loc[en.index]75 d = (em ** 2 - en ** 2).values; n = len(d); dc = d - d.mean()76 lrv = np.mean(dc ** 2) + 2 * sum((1 - k / h) * np.mean(dc[k:] * dc[:-k]) for k in range(1, h)) # Bartlett, h-1 lags77 dm = d.mean() / np.sqrt(lrv / n)78 dm *= np.sqrt((n + 1 - 2 * h + h * (h - 1) / n) / n) # Harvey-Leybourne-Newbold correction79 res[f"dm_stat_h{h}"] = float(dm)80 res[f"dm_pvalue_h{h}"] = float(2 * stats.t.sf(abs(dm), df=n - 1))81for h in H:82 res[f"rel_improvement_rmse_h{h}"] = 1 - res[f"rmse_model_h{h}"] / res[f"rmse_naive_h{h}"]83res.update({"feasibility": "complete", "feasibility_reason": None,84 "data_issues": [{"issue": f"Missing UNRATE value at {d} (empty field)",85 "handling": ("Pairs targeting this month kept in ledger with empty actual and excluded from error stats; "86 "naive at origin 2025-10 uses 2025-09; model forecasts from 2025-09 at that origin and "87 "linearly interpolates the gap (using only data <= origin) at later origins.")} for d in missing]88 + [{"issue": "Final-vintage (revised) data, not real-time vintages",89 "handling": "Errors computed on data as given; no real-time validity claimed."}],90 "n_input_rows": {"UNRATE.csv": int(n_input)}, "n_analysed": last_n, "random_seed": None})91with open(os.path.join(OUT, "results.json"), "w") as f:92 json.dump(res, f, indent=2)93for k in ("n_evaluated_h1", "n_evaluated_h3", "rmse_naive_h1", "rmse_model_h1", "rmse_naive_h3", "rmse_model_h3",94 "rel_improvement_rmse_h1", "rel_improvement_rmse_h3", "dm_stat_h1", "dm_pvalue_h1", "dm_stat_h3", "dm_pvalue_h3"):95 print(k, res[k])96