1"""F3-conf-04: seasonal-naive vs AR(p)+seasonal dummies, expanding window, iterated forecasts.2Run from workspace root: ~/wp29/venv/bin/python out/analysis.py"""3import json, os4import numpy as np, pandas as pd5from scipy import stats67H = [1, 2, 4]; PMAX = 4; MODEL = "ar_seasdum"8q = pd.read_csv("data/cpi_quarterly.csv", parse_dates=["quarter"])9m = pd.read_csv("data/cpi_monthly.csv")10n_rows = {"cpi_quarterly.csv": int(len(q)), "cpi_monthly.csv": int(len(m))}11s = q.dropna(subset=["infl_annualised"]).reset_index(drop=True)12y = s["infl_annualised"].to_numpy(); dates = s["quarter"]; qtr = dates.dt.quarter.to_numpy()13origins = pd.date_range("2015-10-01", "2025-07-01", freq="QS")14assert len(origins) == 4015pos = {d: i for i, d in enumerate(dates)}1617def design(yv, qv, p, start):18 T = len(yv); rows = range(start, T)19 X = [[1.0 * (qv[t] == k) for k in (1, 2, 3, 4)] + [yv[t - j] for j in range(1, p + 1)] for t in rows]20 return np.array(X), yv[start:]2122def fit_forecast(i_end):23 yv, qv = y[: i_end + 1], qtr[: i_end + 1] # only quarters <= origin24 best = None25 for p in range(PMAX + 1): # AIC on common sample (start = PMAX)26 X, Y = design(yv, qv, p, PMAX)27 b, *_ = np.linalg.lstsq(X, Y, rcond=None)28 rss = float(((Y - X @ b) ** 2).sum()); n = len(Y)29 aic = n * np.log(rss / n) + 2 * X.shape[1]30 if best is None or aic < best[0]: best = (aic, p, b, n)31 _, p, b, n = best32 hist = list(yv); qs = list(qv); fc = {}33 for h in range(1, max(H) + 1): # iterated multi-step34 qn = (qs[-1] % 4) + 135 x = [1.0 * (qn == k) for k in (1, 2, 3, 4)] + [hist[-j] for j in range(1, p + 1)]36 f = float(np.dot(x, b)); hist.append(f); qs.append(qn); fc[h] = f37 return fc, p, n3839rows = []; plist = []; n_last = None40for o in origins:41 i = pos[o]; fc, p, n = fit_forecast(i); plist.append(p); n_last = n42 for h in H:43 tgt = i + h; act = float(y[tgt])44 assert dates[tgt] == o + pd.DateOffset(months=3 * h)45 assert tgt - 4 <= i46 rows.append(dict(origin_date=o.date().isoformat(), target_date=dates[tgt].date().isoformat(), horizon=h,47 forecast=float(y[tgt - 4]), actual=act, model="snaive"))48 rows.append(dict(origin_date=o.date().isoformat(), target_date=dates[tgt].date().isoformat(), horizon=h,49 forecast=fc[h], actual=act, model=MODEL))50F = pd.DataFrame(rows); F["err"] = F["actual"] - F["forecast"]51os.makedirs("out", exist_ok=True)52F.drop(columns="err").to_csv("out/forecasts.csv", index=False, float_format="%.6f")5354R = {"task_id": "F3-conf-04", "model_name": MODEL,55 "model_description": ("OLS AR(p) on annualised quarterly inflation with four quarterly seasonal intercepts, "56 "p in {0..4} chosen by AIC at each origin on a common sample, expanding window from 1995Q2 "57 "to the origin, iterated multi-step forecasts."),58 "n_origins": 40, "first_origin": origins[0].date().isoformat(), "last_origin": origins[-1].date().isoformat()}59for h in H:60 for mod, key in (("snaive", "snaive"), (MODEL, "model")):61 e = F[(F.horizon == h) & (F.model == mod)]["err"].to_numpy()62 R[f"rmse_{key}_h{h}"] = float(np.sqrt(np.mean(e ** 2))); R[f"mae_{key}_h{h}"] = float(np.mean(np.abs(e)))63for h in H:64 R[f"rel_improvement_rmse_h{h}"] = 1 - R[f"rmse_model_h{h}"] / R[f"rmse_snaive_h{h}"]6566def hac_var_mean(d, L, kernel="bartlett"):67 d = d - d.mean(); T = len(d); lrv = d @ d / T68 for j in range(1, L + 1):69 w = 1 - j / (L + 1) if kernel == "bartlett" else 1.070 lrv += 2 * w * (d[j:] @ d[:-j]) / T71 return lrv / T7273e_s = F[(F.horizon == 4) & (F.model == "snaive")].sort_values("origin_date")["err"].to_numpy()74e_m = F[(F.horizon == 4) & (F.model == MODEL)].sort_values("origin_date")["err"].to_numpy()75d = e_s ** 2 - e_m ** 2; T = len(d); L = 376dbar = float(d.mean()); se = float(np.sqrt(hac_var_mean(d, L)))77dm = dbar / se78R.update({"loss_diff_mean_h4": dbar, "loss_diff_se_h4": se, "hac_lags_h4": L, "hac_kernel_h4": "bartlett",79 "dm_stat_h4": dm, "dm_pvalue_h4": float(2 * stats.norm.sf(abs(dm)))})80# sensitivity (not primary)81sens = {}82for LL in (4, 6, 8):83 sens[f"bartlett_L{LL}_se"] = float(np.sqrt(hac_var_mean(d, LL)))84 sens[f"bartlett_L{LL}_pvalue"] = float(2 * stats.norm.sf(abs(dbar / np.sqrt(hac_var_mean(d, LL)))))85vr = hac_var_mean(d, 3, "rectangular"); sens["rectangular_L3_se"] = float(np.sqrt(vr)) if vr > 0 else None86hln = np.sqrt((T + 1 - 2 * 4 + 4 * 3 / T) / T) # Harvey-Leybourne-Newbold, h=487sens["hln_stat_bartlett_L3"] = float(dm * hln)88sens["hln_pvalue_t39"] = float(2 * stats.t.sf(abs(dm * hln), T - 1))89sens["iid_se_L0"] = float(np.sqrt(hac_var_mean(d, 0)))90sens["d_autocorr_lag1_to_4"] = [float(np.corrcoef(d[j:], d[:-j])[0, 1]) for j in range(1, 5)]91sens["n_origins_model_better_sqerr_h4"] = int((d > 0).sum())92R["sensitivity_h4"] = sens93R["selected_ar_order_counts"] = {str(k): int(v) for k, v in pd.Series(plist).value_counts().sort_index().items()}94R.update({"feasibility": "complete", "feasibility_reason": None,95 "data_issues": [{"issue": "infl_annualised missing for 1995Q1 (no prior quarter)", "handling": "dropped; series starts 1995Q2"},96 {"issue": "infl_annualised reproduced from cpi_avg only to ~3e-6 due to 6-dp rounding of cpi_avg", "handling": "used infl_annualised as given"},97 {"issue": "COVID-2020 and 2021-2023 inflation-surge outliers within evaluation window", "handling": "kept (genuine data); noted as threat to DM stationarity assumption"}],98 "n_input_rows": n_rows, "n_analysed": int(n_last), "random_seed": None})99with open("out/results.json", "w") as f: json.dump(R, f, indent=2)100print(json.dumps({k: R[k] for k in R if k.startswith(("rmse", "mae", "rel", "dm", "loss"))}, indent=1))101