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).4"""5import json, os6import numpy as np, pandas as pd7from scipy import stats89ROOT = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))10D = lambda *p: os.path.join(ROOT, *p)1112q = pd.read_csv(D("data", "cpi_quarterly.csv"), parse_dates=["quarter"])13m = pd.read_csv(D("data", "cpi_monthly.csv"))14n_q_rows, n_m_rows = len(q), len(m)15q = q.dropna(subset=["infl_annualised"]).reset_index(drop=True) # exclusion step 1 (1995Q1)16y = pd.Series(q["infl_annualised"].values, index=pd.PeriodIndex(q["quarter"], freq="Q"))1718ORIGINS = pd.period_range("2015Q4", "2025Q3", freq="Q")19HORIZONS = [1, 2, 4]20MODEL = "direct_sdum_ybar4"21HAC_LAGS, HAC_KERNEL = 4, "bartlett"222324def snaive(t, h):25 src = t + h - 426 assert src <= t27 return float(y[src])282930def model_forecast(t, h):31 """Direct OLS: y_{s+h} = sum_q g_q*1[quarter(s+h)=q] + b*mean(y_{s-3..s}) + e,32 estimated on all pairs with s+h <= t (only data dated <= t)."""33 ys = y[y.index <= t] # information set at origin t34 ybar4 = ys.rolling(4).mean()35 rows, tgt = [], []36 for s in ys.index:37 if s + h > t or np.isnan(ybar4[s]):38 continue39 rows.append(_x(s + h, ybar4[s])); tgt.append(ys[s + h])40 X, Y = np.array(rows), np.array(tgt)41 beta, *_ = np.linalg.lstsq(X, Y, rcond=None)42 fc = float(_x(t + h, ybar4[t]) @ beta)43 return fc, len(Y)444546def _x(target_period, ybar):47 d = np.zeros(4); d[target_period.quarter - 1] = 1.048 return np.r_[d, ybar]495051recs, n_last = [], {}52for t in ORIGINS:53 for h in HORIZONS:54 act = float(y[t + h])55 recs.append(dict(origin_date=t.start_time.date().isoformat(), target_date=(t + h).start_time.date().isoformat(),56 horizon=h, forecast=snaive(t, h), actual=act, model="snaive"))57 fc, n = model_forecast(t, h)58 if t == ORIGINS[-1]:59 n_last[h] = n60 recs.append(dict(origin_date=t.start_time.date().isoformat(), target_date=(t + h).start_time.date().isoformat(),61 horizon=h, forecast=fc, actual=act, model=MODEL))62fc = pd.DataFrame(recs)63fc["error"] = fc["actual"] - fc["forecast"]64fc.drop(columns="error").to_csv(D("out", "forecasts.csv"), index=False, float_format="%.6f")656667def hac_var_mean(d, L, kernel):68 d = np.asarray(d) - np.mean(d); T = len(d)69 g = lambda k: np.sum(d[k:] * d[:T - k]) / T70 w = (lambda k: 1 - k / (L + 1)) if kernel == "bartlett" else (lambda k: 1.0)71 lrv = g(0) + 2 * sum(w(k) * g(k) for k in range(1, L + 1))72 return lrv / T737475res = dict(task_id="F3-conf-04", model_name=MODEL,76 model_description=("Direct (one regression per horizon) OLS of y_{t+h} on four target-quarter seasonal dummies "77 "and the mean of the last four quarterly inflation rates at the origin (y_t..y_{t-3}), "78 "no intercept, re-estimated at each origin on an expanding window of all pairs with "79 "target quarter <= origin (sample starts 1995Q2); no model selection."),80 n_origins=len(ORIGINS), first_origin=ORIGINS[0].start_time.date().isoformat(),81 last_origin=ORIGINS[-1].start_time.date().isoformat())82err = {}83for h in HORIZONS:84 for mod, key in [("snaive", "snaive"), (MODEL, "model")]:85 e = fc[(fc.horizon == h) & (fc.model == mod)].sort_values("origin_date")["error"].values86 assert len(e) == 4087 err[(key, h)] = e88 res[f"rmse_{key}_h{h}"] = float(np.sqrt(np.mean(e ** 2)))89 res[f"mae_{key}_h{h}"] = float(np.mean(np.abs(e)))90for h in HORIZONS:91 res[f"rel_improvement_rmse_h{h}"] = 1 - res[f"rmse_model_h{h}"] / res[f"rmse_snaive_h{h}"]9293d = err[("snaive", 4)] ** 2 - err[("model", 4)] ** 294dbar = float(np.mean(d)); se = float(np.sqrt(hac_var_mean(d, HAC_LAGS, HAC_KERNEL)))95dm = dbar / se96res.update(loss_diff_mean_h4=dbar, loss_diff_se_h4=se, hac_lags_h4=HAC_LAGS, hac_kernel_h4=HAC_KERNEL,97 dm_stat_h4=dm, dm_pvalue_h4=float(2 * stats.norm.sf(abs(dm))))98# sensitivity analyses (not the primary test)99T = len(d); h = 4100hln = np.sqrt((T + 1 - 2 * h + h * (h - 1) / T) / T)101sens = {"bartlett_lag3": None, "rectangular_lag3": None}102for name, L, k in [("bartlett_lag3", 3, "bartlett"), ("rectangular_lag3", 3, "rectangular"), ("bartlett_lag6", 6, "bartlett")]:103 v = hac_var_mean(d, L, k)104 if v > 0:105 s = dbar / np.sqrt(v)106 sens[name] = dict(se=float(np.sqrt(v)), dm_stat=float(s), p_normal=float(2 * stats.norm.sf(abs(s))))107 else:108 sens[name] = dict(se=None, note="non-positive long-run variance estimate")109s_hln = dm * hln110sens["primary_with_HLN_correction_t39"] = dict(dm_stat=float(s_hln), p_t=float(2 * stats.t.sf(abs(s_hln), T - 1)))111res["sensitivity_dm_h4"] = sens112res.update(feasibility="complete", feasibility_reason=None,113 data_issues=[{"issue": "infl_annualised missing for 1995Q1 (no prior quarter)", "handling": "row dropped (exclusion step 1)"},114 {"issue": "infl_annualised differs from 400*dln(cpi_avg) by up to ~3e-6 due to 6-dp rounding of cpi_avg",115 "handling": "file values used as given"}],116 n_input_rows={"cpi_quarterly.csv": n_q_rows, "cpi_monthly.csv": n_m_rows},117 n_analysed=int(n_last[1]), n_analysed_by_horizon_last_origin={f"h{k}": int(v) for k, v in n_last.items()},118 random_seed=None)119with open(D("out", "results.json"), "w") as f:120 json.dump(res, f, indent=2)121for k in ["rmse_snaive_h1", "rmse_model_h1", "rmse_snaive_h2", "rmse_model_h2", "rmse_snaive_h4", "rmse_model_h4",122 "loss_diff_mean_h4", "loss_diff_se_h4", "dm_stat_h4", "dm_pvalue_h4"]:123 print(f"{k}: {res[k]:.4f}")124