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 json, os5import numpy as np, pandas as pd67DATA = "data/UNRATE.csv"; OUT = "out"8P = 2 # AR lags on first differences (fixed a priori, not selected)9COVID = (pd.Timestamp("2020-04-01"), pd.Timestamp("2020-09-01")) # differences dated here are excluded from estimation10HORIZONS = [1, 3]11MODEL = "AR2_diff"1213raw = pd.read_csv(DATA, dtype={"UNRATE": "float64"})14n_rows = len(raw)15raw["observation_date"] = pd.to_datetime(raw["observation_date"])16s_full = raw.set_index("observation_date")["UNRATE"].asfreq("MS")17origins = pd.date_range("2023-07-01", "2026-06-01", freq="MS")1819def design(u, exclude_covid=True):20 """Regression rows y=du_t on [1, du_{t-1},...,du_{t-P}], using observed (non-interpolated) data only."""21 du = u.diff()22 X = pd.concat({f"l{k}": du.shift(k) for k in range(1, P + 1)}, axis=1)23 df = pd.concat([du.rename("y"), X], axis=1).dropna()24 if exclude_covid:25 bad = pd.Series(False, index=df.index)26 for k in range(0, P + 1): # drop row if y or any lag is dated inside the COVID window27 d = df.index - pd.DateOffset(months=k)28 bad |= (d >= COVID[0]) & (d <= COVID[1])29 df = df[~bad]30 return df3132def ar_forecast(u_obs, origin, hmax, exclude_covid=True):33 """Fit on data <= origin; return dict h->level forecast for h=1..hmax relative to origin."""34 u = u_obs[u_obs.index <= origin]35 df = design(u, exclude_covid)36 X = np.column_stack([np.ones(len(df)), df[[f"l{k}" for k in range(1, P + 1)]].values])37 beta, *_ = np.linalg.lstsq(X, df["y"].values, rcond=None)38 # state: interpolate interior gaps (uses only data <= origin); trailing gaps are forecast39 last_obs = u.last_valid_index()40 ui = u.loc[:last_obs].interpolate(limit_area="inside")41 gap = (origin.year - last_obs.year) * 12 + origin.month - last_obs.month42 du = list(ui.diff().values[-P:])43 level = ui.iloc[-1]; path = {}44 for step in range(1, gap + hmax + 1):45 d = beta[0] + sum(beta[k] * du[-k] for k in range(1, P + 1))46 du.append(d); level = level + d47 if step > gap: path[step - gap] = level48 return path, len(df), beta4950def run(exclude_covid=True):51 rows = []; n_last = None52 for t in origins:53 u = s_full[s_full.index <= t]54 naive = u.loc[u.last_valid_index()]55 path, n_used, beta = ar_forecast(s_full, t, max(HORIZONS), exclude_covid)56 n_last = n_used57 for h in HORIZONS:58 tgt = t + pd.DateOffset(months=h)59 act = s_full.get(tgt, np.nan)60 for m, f in (("naive", naive), (MODEL, path[h])):61 rows.append(dict(origin_date=t.strftime("%Y-%m-%d"), target_date=tgt.strftime("%Y-%m-%d"),62 horizon=h, forecast=round(float(f), 6),63 actual=None if pd.isna(act) else float(act), model=m))64 return pd.DataFrame(rows), n_last6566def stats(fc):67 out = {}68 for h in HORIZONS:69 out[f"n_evaluated_h{h}"] = int(fc[(fc.horizon == h) & (fc.model == "naive")].actual.notna().sum())70 for m, lab in (("naive", "naive"), (MODEL, "model")):71 e = fc[(fc.horizon == h) & (fc.model == m)].dropna(subset=["actual"])72 err = e.forecast - e.actual73 out[f"rmse_{lab}_h{h}"] = float(np.sqrt(np.mean(err ** 2)))74 out[f"mae_{lab}_h{h}"] = float(np.mean(np.abs(err)))75 out[f"rel_improvement_rmse_h{h}"] = 1 - out[f"rmse_model_h{h}"] / out[f"rmse_naive_h{h}"]76 return out7778def dm_test(fc, h):79 """Diebold-Mariano (squared loss) with Newey-West HAC var (h-1 lags) and Harvey-Leybourne-Newbold correction."""80 from scipy import stats as st81 a = fc[(fc.horizon == h) & (fc.model == "naive")].set_index("origin_date")82 b = fc[(fc.horizon == h) & (fc.model == MODEL)].set_index("origin_date")83 ok = a.actual.notna()84 d = ((a.forecast - a.actual) ** 2 - (b.forecast - b.actual) ** 2)[ok].values85 n = len(d); dbar = d.mean(); dc = d - dbar86 v = dc @ dc / n87 for k in range(1, h):88 v += 2 * (1 - k / h) * (dc[k:] @ dc[:-k]) / n89 dmstat = dbar / np.sqrt(v / n)90 hln = dmstat * np.sqrt((n + 1 - 2 * h + h * (h - 1) / n) / n)91 p = 2 * st.t.sf(abs(hln), df=n - 1)92 return float(hln), float(p)9394def boot_rel(fc, h, B=2000, block=6, seed=12345):95 """Moving-block bootstrap (over origins) 90% interval for 1 - RMSE_model/RMSE_naive."""96 rng = np.random.default_rng(seed)97 a = fc[(fc.horizon == h) & (fc.model == "naive")].reset_index(drop=True)98 b = fc[(fc.horizon == h) & (fc.model == MODEL)].reset_index(drop=True)99 ok = a.actual.notna().values100 ea = ((a.forecast - a.actual) ** 2).values[ok]; eb = ((b.forecast - b.actual) ** 2).values[ok]101 n = len(ea); reps = []102 for _ in range(B):103 idx = np.concatenate([np.arange(s, s + block) for s in rng.integers(0, n - block + 1, size=int(np.ceil(n / block)))])[:n]104 reps.append(1 - np.sqrt(eb[idx].mean() / ea[idx].mean()))105 return [float(np.quantile(reps, 0.05)), float(np.quantile(reps, 0.95))]106107if __name__ == "__main__":108 os.makedirs(OUT, exist_ok=True)109 fc, n_last = run(True)110 fc.to_csv(f"{OUT}/forecasts.csv", index=False, na_rep="")111 res = {"task_id": "F3-dev-01", "model_name": MODEL,112 "model_description": ("AR(2) with intercept on first differences of UNRATE, iterated to levels for h=1,3; "113 "lag order fixed a priori; re-estimated by OLS at each origin on an expanding window from 1948-03 "114 "(first available difference rows) to the origin, excluding regression rows whose differences fall in "115 "2020-04..2020-09 (COVID outliers) and rows involving the missing 2025-10 value."),116 "n_origins": len(origins), "first_origin": origins[0].strftime("%Y-%m-%d"),117 "last_origin": origins[-1].strftime("%Y-%m-%d")}118 res.update(stats(fc))119 res.update({"feasibility": "complete", "feasibility_reason": None,120 "data_issues": [{"issue": "UNRATE for 2025-10-01 is missing (empty field).",121 "handling": "Excluded as an evaluation target (origins 2025-09 h=1 and 2025-07 h=3 have empty actual). "122 "Naive at origin 2025-10 uses 2025-09. AR model: regression rows involving it are dropped; for the "123 "forecast state it is linearly interpolated once 2025-11 is observed (only data <= origin), and at "124 "origin 2025-10 the model iterates from 2025-09 (h+1 steps)."},125 {"issue": "Extreme COVID-19 outliers in 2020 (peak 14.8 in 2020-04; monthly changes up to +10.4 pp).",126 "handling": "Regression rows whose differences are dated 2020-04..2020-09 excluded from AR estimation (a fixed rule "127 "decided from pre-2023 data); sensitivity without exclusion reported in report.md."}],128 "n_input_rows": {"UNRATE.csv": n_rows}, "n_analysed": int(n_last), "random_seed": 12345})129 json.dump(res, open(f"{OUT}/results.json", "w"), indent=2)130 # supplementary diagnostics (not required keys)131 fc2, _ = run(False)132 supp = {"sensitivity_no_covid_exclusion": stats(fc2),133 "dm_hln": {f"h{h}": dict(zip(["stat", "p_value"], dm_test(fc, h))) for h in HORIZONS},134 "block_bootstrap_90ci_rel_improvement": {f"h{h}": boot_rel(fc, h) for h in HORIZONS},135 "last_origin_coefficients": list(map(float, ar_forecast(s_full, origins[-1], 3)[2]))}136 json.dump(supp, open(f"{OUT}/supplementary.json", "w"), indent=2)137 print(json.dumps(res, indent=2)); print(json.dumps(supp, indent=2))138