1"""F3-conf-04: AR + seasonal dummies vs seasonal naive, Canadian quarterly CPI inflation."""2import json, os3import numpy as np, pandas as pd4from scipy import stats56ROOT = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))7D = os.path.join(ROOT, "data"); O = os.path.join(ROOT, "out")8q = pd.read_csv(os.path.join(D, "cpi_quarterly.csv"), parse_dates=["quarter"])9m = pd.read_csv(os.path.join(D, "cpi_monthly.csv"), parse_dates=["month"])10issues = []11# consistency checks12mq = m.set_index("month")["cpi"].resample("QS").agg(["mean", "count"])13mq = mq[mq["count"] == 3]14chk = q.set_index("quarter")["cpi_avg"].reindex(mq.index)15if (chk - mq["mean"]).abs().max() > 1e-5:16 issues.append({"issue": "quarterly cpi_avg differs from monthly averages", "handling": "used quarterly file as given"})17rec = 400 * np.log(q["cpi_avg"]).diff()18if (rec - q["infl_annualised"]).abs().max() > 1e-5:19 issues.append({"issue": "infl_annualised inconsistent with cpi_avg", "handling": "used as given"})20issues.append({"issue": "infl_annualised missing for 1995Q1 (first difference)", "handling": "row dropped from target series; 126 inflation observations used"})21y_full = q.dropna(subset=["infl_annualised"]).set_index("quarter")["infl_annualised"]22assert y_full.index.to_series().diff().dropna().dt.days.between(89, 92).all()2324PMAX, HS = 4, [1, 2, 4]25MODEL = "ar_sd_aic"26origins = pd.date_range("2015-10-01", "2025-07-01", freq="QS")27assert len(origins) == 402829def design(y, p, start):30 """rows t=start..n-1: const, 3 seasonal dummies, lags 1..p"""31 n = len(y); qtr = y.index.quarter.values32 X = [np.ones(n - start)]33 for s in (2, 3, 4):34 X.append((qtr[start:] == s).astype(float))35 for l in range(1, p + 1):36 X.append(y.values[start - l:n - l])37 return np.column_stack(X), y.values[start:]3839def fit_forecast(y):40 best = None41 for p in range(PMAX + 1):42 X, Y = design(y, p, PMAX) # common sample for AIC comparability43 b, *_ = np.linalg.lstsq(X, Y, rcond=None)44 rss = ((Y - X @ b) ** 2).sum(); n = len(Y)45 aic = n * np.log(rss / n) + 2 * X.shape[1]46 if best is None or aic < best[0]:47 best = (aic, p, b, n)48 _, p, b, n = best49 # iterate forecasts50 hist = list(y.values); qs = list(y.index)51 fc = {}52 for h in range(1, max(HS) + 1):53 d = qs[-1] + pd.DateOffset(months=3 * h)54 x = [1.0] + [float(d.quarter == s) for s in (2, 3, 4)] + [hist[-l] for l in range(1, p + 1)]55 f = float(np.dot(b, x)); hist.append(f); fc[h] = (d, f)56 return fc, p, n5758rows, plist = [], []59for t in origins:60 y = y_full[y_full.index <= t]61 fc, p, nobs = fit_forecast(y); plist.append(p)62 for h in HS:63 d, f = fc[h]64 act = y_full.get(d, np.nan)65 sn = y_full[d - pd.DateOffset(months=12)]66 for name, val in (("snaive", sn), (MODEL, f)):67 rows.append({"origin_date": t.strftime("%Y-%m-%d"), "target_date": d.strftime("%Y-%m-%d"),68 "horizon": h, "forecast": round(float(val), 6), "actual": float(act), "model": name})69fc = pd.DataFrame(rows)70assert fc["actual"].notna().all() and len(fc) == 24071fc.to_csv(os.path.join(O, "forecasts.csv"), index=False)7273res = {"task_id": "F3-conf-04", "model_name": MODEL,74 "model_description": ("OLS AR(p) on annualised quarterly inflation with intercept and three quarterly seasonal dummies, "75 "p in {0..4} selected by AIC at each origin on a common sample, expanding window from 1995Q2 "76 "(first lag-4-complete obs 1996Q2) to the origin, multi-step forecasts iterated."),77 "n_origins": 40, "first_origin": origins[0].strftime("%Y-%m-%d"), "last_origin": origins[-1].strftime("%Y-%m-%d")}78err = {}79for h in HS:80 for mod, key in (("snaive", "snaive"), (MODEL, "model")):81 s = fc[(fc.horizon == h) & (fc.model == mod)].sort_values("origin_date")82 e = (s.actual - s.forecast).values; err[(h, key)] = e83 res[f"rmse_{key}_h{h}"] = float(np.sqrt(np.mean(e ** 2)))84 res[f"mae_{key}_h{h}"] = float(np.mean(np.abs(e)))85for h in HS:86 res[f"rel_improvement_rmse_h{h}"] = 1 - res[f"rmse_model_h{h}"] / res[f"rmse_snaive_h{h}"]8788def hac_se(d, L, kernel="bartlett"):89 n = len(d); u = d - d.mean()90 lrv = u @ u / n91 for j in range(1, L + 1):92 w = 1 - j / (L + 1) if kernel == "bartlett" else 1.093 lrv += 2 * w * (u[j:] @ u[:-j]) / n94 return np.sqrt(lrv / n), lrv9596d = err[(4, "snaive")] ** 2 - err[(4, "model")] ** 297L = 398se, _ = hac_se(d, L)99dm = d.mean() / se100res.update({"loss_diff_mean_h4": float(d.mean()), "loss_diff_se_h4": float(se), "hac_lags_h4": L,101 "hac_kernel_h4": "bartlett", "dm_stat_h4": float(dm),102 "dm_pvalue_h4": float(2 * stats.norm.sf(abs(dm))),103 "feasibility": "complete", "feasibility_reason": None, "data_issues": issues,104 "n_input_rows": {"cpi_quarterly.csv": int(len(q)), "cpi_monthly.csv": int(len(m))},105 "n_analysed": int(nobs), "random_seed": None})106json.dump(res, open(os.path.join(O, "results.json"), "w"), indent=2)107108# sensitivity (printed only, reported in report.md)109n = len(d)110hln = np.sqrt((n + 1 - 2 * 4 + 4 * 3 / n) / n)111print("selected p per origin:", plist)112print(json.dumps({k: v for k, v in res.items() if isinstance(v, (int, float))}, indent=1))113print("HLN-adjusted stat, t(39) p:", dm * hln, 2 * stats.t.sf(abs(dm * hln), n - 1))114for L2, k in ((3, "rectangular"), (4, "bartlett"), (6, "bartlett"), (8, "bartlett")):115 s2, lrv = hac_se(d, L2, k); print(k, L2, "se", s2, "lrv", lrv, "stat", d.mean() / s2 if lrv > 0 else None)116for h in (1, 2):117 dd = err[(h, "snaive")] ** 2 - err[(h, "model")] ** 2118 s2, _ = hac_se(dd, max(h - 1, 1) + 2); print("h", h, "dbar", dd.mean(), "se", s2, "stat", dd.mean() / s2)119# excluding 2021-2023 targets120s = fc[fc.horizon == 4].pivot(index="origin_date", columns="model", values=["forecast", "actual"])121tgt = fc[(fc.horizon == 4) & (fc.model == "snaive")].sort_values("origin_date").target_date.values122mask = ~((tgt >= "2021-01-01") & (tgt <= "2023-10-01"))123print("h4 excl 2021-2023 targets: n", mask.sum(), "rmse sn", np.sqrt(np.mean(err[(4,'snaive')][mask]**2)),124 "rmse model", np.sqrt(np.mean(err[(4,'model')][mask]**2)))125