1"""F3-conf-04: multi-horizon forecasts of Canadian quarterly CPI inflation.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, hashlib, os5import numpy as np, pandas as pd6from scipy import stats78DATA = "data"; OUT = "out"9HORIZONS = [1, 2, 4]; PMAX = 410MODEL = "sar_aic"1112def sha(p): return hashlib.sha256(open(p, "rb").read()).hexdigest()1314q_raw = pd.read_csv(f"{DATA}/cpi_quarterly.csv", parse_dates=["quarter"])15m_raw = pd.read_csv(f"{DATA}/cpi_monthly.csv", parse_dates=["month"])16n_rows = {"cpi_quarterly.csv": int(len(q_raw)), "cpi_monthly.csv": int(len(m_raw))}17q = q_raw.dropna(subset=["infl_annualised"]).reset_index(drop=True) # drops 1995Q1 only18assert (q.quarter.diff().dropna().dt.days.between(89, 93)).all(), "gap in quarters"19y_all = q.infl_annualised.to_numpy(); dates = q.quarter20qtr_all = dates.dt.quarter.to_numpy()2122def design(y, qtr, p, start):23 """rows t=start..T-1: const, Q2,Q3,Q4 dummies, y_{t-1..t-p}"""24 T = len(y); rows = range(start, T)25 X = [[1.0] + [float(qtr[t] == k) for k in (2, 3, 4)] + [y[t - j] for j in range(1, p + 1)] for t in rows]26 return np.array(X), y[start:]2728def fit(y, qtr, p, start):29 X, Y = design(y, qtr, p, start)30 b, *_ = np.linalg.lstsq(X, Y, rcond=None)31 rss = float(((Y - X @ b) ** 2).sum()); return b, rss, len(Y)3233def forecast_origin(y, qtr):34 """y, qtr: data up to and including the origin only. AIC choice of p on common sample (start=PMAX),35 refit chosen p on its full available sample, iterate forecasts to h=4."""36 aic = {}37 for p in range(PMAX + 1):38 _, rss, n = fit(y, qtr, p, PMAX)39 aic[p] = n * np.log(rss / n) + 2 * (4 + p)40 p = min(aic, key=aic.get)41 b, _, n = fit(y, qtr, p, p)42 hist = list(y); qn = int(qtr[-1]); fc = {}43 for h in range(1, max(HORIZONS) + 1):44 qn = qn % 4 + 145 x = [1.0] + [float(qn == k) for k in (2, 3, 4)] + [hist[-j] for j in range(1, p + 1)]46 f = float(np.dot(b, x)); hist.append(f); fc[h] = f47 return fc, p, n4849origins = pd.date_range("2015-10-01", "2025-07-01", freq="QS")50assert len(origins) == 4051idx = {d: i for i, d in enumerate(dates)}52rows, chosen_p = [], []53for o in origins:54 i = idx[o]55 y, qt = y_all[: i + 1], qtr_all[: i + 1] # information set: quarters <= origin56 assert dates.iloc[i] == o and dates.iloc[: i + 1].max() <= o57 fc, p, n_last = forecast_origin(y, qt); chosen_p.append(p)58 for h in HORIZONS:59 tgt = i + h60 sn = y_all[tgt - 4]; assert tgt - 4 <= i61 for mname, f in (("snaive", sn), (MODEL, fc[h])):62 rows.append(dict(origin_date=o.strftime("%Y-%m-%d"), target_date=dates.iloc[tgt].strftime("%Y-%m-%d"),63 horizon=h, forecast=round(float(f), 6), actual=float(y_all[tgt]), model=mname))64fdf = pd.DataFrame(rows)65fdf.to_csv(f"{OUT}/forecasts.csv", index=False)6667def err(m, h):68 d = fdf[(fdf.model == m) & (fdf.horizon == h)].sort_values("origin_date")69 return (d.actual - d.forecast).to_numpy()7071def lrv(d, L, kernel):72 u = d - d.mean(); T = len(u); s = u @ u / T73 for j in range(1, L + 1):74 w = 1.0 if kernel == "rectangular" else 1 - j / (L + 1)75 s += 2 * w * (u[j:] @ u[:-j]) / T76 return s7778res = {"task_id": "F3-conf-04", "model_name": MODEL,79 "model_description": ("Iterated (recursive) seasonal AR(p): OLS of infl on a constant, three quarterly dummies and p own lags, "80 "p in 0..4 re-selected by AIC at each origin on a common sample, expanding window from 1995Q2 "81 "through the origin quarter; multi-step forecasts by iterating the fitted equation."),82 "n_origins": len(origins), "first_origin": origins[0].strftime("%Y-%m-%d"), "last_origin": origins[-1].strftime("%Y-%m-%d")}83for h in HORIZONS:84 for m, tag in (("snaive", "snaive"), (MODEL, "model")):85 e = err(m, h)86 res[f"rmse_{tag}_h{h}"] = float(np.sqrt(np.mean(e ** 2))); res[f"mae_{tag}_h{h}"] = float(np.mean(np.abs(e)))87for h in HORIZONS:88 res[f"rel_improvement_rmse_h{h}"] = 1 - res[f"rmse_model_h{h}"] / res[f"rmse_snaive_h{h}"]8990d = err("snaive", 4) ** 2 - err(MODEL, 4) ** 2; T = len(d)91L, kern = 3, "rectangular"92v = lrv(d, L, kern)93if v <= 0: # pre-specified fallback94 kern = "bartlett"; v = lrv(d, L, kern)95se = float(np.sqrt(v / T)); dm = float(d.mean() / se)96res.update(loss_diff_mean_h4=float(d.mean()), loss_diff_se_h4=se, hac_lags_h4=L, hac_kernel_h4=kern,97 dm_stat_h4=dm, dm_pvalue_h4=float(2 * stats.norm.sf(abs(dm))))98# sensitivity (not required keys)99sens = {}100for k, l in (("bartlett", 3), ("bartlett", 4), ("bartlett", 6), ("rectangular", 3)):101 vv = lrv(d, l, k); s = float(np.sqrt(vv / T)) if vv > 0 else None102 sens[f"{k}_L{l}"] = {"se": s, "dm": (float(d.mean() / s) if s else None),103 "p_normal": (float(2 * stats.norm.sf(abs(d.mean() / s))) if s else None)}104h = 4; hln = np.sqrt((T + 1 - 2 * h + h * (h - 1) / T) / T)105sens["HLN_rectangular_L3_t39"] = {"dm": dm * hln, "p_t": float(2 * stats.t.sf(abs(dm * hln), T - 1))}106res["sensitivity_h4"] = sens107res["chosen_p_by_origin"] = chosen_p108res.update(feasibility="complete", feasibility_reason=None,109 data_issues=[{"issue": "infl_annualised missing for 1995Q1 (no prior quarter)", "handling": "row dropped (exclusions.csv step 1)"},110 {"issue": "files contain quarters after the last origin (to 2026Q3)", "handling": "estimation truncated at each origin; later quarters used only as realised targets"}],111 n_input_rows=n_rows, n_analysed=int(n_last), random_seed=None,112 input_sha256={"cpi_quarterly.csv": sha(f"{DATA}/cpi_quarterly.csv"), "cpi_monthly.csv": sha(f"{DATA}/cpi_monthly.csv")})113json.dump(res, open(f"{OUT}/results.json", "w"), indent=2)114print(json.dumps({k: v for k, v in res.items() if k.startswith(("rmse", "mae", "rel", "loss", "dm", "hac"))}, indent=1))115