1"""F3-conf-01: pseudo-out-of-sample forecasts of Canadian monthly GDP y/y growth under publication lags.2Run from workspace root: ~/wp29/venv/bin/python out/analysis.py3Reads only data/; writes out/results.json, out/forecasts.csv, out/exclusions.csv. No randomness."""4import json, os5import numpy as np, pandas as pd6from scipy import stats78DATA = "data/canada_gdp_employment.csv"9OUT = "out"10HORIZONS = [1, 3]11MODEL = "dar_emp"12COVID = (pd.Timestamp("2020-03-01"), pd.Timestamp("2021-12-01"))13FIRST_VERIFIABLE = pd.Timestamp("1999-01-01") # 1998 yoy values cannot be checked against levels1415raw = pd.read_csv(DATA)16n_raw = len(raw)17d = raw.copy()18d["month"] = pd.to_datetime(d["month"])19d = d.set_index("month").sort_index()20assert d.index.is_unique21y_all, e_all = d["gdp_yoy"], d["emp_yoy"]22origins = pd.date_range("2022-05-01", "2026-04-01", freq="MS")23M = lambda t, k: t + pd.DateOffset(months=k)24in_covid = lambda m: (m >= COVID[0]) & (m <= COVID[1])252627def train_pairs(t, h, exclude_covid=True):28 """Direct regression pairs y_m ~ y_{m-h-2} + e_{m-h-1}, using only info available at origin t."""29 y = y_all[(y_all.index <= M(t, -2)) & (y_all.index >= FIRST_VERIFIABLE)].dropna()30 e = e_all[(e_all.index <= M(t, -1)) & (e_all.index >= FIRST_VERIFIABLE)].dropna()31 rows = []32 for m in y.index: # targets observed at origin (m <= t-2)33 ml, me = M(m, -h - 2), M(m, -h - 1)34 if ml in y.index and me in e.index:35 if exclude_covid and (in_covid(m) or in_covid(ml) or in_covid(me)):36 continue37 rows.append((m, y[m], y[ml], e[me]))38 return pd.DataFrame(rows, columns=["m", "y", "ylag", "elag"]), y, e394041def fit_forecast(t, h, exclude_covid=True):42 tr, y, e = train_pairs(t, h, exclude_covid)43 X = np.column_stack([np.ones(len(tr)), tr["ylag"], tr["elag"]])44 beta = np.linalg.lstsq(X, tr["y"].to_numpy(), rcond=None)[0]45 yo, eo = M(t, -2), M(t, -1)46 assert yo == y.index.max() and eo == e.index.max()47 f = beta @ np.array([1.0, y[yo], e[eo]])48 return float(f), tr, y.index.max(), e.index.max()495051def run(exclude_covid=True):52 recs = []53 for t in origins:54 yb = y_all[y_all.index <= M(t, -2)].dropna()55 for h in HORIZONS:56 tgt = M(t, h)57 act = float(y_all[tgt]); assert np.isfinite(act)58 recs.append(dict(origin_date=t, target_date=tgt, horizon=h, forecast=float(yb.iloc[-1]), actual=act,59 model="naive", last_target_date_used=yb.index[-1], last_predictor_date_used=pd.NaT))60 f, tr, ly, le = fit_forecast(t, h, exclude_covid)61 recs.append(dict(origin_date=t, target_date=tgt, horizon=h, forecast=f, actual=act, model=MODEL,62 last_target_date_used=max(ly, tr["m"].max()), last_predictor_date_used=le))63 return pd.DataFrame(recs)646566def metrics(fc):67 out = {}68 for h in HORIZONS:69 for mod, lab in [("naive", "naive"), (MODEL, "model")]:70 err = fc.query("horizon==@h and model==@mod")["forecast"] - fc.query("horizon==@h and model==@mod")["actual"]71 out[f"rmse_{lab}_h{h}"] = float(np.sqrt(np.mean(err ** 2)))72 out[f"mae_{lab}_h{h}"] = float(np.mean(np.abs(err)))73 out[f"rel_improvement_rmse_h{h}"] = 1 - out[f"rmse_model_h{h}"] / out[f"rmse_naive_h{h}"]74 return out757677def dm_test(fc, h):78 """Diebold-Mariano (squared loss), Newey-West variance with h-1 lags, HLN small-sample correction, t(n-1)."""79 a = fc.query("horizon==@h and model=='naive'").sort_values("origin_date")80 b = fc.query("horizon==@h and model==@MODEL").sort_values("origin_date")81 dl = (b["forecast"].values - b["actual"].values) ** 2 - (a["forecast"].values - a["actual"].values) ** 282 n = len(dl); dc = dl - dl.mean()83 v = np.sum(dc ** 2) / n84 for k in range(1, h):85 v += 2 * (1 - k / h) * np.sum(dc[k:] * dc[:-k]) / n86 dm = dl.mean() / np.sqrt(v / n)87 dm *= np.sqrt((n + 1 - 2 * h + h * (h - 1) / n) / n)88 return float(dm), float(2 * stats.t.sf(abs(dm), n - 1))899091fc = run(True)92res_main = metrics(fc)93sens = metrics(run(False))9495# chronology checks on the ledger96assert len(fc) == 19297assert (fc["last_target_date_used"] <= fc["origin_date"].apply(lambda t: M(t, -2))).all()98m_ = fc["model"] == MODEL99assert (fc.loc[m_, "last_predictor_date_used"] <= fc.loc[m_, "origin_date"].apply(lambda t: M(t, -1))).all()100101# exclusion ledger: last estimation sample (origin 2026-04-01, h=1)102tL = origins[-1]103steps = []104cur = d.copy(); n0 = len(cur)105def step(name, rule, mask, why):106 global cur107 nb = len(cur); cur = cur[mask(cur)]; steps.append(dict(step=name, rule=rule, n_before=nb, n_after=len(cur), n_removed=nb - len(cur), rationale=why))108step("1_unverifiable", "drop months < 1999-01 as targets", lambda x: x.index >= FIRST_VERIFIABLE,109 "1998 yoy values not verifiable from levels in file (README says they should be empty)")110step("2_availability", f"keep target months <= origin-2 ({M(tL,-2).date()}) at last origin", lambda x: x.index <= M(tL, -2),111 "GDP publication lag: target must be observed at the origin")112step("3_lags", "drop targets whose regressors y_{m-3}, e_{m-2} precede 1999-01 (h=1)", lambda x: x.index >= M(FIRST_VERIFIABLE, 3),113 "regressors must be verifiable observations")114step("4_covid", "drop pairs with target or regressor month in 2020-03..2021-12 (h=1)",115 lambda x: ~(in_covid(x.index) | in_covid(x.index - pd.DateOffset(months=3)) | in_covid(x.index - pd.DateOffset(months=2))),116 "pandemic outliers and base effects; pre-specified, known at first origin")117excl = pd.DataFrame(steps)118n_analysed = len(train_pairs(tL, 1, True)[0])119assert n_analysed == len(cur) == n_raw - excl["n_removed"].sum()120excl.to_csv(f"{OUT}/exclusions.csv", index=False)121122dm1, p1 = dm_test(fc, 1); dm3, p3 = dm_test(fc, 3)123fmt = lambda s: s.dt.strftime("%Y-%m-%d")124led = fc.copy()125for c in ["origin_date", "target_date", "last_target_date_used"]:126 led[c] = fmt(led[c])127led["last_predictor_date_used"] = fmt(led["last_predictor_date_used"]).fillna("")128led["forecast"] = led["forecast"].round(6)129led.to_csv(f"{OUT}/forecasts.csv", index=False)130131results = {132 "task_id": "F3-conf-01",133 "model_name": MODEL,134 "model_description": ("Direct OLS per horizon h: gdp_yoy[m] = a + b*gdp_yoy[m-h-2] + c*emp_yoy[m-h-1], re-estimated at each origin "135 "on an expanding window from 1999-01 using only GDP<=t-2 and employment<=t-1, excluding pairs touching "136 "2020-03..2021-12; forecast = a + b*gdp_yoy[t-2] + c*emp_yoy[t-1]."),137 "n_origins": int(len(origins)),138 "first_origin": str(origins[0].date()),139 "last_origin": str(origins[-1].date()),140 "last_target_date_available": str(M(tL, -2).date()),141 "last_predictor_date_available": str(M(tL, -1).date()),142 **{k: res_main[k] for k in ["rmse_naive_h1", "mae_naive_h1", "rmse_naive_h3", "mae_naive_h3",143 "rmse_model_h1", "mae_model_h1", "rmse_model_h3", "mae_model_h3",144 "rel_improvement_rmse_h1", "rel_improvement_rmse_h3"]},145 "feasibility": "complete",146 "feasibility_reason": None,147 "data_issues": [148 {"issue": "gdp_yoy/emp_yoy populated for 1998 although README says empty; not verifiable from levels",149 "handling": "1998 months excluded from all estimation samples"},150 {"issue": "COVID-19 outliers 2020-03..2021-12 (GDP yoy from -16.2 to +18.5 %)",151 "handling": "training pairs touching this window excluded from model estimation; sensitivity without exclusion reported (sens_* keys)"},152 {"issue": "current-vintage data, not real-time vintages; future months present in file",153 "handling": "per-origin truncation to GDP<=t-2, employment<=t-1; no real-time validity claimed"},154 {"issue": "GDP missing for 2026-08 and 2026-09 (not yet published)",155 "handling": "not needed: latest target is 2026-07"},156 ],157 "n_input_rows": {"canada_gdp_employment.csv": int(n_raw)},158 "n_analysed": int(n_analysed),159 "n_analysed_h3_last_origin": int(len(train_pairs(tL, 3, True)[0])),160 "random_seed": None,161 "dm_stat_h1": dm1, "dm_pvalue_h1": p1, "dm_stat_h3": dm3, "dm_pvalue_h3": p3,162 **{f"sens_nocovidexcl_{k}": v for k, v in sens.items() if k.startswith(("rmse_model", "mae_model", "rel_"))},163}164with open(f"{OUT}/results.json", "w") as fh:165 json.dump(results, fh, indent=2)166for k in ["rmse_naive_h1", "rmse_model_h1", "rel_improvement_rmse_h1", "rmse_naive_h3", "rmse_model_h3", "rel_improvement_rmse_h3"]:167 print(k, round(results[k], 4))168