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 np6import pandas as pd78DATA = "data/canada_gdp_employment.csv"9COVID = (pd.Timestamp("2020-03-01"), pd.Timestamp("2021-12-01")) # training exclusion window10MODEL = "adl_emp_direct"1112raw = pd.read_csv(DATA, parse_dates=["month"])13n_raw = len(raw)14d = raw.set_index("month").sort_index()15months = d.index16pos = {m: i for i, m in enumerate(months)}1718def shift(m, k):19 return m + pd.DateOffset(months=k)2021def in_covid(m):22 return COVID[0] <= m <= COVID[1]2324def training_rows(t, h):25 """Direct regression rows indexed by base month m: target y[m+k], regressors y[m], e[m+1], k=h+2.26 Only data available at end of month t: gdp months <= t-2, emp months <= t-1."""27 k = h + 228 gcut, ecut = shift(t, -2), shift(t, -1)29 log = []30 m_all = [m for m in months if not np.isnan(d.loc[m, "gdp_yoy"])]31 rows = list(m_all)32 log.append(("A", "base month has non-missing gdp_yoy in file", n_raw, len(rows)))33 rows = [m for m in rows if m <= gcut]34 log.append(("B", f"base month gdp available at origin (<= t-2 = {gcut.date()})", log[-1][3], len(rows)))35 rows = [m for m in rows if shift(m, k) <= gcut and shift(m, 1) <= ecut and shift(m, k) in pos]36 log.append(("C", f"training target month m+{k} observed at origin (<= {gcut.date()})", log[-1][3], len(rows)))37 rows = [m for m in rows if not (in_covid(m) or in_covid(shift(m, 1)) or in_covid(shift(m, k)))]38 log.append(("D", "drop rows whose regressor or target month lies in 2020-03..2021-12 (pandemic)", log[-1][3], len(rows)))39 return rows, log4041origins = pd.date_range("2022-05-01", "2026-04-01", freq="MS")42assert len(origins) == 4843recs = []44last_log = {}45for t in origins:46 gcut, ecut = shift(t, -2), shift(t, -1)47 y_now = d.loc[gcut, "gdp_yoy"]; e_now = d.loc[ecut, "emp_yoy"]48 assert not np.isnan(y_now) and not np.isnan(e_now)49 for h in (1, 3):50 k = h + 251 target = shift(t, h)52 actual = float(d.loc[target, "gdp_yoy"])53 recs.append(dict(origin_date=t.date().isoformat(), target_date=target.date().isoformat(), horizon=h,54 forecast=float(y_now), actual=actual, model="naive",55 last_target_date_used=gcut.date().isoformat(), last_predictor_date_used=""))56 rows, log = training_rows(t, h)57 Y = np.array([d.loc[shift(m, k), "gdp_yoy"] for m in rows])58 X = np.column_stack([np.ones(len(rows)), [d.loc[m, "gdp_yoy"] for m in rows],59 [d.loc[shift(m, 1), "emp_yoy"] for m in rows]])60 beta, *_ = np.linalg.lstsq(X, Y, rcond=None)61 fc = float(beta @ np.array([1.0, y_now, e_now]))62 last_tgt = max(max(shift(m, k) for m in rows), gcut)63 last_emp = max(max(shift(m, 1) for m in rows), ecut)64 assert last_tgt <= gcut and last_emp <= ecut65 recs.append(dict(origin_date=t.date().isoformat(), target_date=target.date().isoformat(), horizon=h,66 forecast=fc, actual=actual, model=MODEL,67 last_target_date_used=last_tgt.date().isoformat(),68 last_predictor_date_used=last_emp.date().isoformat()))69 if t == origins[-1]:70 last_log[h] = (log, len(rows))7172fc = pd.DataFrame(recs)73fc.to_csv("out/forecasts.csv", index=False, float_format="%.6f")7475def err(model, h):76 s = fc[(fc.model == model) & (fc.horizon == h)]77 e = s.forecast - s.actual78 return float(np.sqrt((e ** 2).mean())), float(e.abs().mean())7980res = {"task_id": "F3-conf-01", "model_name": MODEL,81 "model_description": ("Direct OLS regression per horizon h of gdp_yoy[m+h+2] on a constant, gdp_yoy[m] and emp_yoy[m+1], "82 "re-estimated at each origin t on an expanding sample from 1998 using only pairs observed at t "83 "(gdp <= t-2, emp <= t-1), excluding pairs touching 2020-03..2021-12; forecast plugs in gdp_yoy[t-2] and emp_yoy[t-1]."),84 "n_origins": int(len(origins)), "first_origin": origins[0].date().isoformat(),85 "last_origin": origins[-1].date().isoformat(),86 "last_target_date_available": shift(origins[-1], -2).date().isoformat(),87 "last_predictor_date_available": shift(origins[-1], -1).date().isoformat()}88for model, tag in (("naive", "naive"), (MODEL, "model")):89 for h in (1, 3):90 r, a = err(model, h)91 res[f"rmse_{tag}_h{h}"] = r; res[f"mae_{tag}_h{h}"] = a92for h in (1, 3):93 res[f"rel_improvement_rmse_h{h}"] = 1 - res[f"rmse_model_h{h}"] / res[f"rmse_naive_h{h}"]94res.update({"feasibility": "complete", "feasibility_reason": None,95 "data_issues": [96 {"issue": "gdp_level/gdp_yoy empty for 2026-08 and 2026-09 (not yet published)",97 "handling": "never needed: latest GDP month used is 2026-02 (t-2 at last origin); latest target is 2026-07."},98 {"issue": "README says gdp_yoy/emp_yoy are empty for the first 12 months, but the file has values for 1998-01..1998-12 that cannot be recomputed from the levels in the file",99 "handling": "used as given (values for 1999 onward match 100*(L_m/L_{m-12}-1) to <5e-5); noted as a documentation inconsistency."},100 {"issue": "pandemic outliers 2020-03..2021-12 (gdp_yoy from -16.2 to +18.5 percent) would dominate OLS estimates",101 "handling": "training pairs whose regressor or target month falls in 2020-03..2021-12 are excluded (rule fixed a priori; window ends before the first origin's information set)."},102 {"issue": "current-vintage data, not historical first releases",103 "handling": "exercise simulates publication timing only; no real-time validity claimed."}],104 "n_input_rows": {"canada_gdp_employment.csv": int(n_raw)},105 "n_analysed": int(last_log[1][1]), "random_seed": None})106# Diebold-Mariano test (squared loss), Newey-West variance with h+1 lags (overlap of k=h+2-step errors)107from scipy import stats108for h in (1, 3):109 a = fc[(fc.model == "naive") & (fc.horizon == h)].reset_index(drop=True)110 b = fc[(fc.model == MODEL) & (fc.horizon == h)].reset_index(drop=True)111 dl = ((a.forecast - a.actual) ** 2 - (b.forecast - b.actual) ** 2).to_numpy()112 n = len(dl); dc = dl - dl.mean(); L = h + 1113 v = dc @ dc / n + 2 * sum((1 - j / (L + 1)) * (dc[j:] @ dc[:-j]) / n for j in range(1, L + 1))114 stat = dl.mean() / np.sqrt(v / n)115 res[f"dm_stat_h{h}"] = float(stat)116 res[f"dm_pvalue_two_sided_h{h}"] = float(2 * (1 - stats.t.cdf(abs(stat), n - 1)))117with open("out/results.json", "w") as f:118 json.dump(res, f, indent=2)119120# exclusion ledger: last origin, h=1 estimation sample (n_analysed); h=3 shown for reference121ex = []122for h in (1,):123 for step, rule, nb, na in last_log[h][0]:124 ex.append(dict(step=f"h{h}_{step}", rule=f"[origin 2026-04-01, h={h}] {rule}", n_before=nb, n_after=na,125 n_removed=nb - na, rationale="information-availability rule" if step in "BC" else126 ("missing value (unpublished GDP)" if step == "A" else "pandemic outliers excluded from estimation")))127pd.DataFrame(ex).to_csv("out/exclusions.csv", index=False)128print(json.dumps({k: v for k, v in res.items() if k.startswith(("rmse", "mae", "rel", "n_"))}, indent=1))129