1"""F3-conf-01: pseudo-real-time forecasts of Canadian monthly GDP YoY growth.2Run from workspace root: ~/wp29/venv/bin/python out/analysis.py3Reads only data/; writes out/results.json, out/forecasts.csv, out/exclusions.csv. Deterministic (no randomness)."""4import json5import numpy as np, pandas as pd6from scipy import stats78DATA = "data/canada_gdp_employment.csv"9OUT = "out"10MODEL = "direct_ols_ar_emp"11COVID = (pd.Timestamp("2020-03-01"), pd.Timestamp("2021-12-01")) # pre-specified exclusion window12H = [1, 3]1314raw = pd.read_csv(DATA)15n_raw = len(raw)16d = raw.copy()17d["month"] = pd.to_datetime(d["month"])18d = d.set_index("month").sort_index()19assert d.index.is_unique and (d.index.to_series().diff().dropna().dt.days <= 31).all()20y = d["gdp_yoy"]; e = d["emp_yoy"]21M = pd.DateOffset(months=1)22origins = pd.date_range("2022-05-01", "2026-04-01", freq="MS")23assert len(origins) == 482425def in_covid(m):26 return (m >= COVID[0]) & (m <= COVID[1])2728def design(t, h):29 """Training rows for the direct model at origin t, horizon h.30 Row indexed by pseudo-origin s: target y(s+h), regressors y(s-2), e(s-1).31 The information set is truncated first: gdp <= t-2, emp <= t-1."""32 gdp_last = t - 2 * M; emp_last = t - 1 * M33 yv = y[y.index <= gdp_last].dropna(); ev = e[e.index <= emp_last].dropna()34 rows = pd.DataFrame({"s": d.index})35 steps = [("0", "candidate pseudo-origins s = every month in file", len(rows), len(rows), "all monthly rows")]36 rows["tgt"] = rows.s + pd.DateOffset(months=h); rows["ylag"] = rows.s - 2 * M; rows["elag"] = rows.s - M37 nb = len(rows); rows = rows[rows.tgt.isin(yv.index)]38 steps.append(("1", f"target y(s+{h}) observed at origin (<= {gdp_last.date()})", nb, len(rows), "training targets must be published at origin"))39 nb = len(rows); rows = rows[rows.ylag.isin(yv.index) & rows.elag.isin(ev.index)]40 steps.append(("2", "lagged regressors y(s-2), e(s-1) exist", nb, len(rows), "start of sample"))41 nb = len(rows); rows = rows[~(in_covid(rows.tgt) | in_covid(rows.ylag) | in_covid(rows.elag))]42 steps.append(("3", "drop rows with target or regressor month in 2020-03..2021-12", nb, len(rows), "pandemic outliers would dominate OLS (pre-specified)"))43 Y = yv.loc[rows.tgt].values44 X = np.column_stack([np.ones(len(rows)), yv.loc[rows.ylag].values, ev.loc[rows.elag].values])45 x_new = np.array([1.0, yv.loc[gdp_last], ev.loc[emp_last]])46 return X, Y, x_new, rows, steps4748recs = []; last_info = {}49for t in origins:50 for h in H:51 tgt = t + pd.DateOffset(months=h)52 actual = float(y.loc[tgt])53 gdp_last = t - 2 * M54 recs.append(dict(origin_date=t.date().isoformat(), target_date=tgt.date().isoformat(), horizon=h,55 forecast=float(y.loc[gdp_last]), actual=actual, model="naive",56 last_target_date_used=gdp_last.date().isoformat(), last_predictor_date_used=""))57 X, Y, xn, rows, steps = design(t, h)58 beta, *_ = np.linalg.lstsq(X, Y, rcond=None)59 last_t = max(rows.tgt.max(), rows.ylag.max(), gdp_last)60 last_e = max(rows.elag.max(), t - M)61 assert last_t <= t - 2 * M and last_e <= t - M62 recs.append(dict(origin_date=t.date().isoformat(), target_date=tgt.date().isoformat(), horizon=h,63 forecast=float(xn @ beta), actual=actual, model=MODEL,64 last_target_date_used=last_t.date().isoformat(), last_predictor_date_used=last_e.date().isoformat()))65 last_info[h] = (len(Y), steps, beta)6667fc = pd.DataFrame(recs)68fc["forecast"] = fc.forecast.round(6)69fc.to_csv(f"{OUT}/forecasts.csv", index=False)7071def err(m, h):72 s = fc[(fc.model == m) & (fc.horizon == h)]73 u = s.forecast - s.actual74 return float(np.sqrt((u ** 2).mean())), float(u.abs().mean())7576res = {"task_id": "F3-conf-01", "model_name": MODEL,77 "model_description": ("Direct OLS per horizon h re-estimated at each origin t on an expanding window: "78 "gdp_yoy(s+h) = a + b*gdp_yoy(s-2) + c*emp_yoy(s-1), using all pseudo-origins s with s+h <= t-2 "79 "(1998 onward, excluding rows whose target or regressor month lies in 2020-03..2021-12); "80 "the forecast plugs in gdp_yoy(t-2) and emp_yoy(t-1)."),81 "n_origins": len(origins), "first_origin": origins[0].date().isoformat(), "last_origin": origins[-1].date().isoformat(),82 "last_target_date_available": (origins[-1] - 2 * M).date().isoformat(),83 "last_predictor_date_available": (origins[-1] - M).date().isoformat()}84for m, k in [("naive", "naive"), (MODEL, "model")]:85 for h in H:86 r, a = err(m, h); res[f"rmse_{k}_h{h}"] = round(r, 6); res[f"mae_{k}_h{h}"] = round(a, 6)87for h in H:88 res[f"rel_improvement_rmse_h{h}"] = round(1 - res[f"rmse_model_h{h}"] / res[f"rmse_naive_h{h}"], 6)89res["feasibility"] = "complete"; res["feasibility_reason"] = None90res["data_issues"] = [91 {"issue": "gdp_yoy and emp_yoy are populated for 1998-01..1998-12 although data/README.md says they are empty for the first 12 months; the 1997 levels needed to verify them are not in the file.",92 "handling": "Used as given (values from 1999-01 on were verified against the levels to within 5e-5 rounding error)."},93 {"issue": "GDP columns missing for 2026-08 and 2026-09 (not yet published).",94 "handling": "Not needed: latest target is 2026-07 (origin 2026-04, h=3); no rows dropped."},95 {"issue": "Pandemic outliers (GDP YoY between -16.2% and +18.5% in 2020-03..2021-12).",96 "handling": "Pre-specified exclusion of estimation rows whose target or regressor month is in 2020-03..2021-12; evaluation targets (2022-06 onward) are unaffected."},97 {"issue": "Data are the current vintage, not first releases.",98 "handling": "Exercise simulates publication timing only; no real-time validity claimed."}]99res["n_input_rows"] = {"canada_gdp_employment.csv": n_raw}100res["n_analysed"] = last_info[1][0]101res["n_analysed_h3"] = last_info[3][0]102res["last_origin_coefficients"] = {f"h{h}": {"const": round(float(b[0]), 6), "gdp_yoy_lag2": round(float(b[1]), 6), "emp_yoy_lag1": round(float(b[2]), 6)} for h, (_, _, b) in last_info.items()}103res["random_seed"] = None104# Diebold-Mariano (squared loss, rectangular HAC with h-1 lags, Harvey-Leybourne-Newbold adjustment, t(n-1))105dm = {}106for h in H:107 a = fc[(fc.model == "naive") & (fc.horizon == h)].sort_values("origin_date")108 b = fc[(fc.model == MODEL) & (fc.horizon == h)].sort_values("origin_date")109 dl = (a.forecast.values - a.actual.values) ** 2 - (b.forecast.values - b.actual.values) ** 2110 n = len(dl); dc = dl - dl.mean()111 v = dc @ dc / n112 for k in range(1, h):113 v += 2 * (dc[k:] @ dc[:-k]) / n114 st = float(dl.mean() / np.sqrt(v / n) * np.sqrt((n + 1 - 2 * h + h * (h - 1) / n) / n))115 dm[f"h{h}"] = {"dm_stat": round(st, 4), "p_value_two_sided": round(float(2 * stats.t.sf(abs(st), n - 1)), 4)}116res["dm_test_naive_vs_model"] = dm117with open(f"{OUT}/results.json", "w") as f:118 json.dump(res, f, indent=2)119120ex = []121for h in H:122 for st in last_info[h][1]:123 ex.append(dict(step=f"h{h}_last_origin_{st[0]}", rule=st[1], n_before=st[2], n_after=st[3], n_removed=st[2] - st[3], rationale=st[4]))124pd.DataFrame(ex).to_csv(f"{OUT}/exclusions.csv", index=False)125print(json.dumps({k: res[k] for k in res if k.startswith(("rmse", "mae", "rel", "n_analysed"))}, indent=1))126