1"""F3-conf-04: multi-horizon forecasts of Canadian quarterly CPI inflation.2Run from workspace root: ~/wp29/venv/bin/python out/analysis.py3Reads only data/, writes out/results.json and out/forecasts.csv. Deterministic (no randomness)."""4import json, os5import numpy as np, pandas as pd6from scipy import stats78ROOT = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))9D = os.path.join(ROOT, "data"); O = os.path.join(ROOT, "out")10HORIZONS = [1, 2, 4]; PMAX = 4; HAC_LAGS = 4; KERNEL = "bartlett"11MODEL = "sar_dummy_aic"1213q_raw = pd.read_csv(os.path.join(D, "cpi_quarterly.csv"), parse_dates=["quarter"])14m_raw = pd.read_csv(os.path.join(D, "cpi_monthly.csv"), parse_dates=["month"])15q = q_raw.dropna(subset=["infl_annualised"]).reset_index(drop=True) # exclusion step 116y_all = q["infl_annualised"].to_numpy(); dates = q["quarter"]; qtr = dates.dt.quarter.to_numpy()1718# data-issue check: quarterly series reproducible from monthly file19m = m_raw.copy(); m["q"] = m["month"].dt.to_period("Q")20agg = m.groupby("q")["cpi"].agg(["mean", "count"]); agg = agg[agg["count"] == 3]21rec = 400 * np.log(agg["mean"]).diff().dropna().to_numpy()22max_dev = float(np.max(np.abs(rec - y_all)))23data_issues = [{"issue": "infl_annualised missing for 1995Q1 (first difference undefined)",24 "handling": "row excluded from modelling (exclusions.csv step 1)"},25 {"issue": "large genuine shocks: 2008Q4 (-6.0), 2020Q2 COVID, 2022Q2 (+12.1)",26 "handling": "retained; no outlier treatment, since any such treatment would require post-origin judgement"}]27assert max_dev < 1e-5, "quarterly file not reproducible from monthly file"2829def design(y, s, p):30 """rows t=p..n-1: 4 quarter dummies + p own lags."""31 n = len(y); X = []32 for t in range(p, n):33 X.append([1.0 * (s[t] == k) for k in (1, 2, 3, 4)] + [y[t - i] for i in range(1, p + 1)])34 return np.array(X), y[p:]3536def fit_forecast(y, s, hmax):37 """Expanding-window fit on data <= origin; AIC selects p in 0..PMAX on common sample; iterated forecasts."""38 best = None39 for p in range(PMAX + 1):40 X, Y = design(y, s, p); X, Y = X[PMAX - p:], Y[PMAX - p:] # common sample41 b, *_ = np.linalg.lstsq(X, Y, rcond=None); r = Y - X @ b42 aic = len(Y) * np.log(r @ r / len(Y)) + 2 * X.shape[1]43 if best is None or aic < best[0]: best = (aic, p)44 p = best[1]45 X, Y = design(y, s, p); b, *_ = np.linalg.lstsq(X, Y, rcond=None)46 hist = list(y); fc = []; s_last = s[-1]47 for h in range(1, hmax + 1):48 sq = (s_last - 1 + h) % 4 + 149 x = [1.0 * (sq == k) for k in (1, 2, 3, 4)] + [hist[-i] for i in range(1, p + 1)]50 f = float(np.dot(x, b)); fc.append(f); hist.append(f)51 return fc, p, len(Y)5253origins = pd.date_range("2015-10-01", "2025-07-01", freq="QS")54rows = []; p_used = []; n_last = None55for o in origins:56 i0 = int(np.where(dates == o)[0][0])57 y, s = y_all[: i0 + 1], qtr[: i0 + 1] # information set: quarters <= origin58 fc, p, nobs = fit_forecast(y, s, max(HORIZONS)); p_used.append(p); n_last = nobs59 for h in HORIZONS:60 tgt = i0 + h; actual = float(y_all[tgt])61 rows.append(dict(origin_date=o.date().isoformat(), target_date=dates[tgt].date().isoformat(), horizon=h,62 forecast=float(y_all[tgt - 4]), actual=actual, model="snaive"))63 assert tgt - 4 <= i064 rows.append(dict(origin_date=o.date().isoformat(), target_date=dates[tgt].date().isoformat(), horizon=h,65 forecast=fc[h - 1], actual=actual, model=MODEL))66F = pd.DataFrame(rows); F.to_csv(os.path.join(O, "forecasts.csv"), index=False, float_format="%.6f")6768def hac_se(d, L, kernel):69 T = len(d); u = d - d.mean(); lrv = u @ u / T70 for j in range(1, L + 1):71 w = 1 - j / (L + 1) if kernel == "bartlett" else 1.072 lrv += 2 * w * (u[j:] @ u[:-j]) / T73 return float(np.sqrt(lrv / T)) if lrv > 0 else None, float(lrv)7475res = {"task_id": "F3-conf-04", "model_name": MODEL,76 "model_description": ("Seasonal autoregression: OLS of quarterly annualised inflation on four quarter dummies and p own lags, "77 "p in {0,...,4} chosen by AIC on a common sample at each origin, re-estimated by OLS on the expanding window "78 "from 1995Q2 to the origin, iterated multi-step forecasts."),79 "n_origins": len(origins), "first_origin": origins[0].date().isoformat(), "last_origin": origins[-1].date().isoformat()}80E = {}81for h in HORIZONS:82 for mod, key in (("snaive", "snaive"), (MODEL, "model")):83 g = F[(F.horizon == h) & (F.model == mod)].sort_values("origin_date")84 e = (g.actual - g.forecast).to_numpy(); E[(h, key)] = e85 res[f"rmse_{key}_h{h}"] = float(np.sqrt(np.mean(e ** 2))); res[f"mae_{key}_h{h}"] = float(np.mean(np.abs(e)))86for h in HORIZONS:87 res[f"rel_improvement_rmse_h{h}"] = 1 - res[f"rmse_model_h{h}"] / res[f"rmse_snaive_h{h}"]88d = E[(4, "snaive")] ** 2 - E[(4, "model")] ** 2; T = len(d)89se, _ = hac_se(d, HAC_LAGS, KERNEL)90res.update({"loss_diff_mean_h4": float(d.mean()), "loss_diff_se_h4": se, "hac_lags_h4": HAC_LAGS, "hac_kernel_h4": KERNEL,91 "dm_stat_h4": float(d.mean() / se), "dm_pvalue_h4": float(2 * stats.norm.sf(abs(d.mean() / se)))})92# sensitivity of the h=4 test (supplementary keys)93sens = {}94for L, k in ((3, "bartlett"), (3, "rectangular"), (6, "bartlett")):95 s_, lrv = hac_se(d, L, k)96 sens[f"{k}_L{L}"] = {"se": s_, "lrv": lrv, "dm_stat": (float(d.mean() / s_) if s_ else None),97 "pvalue_normal": (float(2 * stats.norm.sf(abs(d.mean() / s_))) if s_ else None)}98h = 4; hln = np.sqrt((T + 1 - 2 * h + h * (h - 1) / T) / T) # Harvey-Leybourne-Newbold correction99dm_hln = res["dm_stat_h4"] * hln100sens["hln_corrected_main"] = {"dm_stat": float(dm_hln), "pvalue_t_T-1": float(2 * stats.t.sf(abs(dm_hln), T - 1))}101res["sensitivity_h4"] = sens102res["ar_lag_order_selected"] = {"min": int(min(p_used)), "max": int(max(p_used)),103 "counts": {str(p): int(p_used.count(p)) for p in sorted(set(p_used))}}104res.update({"feasibility": "complete", "feasibility_reason": None, "data_issues": data_issues,105 "n_input_rows": {"cpi_quarterly.csv": int(len(q_raw)), "cpi_monthly.csv": int(len(m_raw))},106 "n_analysed": int(n_last), "random_seed": None})107with open(os.path.join(O, "results.json"), "w") as f: json.dump(res, f, indent=2)108print(json.dumps({k: v for k, v in res.items() if k not in ("data_issues", "model_description")}, indent=1))109