1"""F3-conf-04: multi-horizon forecasts of Canadian quarterly CPI inflation vs seasonal naive.2Run from workspace root: ~/wp29/venv/bin/python out/analysis.py3Reads only data/; writes out/results.json and out/forecasts.csv. Deterministic (no randomness)."""4import json, os5import numpy as np6import pandas as pd7from scipy import stats89ROOT = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))10DATA = os.path.join(ROOT, "data"); OUT = os.path.join(ROOT, "out")11HORIZONS = [1, 2, 4]; PMAX = 8; HAC_LAGS = 3; MODEL = "ar_seas_aic"1213q = pd.read_csv(os.path.join(DATA, "cpi_quarterly.csv"), parse_dates=["quarter"])14m = pd.read_csv(os.path.join(DATA, "cpi_monthly.csv"))15n_q_raw, n_m_raw = len(q), len(m)16q = q.dropna(subset=["infl_annualised"]).reset_index(drop=True) # drops 1995Q1 only17y_all = q.set_index("quarter")["infl_annualised"]18dates = list(y_all.index)19origins = pd.date_range("2015-10-01", "2025-07-01", freq="QS-JAN")20assert len(origins) == 402122def design(y, qtr, p, start):23 """Rows t=start..T-1: [Q1..Q4 dummies, y_{t-1..t-p}]; target y_t."""24 T = len(y)25 X = [np.eye(4)[qtr[start:T]]]26 for j in range(1, p + 1):27 X.append(y[start - j:T - j, None])28 return np.hstack(X), y[start:T]2930def fit_select(y, qtr):31 """AIC lag selection p in 0..PMAX on common sample (first PMAX obs held out), then OLS refit on max sample."""32 best = None33 for p in range(PMAX + 1):34 X, Y = design(y, qtr, p, PMAX)35 b, *_ = np.linalg.lstsq(X, Y, rcond=None)36 ssr = np.sum((Y - X @ b) ** 2); n = len(Y)37 aic = n * np.log(ssr / n) + 2 * X.shape[1]38 if best is None or aic < best[0]:39 best = (aic, p)40 p = best[1]41 X, Y = design(y, qtr, p, p)42 b, *_ = np.linalg.lstsq(X, Y, rcond=None)43 return p, b, len(Y)4445def iterate(y, qtr_last, p, b, H):46 hist = list(y); out = []47 for h in range(1, H + 1):48 qq = (qtr_last + h) % 449 f = b[qq] + sum(b[4 + j - 1] * hist[-j] for j in range(1, p + 1))50 hist.append(f); out.append(f)51 return out5253rows = []; lag_choice = {}; n_last = None54for t in origins:55 idx = dates.index(t)56 hist = y_all.iloc[: idx + 1] # quarters <= t only57 assert hist.index.max() == t58 y = hist.values; qtr = ((hist.index.month - 1) // 3).values59 p, b, n_est = fit_select(y, qtr)60 lag_choice[str(t.date())] = p; n_last = n_est61 fc = iterate(y, qtr[-1], p, b, max(HORIZONS))62 for h in HORIZONS:63 tgt = t + pd.DateOffset(months=3 * h)64 actual = float(y_all.loc[tgt])65 sn_date = tgt - pd.DateOffset(months=12)66 assert sn_date <= t67 sn = float(y_all.loc[sn_date])68 rows.append(dict(origin_date=str(t.date()), target_date=str(tgt.date()), horizon=h, forecast=sn, actual=actual, model="snaive"))69 rows.append(dict(origin_date=str(t.date()), target_date=str(tgt.date()), horizon=h, forecast=float(fc[h - 1]), actual=actual, model=MODEL))70fc_df = pd.DataFrame(rows).sort_values(["model", "horizon", "origin_date"]).reset_index(drop=True)71fc_df.to_csv(os.path.join(OUT, "forecasts.csv"), index=False, float_format="%.6f")72fc_df = pd.read_csv(os.path.join(OUT, "forecasts.csv")) # metrics computed from the written file7374def hac_var_mean(d, L, kernel):75 d = np.asarray(d) - np.mean(d); n = len(d)76 g = [np.sum(d[j:] * d[:n - j]) / n for j in range(L + 1)]77 w = [1 - j / (L + 1) if kernel == "bartlett" else 1.0 for j in range(L + 1)]78 lrv = g[0] + 2 * sum(w[j] * g[j] for j in range(1, L + 1))79 return lrv / n8081res = {"task_id": "F3-conf-04", "model_name": MODEL,82 "model_description": ("Iterated (recursive) AR(p) for annualised quarterly inflation with four quarter-of-year intercepts, "83 "p in 0..8 selected by AIC on a common sample at each origin, estimated by OLS on an expanding window of all quarters <= origin (from 1995Q2)."),84 "n_origins": len(origins), "first_origin": str(origins[0].date()), "last_origin": str(origins[-1].date())}85err = {}86for mod, key in [("snaive", "snaive"), (MODEL, "model")]:87 for h in HORIZONS:88 s = fc_df[(fc_df.model == mod) & (fc_df.horizon == h)].sort_values("origin_date")89 e = (s.actual - s.forecast).values; err[(key, h)] = e90 res[f"rmse_{key}_h{h}"] = float(np.sqrt(np.mean(e ** 2)))91 res[f"mae_{key}_h{h}"] = float(np.mean(np.abs(e)))92for h in HORIZONS:93 res[f"rel_improvement_rmse_h{h}"] = 1 - res[f"rmse_model_h{h}"] / res[f"rmse_snaive_h{h}"]94d = err[("snaive", 4)] ** 2 - err[("model", 4)] ** 295n = len(d); dbar = float(np.mean(d))96se = float(np.sqrt(hac_var_mean(d, HAC_LAGS, "bartlett")))97dm = dbar / se98res.update(loss_diff_mean_h4=dbar, loss_diff_se_h4=se, hac_lags_h4=HAC_LAGS, hac_kernel_h4="bartlett",99 dm_stat_h4=dm, dm_pvalue_h4=float(2 * stats.norm.sf(abs(dm))))100# sensitivity (supplementary keys)101sens = {}102for L, k in [(4, "bartlett"), (6, "bartlett"), (3, "rectangular")]:103 v = hac_var_mean(d, L, k)104 if v > 0:105 s_ = np.sqrt(v); sens[f"{k}_L{L}"] = {"se": float(s_), "dm": float(dbar / s_), "p_normal": float(2 * stats.norm.sf(abs(dbar / s_)))}106 else:107 sens[f"{k}_L{L}"] = {"se": None, "dm": None, "p_normal": None, "note": "non-positive long-run variance"}108hln = np.sqrt((n + 1 - 2 * 4 + 4 * 3 / n) / n) # Harvey-Leybourne-Newbold small-sample factor, h=4109sens["hln_bartlett_L3_t39"] = {"dm": float(dm * hln), "p_t": float(2 * stats.t.sf(abs(dm * hln), n - 1))}110res["sensitivity_h4"] = sens111res["n_loss_diff_positive_h4"] = int(np.sum(d > 0))112res["selected_lags_by_origin"] = lag_choice113res.update(feasibility="complete", feasibility_reason=None,114 data_issues=[{"issue": "infl_annualised missing for 1995Q1 (no preceding quarter)", "handling": "row dropped; estimation starts 1995Q2"},115 {"issue": "extreme value 12.15 in 2022Q2 (inflation surge)", "handling": "genuine observation, kept"}],116 n_input_rows={"cpi_quarterly.csv": n_q_raw, "cpi_monthly.csv": n_m_raw},117 n_analysed=int(n_last), random_seed=None)118with open(os.path.join(OUT, "results.json"), "w") as f:119 json.dump(res, f, indent=2)120print(json.dumps({k: v for k, v in res.items() if k != "selected_lags_by_origin"}, indent=2))121print("selected lags:", sorted(set(lag_choice.values())))122