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 np6import pandas as pd78DATA = "data/UNRATE.csv"9OUT = "out"10MODEL = "AR2_diff"11P = 212HORIZONS = [1, 3]13ORIGINS = pd.date_range("2023-07-01", "2026-06-01", freq="MS")1415raw = pd.read_csv(DATA, dtype={"observation_date": str, "UNRATE": str}, keep_default_na=False)16n_rows = len(raw)17df = pd.DataFrame({"date": pd.to_datetime(raw["observation_date"], format="%Y-%m-%d"),18 "y": pd.to_numeric(raw["UNRATE"].replace("", np.nan), errors="coerce")})19issues = []20miss = df.loc[df.y.isna(), "date"].dt.strftime("%Y-%m-%d").tolist()21assert df.date.is_unique and df.date.is_monotonic_increasing22full = pd.date_range(df.date.min(), df.date.max(), freq="MS")23assert len(full) == len(df), "calendar gaps present"24s = df.set_index("date")["y"]25if miss:26 issues.append({"issue": f"Missing value(s) (empty field) at {', '.join(miss)}",27 "handling": "Excluded from estimation (any regression row whose dependent or lagged difference touches the missing month is dropped); "28 "naive uses latest earlier observation; when missing is the latest month at an origin the AR model iterates forward from the last observed month "29 "(extra steps); when it is interior to the window it is linearly interpolated from neighbours dated <= origin only to build lag inputs; "30 "(origin,horizon) pairs whose target is missing are kept in the ledger with empty actual and excluded from error statistics."})31covid = s["2020-03-01":"2020-06-01"]32issues.append({"issue": "COVID-19 outlier: unemployment jumped from %.1f (2020-03) to %.1f (2020-04), extreme first differences in training sample" % (s["2020-03-01"], s["2020-04-01"]),33 "handling": "Retained in the primary specification (pre-specified, no outlier treatment); a sensitivity variant excluding 2020-03..2020-12 differences from estimation is reported in report.md only."})343536def fit_ar_diff(y, p=P, exclude=None):37 """OLS AR(p) with intercept on dy, using rows with no missing values. y: level Series up to origin."""38 d = y.diff()39 X = pd.concat({f"l{k}": d.shift(k) for k in range(1, p + 1)}, axis=1)40 data = pd.concat([d.rename("dy"), X], axis=1)41 if exclude is not None:42 data = data[~((data.index >= exclude[0]) & (data.index <= exclude[1]))]43 data = data.dropna()44 A = np.column_stack([np.ones(len(data)), data[[f"l{k}" for k in range(1, p + 1)]].values])45 beta, *_ = np.linalg.lstsq(A, data["dy"].values, rcond=None)46 return beta, len(data)474849def forecast_ar_diff(y, t, horizons, exclude=None):50 y = y[:t]51 beta, nobs = fit_ar_diff(y, exclude=exclude)52 last_obs = y.last_valid_index()53 yy = y[:last_obs].interpolate(method="linear", limit_area="inside") # interior gaps only, data <= t54 lev = list(yy.values)55 dif = list(np.diff(yy.values))56 steps_extra = (t.year - last_obs.year) * 12 + (t.month - last_obs.month)57 out = {}58 H = max(horizons) + steps_extra59 for step in range(1, H + 1):60 dnew = beta[0] + sum(beta[k] * dif[-k] for k in range(1, P + 1))61 dif.append(dnew)62 lev.append(lev[-1] + dnew)63 hh = step - steps_extra64 if hh in horizons:65 out[hh] = lev[-1]66 return out, nobs676869def run(exclude=None):70 rows, nobs_last = [], None71 for t in ORIGINS:72 yt = s[:t]73 naive = yt.dropna().iloc[-1]74 fc, nobs = forecast_ar_diff(s, t, HORIZONS, exclude)75 nobs_last = nobs76 for h in HORIZONS:77 tgt = t + pd.DateOffset(months=h)78 act = s.get(tgt, np.nan)79 for m, f in (("naive", naive), (MODEL, fc[h])):80 rows.append({"origin_date": t.strftime("%Y-%m-%d"), "target_date": tgt.strftime("%Y-%m-%d"),81 "horizon": h, "forecast": round(float(f), 6),82 "actual": act, "model": m})83 return pd.DataFrame(rows), nobs_last848586def stats(led):87 r = {}88 for h in HORIZONS:89 for m, tag in (("naive", "naive"), (MODEL, "model")):90 sub = led[(led.horizon == h) & (led.model == m) & led.actual.notna()]91 e = sub.actual - sub.forecast92 r[f"rmse_{tag}_h{h}"] = float(np.sqrt(np.mean(e ** 2)))93 r[f"mae_{tag}_h{h}"] = float(np.mean(np.abs(e)))94 r[f"n_evaluated_h{h}"] = int(len(sub))95 r[f"rel_improvement_rmse_h{h}"] = 1 - r[f"rmse_model_h{h}"] / r[f"rmse_naive_h{h}"]96 return r979899if __name__ == "__main__":100 led, nobs_last = run()101 st = stats(led)102 os.makedirs(OUT, exist_ok=True)103 led.to_csv(os.path.join(OUT, "forecasts.csv"), index=False, float_format="%.6g")104 res = {105 "task_id": "F3-dev-01",106 "model_name": MODEL,107 "model_description": "AR(2) with intercept on monthly first differences of UNRATE, estimated by OLS on an expanding window starting with the first available difference (1948-04 as dependent observation) up to the origin, iterated recursively and cumulated to levels for h-step forecasts.",108 "n_origins": len(ORIGINS),109 "first_origin": ORIGINS[0].strftime("%Y-%m-%d"),110 "last_origin": ORIGINS[-1].strftime("%Y-%m-%d"),111 }112 for k in ["n_evaluated_h1", "n_evaluated_h3", "rmse_naive_h1", "mae_naive_h1", "rmse_naive_h3", "mae_naive_h3",113 "rmse_model_h1", "mae_model_h1", "rmse_model_h3", "mae_model_h3",114 "rel_improvement_rmse_h1", "rel_improvement_rmse_h3"]:115 res[k] = st[k]116 res.update({"feasibility": "complete", "feasibility_reason": None, "data_issues": issues,117 "n_input_rows": {"UNRATE.csv": n_rows}, "n_analysed": int(nobs_last), "random_seed": None})118 with open(os.path.join(OUT, "results.json"), "w") as f:119 json.dump(res, f, indent=2)120 print(json.dumps(res, indent=2))121122 # ---------------- supplementary (report only) ----------------123 from scipy import stats as sps124 sup = {}125 beta_last, _ = fit_ar_diff(s[:ORIGINS[-1]])126 beta_first, _ = fit_ar_diff(s[:ORIGINS[0]])127 sup["coef_first_origin"] = [float(b) for b in beta_first]128 sup["coef_last_origin"] = [float(b) for b in beta_last]129130 def dm(led, h):131 a = led[(led.horizon == h) & (led.model == "naive") & led.actual.notna()].reset_index(drop=True)132 b = led[(led.horizon == h) & (led.model == MODEL) & led.actual.notna()].reset_index(drop=True)133 d = (a.actual - b.forecast) ** 2 - (a.actual - a.forecast) ** 2 # model loss - naive loss134 n = len(d); dbar = d.mean(); dc = d - dbar135 lrv = np.sum(dc ** 2) / n136 for k in range(1, h): # rectangular HAC with h-1 lags137 lrv += 2 * np.sum(dc[k:].values * dc[:-k].values) / n138 stat = dbar / np.sqrt(lrv / n)139 hln = np.sqrt((n + 1 - 2 * h + h * (h - 1) / n) / n) # Harvey-Leybourne-Newbold correction140 stat *= hln141 p = 2 * sps.t.sf(abs(stat), df=n - 1)142 return {"mean_loss_diff_model_minus_naive": float(dbar), "dm_hln_stat": float(stat), "p_value_two_sided": float(p), "n": int(n)}143 sup["dm_h1"] = dm(led, 1); sup["dm_h3"] = dm(led, 3)144145 # bootstrap-free descriptive: share of origins where model abs error < naive abs error146 for h in HORIZONS:147 a = led[(led.horizon == h) & (led.model == "naive") & led.actual.notna()].reset_index(drop=True)148 b = led[(led.horizon == h) & (led.model == MODEL) & led.actual.notna()].reset_index(drop=True)149 sup[f"share_model_better_abs_h{h}"] = float(np.mean(np.abs(b.actual - b.forecast) < np.abs(a.actual - a.forecast)))150 sup[f"share_ties_abs_h{h}"] = float(np.mean(np.isclose(np.abs(b.actual - b.forecast), np.abs(a.actual - a.forecast))))151152 led2, _ = run(exclude=(pd.Timestamp("2020-03-01"), pd.Timestamp("2020-12-01")))153 st2 = stats(led2)154 sup["sensitivity_exclude_covid_2020"] = {k: st2[k] for k in st2 if k.startswith(("rmse_model", "mae_model", "rel_"))}155 sup["sensitivity_exclude_covid_coef_last_origin"] = [float(b) for b in fit_ar_diff(s[:ORIGINS[-1]], exclude=(pd.Timestamp("2020-03-01"), pd.Timestamp("2020-12-01")))[0]]156 with open(os.path.join(OUT, "supplementary.json"), "w") as f:157 json.dump(sup, f, indent=2)158 print(json.dumps(sup, indent=2))159