"""F3-conf-01: pseudo-out-of-sample forecasts of Canadian monthly GDP y/y growth under publication lags. Run from workspace root: ~/wp29/venv/bin/python out/analysis.py Reads only data/; writes out/results.json, out/forecasts.csv, out/exclusions.csv. No randomness.""" import json, os import numpy as np, pandas as pd from scipy import stats DATA = "data/canada_gdp_employment.csv" OUT = "out" HORIZONS = [1, 3] MODEL = "dar_emp" COVID = (pd.Timestamp("2020-03-01"), pd.Timestamp("2021-12-01")) FIRST_VERIFIABLE = pd.Timestamp("1999-01-01") # 1998 yoy values cannot be checked against levels raw = pd.read_csv(DATA) n_raw = len(raw) d = raw.copy() d["month"] = pd.to_datetime(d["month"]) d = d.set_index("month").sort_index() assert d.index.is_unique y_all, e_all = d["gdp_yoy"], d["emp_yoy"] origins = pd.date_range("2022-05-01", "2026-04-01", freq="MS") M = lambda t, k: t + pd.DateOffset(months=k) in_covid = lambda m: (m >= COVID[0]) & (m <= COVID[1]) def train_pairs(t, h, exclude_covid=True): """Direct regression pairs y_m ~ y_{m-h-2} + e_{m-h-1}, using only info available at origin t.""" y = y_all[(y_all.index <= M(t, -2)) & (y_all.index >= FIRST_VERIFIABLE)].dropna() e = e_all[(e_all.index <= M(t, -1)) & (e_all.index >= FIRST_VERIFIABLE)].dropna() rows = [] for m in y.index: # targets observed at origin (m <= t-2) ml, me = M(m, -h - 2), M(m, -h - 1) if ml in y.index and me in e.index: if exclude_covid and (in_covid(m) or in_covid(ml) or in_covid(me)): continue rows.append((m, y[m], y[ml], e[me])) return pd.DataFrame(rows, columns=["m", "y", "ylag", "elag"]), y, e def fit_forecast(t, h, exclude_covid=True): tr, y, e = train_pairs(t, h, exclude_covid) X = np.column_stack([np.ones(len(tr)), tr["ylag"], tr["elag"]]) beta = np.linalg.lstsq(X, tr["y"].to_numpy(), rcond=None)[0] yo, eo = M(t, -2), M(t, -1) assert yo == y.index.max() and eo == e.index.max() f = beta @ np.array([1.0, y[yo], e[eo]]) return float(f), tr, y.index.max(), e.index.max() def run(exclude_covid=True): recs = [] for t in origins: yb = y_all[y_all.index <= M(t, -2)].dropna() for h in HORIZONS: tgt = M(t, h) act = float(y_all[tgt]); assert np.isfinite(act) recs.append(dict(origin_date=t, target_date=tgt, horizon=h, forecast=float(yb.iloc[-1]), actual=act, model="naive", last_target_date_used=yb.index[-1], last_predictor_date_used=pd.NaT)) f, tr, ly, le = fit_forecast(t, h, exclude_covid) recs.append(dict(origin_date=t, target_date=tgt, horizon=h, forecast=f, actual=act, model=MODEL, last_target_date_used=max(ly, tr["m"].max()), last_predictor_date_used=le)) return pd.DataFrame(recs) def metrics(fc): out = {} for h in HORIZONS: for mod, lab in [("naive", "naive"), (MODEL, "model")]: err = fc.query("horizon==@h and model==@mod")["forecast"] - fc.query("horizon==@h and model==@mod")["actual"] out[f"rmse_{lab}_h{h}"] = float(np.sqrt(np.mean(err ** 2))) out[f"mae_{lab}_h{h}"] = float(np.mean(np.abs(err))) out[f"rel_improvement_rmse_h{h}"] = 1 - out[f"rmse_model_h{h}"] / out[f"rmse_naive_h{h}"] return out def dm_test(fc, h): """Diebold-Mariano (squared loss), Newey-West variance with h-1 lags, HLN small-sample correction, t(n-1).""" a = fc.query("horizon==@h and model=='naive'").sort_values("origin_date") b = fc.query("horizon==@h and model==@MODEL").sort_values("origin_date") dl = (b["forecast"].values - b["actual"].values) ** 2 - (a["forecast"].values - a["actual"].values) ** 2 n = len(dl); dc = dl - dl.mean() v = np.sum(dc ** 2) / n for k in range(1, h): v += 2 * (1 - k / h) * np.sum(dc[k:] * dc[:-k]) / n dm = dl.mean() / np.sqrt(v / n) dm *= np.sqrt((n + 1 - 2 * h + h * (h - 1) / n) / n) return float(dm), float(2 * stats.t.sf(abs(dm), n - 1)) fc = run(True) res_main = metrics(fc) sens = metrics(run(False)) # chronology checks on the ledger assert len(fc) == 192 assert (fc["last_target_date_used"] <= fc["origin_date"].apply(lambda t: M(t, -2))).all() m_ = fc["model"] == MODEL assert (fc.loc[m_, "last_predictor_date_used"] <= fc.loc[m_, "origin_date"].apply(lambda t: M(t, -1))).all() # exclusion ledger: last estimation sample (origin 2026-04-01, h=1) tL = origins[-1] steps = [] cur = d.copy(); n0 = len(cur) def step(name, rule, mask, why): global cur 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)) step("1_unverifiable", "drop months < 1999-01 as targets", lambda x: x.index >= FIRST_VERIFIABLE, "1998 yoy values not verifiable from levels in file (README says they should be empty)") step("2_availability", f"keep target months <= origin-2 ({M(tL,-2).date()}) at last origin", lambda x: x.index <= M(tL, -2), "GDP publication lag: target must be observed at the origin") step("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), "regressors must be verifiable observations") step("4_covid", "drop pairs with target or regressor month in 2020-03..2021-12 (h=1)", lambda x: ~(in_covid(x.index) | in_covid(x.index - pd.DateOffset(months=3)) | in_covid(x.index - pd.DateOffset(months=2))), "pandemic outliers and base effects; pre-specified, known at first origin") excl = pd.DataFrame(steps) n_analysed = len(train_pairs(tL, 1, True)[0]) assert n_analysed == len(cur) == n_raw - excl["n_removed"].sum() excl.to_csv(f"{OUT}/exclusions.csv", index=False) dm1, p1 = dm_test(fc, 1); dm3, p3 = dm_test(fc, 3) fmt = lambda s: s.dt.strftime("%Y-%m-%d") led = fc.copy() for c in ["origin_date", "target_date", "last_target_date_used"]: led[c] = fmt(led[c]) led["last_predictor_date_used"] = fmt(led["last_predictor_date_used"]).fillna("") led["forecast"] = led["forecast"].round(6) led.to_csv(f"{OUT}/forecasts.csv", index=False) results = { "task_id": "F3-conf-01", "model_name": MODEL, "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 " "on an expanding window from 1999-01 using only GDP<=t-2 and employment<=t-1, excluding pairs touching " "2020-03..2021-12; forecast = a + b*gdp_yoy[t-2] + c*emp_yoy[t-1]."), "n_origins": int(len(origins)), "first_origin": str(origins[0].date()), "last_origin": str(origins[-1].date()), "last_target_date_available": str(M(tL, -2).date()), "last_predictor_date_available": str(M(tL, -1).date()), **{k: res_main[k] for k in ["rmse_naive_h1", "mae_naive_h1", "rmse_naive_h3", "mae_naive_h3", "rmse_model_h1", "mae_model_h1", "rmse_model_h3", "mae_model_h3", "rel_improvement_rmse_h1", "rel_improvement_rmse_h3"]}, "feasibility": "complete", "feasibility_reason": None, "data_issues": [ {"issue": "gdp_yoy/emp_yoy populated for 1998 although README says empty; not verifiable from levels", "handling": "1998 months excluded from all estimation samples"}, {"issue": "COVID-19 outliers 2020-03..2021-12 (GDP yoy from -16.2 to +18.5 %)", "handling": "training pairs touching this window excluded from model estimation; sensitivity without exclusion reported (sens_* keys)"}, {"issue": "current-vintage data, not real-time vintages; future months present in file", "handling": "per-origin truncation to GDP<=t-2, employment<=t-1; no real-time validity claimed"}, {"issue": "GDP missing for 2026-08 and 2026-09 (not yet published)", "handling": "not needed: latest target is 2026-07"}, ], "n_input_rows": {"canada_gdp_employment.csv": int(n_raw)}, "n_analysed": int(n_analysed), "n_analysed_h3_last_origin": int(len(train_pairs(tL, 3, True)[0])), "random_seed": None, "dm_stat_h1": dm1, "dm_pvalue_h1": p1, "dm_stat_h3": dm3, "dm_pvalue_h3": p3, **{f"sens_nocovidexcl_{k}": v for k, v in sens.items() if k.startswith(("rmse_model", "mae_model", "rel_"))}, } with open(f"{OUT}/results.json", "w") as fh: json.dump(results, fh, indent=2) for k in ["rmse_naive_h1", "rmse_model_h1", "rel_improvement_rmse_h1", "rmse_naive_h3", "rmse_model_h3", "rel_improvement_rmse_h3"]: print(k, round(results[k], 4))