1"""F3-conf-04: multi-horizon forecasts of Canadian quarterly CPI inflation.2Run from workspace root: ~/wp29/venv/bin/python out/analysis.py3Model: direct h-step OLS regression of y_{s+h} on 4 quarter-of-target dummies and4p lags y_s..y_{s-p+1}, p in {0..4} chosen by BIC at each origin and horizon on a5common sample, expanding window using only quarters <= origin. No randomness."""6import json, os7import numpy as np, pandas as pd8from scipy import stats910ROOT = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))11DATA, OUT = os.path.join(ROOT, "data"), os.path.join(ROOT, "out")12H = [1, 2, 4]; PMAX = 4; MODEL = "dirSAR_bic"1314q = pd.read_csv(os.path.join(DATA, "cpi_quarterly.csv"), parse_dates=["quarter"])15m = pd.read_csv(os.path.join(DATA, "cpi_monthly.csv"))16n_rows = {"cpi_quarterly.csv": int(len(q)), "cpi_monthly.csv": int(len(m))}17q = q.dropna(subset=["infl_annualised"]).reset_index(drop=True) # drops 1995Q1 only18dates = q["quarter"]; y = q["infl_annualised"].to_numpy(); qtr = dates.dt.quarter.to_numpy()19origins = pd.date_range("2015-10-01", "2025-07-01", freq="QS")20assert len(origins) == 402122def design(yy, qq, h, p, start):23 """Rows s (index of origin-in-sample) with s>=start and s+h<len(yy)."""24 S = np.arange(start, len(yy) - h)25 X = [(qq[S + h] == k).astype(float) for k in (1, 2, 3, 4)]26 X += [yy[S - i] for i in range(p)]27 return np.column_stack(X), yy[S + h]2829def fit_forecast(yy, qq, h):30 """yy, qq: data up to and including origin (last element = origin)."""31 start = PMAX - 1 # common sample for all p32 best = None33 for p in range(PMAX + 1):34 X, Y = design(yy, qq, h, p, start)35 b, *_ = np.linalg.lstsq(X, Y, rcond=None)36 n = len(Y); rss = np.sum((Y - X @ b) ** 2)37 bic = n * np.log(rss / n) + X.shape[1] * np.log(n)38 if best is None or bic < best[0]:39 best = (bic, p)40 p = best[1]41 X, Y = design(yy, qq, h, p, max(p - 1, 0)) # re-estimate on all usable rows42 b, *_ = np.linalg.lstsq(X, Y, rcond=None)43 tq = (qq[-1] - 1 + h) % 4 + 144 x = [float(tq == k) for k in (1, 2, 3, 4)] + [yy[-1 - i] for i in range(p)]45 return float(np.dot(b, x)), p, len(Y)4647rows, chosen_p, nobs_last = [], {h: [] for h in H}, {}48for o in origins:49 io = int(np.where(dates == o)[0][0])50 yy, qq = y[: io + 1], qtr[: io + 1] # information set: quarters <= origin51 assert dates.iloc[io] == o and dates.iloc[: io + 1].max() <= o52 for h in H:53 it = io + h; tdate = dates.iloc[it]; act = y[it]54 sn = y[it - 4]; assert it - 4 <= io55 f, p, n = fit_forecast(yy, qq, h)56 chosen_p[h].append(p)57 if o == origins[-1]: nobs_last[h] = n58 for mod, val in (("snaive", sn), (MODEL, f)):59 rows.append(dict(origin_date=o.strftime("%Y-%m-%d"), target_date=tdate.strftime("%Y-%m-%d"),60 horizon=h, forecast=round(val, 6), actual=round(act, 6), model=mod))61fc = pd.DataFrame(rows)62fc.to_csv(os.path.join(OUT, "forecasts.csv"), index=False)6364def hac_se(d, L, kernel="bartlett"):65 d = np.asarray(d) - np.mean(d); n = len(d)66 lrv = np.sum(d * d) / n67 for j in range(1, L + 1):68 w = 1 - j / (L + 1) if kernel == "bartlett" else 1.069 lrv += 2 * w * np.sum(d[j:] * d[:-j]) / n70 return float(np.sqrt(lrv / n)) if lrv > 0 else None7172res = {"task_id": "F3-conf-04", "model_name": MODEL,73 "model_description": ("Direct h-step OLS regression of inflation at t+h on four target-quarter dummies and p lags "74 "of inflation (y_t..y_{t-p+1}), p in 0..4 selected by BIC per origin and horizon, "75 "expanding window from 1995Q2 to the origin."),76 "n_origins": len(origins), "first_origin": origins[0].strftime("%Y-%m-%d"),77 "last_origin": origins[-1].strftime("%Y-%m-%d")}78err = {}79for mod, tag in (("snaive", "snaive"), (MODEL, "model")):80 for h in H:81 s = fc[(fc.model == mod) & (fc.horizon == h)]82 e = (s.actual - s.forecast).to_numpy(); err[(tag, h)] = e83 res[f"rmse_{tag}_h{h}"] = float(np.sqrt(np.mean(e ** 2)))84 res[f"mae_{tag}_h{h}"] = float(np.mean(np.abs(e)))85for h in H:86 res[f"rel_improvement_rmse_h{h}"] = 1 - res[f"rmse_model_h{h}"] / res[f"rmse_snaive_h{h}"]87d = err[("snaive", 4)] ** 2 - err[("model", 4)] ** 288L = 4; se = hac_se(d, L, "bartlett"); n = len(d)89res.update({"loss_diff_mean_h4": float(np.mean(d)), "loss_diff_se_h4": se, "hac_lags_h4": L,90 "hac_kernel_h4": "bartlett", "dm_stat_h4": float(np.mean(d) / se)})91res["dm_pvalue_h4"] = float(2 * stats.norm.sf(abs(res["dm_stat_h4"])))92# sensitivity (not required): other HAC choices and small-sample t reference93sens = {}94for k, LL in (("bartlett", 3), ("bartlett", 8), ("rectangular", 3)):95 s2 = hac_se(d, LL, k)96 sens[f"{k}_L{LL}"] = None if s2 is None else {"se": s2, "stat": float(np.mean(d) / s2),97 "p_normal": float(2 * stats.norm.sf(abs(np.mean(d) / s2))),98 "p_t39": float(2 * stats.t.sf(abs(np.mean(d) / s2), n - 1))}99sens["main_p_t39"] = float(2 * stats.t.sf(abs(res["dm_stat_h4"]), n - 1))100res["sensitivity_dm_h4"] = sens101res["dm_tests_h1_h2_bartlett_L_h"] = {}102for h in (1, 2):103 dh = err[("snaive", h)] ** 2 - err[("model", h)] ** 2; sh = hac_se(dh, h)104 res["dm_tests_h1_h2_bartlett_L_h"][f"h{h}"] = {"mean": float(np.mean(dh)), "se": sh,105 "stat": float(np.mean(dh) / sh), "p_normal": float(2 * stats.norm.sf(abs(np.mean(dh) / sh)))}106res.update({"feasibility": "complete", "feasibility_reason": None,107 "data_issues": [{"issue": "infl_annualised missing for 1995Q1 (no previous quarter)",108 "handling": "row dropped; inflation sample starts 1995Q2"},109 {"issue": "infl_annualised recomputed from rounded cpi_avg differs by up to ~3e-6",110 "handling": "file values used as given (negligible)"}],111 "n_input_rows": n_rows,112 "n_analysed": int(np.sum(dates <= origins[-1])),113 "n_regression_rows_last_origin": {f"h{h}": nobs_last[h] for h in H},114 "selected_p_counts": {f"h{h}": {str(k): chosen_p[h].count(k) for k in range(PMAX + 1)} for h in H},115 "baseline_note": "seasonal-naive baseline = keys *_snaive_*",116 "random_seed": None})117with open(os.path.join(OUT, "results.json"), "w") as fh:118 json.dump(res, fh, indent=2)119print(json.dumps({k: v for k, v in res.items() if k.startswith(("rmse", "mae", "rel", "loss", "dm_stat", "dm_p"))}, indent=1))120