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.py3Reads only data/UNRATE.csv; writes out/results.json, out/forecasts.csv, out/exclusions.csv. Deterministic (no randomness)."""4import json, os5import numpy as np6import pandas as pd7from scipy import stats89DATA = "data/UNRATE.csv"10OUT = "out"11P = 2 # AR order on first differences12HORIZONS = (1, 3)13MODEL_NAME = "AR2_diff"14ORIGINS = pd.date_range("2023-07-01", "2026-06-01", freq="MS")1516raw = pd.read_csv(DATA, dtype={"observation_date": str})17n_raw = len(raw)18s = pd.Series(raw["UNRATE"].values, index=pd.to_datetime(raw["observation_date"]), dtype=float)19assert s.index.is_monotonic_increasing and not s.index.duplicated().any()20assert len(pd.date_range(s.index[0], s.index[-1], freq="MS")) == len(s)212223def fit_ar_diff(y):24 """OLS of d_t on const, d_{t-1..t-P}; rows with any NaN dropped. y: level series up to origin."""25 d = y.diff()26 X = pd.concat([d.shift(k) for k in range(1, P + 1)], axis=1)27 df = pd.concat([d, X], axis=1).dropna()28 Y = df.iloc[:, 0].values29 Xm = np.column_stack([np.ones(len(df)), df.iloc[:, 1:].values])30 beta, *_ = np.linalg.lstsq(Xm, Y, rcond=None)31 return beta, len(df)323334def forecast_model(y, origin, h):35 """Forecast level at origin+h using only y (already truncated to <= origin)."""36 beta, nreg = fit_ar_diff(y)37 last_obs = y.last_valid_index()38 # state: interior missing values interpolated from data <= origin (state only, not training)39 ys = y.loc[:last_obs].interpolate(method="linear", limit_area="inside")40 steps = h + (len(pd.date_range(last_obs, origin, freq="MS")) - 1) # extra steps if origin missing41 level = ys.iloc[-1]42 dl = list(ys.diff().iloc[-P:].values[::-1]) # most recent first43 for _ in range(steps):44 dn = beta[0] + sum(beta[k + 1] * dl[k] for k in range(P))45 level += dn46 dl = [dn] + dl[:-1]47 return level, nreg484950rows = []51for t in ORIGINS:52 y = s.loc[:t] # information set: observations dated <= t53 naive = y.dropna().iloc[-1]54 for h in HORIZONS:55 tgt = t + pd.DateOffset(months=h)56 actual = s.get(tgt, np.nan)57 fm, nreg = forecast_model(y, t, h)58 for m, f in (("naive", naive), (MODEL_NAME, fm)):59 rows.append(dict(origin_date=t.strftime("%Y-%m-%d"), target_date=tgt.strftime("%Y-%m-%d"),60 horizon=h, forecast=round(float(f), 10),61 actual=(None if pd.isna(actual) else float(actual)), model=m))62fc = pd.DataFrame(rows)63os.makedirs(OUT, exist_ok=True)64fc.to_csv(os.path.join(OUT, "forecasts.csv"), index=False)6566res = {"task_id": "F3-dev-01", "model_name": MODEL_NAME,67 "model_description": ("AR(2) with intercept on first differences of UNRATE, estimated by OLS on all available "68 "observations from 1948-01 up to each origin (expanding window; rows touching missing "69 "values dropped), iterated h-step forecasts of differences cumulated onto the last "70 "observed level."),71 "n_origins": int(len(ORIGINS)), "first_origin": ORIGINS[0].strftime("%Y-%m-%d"),72 "last_origin": ORIGINS[-1].strftime("%Y-%m-%d")}73extra = {}74for h in HORIZONS:75 sub = fc[(fc.horizon == h) & fc.actual.notna()]76 en = (sub[sub.model == "naive"].set_index("origin_date").pipe(lambda d: d.actual - d.forecast))77 em = (sub[sub.model == MODEL_NAME].set_index("origin_date").pipe(lambda d: d.actual - d.forecast))78 assert (en.index == em.index).all()79 n = len(en)80 res[f"n_evaluated_h{h}"] = int(n)81 for name, e in (("naive", en), ("model", em)):82 res[f"rmse_{name}_h{h}"] = float(np.sqrt(np.mean(e ** 2)))83 res[f"mae_{name}_h{h}"] = float(np.mean(np.abs(e)))84 res[f"rel_improvement_rmse_h{h}"] = 1 - res[f"rmse_model_h{h}"] / res[f"rmse_naive_h{h}"]85 # Diebold-Mariano (squared-error loss), rectangular HAC with h-1 lags, Harvey-Leybourne-Newbold correction.86 # Note: with a missing pair the series has a gap; treated as consecutive (approximation, documented).87 dlt = (em.values ** 2 - en.values ** 2)88 dbar = dlt.mean(); dc = dlt - dbar89 gam = [np.sum(dc[k:] * dc[:n - k]) / n for k in range(h)]90 var = (gam[0] + 2 * sum(gam[1:])) / n91 if var <= 0:92 var = gam[0] / n93 dm = dbar / np.sqrt(var)94 hln = dm * np.sqrt((n + 1 - 2 * h + h * (h - 1) / n) / n)95 extra[f"dm_stat_h{h}"] = float(hln)96 extra[f"dm_pvalue_two_sided_h{h}"] = float(2 * stats.t.sf(abs(hln), df=n - 1))97 extra[f"mean_error_naive_h{h}"] = float(en.mean())98 extra[f"mean_error_model_h{h}"] = float(em.mean())99100# exclusion ledger (observation level, last training window)101n_nonmiss = int(s.notna().sum())102last_window = s.loc[:ORIGINS[-1]]103n_analysed = int(last_window.notna().sum())104_, nreg_last = fit_ar_diff(last_window)105ex = pd.DataFrame([106 dict(step=1, rule="drop rows with missing UNRATE value", n_before=n_raw, n_after=n_nonmiss,107 n_removed=n_raw - n_nonmiss, rationale="empty value field (2025-10-01) cannot be used for estimation"),108 dict(step=2, rule="drop observations dated after last origin 2026-06-01", n_before=n_nonmiss,109 n_after=n_analysed, n_removed=n_nonmiss - n_analysed,110 rationale="information rule: last training window uses only data <= last origin (later months are evaluation targets only)"),111])112ex.to_csv(os.path.join(OUT, "exclusions.csv"), index=False)113114missing_dates = [d.strftime("%Y-%m-%d") for d in s.index[s.isna()]]115res.update({116 "feasibility": "complete", "feasibility_reason": None,117 "data_issues": [118 {"issue": f"missing observation(s) at {missing_dates}",119 "handling": ("target pairs with missing actual kept in ledger with empty actual and excluded from error "120 "statistics; origin 2025-10 uses latest earlier observation (naive) / iterates one extra step "121 "(model); for later origins the missing level is linearly interpolated from data <= origin for "122 "the AR state only; regression rows touching it are dropped")},123 {"issue": "COVID-19 outliers (2020-04 = 14.8%, large swings 2020-04..2020-08) in the estimation sample",124 "handling": "kept (no outlier treatment, decided a priori); they influence OLS coefficients - see report"},125 {"issue": "final-vintage (revised) data, not real-time vintages",126 "handling": "results describe pseudo-out-of-sample accuracy on revised data; no real-time claim"},127 ],128 "n_input_rows": {"UNRATE.csv": int(n_raw)},129 "n_analysed": n_analysed,130 "n_regression_rows_last_window": int(nreg_last),131 "random_seed": None,132})133res.update(extra)134beta_last, _ = fit_ar_diff(last_window)135res["ar_coefficients_last_window"] = {"const": float(beta_last[0]), "phi1": float(beta_last[1]), "phi2": float(beta_last[2])}136with open(os.path.join(OUT, "results.json"), "w") as fh:137 json.dump(res, fh, indent=2)138for h in HORIZONS:139 print(f"h={h}: n={res[f'n_evaluated_h{h}']} RMSE naive={res[f'rmse_naive_h{h}']:.4f} model={res[f'rmse_model_h{h}']:.4f} "140 f"rel_impr={res[f'rel_improvement_rmse_h{h}']:.4f} DM={res[f'dm_stat_h{h}']:.3f} p={res[f'dm_pvalue_two_sided_h{h}']:.3f}")141