1"""F3-conf-01: pseudo-out-of-sample forecasts of Canadian monthly GDP y/y growth under publication lags.2Run from workspace root: ~/wp29/venv/bin/python out/analysis.py3"""4import json5import numpy as np6import pandas as pd78DATA = "data/canada_gdp_employment.csv"9MODEL = "direct_ardl_emp"10COVID = (pd.Timestamp("2020-03-01"), pd.Timestamp("2021-12-01")) # excluded from training (pre-specified)111213def fit_predict(d, t, h, exclude_covid=True):14 """Direct regression y_{m+k} = a + b*y_m + c*e_{m+1} + u, k = h+2, using data available at end of t."""15 k = h + 216 T = t - pd.DateOffset(months=2) # last GDP month available17 E = t - pd.DateOffset(months=1) # last employment month available18 y = d["gdp_yoy"].where(d.index <= T)19 e = d["emp_yoy"].where(d.index <= E)20 idx = d.index21 rows = []22 for i in range(len(idx) - k):23 m = idx[i]24 tgt = idx[i + k]25 if tgt > T:26 break27 ym, em1, yk = y.iloc[i], e.iloc[i + 1], y.iloc[i + k]28 if np.isnan(ym) or np.isnan(em1) or np.isnan(yk):29 continue30 if exclude_covid:31 months = idx[i:i + k + 1]32 if ((months >= COVID[0]) & (months <= COVID[1])).any():33 continue34 rows.append((1.0, ym, em1, yk, tgt))35 A = np.array([r[:3] for r in rows]); b = np.array([r[3] for r in rows])36 coef, *_ = np.linalg.lstsq(A, b, rcond=None)37 x = np.array([1.0, y.loc[T], e.loc[E]])38 last_tgt = max(max(r[4] for r in rows), T)39 return float(x @ coef), last_tgt, E, len(rows)404142def main():43 raw = pd.read_csv(DATA)44 n_rows = len(raw)45 d = raw.copy(); d["month"] = pd.to_datetime(d["month"]); d = d.set_index("month").sort_index()46 issues = []47 # checks48 assert d.index.is_unique and d.index.is_monotonic_increasing49 full = pd.date_range(d.index[0], d.index[-1], freq="MS")50 assert len(full) == len(d), "gaps in monthly index"51 chk = (100 * (d.gdp_level / d.gdp_level.shift(12) - 1)).round(4) - d.gdp_yoy52 chk_e = (100 * (d.emp_level / d.emp_level.shift(12) - 1)).round(4) - d.emp_yoy53 issues.append({"issue": "gdp_yoy/emp_yoy are populated for 1998-01..1998-12 although the file holds no 1997 levels (README says first 12 months empty); they could not be re-derived from the file.",54 "handling": f"Kept as given; for later months the yoy columns were verified against the levels (max abs diff GDP {np.nanmax(np.abs(chk)):.4g}, EMP {np.nanmax(np.abs(chk_e)):.4g})."})55 issues.append({"issue": "GDP columns empty for 2026-08 and 2026-09 (not yet published).",56 "handling": "Expected; not needed (last target month used is 2026-07)."})57 issues.append({"issue": "COVID-19 outliers in 2020-2021 (gdp_yoy range %.1f to %.1f %%)." % (d.gdp_yoy.min(), d.gdp_yoy.max()),58 "handling": "Model training pairs touching any month in 2020-03..2021-12 are excluded (pre-specified choice); naive unaffected. Sensitivity without exclusion reported in report.md."})5960 origins = pd.date_range("2022-05-01", "2026-04-01", freq="MS")61 assert len(origins) == 4862 ledger = []; sens = []63 n_last = None64 for t in origins:65 T = t - pd.DateOffset(months=2)66 for h in (1, 3):67 tgt = t + pd.DateOffset(months=h)68 act = float(d.loc[tgt, "gdp_yoy"])69 ledger.append(dict(origin_date=t.strftime("%Y-%m-%d"), target_date=tgt.strftime("%Y-%m-%d"), horizon=h,70 forecast=float(d.loc[T, "gdp_yoy"]), actual=act, model="naive",71 last_target_date_used=T.strftime("%Y-%m-%d"), last_predictor_date_used=""))72 f, lt, le, n = fit_predict(d, t, h)73 n_last = n74 ledger.append(dict(origin_date=t.strftime("%Y-%m-%d"), target_date=tgt.strftime("%Y-%m-%d"), horizon=h,75 forecast=f, actual=act, model=MODEL,76 last_target_date_used=lt.strftime("%Y-%m-%d"), last_predictor_date_used=le.strftime("%Y-%m-%d")))77 f2, *_ = fit_predict(d, t, h, exclude_covid=False)78 sens.append(dict(h=h, err=f2 - act))79 L = pd.DataFrame(ledger)80 L.to_csv("out/forecasts.csv", index=False)8182 res = {"task_id": "F3-conf-01", "model_name": MODEL,83 "model_description": ("Direct OLS per horizon h: gdp_yoy(m+h+2) = a + b*gdp_yoy(m) + c*emp_yoy(m+1), expanding window from 1998 "84 "using only pairs observed at the origin, excluding pairs touching 2020-03..2021-12; forecast plugs in gdp_yoy(t-2) and emp_yoy(t-1)."),85 "n_origins": 48, "first_origin": "2022-05-01", "last_origin": "2026-04-01",86 "last_target_date_available": "2026-02-01", "last_predictor_date_available": "2026-03-01"}87 for mod, tag in (("naive", "naive"), (MODEL, "model")):88 for h in (1, 3):89 s = L[(L.model == mod) & (L.horizon == h)]90 err = s.forecast - s.actual91 res[f"rmse_{tag}_h{h}"] = float(np.sqrt((err ** 2).mean()))92 res[f"mae_{tag}_h{h}"] = float(err.abs().mean())93 for h in (1, 3):94 res[f"rel_improvement_rmse_h{h}"] = 1 - res[f"rmse_model_h{h}"] / res[f"rmse_naive_h{h}"]95 res.update({"feasibility": "complete", "feasibility_reason": None, "data_issues": issues,96 "n_input_rows": {"canada_gdp_employment.csv": n_rows}, "n_analysed": n_last, "random_seed": None})97 # supplementary diagnostics (not required keys)98 S = pd.DataFrame(sens)99 supp = {}100 from scipy import stats101 for h in (1, 3):102 supp[f"rmse_model_nocovidexcl_h{h}"] = float(np.sqrt((S[S.h == h].err ** 2).mean()))103 a = L[(L.model == "naive") & (L.horizon == h)].reset_index(drop=True)104 b = L[(L.model == MODEL) & (L.horizon == h)].reset_index(drop=True)105 dl = (a.forecast - a.actual) ** 2 - (b.forecast - b.actual) ** 2106 # Diebold-Mariano with Newey-West HAC variance, lag h-1 (+ HLN small-sample correction)107 n = len(dl); dbar = dl.mean(); u = dl - dbar108 lag = h - 1109 v = (u @ u) / n + sum(2 * (1 - j / (lag + 1)) * (u[j:].values @ u[:-j].values) / n for j in range(1, lag + 1))110 dm = dbar / np.sqrt(v / n)111 hln = dm * np.sqrt((n + 1 - 2 * h + h * (h - 1) / n) / n)112 supp[f"dm_stat_h{h}"] = float(hln)113 supp[f"dm_pvalue_two_sided_h{h}"] = float(2 * stats.t.sf(abs(hln), df=n - 1))114 res["supplementary"] = supp115 with open("out/results.json", "w") as fh:116 json.dump(res, fh, indent=2)117 print(json.dumps(res, indent=2))118119120if __name__ == "__main__":121 main()122