1"""F3-conf-04: seasonal AR (iterated) vs seasonal-naive forecasts of Canadian quarterly CPI inflation.2Run from workspace root: ~/wp29/venv/bin/python out/analysis.py3Reads data/ only; writes out/results.json, out/forecasts.csv, out/exclusions.csv. Deterministic (no randomness)."""4import json, os5import numpy as np, pandas as pd6from scipy import stats78OUT = "out"9os.makedirs(OUT, exist_ok=True)10q = pd.read_csv("data/cpi_quarterly.csv", parse_dates=["quarter"])11m = pd.read_csv("data/cpi_monthly.csv")12n_q_raw, n_m_raw = len(q), len(m)13q = q.sort_values("quarter").reset_index(drop=True)1415data_issues = []16# consistency checks against monthly file17mm = m.copy(); mm["q"] = pd.to_datetime(mm["month"]).dt.to_period("Q")18agg = mm.groupby("q")["cpi"].agg(["mean", "size"])19assert (agg["size"] == 3).all() and len(agg) == n_q_raw20assert np.abs(agg["mean"].values - q["cpi_avg"].values).max() < 1e-521recomp = 400 * np.diff(np.log(q["cpi_avg"].values))22assert np.abs(recomp - q["infl_annualised"].values[1:]).max() < 1e-423assert (q["quarter"].diff().dropna().dt.days.between(89, 92)).all()24data_issues.append({"issue": "infl_annualised is empty for 1995Q1 (no previous quarter)",25 "handling": "row excluded from all estimation samples (exclusion ledger step 1)"})2627excl = []28df = q.dropna(subset=["infl_annualised"]).reset_index(drop=True)29excl.append(("1", "drop rows with missing infl_annualised", n_q_raw, len(df),30 "first difference undefined for 1995Q1"))31y_all = df["infl_annualised"].values32dates = df["quarter"]33qtr = dates.dt.quarter.values # 1..434pos = {d: i for i, d in enumerate(dates)}3536PMAX = 437H = [1, 2, 4]38origins = pd.date_range("2015-10-01", "2025-07-01", freq="QS")39assert len(origins) == 40404142def design(y, qq, p, start):43 """rows t = start..len(y)-1; regressors: 4 seasonal dummies + p lags."""44 T = len(y)45 X = [np.eye(4)[qq[start:T] - 1]]46 for i in range(1, p + 1):47 X.append(y[start - i:T - i][:, None])48 return np.hstack(X), y[start:T]495051def fit_select(y, qq):52 best = None53 for p in range(PMAX + 1):54 X, Y = design(y, qq, p, PMAX) # common sample for AIC comparison55 b, *_ = np.linalg.lstsq(X, Y, rcond=None)56 rss = float(((Y - X @ b) ** 2).sum()); n = len(Y)57 aic = n * np.log(rss / n) + 2 * X.shape[1]58 if best is None or aic < best[0]:59 best = (aic, p, b, n)60 return best[1], best[2], best[3]616263rows = []; sel = []64for t in origins:65 i = pos[t]66 y = y_all[: i + 1]; qq = qtr[: i + 1] # information set: quarters <= t67 p, b, n_est = fit_select(y, qq)68 sel.append(p)69 hist = list(y); fc = {}70 for h in range(1, max(H) + 1): # iterated multi-step71 qt = ((qtr[i] - 1 + h) % 4) + 172 f = b[qt - 1] + sum(b[4 + k] * hist[-1 - k] for k in range(p))73 hist.append(f); fc[h] = f74 for h in H:75 tgt = t + pd.DateOffset(months=3 * h)76 j = pos[tgt]; assert j == i + h77 actual = y_all[j]78 sn = y_all[j - 4]; assert j - 4 <= i79 rows.append((t.date().isoformat(), tgt.date().isoformat(), h, sn, actual, "snaive"))80 rows.append((t.date().isoformat(), tgt.date().isoformat(), h, fc[h], actual, "sar_aic"))81 last_n = n_est; last_p = p8283fc_df = pd.DataFrame(rows, columns=["origin_date", "target_date", "horizon", "forecast", "actual", "model"])84fc_df.to_csv(f"{OUT}/forecasts.csv", index=False, float_format="%.6f")8586n_last_avail = pos[origins[-1]] + 187excl.append(("2", "last estimation sample: quarters after last origin 2025-07-01 not used",88 len(df), n_last_avail, "information-availability rule (quarters <= origin only)"))89excl.append(("3", f"drop first {PMAX} quarters to form lags (common AIC sample, PMAX={PMAX})",90 n_last_avail, last_n, "lagged regressors undefined at start of sample"))91pd.DataFrame([dict(step=s, rule=r, n_before=a, n_after=b_, n_removed=a - b_, rationale=ra)92 for s, r, a, b_, ra in excl]).to_csv(f"{OUT}/exclusions.csv", index=False)9394res = {"task_id": "F3-conf-04", "model_name": "sar_aic",95 "model_description": ("Iterated AR(p) with four quarterly seasonal intercepts for annualised quarterly inflation, "96 "p in {0..4} chosen by AIC at each origin on a common sample (first 4 obs dropped), "97 "OLS on an expanding window from 1995Q2 to the origin, multi-step forecasts by iteration."),98 "n_origins": len(origins), "first_origin": origins[0].date().isoformat(),99 "last_origin": origins[-1].date().isoformat()}100err = {}101for mod, key in [("snaive", "snaive"), ("sar_aic", "model")]:102 for h in H:103 s = fc_df[(fc_df.model == mod) & (fc_df.horizon == h)].sort_values("origin_date")104 e = (s.actual - s.forecast).values; err[(key, h)] = e105 res[f"rmse_{key}_h{h}"] = float(np.sqrt(np.mean(e ** 2)))106 res[f"mae_{key}_h{h}"] = float(np.mean(np.abs(e)))107for h in H:108 res[f"rel_improvement_rmse_h{h}"] = 1 - res[f"rmse_model_h{h}"] / res[f"rmse_snaive_h{h}"]109110# DM test at h=4: Bartlett (Newey-West) long-run variance, lag 3 (= h-1), no small-sample correction111d = err[("snaive", 4)] ** 2 - err[("model", 4)] ** 2112T = len(d); L = 3113dc = d - d.mean()114lrv = dc @ dc / T115for k in range(1, L + 1):116 lrv += 2 * (1 - k / (L + 1)) * (dc[k:] @ dc[:-k]) / T117se = float(np.sqrt(lrv / T))118dm = float(d.mean() / se)119res.update({"loss_diff_mean_h4": float(d.mean()), "loss_diff_se_h4": se, "hac_lags_h4": L,120 "hac_kernel_h4": "bartlett", "dm_stat_h4": dm,121 "dm_pvalue_h4": float(2 * stats.norm.sf(abs(dm))),122 "feasibility": "complete", "feasibility_reason": None, "data_issues": data_issues,123 "n_input_rows": {"cpi_quarterly.csv": n_q_raw, "cpi_monthly.csv": n_m_raw},124 "n_analysed": int(last_n), "random_seed": None})125# supplementary (not required): sensitivity of the DM SE to lag/kernel, and lag orders selected126sens = {}127for LL in [3, 4, 6]:128 v = dc @ dc / T129 for k in range(1, LL + 1):130 v += 2 * (1 - k / (LL + 1)) * (dc[k:] @ dc[:-k]) / T131 sens[f"bartlett_L{LL}"] = float(np.sqrt(v / T))132v = dc @ dc / T + 2 * sum((dc[k:] @ dc[:-k]) / T for k in range(1, 4))133sens["rectangular_L3"] = float(np.sqrt(v / T)) if v > 0 else None134res["supp_dm_se_sensitivity_h4"] = sens135res["supp_hln_pvalue_h4"] = float(2 * stats.t.sf(abs(dm) * np.sqrt((T + 1 - 2 * 4 + 4 * 3 / T) / T), T - 1))136res["supp_selected_p_counts"] = {str(k): int(v) for k, v in zip(*np.unique(sel, return_counts=True))}137res["supp_selected_p_last_origin"] = int(last_p)138with open(f"{OUT}/results.json", "w") as f:139 json.dump(res, f, indent=2)140print(json.dumps({k: v for k, v in res.items() if k != "data_issues"}, indent=1))141