1"""F3-conf-01: pseudo-out-of-sample forecasts of Canadian monthly GDP y/y growth2under simulated publication lags. Run from workspace root:3 ~/wp29/venv/bin/python out/analysis.py4"""5import json6import numpy as np7import pandas as pd89DATA = "data/canada_gdp_employment.csv"10MODEL = "ardl_direct_emp"11HORIZONS = [1, 3]12# Pre-specified (not tuned on the evaluation sample): drop training rows touching13# the COVID collapse/rebound months, whose y/y values (-16% .. +18.5%) dominate OLS.14COVID_START, COVID_END = pd.Timestamp("2020-03-01"), pd.Timestamp("2021-06-01")1516raw = pd.read_csv(DATA)17n_rows = len(raw)18d = raw.copy()19d["month"] = pd.to_datetime(d["month"])20d = d.set_index("month").sort_index()21assert d.index.is_unique22assert (d.index.to_series().diff().dropna().dt.days.between(28, 31)).all()2324# data checks25issues = []26recomp_g = 100 * (d.gdp_level / d.gdp_level.shift(12) - 1)27recomp_e = 100 * (d.emp_level / d.emp_level.shift(12) - 1)28mg = (recomp_g - d.gdp_yoy).abs().max()29me = (recomp_e - d.emp_yoy).abs().max()30issues.append({"issue": f"yoy columns are rounded to 4 decimals; recomputation from levels (where 12-month lags exist in file) differs by at most {mg:.2e} (GDP) / {me:.2e} (employment) pp",31 "handling": "used the provided gdp_yoy/emp_yoy columns as given"})32issues.append({"issue": "gdp_level/gdp_yoy missing for 2026-08 and 2026-09 (not yet published)",33 "handling": "not needed: last target month used is 2026-07 (origin 2026-04, h=3); availability rule applied anyway"})34issues.append({"issue": "COVID-19 outliers: gdp_yoy ranges from -16.2% (2020-04) to +18.5% (2021-04), emp_yoy similarly",35 "handling": "training rows whose target or predictor month lies in 2020-03..2021-06 are excluded from model estimation (pre-specified); naive forecasts unaffected; no evaluation target falls in that window"})36issues.append({"issue": "file is the current (revised) vintage, not real-time first releases",37 "handling": "publication timing simulated only; no real-time validity claimed"})3839y = d["gdp_yoy"]40e = d["emp_yoy"]41months = d.index4243def shift(m, k):44 return m + pd.DateOffset(months=k)4546def in_covid(m):47 return COVID_START <= m <= COVID_END4849origins = pd.date_range("2022-05-01", "2026-04-01", freq="MS")50assert len(origins) == 485152rows = []53n_last = {}54for t in origins:55 g_last = shift(t, -2) # latest GDP month available56 e_last = shift(t, -1) # latest employment month available57 y_av = y[(y.index <= g_last)].dropna()58 e_av = e[(e.index <= e_last)].dropna()59 assert y_av.index.max() == g_last and e_av.index.max() == e_last60 for h in HORIZONS:61 target = shift(t, h)62 actual = float(y.loc[target])63 # naive64 rows.append(dict(origin_date=t, target_date=target, horizon=h,65 forecast=float(y_av.loc[g_last]), actual=actual, model="naive",66 last_target_date_used=g_last, last_predictor_date_used=None))67 # direct regression: y_{s+h} = a + b*y_{s-2} + c*e_{s-1}, s = pseudo-origin68 X, Y = [], []69 for s in months:70 ty, py, pe = shift(s, h), shift(s, -2), shift(s, -1)71 if ty > g_last: # training target must be observed at origin t72 continue73 if py not in y_av.index or pe not in e_av.index or ty not in y_av.index:74 continue75 if any(in_covid(m) for m in (ty, py, pe)):76 continue77 X.append([1.0, y_av.loc[py], e_av.loc[pe]])78 Y.append(y_av.loc[ty])79 X, Y = np.array(X), np.array(Y)80 beta, *_ = np.linalg.lstsq(X, Y, rcond=None)81 xf = np.array([1.0, y_av.loc[g_last], e_av.loc[e_last]])82 n_last[h] = len(Y)83 rows.append(dict(origin_date=t, target_date=target, horizon=h,84 forecast=float(xf @ beta), actual=actual, model=MODEL,85 last_target_date_used=g_last, last_predictor_date_used=e_last))8687fc = pd.DataFrame(rows)88fc["error"] = fc["forecast"] - fc["actual"]89out = fc.drop(columns="error").copy()90for c in ["origin_date", "target_date", "last_target_date_used"]:91 out[c] = out[c].dt.strftime("%Y-%m-%d")92out["last_predictor_date_used"] = out["last_predictor_date_used"].apply(93 lambda v: "" if v is None or pd.isna(v) else pd.Timestamp(v).strftime("%Y-%m-%d"))94out.to_csv("out/forecasts.csv", index=False, float_format="%.6f")95assert len(out) == 1929697def metrics(m, h):98 err = fc[(fc.model == m) & (fc.horizon == h)]["error"].values99 return float(np.sqrt(np.mean(err ** 2))), float(np.mean(np.abs(err)))100101def dm_test(h):102 """Diebold-Mariano (squared loss). Effective forecast distance from last observed GDP month is h+2,103 so Newey-West HAC uses h+1 lags and the HLN correction uses h+2."""104 from scipy import stats105 a = fc[(fc.model == MODEL) & (fc.horizon == h)].sort_values("origin_date")["error"].values106 b = fc[(fc.model == "naive") & (fc.horizon == h)].sort_values("origin_date")["error"].values107 dl = a ** 2 - b ** 2108 T = len(dl); L = h + 1; he = h + 2109 dc = dl - dl.mean()110 lrv = dc @ dc / T111 for k in range(1, L + 1):112 lrv += 2 * (1 - k / (L + 1)) * (dc[k:] @ dc[:-k]) / T113 stat = dl.mean() / np.sqrt(lrv / T)114 hln = np.sqrt((T + 1 - 2 * he + he * (he - 1) / T) / T)115 stat *= hln116 p = 2 * stats.t.sf(abs(stat), df=T - 1)117 return float(dl.mean()), float(stat), float(p)118119res = {"task_id": "F3-conf-01", "model_name": MODEL,120 "model_description": ("Direct h-step OLS of gdp_yoy(t+h) on a constant, gdp_yoy(t-2) and emp_yoy(t-1), "121 "re-estimated at each origin on an expanding window of all pseudo-origins since 1998 whose "122 "target month is <= t-2, excluding rows touching 2020-03..2021-06 (COVID)."),123 "n_origins": 48, "first_origin": "2022-05-01", "last_origin": "2026-04-01",124 "last_target_date_available": "2026-02-01", "last_predictor_date_available": "2026-03-01"}125for h in HORIZONS:126 r, a = metrics("naive", h); res[f"rmse_naive_h{h}"] = r; res[f"mae_naive_h{h}"] = a127 r, a = metrics(MODEL, h); res[f"rmse_model_h{h}"] = r; res[f"mae_model_h{h}"] = a128for h in HORIZONS:129 res[f"rel_improvement_rmse_h{h}"] = 1 - res[f"rmse_model_h{h}"] / res[f"rmse_naive_h{h}"]130res.update({"feasibility": "complete", "feasibility_reason": None, "data_issues": issues,131 "n_input_rows": {"canada_gdp_employment.csv": n_rows},132 "n_analysed": int(n_last[1]), "random_seed": None})133# supplementary diagnostics (extra keys)134res["n_analysed_h3_last_origin"] = int(n_last[3])135for h in HORIZONS:136 md, s, p = dm_test(h)137 res[f"dm_mean_lossdiff_h{h}"] = md; res[f"dm_stat_h{h}"] = s; res[f"dm_pvalue_h{h}"] = p138 sub = fc[(fc.horizon == h)]139 res[f"mean_error_naive_h{h}"] = float(sub[sub.model == "naive"].error.mean())140 res[f"mean_error_model_h{h}"] = float(sub[sub.model == MODEL].error.mean())141with open("out/results.json", "w") as f:142 json.dump(res, f, indent=2)143print(json.dumps({k: v for k, v in res.items() if k != "data_issues"}, indent=2))144