1"""F3-conf-04: multi-horizon forecasts of Canadian quarterly CPI inflation.2Run from workspace root: ~/wp29/venv/bin/python out/analysis.py3Reads data/ only; writes out/results.json, out/forecasts.csv, out/exclusions.csv,4out/sensitivity_dm.csv. Deterministic (no randomness)."""5import json, os6import numpy as np, pandas as pd7from scipy import stats89OUT = "out"10HORIZONS = [1, 2, 4]11PMAX = 412MODEL = "dsar_bic"1314q = pd.read_csv("data/cpi_quarterly.csv", parse_dates=["quarter"])15m = pd.read_csv("data/cpi_monthly.csv")16n_q, n_m = len(q), len(m)17q = q.sort_values("quarter").reset_index(drop=True)1819excl = []20n0 = len(q)21y_all = q.dropna(subset=["infl_annualised"]).reset_index(drop=True)22excl.append(["1", "drop quarters with missing infl_annualised", n0, len(y_all), n0 - len(y_all),23 "1995Q1 has no previous quarter, so its inflation is undefined"])2425origins = pd.date_range("2015-10-01", "2025-07-01", freq="QS")26assert len(origins) == 4027dates = y_all["quarter"].reset_index(drop=True)28y = y_all["infl_annualised"].to_numpy()29qtr = dates.dt.quarter.to_numpy()30pos = {d: i for i, d in enumerate(dates)}313233def design(yv, qv, idx, h, p):34 """Rows for direct regression of y[s+h] on target-quarter dummies and y[s],...,y[s-p+1]."""35 X = []36 for s in idx:37 d = [1.0 if qv[s + h] == k else 0.0 for k in (1, 2, 3, 4)] if s + h < len(qv) else None38 X.append((d, [yv[s - j] for j in range(p)]))39 return X404142def fit_forecast(t_idx, h):43 """Direct model using only observations with index <= t_idx (quarters <= origin)."""44 yv = y[: t_idx + 1]45 qv = qtr[: t_idx + 1]46 # target quarter of the forecast (deterministic calendar; known at origin)47 tq = ((qv[t_idx] - 1 + h) % 4) + 148 # common estimation sample for all p: s from PMAX-1 to t_idx-h49 S = np.arange(PMAX - 1, t_idx - h + 1)50 Y = yv[S + h]51 D = np.array([[1.0 if qv[s + h] == k else 0.0 for k in (1, 2, 3, 4)] for s in S])52 best = None53 for p in range(PMAX + 1):54 L = np.column_stack([yv[S - j] for j in range(p)]) if p else np.empty((len(S), 0))55 X = np.hstack([D, L])56 beta, *_ = np.linalg.lstsq(X, Y, rcond=None)57 rss = float(np.sum((Y - X @ beta) ** 2))58 n, k = X.shape59 bic = n * np.log(rss / n) + k * np.log(n)60 if best is None or bic < best[0]:61 best = (bic, p, beta, n)62 _, p, beta, n = best63 x = np.r_[[1.0 if tq == k else 0.0 for k in (1, 2, 3, 4)], [yv[t_idx - j] for j in range(p)]]64 return float(x @ beta), p, n656667rows, psel = [], []68for t in origins:69 ti = pos[t]70 for h in HORIZONS:71 tgt = dates[ti + h]72 assert tgt == t + pd.DateOffset(months=3 * h)73 actual = float(y[ti + h])74 sn = float(y[ti + h - 4]) # quarter t+h-4 <= t75 assert dates[ti + h - 4] <= t76 f, p, n = fit_forecast(ti, h)77 psel.append({"origin": t, "h": h, "p": p, "n_reg": n})78 for mod, fc in (("snaive", sn), (MODEL, f)):79 rows.append({"origin_date": t.strftime("%Y-%m-%d"), "target_date": tgt.strftime("%Y-%m-%d"),80 "horizon": h, "forecast": fc, "actual": actual, "model": mod})81fc = pd.DataFrame(rows)82fc.to_csv(os.path.join(OUT, "forecasts.csv"), index=False, float_format="%.6f")83psel = pd.DataFrame(psel)8485last_ti = pos[origins[-1]]86n_last = last_ti + 187excl.append(["2", "last estimation sample: quarters <= last origin 2025-07-01", len(y_all), n_last,88 len(y_all) - n_last, "information rule: forecasts at origin t use only quarters <= t; "89 "later quarters serve only as targets"])90pd.DataFrame(excl, columns=["step", "rule", "n_before", "n_after", "n_removed", "rationale"]).to_csv(91 os.path.join(OUT, "exclusions.csv"), index=False)9293res = {"task_id": "F3-conf-04", "model_name": MODEL,94 "model_description": ("Direct (one regression per horizon) OLS model: infl(s+h) on four target-quarter "95 "seasonal dummies plus p lags infl(s),...,infl(s-p+1), p in 0..4 chosen by BIC on a "96 "common sample, re-estimated and re-selected at each origin on an expanding window "97 "from 1995Q2 using only quarters <= origin."),98 "n_origins": len(origins), "first_origin": origins[0].strftime("%Y-%m-%d"),99 "last_origin": origins[-1].strftime("%Y-%m-%d")}100err = {}101for h in HORIZONS:102 for mod, key in (("snaive", "snaive"), (MODEL, "model")):103 s = fc[(fc.horizon == h) & (fc.model == mod)].sort_values("origin_date")104 e = (s.actual - s.forecast).to_numpy()105 err[(key, h)] = e106 res[f"rmse_{key}_h{h}"] = float(np.sqrt(np.mean(e ** 2)))107 res[f"mae_{key}_h{h}"] = float(np.mean(np.abs(e)))108for h in HORIZONS:109 res[f"rel_improvement_rmse_h{h}"] = 1 - res[f"rmse_model_h{h}"] / res[f"rmse_snaive_h{h}"]110111112def lrv(d, lags, kernel):113 d = d - d.mean()114 T = len(d)115 g = [np.sum(d[j:] * d[: T - j]) / T for j in range(lags + 1)]116 w = [1.0] + [(1 - j / (lags + 1)) if kernel == "bartlett" else 1.0 for j in range(1, lags + 1)]117 return g[0] + 2 * sum(w[j] * g[j] for j in range(1, lags + 1))118119120d = err[("snaive", 4)] ** 2 - err[("model", 4)] ** 2121T = len(d)122LAGS, KERNEL = 3, "bartlett"123dbar = float(d.mean())124se = float(np.sqrt(lrv(d, LAGS, KERNEL) / T))125dm = dbar / se126res.update({"loss_diff_mean_h4": dbar, "loss_diff_se_h4": se, "hac_lags_h4": LAGS, "hac_kernel_h4": KERNEL,127 "dm_stat_h4": dm, "dm_pvalue_h4": float(2 * stats.norm.sf(abs(dm)))})128129# sensitivity of the h=4 test (reported, not primary)130sens = []131for kern, L in (("bartlett", 3), ("bartlett", 4), ("bartlett", 6), ("rectangular", 3)):132 v = lrv(d, L, kern)133 if v > 0:134 s_ = np.sqrt(v / T); st = dbar / s_135 hln = st * np.sqrt((T + 1 - 2 * 4 + 4 * 3 / T) / T) if kern == "rectangular" else np.nan136 sens.append({"kernel": kern, "lags": L, "se": s_, "dm_stat": st, "p_normal": 2 * stats.norm.sf(abs(st)),137 "hln_stat": hln, "p_hln_t39": 2 * stats.t.sf(abs(hln), T - 1) if kern == "rectangular" else np.nan})138 else:139 sens.append({"kernel": kern, "lags": L, "se": np.nan, "dm_stat": np.nan, "p_normal": np.nan,140 "hln_stat": np.nan, "p_hln_t39": np.nan})141pd.DataFrame(sens).to_csv(os.path.join(OUT, "sensitivity_dm.csv"), index=False, float_format="%.6f")142143res.update({"feasibility": "complete", "feasibility_reason": None,144 "data_issues": [{"issue": "infl_annualised missing for 1995Q1 (no prior quarter)",145 "handling": "row excluded from all estimation samples"},146 {"issue": "large genuine outliers (2008Q4, 2020Q2, 2022Q1-Q2) in the target",147 "handling": "kept as observed; no winsorising; noted in report"}],148 "n_input_rows": {"cpi_quarterly.csv": n_q, "cpi_monthly.csv": n_m},149 "n_analysed": int(n_last), "random_seed": None,150 "supplementary": {151 "n_analysed_definition": "quarters with observed inflation <= last origin (1995Q2-2025Q3); regression obs in last fit are fewer (lags/horizon)",152 "n_reg_last_origin": {f"h{h}": int(psel[(psel.origin == origins[-1]) & (psel.h == h)].n_reg.iloc[0]) for h in HORIZONS},153 "selected_p_counts": {f"h{h}": {str(k): int(v) for k, v in psel[psel.h == h].p.value_counts().sort_index().items()} for h in HORIZONS},154 "mean_loss_diff_h1": float(np.mean(err[("snaive", 1)] ** 2 - err[("model", 1)] ** 2)),155 "mean_loss_diff_h2": float(np.mean(err[("snaive", 2)] ** 2 - err[("model", 2)] ** 2)),156 "dm_sensitivity_h4": [{k: (None if (isinstance(v, float) and np.isnan(v)) else (float(v) if not isinstance(v, str) and not isinstance(v, int) else v)) for k, v in r.items()} for r in sens],157 "loss_diff_mean_h4_excl_2021_2022_targets": float(np.mean(d[~pd.to_datetime(fc[(fc.horizon == 4) & (fc.model == 'snaive')].sort_values('origin_date').target_date).dt.year.isin([2021, 2022]).to_numpy()])),158 }})159with open(os.path.join(OUT, "results.json"), "w") as f:160 json.dump(res, f, indent=2)161for k in [k for k in res if k.startswith(("rmse", "mae", "rel_", "loss", "dm_"))]:162 print(f"{k}: {res[k]:.4f}")163