1"""F3-conf-04: multi-horizon forecasts of Canadian quarterly CPI inflation (NSA, annualised).2Model: AR(p) with quarterly seasonal dummies, OLS, iterated multi-step, p in 0..4 chosen by AIC3on each expanding window (common sample). Baseline: seasonal naive. DM test at h=4 with Bartlett HAC."""4import json, numpy as np, pandas as pd5from scipy import stats6from pathlib import Path78ROOT = Path(__file__).resolve().parent.parent9OUT = ROOT / "out"10q = pd.read_csv(ROOT / "data/cpi_quarterly.csv", parse_dates=["quarter"])11m = pd.read_csv(ROOT / "data/cpi_monthly.csv", parse_dates=["month"])12issues = []13# --- data checks14m = m.sort_values("month").reset_index(drop=True)15if m["cpi"].isna().any(): issues.append({"issue": "missing monthly CPI", "handling": "none needed for quarterly input"})16q = q.sort_values("quarter").reset_index(drop=True)17exp_q = pd.date_range(q.quarter.iloc[0], q.quarter.iloc[-1], freq="QS")18assert len(exp_q) == len(q) and (exp_q == q.quarter).all(), "quarterly gaps"19mq = m.assign(quarter=m.month.dt.to_period("Q").dt.start_time).groupby("quarter")["cpi"].agg(["mean", "size"])20mq = mq[mq["size"] == 3]["mean"]21chk = q.set_index("quarter")["cpi_avg"] - mq.reindex(q.quarter).values22infl_chk = 400 * np.log(q.cpi_avg).diff()23maxdev_avg = float(np.nanmax(np.abs(chk.values))); maxdev_inf = float(np.nanmax(np.abs(infl_chk - q.infl_annualised)))24issues.append({"issue": "infl_annualised is empty for 1995Q1 (no previous quarter)",25 "handling": "1995Q1 dropped from the inflation series; estimation starts 1995Q2"})26if maxdev_avg > 1e-5 or maxdev_inf > 1e-4:27 issues.append({"issue": f"quarterly file deviates from monthly recomputation (max cpi_avg dev {maxdev_avg:.2e}, infl dev {maxdev_inf:.2e})",28 "handling": "quarterly file used as given (it is the designated input)"})2930y = q.set_index("quarter")["infl_annualised"].dropna()31dates = list(y.index)32H = [1, 2, 4]; PMAX = 4; MODEL = "ar_seas_aic"33origins = pd.date_range("2015-10-01", "2025-07-01", freq="QS")34assert len(origins) == 403536def design(yv, qtr, p, start):37 # rows t=start..len-1 ; regressors: 4 seasonal dummies + p lags38 X = [np.eye(4)[qtr[start:]]]39 for l in range(1, p + 1):40 X.append(yv[start - l:len(yv) - l][:, None])41 return np.hstack(X), yv[start:]4243def fit_forecast(ys, hmax):44 yv = ys.values; qtr = np.array([d.month // 3 for d in ys.index]) # 0..345 best = None46 for p in range(PMAX + 1):47 X, Y = design(yv, qtr, p, PMAX) # common sample for AIC48 b, *_ = np.linalg.lstsq(X, Y, rcond=None)49 n = len(Y); s2 = np.sum((Y - X @ b) ** 2) / n50 aic = n * np.log(s2) + 2 * X.shape[1]51 if best is None or aic < best[0]: best = (aic, p)52 p = best[1]53 X, Y = design(yv, qtr, p, p) # re-estimate on max sample for chosen p54 b, *_ = np.linalg.lstsq(X, Y, rcond=None)55 hist = list(yv); lastq = qtr[-1]; fc = []56 for h in range(1, hmax + 1):57 qq = (lastq + h) % 458 f = b[qq] + sum(b[4 + l - 1] * hist[-l] for l in range(1, p + 1))59 hist.append(f); fc.append(f)60 return fc, p, len(Y)6162rows = []; plist = []63for o in origins:64 ys = y[y.index <= o]65 fc, p, nobs = fit_forecast(ys, max(H)); plist.append(p)66 for h in H:67 tgt = o + pd.DateOffset(months=3 * h)68 act = float(y.loc[tgt])69 sn = float(y.loc[tgt - pd.DateOffset(months=12)])70 assert tgt - pd.DateOffset(months=12) <= o71 rows.append(dict(origin_date=o.date(), target_date=tgt.date(), horizon=h, forecast=sn, actual=act, model="snaive"))72 rows.append(dict(origin_date=o.date(), target_date=tgt.date(), horizon=h, forecast=fc[h - 1], actual=act, model=MODEL))73n_last = nobs74fc = pd.DataFrame(rows); fc.to_csv(OUT / "forecasts.csv", index=False, float_format="%.6f")75fc["e"] = fc.actual - fc.forecast7677res = {"task_id": "F3-conf-04", "model_name": MODEL,78 "model_description": ("Iterated multi-step AR(p) for annualised NSA quarterly inflation with four quarterly seasonal intercepts, "79 "p in {0,...,4} chosen by AIC on a common sample at each origin, estimated by OLS on an expanding window from 1995Q2 to the origin."),80 "n_origins": 40, "first_origin": str(origins[0].date()), "last_origin": str(origins[-1].date())}81for h in H:82 for mod, tag in [("snaive", "snaive"), (MODEL, "model")]:83 e = fc[(fc.horizon == h) & (fc.model == mod)].e.values84 res[f"rmse_{tag}_h{h}"] = float(np.sqrt(np.mean(e ** 2)))85 res[f"mae_{tag}_h{h}"] = float(np.mean(np.abs(e)))86for h in H:87 res[f"rel_improvement_rmse_h{h}"] = 1 - res[f"rmse_model_h{h}"] / res[f"rmse_snaive_h{h}"]8889def hac_se(d, L, kernel="bartlett"):90 d = np.asarray(d); n = len(d); u = d - d.mean()91 lrv = u @ u / n92 for j in range(1, L + 1):93 w = 1 - j / (L + 1) if kernel == "bartlett" else 1.094 lrv += 2 * w * (u[j:] @ u[:-j]) / n95 return np.sqrt(lrv / n), lrv9697def dloss(h):98 es = fc[(fc.horizon == h) & (fc.model == "snaive")].sort_values("origin_date").e.values99 em = fc[(fc.horizon == h) & (fc.model == MODEL)].sort_values("origin_date").e.values100 return es ** 2 - em ** 2101d4 = dloss(4); L = 4102se, _ = hac_se(d4, L)103dm = d4.mean() / se104res.update({"loss_diff_mean_h4": float(d4.mean()), "loss_diff_se_h4": float(se), "hac_lags_h4": L,105 "hac_kernel_h4": "bartlett", "dm_stat_h4": float(dm),106 "dm_pvalue_h4": float(2 * stats.norm.sf(abs(dm)))})107# sensitivity (reported in report, not required keys)108sens = {}109for k, LL in [("bartlett", 3), ("bartlett", 4), ("bartlett", 6), ("rectangular", 3)]:110 s, lrv = hac_se(d4, LL, k)111 t = d4.mean() / s if lrv > 0 else float("nan")112 n = len(d4); hh = 4113 hln = np.sqrt((n + 1 - 2 * hh + hh * (hh - 1) / n) / n)114 sens[f"{k}_L{LL}"] = {"se": float(s), "stat": float(t), "p_normal": float(2 * stats.norm.sf(abs(t))),115 "p_t39": float(2 * stats.t.sf(abs(t), n - 1)),116 "hln_stat": float(t * hln), "hln_p_t39": float(2 * stats.t.sf(abs(t * hln), n - 1))}117for h in [1, 2]:118 d = dloss(h); s, _ = hac_se(d, max(h - 1, 1))119 sens[f"h{h}_bartlett_L{max(h-1,1)}"] = {"mean": float(d.mean()), "se": float(s), "stat": float(d.mean() / s),120 "p_normal": float(2 * stats.norm.sf(abs(d.mean() / s)))}121res.update({"feasibility": "complete", "feasibility_reason": None, "data_issues": issues,122 "n_input_rows": {"cpi_quarterly.csv": int(len(q)), "cpi_monthly.csv": int(len(m))},123 "n_analysed": int(n_last), "random_seed": None})124json.dump(res, open(OUT / "results.json", "w"), indent=2)125json.dump({"sensitivity_h4_and_h12": sens, "ar_order_by_origin": dict(zip([str(o.date()) for o in origins], plist)),126 "check_max_dev_cpi_avg": maxdev_avg, "check_max_dev_infl": maxdev_inf},127 open(OUT / "diagnostics.json", "w"), indent=2)128print(json.dumps(res, indent=2))129