"""F3-conf-01: pseudo-out-of-sample forecasts of Canadian monthly GDP y/y growth under publication lags. Run from workspace root: ~/wp29/venv/bin/python out/analysis.py """ import json import numpy as np import pandas as pd DATA = "data/canada_gdp_employment.csv" MODEL = "direct_ardl_emp" COVID = (pd.Timestamp("2020-03-01"), pd.Timestamp("2021-12-01")) # excluded from training (pre-specified) def fit_predict(d, t, h, exclude_covid=True): """Direct regression y_{m+k} = a + b*y_m + c*e_{m+1} + u, k = h+2, using data available at end of t.""" k = h + 2 T = t - pd.DateOffset(months=2) # last GDP month available E = t - pd.DateOffset(months=1) # last employment month available y = d["gdp_yoy"].where(d.index <= T) e = d["emp_yoy"].where(d.index <= E) idx = d.index rows = [] for i in range(len(idx) - k): m = idx[i] tgt = idx[i + k] if tgt > T: break ym, em1, yk = y.iloc[i], e.iloc[i + 1], y.iloc[i + k] if np.isnan(ym) or np.isnan(em1) or np.isnan(yk): continue if exclude_covid: months = idx[i:i + k + 1] if ((months >= COVID[0]) & (months <= COVID[1])).any(): continue rows.append((1.0, ym, em1, yk, tgt)) A = np.array([r[:3] for r in rows]); b = np.array([r[3] for r in rows]) coef, *_ = np.linalg.lstsq(A, b, rcond=None) x = np.array([1.0, y.loc[T], e.loc[E]]) last_tgt = max(max(r[4] for r in rows), T) return float(x @ coef), last_tgt, E, len(rows) def main(): raw = pd.read_csv(DATA) n_rows = len(raw) d = raw.copy(); d["month"] = pd.to_datetime(d["month"]); d = d.set_index("month").sort_index() issues = [] # checks assert d.index.is_unique and d.index.is_monotonic_increasing full = pd.date_range(d.index[0], d.index[-1], freq="MS") assert len(full) == len(d), "gaps in monthly index" chk = (100 * (d.gdp_level / d.gdp_level.shift(12) - 1)).round(4) - d.gdp_yoy chk_e = (100 * (d.emp_level / d.emp_level.shift(12) - 1)).round(4) - d.emp_yoy 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.", "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})."}) issues.append({"issue": "GDP columns empty for 2026-08 and 2026-09 (not yet published).", "handling": "Expected; not needed (last target month used is 2026-07)."}) issues.append({"issue": "COVID-19 outliers in 2020-2021 (gdp_yoy range %.1f to %.1f %%)." % (d.gdp_yoy.min(), d.gdp_yoy.max()), "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."}) origins = pd.date_range("2022-05-01", "2026-04-01", freq="MS") assert len(origins) == 48 ledger = []; sens = [] n_last = None for t in origins: T = t - pd.DateOffset(months=2) for h in (1, 3): tgt = t + pd.DateOffset(months=h) act = float(d.loc[tgt, "gdp_yoy"]) ledger.append(dict(origin_date=t.strftime("%Y-%m-%d"), target_date=tgt.strftime("%Y-%m-%d"), horizon=h, forecast=float(d.loc[T, "gdp_yoy"]), actual=act, model="naive", last_target_date_used=T.strftime("%Y-%m-%d"), last_predictor_date_used="")) f, lt, le, n = fit_predict(d, t, h) n_last = n ledger.append(dict(origin_date=t.strftime("%Y-%m-%d"), target_date=tgt.strftime("%Y-%m-%d"), horizon=h, forecast=f, actual=act, model=MODEL, last_target_date_used=lt.strftime("%Y-%m-%d"), last_predictor_date_used=le.strftime("%Y-%m-%d"))) f2, *_ = fit_predict(d, t, h, exclude_covid=False) sens.append(dict(h=h, err=f2 - act)) L = pd.DataFrame(ledger) L.to_csv("out/forecasts.csv", index=False) res = {"task_id": "F3-conf-01", "model_name": MODEL, "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 " "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)."), "n_origins": 48, "first_origin": "2022-05-01", "last_origin": "2026-04-01", "last_target_date_available": "2026-02-01", "last_predictor_date_available": "2026-03-01"} for mod, tag in (("naive", "naive"), (MODEL, "model")): for h in (1, 3): s = L[(L.model == mod) & (L.horizon == h)] err = s.forecast - s.actual res[f"rmse_{tag}_h{h}"] = float(np.sqrt((err ** 2).mean())) res[f"mae_{tag}_h{h}"] = float(err.abs().mean()) for h in (1, 3): res[f"rel_improvement_rmse_h{h}"] = 1 - res[f"rmse_model_h{h}"] / res[f"rmse_naive_h{h}"] res.update({"feasibility": "complete", "feasibility_reason": None, "data_issues": issues, "n_input_rows": {"canada_gdp_employment.csv": n_rows}, "n_analysed": n_last, "random_seed": None}) # supplementary diagnostics (not required keys) S = pd.DataFrame(sens) supp = {} from scipy import stats for h in (1, 3): supp[f"rmse_model_nocovidexcl_h{h}"] = float(np.sqrt((S[S.h == h].err ** 2).mean())) a = L[(L.model == "naive") & (L.horizon == h)].reset_index(drop=True) b = L[(L.model == MODEL) & (L.horizon == h)].reset_index(drop=True) dl = (a.forecast - a.actual) ** 2 - (b.forecast - b.actual) ** 2 # Diebold-Mariano with Newey-West HAC variance, lag h-1 (+ HLN small-sample correction) n = len(dl); dbar = dl.mean(); u = dl - dbar lag = h - 1 v = (u @ u) / n + sum(2 * (1 - j / (lag + 1)) * (u[j:].values @ u[:-j].values) / n for j in range(1, lag + 1)) dm = dbar / np.sqrt(v / n) hln = dm * np.sqrt((n + 1 - 2 * h + h * (h - 1) / n) / n) supp[f"dm_stat_h{h}"] = float(hln) supp[f"dm_pvalue_two_sided_h{h}"] = float(2 * stats.t.sf(abs(hln), df=n - 1)) res["supplementary"] = supp with open("out/results.json", "w") as fh: json.dump(res, fh, indent=2) print(json.dumps(res, indent=2)) if __name__ == "__main__": main()