1"""F3-conf-04: multi-horizon forecasts of Canadian quarterly NSA CPI inflation.2Run from workspace root: ~/wp29/venv/bin/python out/analysis.py3Model: SAR-AIC = OLS AR(p) with 4 quarterly intercepts (seasonal dummies), p in {0..4}4chosen by AIC at each origin on a common sample, expanding window, iterated multi-step forecasts.5No randomness is used.6"""7import json, os8import numpy as np, pandas as pd9from scipy import stats1011ROOT = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))12DATA, OUT = os.path.join(ROOT, "data"), os.path.join(ROOT, "out")13MODEL = "sar_aic"14PMAX, HS, HAC_LAG, KERNEL = 4, (1, 2, 4), 4, "bartlett"1516q_raw = pd.read_csv(os.path.join(DATA, "cpi_quarterly.csv"), parse_dates=["quarter"])17m_raw = pd.read_csv(os.path.join(DATA, "cpi_monthly.csv"))18n_q, n_m = len(q_raw), len(m_raw)1920excl = []21q = q_raw.dropna(subset=["infl_annualised"]).sort_values("quarter").reset_index(drop=True)22excl.append(dict(step=1, rule="drop rows with missing infl_annualised (1995Q1, no previous quarter)",23 n_before=n_q, n_after=len(q), n_removed=n_q - len(q),24 rationale="target undefined for the first quarter"))25pd.DataFrame(excl).to_csv(os.path.join(OUT, "exclusions.csv"), index=False)2627y = q["infl_annualised"].to_numpy()28dates = q["quarter"]29qtr = dates.dt.quarter.to_numpy()30idx = {d: i for i, d in enumerate(dates)}31origins = pd.date_range("2015-10-01", "2025-07-01", freq="QS-JAN")32assert len(origins) == 40333435def fit_sar(yw, qw):36 """Fit AR(p)+seasonal dummies for p=0..PMAX on common sample; return best by AIC."""37 T = len(yw)38 best = None39 for p in range(PMAX + 1):40 rows = np.arange(PMAX, T)41 D = (qw[rows][:, None] == np.arange(1, 5)[None, :]).astype(float)42 L = np.column_stack([yw[rows - k] for k in range(1, p + 1)]) if p else np.empty((len(rows), 0))43 X = np.hstack([D, L]); Y = yw[rows]44 b, *_ = np.linalg.lstsq(X, Y, rcond=None)45 rss = np.sum((Y - X @ b) ** 2); n = len(Y)46 aic = n * np.log(rss / n) + 2 * X.shape[1]47 if best is None or aic < best[0]:48 best = (aic, p, b, n)49 return best505152def iterate(yw, qw, p, b, hmax):53 hist = list(yw); qs = list(qw); out = []54 for h in range(1, hmax + 1):55 qn = qs[-1] % 4 + 156 f = b[qn - 1] + sum(b[4 + k - 1] * hist[-k] for k in range(1, p + 1))57 out.append(f); hist.append(f); qs.append(qn)58 return out596061rows, chosen_p, n_last = [], [], None62for t in origins:63 i = idx[t]64 yw, qw = y[: i + 1], qtr[: i + 1] # information set: quarters <= t65 assert dates.iloc[i] == t and dates.iloc[: i + 1].max() <= t66 aic, p, b, n = fit_sar(yw, qw)67 chosen_p.append(p); n_last = n68 fc = iterate(yw, qw, p, b, max(HS))69 for h in HS:70 tgt = i + h71 assert tgt - 4 <= i72 common = dict(origin_date=t.strftime("%Y-%m-%d"), target_date=dates.iloc[tgt].strftime("%Y-%m-%d"),73 horizon=h, actual=y[tgt])74 rows.append({**common, "forecast": y[tgt - 4], "model": "snaive"})75 rows.append({**common, "forecast": fc[h - 1], "model": MODEL})7677F = pd.DataFrame(rows)[["origin_date", "target_date", "horizon", "forecast", "actual", "model"]]78F = F.sort_values(["model", "horizon", "origin_date"]).reset_index(drop=True)79F.to_csv(os.path.join(OUT, "forecasts.csv"), index=False, float_format="%.6f")80F["e"] = F["actual"] - F["forecast"]818283def lrv(d, L, kernel):84 d = d - d.mean(); T = len(d)85 g = lambda j: np.sum(d[j:] * d[: T - j]) / T86 w = (lambda j: 1 - j / (L + 1)) if kernel == "bartlett" else (lambda j: 1.0)87 return g(0) + 2 * sum(w(j) * g(j) for j in range(1, L + 1))888990res = dict(task_id="F3-conf-04", model_name=MODEL,91 model_description=("OLS regression of quarterly annualised inflation on four quarterly seasonal intercepts "92 "plus p own lags, p in {0,...,4} selected by AIC at each origin on a common sample "93 "(first 4 obs held out for lags), expanding window from 1995Q2 to the origin, "94 "iterated (recursive) multi-step forecasts for h=1,2,4."),95 n_origins=len(origins), first_origin=origins[0].strftime("%Y-%m-%d"),96 last_origin=origins[-1].strftime("%Y-%m-%d"))97for h in HS:98 for mname, key in (("snaive", "snaive"), (MODEL, "model")):99 e = F[(F.model == mname) & (F.horizon == h)]["e"].to_numpy()100 assert len(e) == 40101 res[f"rmse_{key}_h{h}"] = float(np.sqrt(np.mean(e ** 2)))102 res[f"mae_{key}_h{h}"] = float(np.mean(np.abs(e)))103for h in HS:104 res[f"rel_improvement_rmse_h{h}"] = 1 - res[f"rmse_model_h{h}"] / res[f"rmse_snaive_h{h}"]105106e_s = F[(F.model == "snaive") & (F.horizon == 4)].sort_values("origin_date")["e"].to_numpy()107e_m = F[(F.model == MODEL) & (F.horizon == 4)].sort_values("origin_date")["e"].to_numpy()108d = e_s ** 2 - e_m ** 2; T = len(d)109se = float(np.sqrt(lrv(d, HAC_LAG, KERNEL) / T))110dm = float(d.mean() / se)111res.update(loss_diff_mean_h4=float(d.mean()), loss_diff_se_h4=se, hac_lags_h4=HAC_LAG, hac_kernel_h4=KERNEL,112 dm_stat_h4=dm, dm_pvalue_h4=float(2 * stats.norm.sf(abs(dm))))113# sensitivity (extra keys)114v_rect = lrv(d, 3, "rectangular")115res["sens_rect_lag3_lrv_positive"] = bool(v_rect > 0)116res["sens_rect_lag3_dm_stat_h4"] = float(d.mean() / np.sqrt(v_rect / T)) if v_rect > 0 else None117h = 4118hln = np.sqrt((T + 1 - 2 * h + h * (h - 1) / T) / T)119se_r3 = np.sqrt(v_rect / T) if v_rect > 0 else se120res["sens_hln_dm_stat_h4"] = float(hln * d.mean() / se_r3)121res["sens_hln_dm_pvalue_h4_t39"] = float(2 * stats.t.sf(abs(res["sens_hln_dm_stat_h4"]), T - 1))122for L in (3, 6, 8):123 res[f"sens_bartlett_lag{L}_dm_stat_h4"] = float(d.mean() / np.sqrt(lrv(d, L, "bartlett") / T))124res["selected_p_counts"] = {str(k): int(v) for k, v in zip(*np.unique(chosen_p, return_counts=True))}125res.update(feasibility="complete", feasibility_reason=None,126 data_issues=[{"issue": "infl_annualised missing for 1995Q1 (no previous quarter)",127 "handling": "row dropped (exclusions.csv step 1)"},128 {"issue": "recomputing 400*dln(cpi_avg) differs from infl_annualised by up to ~3e-6 due to 6-decimal rounding of cpi_avg",129 "handling": "used infl_annualised as given (immaterial)"}],130 n_input_rows={"cpi_quarterly.csv": n_q, "cpi_monthly.csv": n_m},131 n_analysed=int(n_last), random_seed=None)132with open(os.path.join(OUT, "results.json"), "w") as fh:133 json.dump(res, fh, indent=2)134print(json.dumps(res, indent=2))135