1"""F3-conf-04: AR(p)+seasonal dummies vs seasonal naive, Canadian quarterly CPI inflation.2Run from workspace root: ~/wp29/venv/bin/python out/analysis.py"""3import json4import numpy as np, pandas as pd5from scipy import stats67D = "data/"; O = "out/"8q = pd.read_csv(D + "cpi_quarterly.csv", parse_dates=["quarter"])9m = pd.read_csv(D + "cpi_monthly.csv", parse_dates=["month"])10n_rows = {"cpi_quarterly.csv": int(len(q)), "cpi_monthly.csv": int(len(m))}1112# ---- data checks ----13issues = []14assert q["quarter"].is_monotonic_increasing15gaps = (q["quarter"].dt.to_period("Q").astype(int).diff().dropna() != 1).sum()16assert gaps == 0, "non-contiguous quarters"17# recompute quarterly series from monthly file to verify derivation18m["q"] = m["month"].dt.to_period("Q")19mq = m.groupby("q")["cpi"].agg(["mean", "size"])20mq = mq[mq["size"] == 3]21rec = 400 * np.log(mq["mean"]).diff()22chk = pd.Series(rec.values, index=mq.index.to_timestamp()).reindex(q["quarter"])23maxdiff = float(np.nanmax(np.abs(chk.values - q["infl_annualised"].values)))24assert maxdiff < 1e-4, maxdiff25issues.append({"issue": "infl_annualised is empty for 1995Q1 (no previous quarter)",26 "handling": "first quarter dropped; series starts 1995Q2 (126 obs)"})27y = q.set_index("quarter")["infl_annualised"].dropna()28dates = list(y.index)29yv = y.values30qtr = np.array([d.quarter for d in dates])3132PMAX = 433H = [1, 2, 4]34origins = pd.date_range("2015-10-01", "2025-07-01", freq="QS")35assert len(origins) == 403637def design(yv, qtr, p, start):38 """rows i=start..end: y_i on const, Q2-Q4 dummies, y_{i-1..i-p}"""39 idx = np.arange(start, len(yv))40 X = [np.ones(len(idx))] + [(qtr[idx] == s).astype(float) for s in (2, 3, 4)]41 X += [yv[idx - k] for k in range(1, p + 1)]42 return np.column_stack(X), yv[idx]4344def fit_forecast(ys, qs, hmax):45 # AIC lag selection p in 0..PMAX on common sample (first PMAX obs held out)46 best = None47 for p in range(PMAX + 1):48 X, Y = design(ys, qs, p, PMAX)49 b, *_ = np.linalg.lstsq(X, Y, rcond=None)50 rss = float(((Y - X @ b) ** 2).sum()); n = len(Y)51 aic = n * np.log(rss / n) + 2 * X.shape[1]52 if best is None or aic < best[0]:53 best = (aic, p)54 p = best[1]55 X, Y = design(ys, qs, p, p) # re-estimate on all usable obs56 b, *_ = np.linalg.lstsq(X, Y, rcond=None)57 hist = list(ys); qcur = qs[-1]; fc = []58 for h in range(1, hmax + 1): # iterated forecasts59 qn = (qcur + h - 1) % 4 + 160 x = [1.0] + [float(qn == s) for s in (2, 3, 4)] + [hist[-k] for k in range(1, p + 1)]61 f = float(np.dot(b, x)); fc.append(f); hist.append(f)62 return fc, p, len(Y)6364rows = []; lags_used = []65for t in origins:66 i = dates.index(t)67 ys, qs = yv[: i + 1], qtr[: i + 1] # only quarters <= t68 fc, p, nobs = fit_forecast(ys, qs, max(H))69 lags_used.append(p)70 for h in H:71 tgt = dates[i + h]72 rows.append(dict(origin_date=t.date(), target_date=tgt.date(), horizon=h,73 forecast=fc[h - 1], actual=yv[i + h], model="ar_seas_aic"))74 rows.append(dict(origin_date=t.date(), target_date=tgt.date(), horizon=h,75 forecast=yv[i + h - 4], actual=yv[i + h], model="snaive"))76fcst = pd.DataFrame(rows).sort_values(["model", "horizon", "origin_date"])77fcst.to_csv(O + "forecasts.csv", index=False)78n_analysed = nobs # obs in last (2025Q3 origin) estimation regression7980def err(mod, h):81 s = fcst[(fcst.model == mod) & (fcst.horizon == h)].sort_values("origin_date")82 return (s.actual - s.forecast).values8384def hac_se(d, L, kernel="bartlett"):85 d = d - d.mean(); n = len(d)86 lrv = d @ d / n87 for k in range(1, L + 1):88 w = 1 - k / (L + 1) if kernel == "bartlett" else 1.089 lrv += 2 * w * (d[k:] @ d[:-k]) / n90 return np.sqrt(lrv / n), lrv9192res = {"task_id": "F3-conf-04", "model_name": "ar_seas_aic",93 "model_description": ("OLS autoregression of quarterly annualised inflation on a constant, "94 "three quarterly seasonal dummies and p own lags (p in 0..4 chosen by AIC at each origin on a "95 "common sample), estimated on an expanding window from 1995Q2 to the origin, with iterated "96 "multi-step forecasts (h=1,2,4) and known future seasonal dummies."),97 "n_origins": 40, "first_origin": str(origins[0].date()), "last_origin": str(origins[-1].date())}98for h in H:99 es, em = err("snaive", h), err("ar_seas_aic", h)100 res[f"rmse_snaive_h{h}"] = float(np.sqrt(np.mean(es ** 2)))101 res[f"mae_snaive_h{h}"] = float(np.mean(np.abs(es)))102 res[f"rmse_model_h{h}"] = float(np.sqrt(np.mean(em ** 2)))103 res[f"mae_model_h{h}"] = float(np.mean(np.abs(em)))104for h in H:105 res[f"rel_improvement_rmse_h{h}"] = 1 - res[f"rmse_model_h{h}"] / res[f"rmse_snaive_h{h}"]106107d = err("snaive", 4) ** 2 - err("ar_seas_aic", 4) ** 2108L = 3109se, _ = hac_se(d, L, "bartlett")110dm = d.mean() / se111res.update({"loss_diff_mean_h4": float(d.mean()), "loss_diff_se_h4": float(se),112 "hac_lags_h4": L, "hac_kernel_h4": "bartlett", "dm_stat_h4": float(dm),113 "dm_pvalue_h4": float(2 * stats.norm.sf(abs(dm))),114 "feasibility": "complete", "feasibility_reason": None, "data_issues": issues,115 "n_input_rows": n_rows, "n_analysed": int(n_analysed), "random_seed": None})116117# sensitivity (not required keys; stored separately)118sens = {"selected_lags_by_origin": lags_used,119 "derivation_check_max_abs_diff": maxdiff}120n = len(d)121for kern, lag in [("bartlett", 4), ("bartlett", 6), ("rectangular", 3)]:122 s_, lrv = hac_se(d, lag, kern)123 sens[f"{kern}_L{lag}"] = {"se": float(s_) if lrv > 0 else None,124 "dm": float(d.mean() / s_) if lrv > 0 else None,125 "p_norm": float(2 * stats.norm.sf(abs(d.mean() / s_))) if lrv > 0 else None}126hln = np.sqrt((n + 1 - 2 * 4 + 4 * 3 / n) / n)127sens["hln_bartlett_L3"] = {"dm": float(dm * hln), "p_t39": float(2 * stats.t.sf(abs(dm * hln), n - 1))}128for h in (1, 2):129 dh = err("snaive", h) ** 2 - err("ar_seas_aic", h) ** 2130 s_, _ = hac_se(dh, max(h - 1, 1))131 sens[f"dm_h{h}_bartlett_L{max(h-1,1)}"] = {"mean": float(dh.mean()), "se": float(s_),132 "dm": float(dh.mean() / s_), "p_norm": float(2 * stats.norm.sf(abs(dh.mean() / s_)))}133# excluding 2021-2023 inflation surge targets134fm = fcst[fcst.horizon == 4].copy()135fm["td"] = pd.to_datetime(fm.target_date)136keep = ~fm.td.between("2021-01-01", "2023-12-31")137for mod in ("snaive", "ar_seas_aic"):138 s = fm[keep & (fm.model == mod)]139 sens[f"rmse_h4_excl_2021_2023_{mod}"] = float(np.sqrt(((s.actual - s.forecast) ** 2).mean()))140sens["n_h4_excl_2021_2023"] = int((keep & (fm.model == "snaive")).sum())141142json.dump(res, open(O + "results.json", "w"), indent=2)143json.dump(sens, open(O + "sensitivity.json", "w"), indent=2)144print(json.dumps(res, indent=2)); print(json.dumps(sens, indent=2))145