1"""F3-conf-04: multi-horizon forecasts of Canadian quarterly CPI inflation (NSA, annualised).2Model: AR(p) + quarter-of-year dummies, OLS, iterated multi-step, p in 0..4 chosen by AIC3on each expanding window (data dated <= origin only). Baseline: seasonal naive.4Run from workspace root: ~/wp29/venv/bin/python out/analysis.py5"""6import json7import numpy as np8import pandas as pd9from scipy import stats1011DATA = "data"12OUT = "out"13MODEL = "ar_sdum_aic"14PMAX = 415H = [1, 2, 4]1617q = pd.read_csv(f"{DATA}/cpi_quarterly.csv")18m = pd.read_csv(f"{DATA}/cpi_monthly.csv")19n_rows = {"cpi_quarterly.csv": int(len(q)), "cpi_monthly.csv": int(len(m))}20q["quarter"] = pd.to_datetime(q["quarter"])21q = q.sort_values("quarter").reset_index(drop=True)2223# ---- data checks ----24issues = []25assert not q["quarter"].duplicated().any()26assert (q["quarter"].dt.to_period("Q").diff().iloc[1:].apply(lambda x: x.n) == 1).all()27m["q"] = pd.to_datetime(m["month"]).dt.to_period("Q")28agg = m.groupby("q")["cpi"].agg(["mean", "count"])29assert (agg["count"] == 3).all() and len(agg) == len(q)30assert np.max(np.abs(agg["mean"].values - q["cpi_avg"].values)) < 1e-531recomp = 400 * np.log(q["cpi_avg"]).diff()32assert np.nanmax(np.abs(recomp - q["infl_annualised"])) < 1e-433issues.append({"issue": "infl_annualised is missing for 1995Q1 (first difference undefined)",34 "handling": "1995Q1 excluded from all estimation samples; consistency of cpi_avg with the monthly file and of infl_annualised with 400*dln(cpi_avg) verified (max abs diff < 1e-5)"})35issues.append({"issue": "Evaluation window contains the 2020 COVID shock and the 2021-2023 inflation surge (outliers relative to the pre-2015 estimation sample)",36 "handling": "No winsorising or exclusion; errors computed on the data as given, noted as a limitation"})3738y_all = q["infl_annualised"].values39dates = q["quarter"].values40qtr = q["quarter"].dt.quarter.values414243def fit_ar_sdum(y, qq, p, start):44 """OLS of y_t on 4 quarter dummies and p lags, t = start..T-1 (common sample)."""45 T = len(y)46 rows = range(start, T)47 X = []48 for t in rows:49 d = [1.0 if qq[t] == k else 0.0 for k in (1, 2, 3, 4)]50 X.append(d + [y[t - j] for j in range(1, p + 1)])51 X = np.array(X)52 Y = y[start:]53 beta, *_ = np.linalg.lstsq(X, Y, rcond=None)54 resid = Y - X @ beta55 n = len(Y)56 k = X.shape[1]57 aic = n * np.log(resid @ resid / n) + 2 * k58 return beta, aic, n596061def forecast_iter(y, qq, beta, p, hmax):62 path = list(y)63 qpath = list(qq)64 out = []65 for h in range(1, hmax + 1):66 qn = (qpath[-1] % 4) + 167 d = [1.0 if qn == k else 0.0 for k in (1, 2, 3, 4)]68 x = np.array(d + [path[-j] for j in range(1, p + 1)])69 f = float(x @ beta)70 out.append(f)71 path.append(f)72 qpath.append(qn)73 return out747576origins = pd.date_range("2015-10-01", "2025-07-01", freq="QS")77assert len(origins) == 4078rows = []79chosen_p = []80n_last = None81for o in origins:82 io = int(np.where(dates == np.datetime64(o))[0][0])83 y = y_all[1:io + 1] # drop 1995Q1 (NaN); data dated <= origin only84 qq = qtr[1:io + 1]85 start = PMAX # common sample across p for AIC comparability86 best = None87 for p in range(0, PMAX + 1):88 beta, aic, n = fit_ar_sdum(y, qq, p, start)89 if best is None or aic < best[1]:90 best = (p, aic, beta)91 p, _, beta = best92 beta, _, n = fit_ar_sdum(y, qq, p, p) # re-estimate on the longest sample for chosen p93 n_last = n94 chosen_p.append(p)95 fc = forecast_iter(y, qq, beta, p, max(H))96 for h in H:97 it = io + h98 act = float(y_all[it])99 tdate = pd.Timestamp(dates[it])100 rows.append(dict(origin_date=o.date().isoformat(), target_date=tdate.date().isoformat(),101 horizon=h, forecast=fc[h - 1], actual=act, model=MODEL))102 rows.append(dict(origin_date=o.date().isoformat(), target_date=tdate.date().isoformat(),103 horizon=h, forecast=float(y_all[it - 4]), actual=act, model="snaive"))104105fc = pd.DataFrame(rows)106fc.to_csv(f"{OUT}/forecasts.csv", index=False, float_format="%.6f")107fc["e"] = fc["actual"] - fc["forecast"]108109res = {"task_id": "F3-conf-04", "model_name": MODEL,110 "model_description": ("Iterated multi-step AR(p) with four quarter-of-year intercept dummies, "111 "estimated by OLS on an expanding window (1995Q2 to origin), lag order p in 0..4 "112 "re-selected by AIC at every origin on a common sample, then re-estimated on the full window."),113 "n_origins": 40, "first_origin": origins[0].date().isoformat(),114 "last_origin": origins[-1].date().isoformat()}115for h in H:116 for mod, tag in (("snaive", "snaive"), (MODEL, "model")):117 e = fc[(fc.horizon == h) & (fc.model == mod)]["e"].values118 assert len(e) == 40119 res[f"rmse_{tag}_h{h}"] = float(np.sqrt(np.mean(e ** 2)))120 res[f"mae_{tag}_h{h}"] = float(np.mean(np.abs(e)))121for h in H:122 res[f"rel_improvement_rmse_h{h}"] = 1 - res[f"rmse_model_h{h}"] / res[f"rmse_snaive_h{h}"]123124125def hac_se(d, L, kernel):126 d = np.asarray(d) - np.mean(d)127 n = len(d)128 lrv = d @ d / n129 for j in range(1, L + 1):130 w = 1 - j / (L + 1) if kernel == "bartlett" else 1.0131 lrv += 2 * w * (d[j:] @ d[:-j]) / n132 return np.sqrt(lrv / n) if lrv > 0 else np.nan133134135def dm_d(h):136 es = fc[(fc.horizon == h) & (fc.model == "snaive")].sort_values("origin_date")["e"].values137 em = fc[(fc.horizon == h) & (fc.model == MODEL)].sort_values("origin_date")["e"].values138 return es ** 2 - em ** 2139140141d4 = dm_d(4)142L, K = 4, "bartlett" # Bartlett with L=4: non-zero weights on autocovariances 1..4 (>= h-1 = 3)143se = hac_se(d4, L, K)144res["loss_diff_mean_h4"] = float(np.mean(d4))145res["loss_diff_se_h4"] = float(se)146res["hac_lags_h4"] = L147res["hac_kernel_h4"] = K148res["dm_stat_h4"] = float(np.mean(d4) / se)149# two-sided p-value from Student t with n-1 df (Harvey-Leybourne-Newbold style small-sample reference)150res["dm_pvalue_h4"] = float(2 * stats.t.sf(abs(res["dm_stat_h4"]), df=len(d4) - 1))151152# sensitivity (not required keys)153sens = {}154for (LL, KK) in [(3, "bartlett"), (4, "bartlett"), (6, "bartlett"), (8, "bartlett"), (3, "rectangular")]:155 s = hac_se(d4, LL, KK)156 st = np.mean(d4) / s157 sens[f"{KK}_L{LL}"] = {"se": float(s), "stat": float(st),158 "p_t": float(2 * stats.t.sf(abs(st), len(d4) - 1)),159 "p_normal": float(2 * stats.norm.sf(abs(st)))}160# HLN correction factor for h=4 with rectangular L=3161n = len(d4); h = 4162hln = np.sqrt((n + 1 - 2 * h + h * (h - 1) / n) / n)163sens["hln_rectangular_L3"] = {"stat": float(sens["rectangular_L3"]["stat"] * hln),164 "p_t": float(2 * stats.t.sf(abs(sens["rectangular_L3"]["stat"] * hln), n - 1))}165for hh in (1, 2):166 dd = dm_d(hh)167 s = hac_se(dd, max(hh, 1), "bartlett")168 sens[f"dm_h{hh}_bartlett_L{max(hh,1)}"] = {"mean": float(np.mean(dd)), "se": float(s),169 "stat": float(np.mean(dd) / s),170 "p_t": float(2 * stats.t.sf(abs(np.mean(dd) / s), n - 1))}171# excluding 2020-2023 targets at h=4172mask = ~fc[(fc.horizon == 4) & (fc.model == "snaive")].sort_values("origin_date")["target_date"].str[:4].isin(["2020", "2021", "2022", "2023"]).values173sens["h4_excl_targets_2020_2023"] = {"n": int(mask.sum()), "loss_diff_mean": float(np.mean(d4[mask]))}174175res["feasibility"] = "complete"176res["feasibility_reason"] = None177res["data_issues"] = issues178res["n_input_rows"] = n_rows179res["n_analysed"] = int(n_last)180res["random_seed"] = None181res["supplementary"] = {"chosen_p_by_origin": dict(zip([o.date().isoformat() for o in origins], chosen_p)),182 "dm_sensitivity": sens,183 "dm_pvalue_reference": "Student t with n-1 = 39 df"}184with open(f"{OUT}/results.json", "w") as f:185 json.dump(res, f, indent=2)186print(json.dumps(res, indent=2))187