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.py3"""4import json5import numpy as np6import pandas as pd78DATA = "data/UNRATE.csv"9MODEL = "AR2_diff"10P = 211HORIZONS = [1, 3]1213raw = pd.read_csv(DATA, dtype={"observation_date": str, "UNRATE": str}, keep_default_na=False)14n_rows = len(raw)15raw["date"] = pd.to_datetime(raw["observation_date"], format="%Y-%m-%d")16raw["y"] = pd.to_numeric(raw["UNRATE"].replace("", np.nan))17s = raw.set_index("date")["y"]18assert s.index.is_monotonic_increasing and not s.index.duplicated().any()19full_idx = pd.date_range(s.index[0], s.index[-1], freq="MS")20assert len(full_idx) == len(s), "non-contiguous months"21missing = [d.strftime("%Y-%m-%d") for d in s.index[s.isna()]]2223origins = pd.date_range("2023-07-01", "2026-06-01", freq="MS")24assert len(origins) == 36252627def ar_forecast(train, steps):28 """train: series of levels (no trailing NaN, interior NaN interpolated). AR(P)+const on diffs, OLS."""29 d = train.diff().dropna().values30 Y = d[P:]31 X = np.column_stack([np.ones(len(Y))] + [d[P - k:len(d) - k] for k in range(1, P + 1)])32 beta, *_ = np.linalg.lstsq(X, Y, rcond=None)33 hist = list(d[-P:])34 level = train.values[-1]35 path = []36 for _ in range(steps):37 dn = beta[0] + sum(beta[k] * hist[-k] for k in range(1, P + 1))38 hist.append(dn)39 level += dn40 path.append(level)41 return np.array(path), len(train)424344rows = []45n_last = None46for t in origins:47 avail = s[s.index <= t]48 last_obs_date = avail.last_valid_index()49 train = avail.loc[:last_obs_date] # drop trailing missing50 # interior gaps: linear interpolation using only data <= last_obs_date (<= t)51 train = train.interpolate(method="linear", limit_area="inside")52 naive = train.iloc[-1]53 gap = (t.year - last_obs_date.year) * 12 + (t.month - last_obs_date.month)54 path, n_used = ar_forecast(train, gap + max(HORIZONS))55 n_last = n_used56 for h in HORIZONS:57 target = t + pd.DateOffset(months=h)58 actual = s.get(target, np.nan)59 for model, fc in (("naive", naive), (MODEL, path[gap + h - 1])):60 rows.append({"origin_date": t.strftime("%Y-%m-%d"), "target_date": target.strftime("%Y-%m-%d"),61 "horizon": h, "forecast": round(float(fc), 6),62 "actual": None if pd.isna(actual) else float(actual), "model": model})6364fc = pd.DataFrame(rows)65fc.to_csv("out/forecasts.csv", index=False)6667res = {"task_id": "F3-dev-01", "model_name": MODEL,68 "model_description": ("AR(2) with intercept on first differences of the unemployment rate, estimated by OLS on an "69 "expanding window starting 1948-03 (first usable difference lag) up to each origin, interior "70 "missing months linearly interpolated within the window, iterated multi-step forecasts "71 "cumulated back to levels."),72 "n_origins": len(origins), "first_origin": origins[0].strftime("%Y-%m-%d"),73 "last_origin": origins[-1].strftime("%Y-%m-%d")}74for h in HORIZONS:75 ev = fc[(fc.horizon == h) & fc.actual.notna()]76 res[f"n_evaluated_h{h}"] = int((ev.model == "naive").sum())77for h in HORIZONS:78 for m, tag in (("naive", "naive"), (MODEL, "model")):79 ev = fc[(fc.horizon == h) & (fc.model == m) & fc.actual.notna()]80 e = ev.forecast - ev.actual81 res[f"rmse_{tag}_h{h}"] = float(np.sqrt(np.mean(e ** 2)))82 res[f"mae_{tag}_h{h}"] = float(np.mean(np.abs(e)))83for h in HORIZONS:84 res[f"rel_improvement_rmse_h{h}"] = 1 - res[f"rmse_model_h{h}"] / res[f"rmse_naive_h{h}"]8586# Diebold-Mariano style check (squared loss, HAC variance with h-1 lags), reported in report only87dm = {}88for h in HORIZONS:89 a = fc[(fc.horizon == h) & (fc.model == "naive") & fc.actual.notna()].set_index("origin_date")90 b = fc[(fc.horizon == h) & (fc.model == MODEL) & fc.actual.notna()].set_index("origin_date")91 dl = ((b.forecast - b.actual) ** 2 - (a.forecast - a.actual) ** 2).values92 n = len(dl); dbar = dl.mean(); dc = dl - dbar93 v = dc @ dc / n94 for k in range(1, h):95 v += 2 * (1 - k / h) * (dc[k:] @ dc[:-k]) / n96 stat = dbar / np.sqrt(v / n)97 from scipy import stats98 dm[f"h{h}"] = {"mean_loss_diff_model_minus_naive": float(dbar), "dm_stat": float(stat),99 "p_two_sided_t": float(2 * stats.t.sf(abs(stat), n - 1)), "n": n}100101res.update({"feasibility": "complete", "feasibility_reason": None,102 "data_issues": [{"issue": f"Missing observation(s) (empty value) at {', '.join(missing)}",103 "handling": ("Targets at missing months kept in ledger with empty actual and excluded from "104 "error statistics (affects origin 2025-09 h=1 and origin 2025-07 h=3); at origin "105 "2025-10 the naive forecast uses 2025-09 and the AR forecast iterates from 2025-09; "106 "for later origins the interior gap is linearly interpolated inside the training "107 "window (uses only data <= origin) for AR estimation.")}],108 "n_input_rows": {"UNRATE.csv": n_rows}, "n_analysed": int(n_last), "random_seed": None})109with open("out/results.json", "w") as f:110 json.dump(res, f, indent=2)111with open("out/dm_tests.json", "w") as f:112 json.dump(dm, f, indent=2)113print(json.dumps(res, indent=2)); print(json.dumps(dm, indent=2))114