1"""F3-conf-01: direct forecasts of Canadian monthly GDP yoy growth under publication lags.2Run from workspace root: ~/wp29/venv/bin/python out/analysis.py"""3import json, numpy as np, pandas as pd45DATA = "data/canada_gdp_employment.csv"6d = pd.read_csv(DATA, parse_dates=["month"])7n_rows = len(d)8d = d.set_index("month").asfreq("MS")9assert len(d) == n_rows # no gaps10MODEL = "ardl_emp_direct"11COVID = (pd.Timestamp("2020-03-01"), pd.Timestamp("2021-12-01"))12HS = [1, 3]13origins = pd.date_range("2022-05-01", "2026-04-01", freq="MS")14assert len(origins) == 481516def in_covid(m):17 return (m >= COVID[0]) & (m <= COVID[1])1819rows, n_last = [], {}20for t in origins:21 g_last = t - pd.DateOffset(months=2) # latest GDP month published at end of t22 e_last = t - pd.DateOffset(months=1) # latest LFS month published at end of t23 y = d["gdp_yoy"].where(d.index <= g_last) # mask unavailable info24 e = d["emp_yoy"].where(d.index <= e_last)25 for h in HS:26 target = t + pd.DateOffset(months=h)27 actual = d["gdp_yoy"].get(target, np.nan)28 # naive29 rows.append(dict(origin_date=t, target_date=target, horizon=h, forecast=y[g_last],30 actual=actual, model="naive", last_target_date_used=g_last,31 last_predictor_date_used=pd.NaT))32 # direct regression: y_s = a + b*y_{s-h-2} + c*e_{s-h-1}, all with s <= g_last33 tr = pd.DataFrame({"y": y, "ylag": y.shift(h + 2), "elag": e.shift(h + 1)})34 idx = tr.index35 bad = in_covid(idx) | in_covid(idx - pd.DateOffset(months=h + 2)) | in_covid(idx - pd.DateOffset(months=h + 1))36 tr = tr[~bad].dropna()37 assert tr.index.max() <= g_last38 X = np.column_stack([np.ones(len(tr)), tr["ylag"], tr["elag"]])39 beta, *_ = np.linalg.lstsq(X, tr["y"].values, rcond=None)40 assert not np.isnan(e[e_last]) and not np.isnan(y[g_last])41 fc = beta @ np.array([1.0, y[g_last], e[e_last]])42 rows.append(dict(origin_date=t, target_date=target, horizon=h, forecast=fc,43 actual=actual, model=MODEL, last_target_date_used=g_last,44 last_predictor_date_used=e_last))45 if t == origins[-1]:46 n_last[h] = len(tr)4748f = pd.DataFrame(rows)49for c in ["origin_date", "target_date", "last_target_date_used", "last_predictor_date_used"]:50 f[c] = f[c].dt.strftime("%Y-%m-%d")51f = f.fillna({"last_predictor_date_used": ""})52f["forecast"] = f["forecast"].round(6)53f.to_csv("out/forecasts.csv", index=False)5455res = {"task_id": "F3-conf-01", "model_name": MODEL,56 "model_description": ("Direct OLS per horizon h: gdp_yoy(s) = a + b*gdp_yoy(s-h-2) + c*emp_yoy(s-h-1), "57 "expanding window from 1998 using only targets observed at the origin (s <= t-2), "58 "dropping training rows whose target or regressor month lies in 2020-03..2021-12 (COVID); "59 "forecast plugs in gdp_yoy(t-2) and emp_yoy(t-1)."),60 "n_origins": len(origins), "first_origin": "2022-05-01", "last_origin": "2026-04-01",61 "last_target_date_available": (origins[-1] - pd.DateOffset(months=2)).strftime("%Y-%m-%d"),62 "last_predictor_date_available": (origins[-1] - pd.DateOffset(months=1)).strftime("%Y-%m-%d")}63stats = {}64for m, tag in [("naive", "naive"), (MODEL, "model")]:65 for h in HS:66 s = f[(f.model == m) & (f.horizon == h)]67 err = s.forecast - s.actual68 assert err.notna().sum() == 4869 res[f"rmse_{tag}_h{h}"] = float(np.sqrt(np.mean(err ** 2)))70 res[f"mae_{tag}_h{h}"] = float(np.mean(np.abs(err)))71for h in HS:72 res[f"rel_improvement_rmse_h{h}"] = 1 - res[f"rmse_model_h{h}"] / res[f"rmse_naive_h{h}"]73# Diebold-Mariano (squared loss, HAC with h-1 lags, Harvey small-sample correction) - supplementary74from scipy import stats as st75dm = {}76for h in HS:77 a = f[(f.model == "naive") & (f.horizon == h)].reset_index(drop=True)78 b = f[(f.model == MODEL) & (f.horizon == h)].reset_index(drop=True)79 dl = ((b.forecast - b.actual) ** 2 - (a.forecast - a.actual) ** 2).values80 n = len(dl); dbar = dl.mean(); u = dl - dbar81 lr = u @ u / n + 2 * sum((1 - k / h) * (u[k:] @ u[:-k]) / n for k in range(1, h))82 stat = dbar / np.sqrt(lr / n) * np.sqrt((n + 1 - 2 * h + h * (h - 1) / n) / n)83 dm[f"h{h}"] = {"dm_stat": float(stat), "p_two_sided": float(2 * st.t.sf(abs(stat), n - 1))}84res.update({"feasibility": "complete", "feasibility_reason": None,85 "data_issues": [86 {"issue": "COVID-19 outliers in 2020-2021 (gdp_yoy from -16.2% to +18.5%, emp_yoy similar) would dominate OLS fits.",87 "handling": "Training rows with target or regressor month in 2020-03..2021-12 excluded from model estimation; evaluation sample (targets 2022-06..2026-07) unaffected; naive unaffected."},88 {"issue": "Data are the current (final) vintage, not first releases; yoy columns match levels to rounding (max abs diff ~5e-5).",89 "handling": "Used as given; results simulate publication timing only, no real-time validity claimed."},90 {"issue": "GDP columns empty for 2026-08 and 2026-09.",91 "handling": "Not needed: latest target month is 2026-07; GDP never used beyond t-2."}],92 "n_input_rows": {"canada_gdp_employment.csv": n_rows},93 "n_analysed": int(n_last[1]),94 "n_analysed_note": f"training rows at last origin: h=1 {n_last[1]}, h=3 {n_last[3]} (n_analysed reports h=1)", "random_seed": None,95 "supplementary_dm_test": dm})96json.dump(res, open("out/results.json", "w"), indent=2)97print(json.dumps(res, indent=2))98