1"""F3-conf-01: Canadian monthly GDP y/y growth forecasts 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. Deterministic (no randomness)."""4import json, os5import numpy as np, pandas as pd6from scipy import stats78DATA = "data/canada_gdp_employment.csv"9OUT = "out"10MODEL = "ardl_emp_direct"11COVID = (pd.Timestamp("2020-03-01"), pd.Timestamp("2021-12-01")) # pre-specified training exclusion window12HORIZONS = (1, 3)1314raw = pd.read_csv(DATA, parse_dates=["month"])15n_raw = len(raw)16d = raw.set_index("month").sort_index()17assert d.index.is_unique18origins = pd.date_range("2022-05-01", "2026-04-01", freq="MS")19assert len(origins) == 4820M = lambda ts, k: ts + pd.DateOffset(months=k)2122def design(gdp, emp, h, exclude_covid):23 """Direct regression rows: target y(m) on [1, y(m-h-2), emp(m-h-1)], using only the masked series."""24 rows = []25 for m in gdp.dropna().index:26 ml_y, ml_e = M(m, -(h + 2)), M(m, -(h + 1))27 if ml_y not in gdp.index or ml_e not in emp.index or pd.isna(gdp.get(ml_y)) or pd.isna(emp.get(ml_e)):28 continue29 if exclude_covid and any(COVID[0] <= x <= COVID[1] for x in (m, ml_y, ml_e)):30 continue31 rows.append((m, gdp[m], gdp[ml_y], emp[ml_e]))32 return pd.DataFrame(rows, columns=["m", "y", "ylag", "elag"]).set_index("m")3334def fit_forecast(t, h, exclude_covid=True):35 gdp = d["gdp_yoy"][d.index <= M(t, -2)] # GDP available: months <= t-236 emp = d["emp_yoy"][d.index <= M(t, -1)] # LFS available: months <= t-137 X = design(gdp, emp, h, exclude_covid)38 A = np.column_stack([np.ones(len(X)), X["ylag"], X["elag"]])39 beta, *_ = np.linalg.lstsq(A, X["y"].values, rcond=None)40 x_new = np.array([1.0, gdp[M(t, -2)], emp[M(t, -1)]])41 return float(x_new @ beta), X, gdp.index.max(), emp.index.max()4243recs, sens, n_last = [], [], {}44for t in origins:45 for h in HORIZONS:46 tgt = M(t, h)47 actual = float(d.loc[tgt, "gdp_yoy"])48 assert not np.isnan(actual)49 last_g = M(t, -2)50 recs.append(dict(origin_date=t, target_date=tgt, horizon=h, forecast=float(d.loc[last_g, "gdp_yoy"]),51 actual=actual, model="naive", last_target_date_used=last_g, last_predictor_date_used=pd.NaT))52 f, X, lg, le = fit_forecast(t, h, True)53 assert lg <= M(t, -2) and le <= M(t, -1) and X.index.max() <= M(t, -2)54 recs.append(dict(origin_date=t, target_date=tgt, horizon=h, forecast=f, actual=actual, model=MODEL,55 last_target_date_used=lg, last_predictor_date_used=le))56 fs, *_ = fit_forecast(t, h, False)57 sens.append(dict(h=h, e=fs - actual))58 if t == origins[-1]:59 n_last[h] = len(X)60fc = pd.DataFrame(recs)61fc["error"] = fc["forecast"] - fc["actual"]6263def rmse(e): return float(np.sqrt(np.mean(np.square(e))))64def mae(e): return float(np.mean(np.abs(e)))6566def dm_test(e1, e2, h):67 """Diebold-Mariano (squared loss), HAC variance with h-1 Bartlett-free rectangular lags, Harvey et al. correction."""68 dl = e1 ** 2 - e2 ** 269 n = len(dl); dbar = dl.mean(); dc = dl - dbar70 g = [np.sum(dc[k:] * dc[:n - k]) / n for k in range(h)]71 v = (g[0] + 2 * sum(g[1:])) / n72 if v <= 0:73 v = g[0] / n74 stat = dbar / np.sqrt(v) * np.sqrt((n + 1 - 2 * h + h * (h - 1) / n) / n)75 return float(stat), float(2 * stats.t.sf(abs(stat), df=n - 1))7677res = {"task_id": "F3-conf-01", "model_name": MODEL,78 "model_description": ("Direct h-step OLS regression per horizon and origin of gdp_yoy(m) on a constant, gdp_yoy(m-h-2) "79 "and emp_yoy(m-h-1), estimated on an expanding window from 1998 using only training pairs whose target "80 "month is <= t-2 (and employment <= t-1), excluding pairs touching 2020-03..2021-12 (COVID); "81 "forecast plugs in gdp_yoy(t-2) and emp_yoy(t-1)."),82 "n_origins": len(origins), "first_origin": str(origins[0].date()), "last_origin": str(origins[-1].date()),83 "last_target_date_available": str(M(origins[-1], -2).date()),84 "last_predictor_date_available": str(M(origins[-1], -1).date())}85sub = {}86for h in HORIZONS:87 en = fc[(fc.horizon == h) & (fc.model == "naive")].sort_values("origin_date")["error"].values88 em = fc[(fc.horizon == h) & (fc.model == MODEL)].sort_values("origin_date")["error"].values89 res[f"rmse_naive_h{h}"], res[f"mae_naive_h{h}"] = rmse(en), mae(en)90 res[f"rmse_model_h{h}"], res[f"mae_model_h{h}"] = rmse(em), mae(em)91 sub[h] = (en, em)92for h in HORIZONS: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": "README states gdp_yoy is empty for the first 12 months, but gdp_yoy and emp_yoy are populated from 1998-01 (presumably from pre-1998 levels not in the file).",97 "handling": "Used the provided yoy values as given; they match recomputation from levels to <=5e-5 where verifiable."},98 {"issue": "gdp_level/gdp_yoy missing for 2026-08 and 2026-09 (not yet published).",99 "handling": "Expected; never needed (last target month 2026-07, last GDP used 2026-02)."},100 {"issue": "Extreme COVID-19 swings 2020-03..2021-06 and base effects through 2021 in both series.",101 "handling": "Pre-specified exclusion of training pairs touching 2020-03..2021-12 from the model; sensitivity without exclusion reported (sens_* keys). Naive baseline unaffected."},102 {"issue": "Data are the current vintage (revised), not first releases.",103 "handling": "Exercise simulates publication timing only; no real-time validity claimed."}],104 "n_input_rows": {"canada_gdp_employment.csv": n_raw},105 "n_analysed": int(n_last[1]), "n_analysed_last_origin_h3": int(n_last[3]),106 "random_seed": None})107for h in HORIZONS:108 en, em = sub[h]109 s, p = dm_test(em, en, h)110 res[f"dm_stat_h{h}"], res[f"dm_pvalue_h{h}"] = s, p111 se = np.array([r["e"] for r in sens if r["h"] == h])112 res[f"sens_rmse_model_nocovidexcl_h{h}"] = rmse(se)113 res[f"sens_rel_improvement_rmse_nocovidexcl_h{h}"] = 1 - rmse(se) / res[f"rmse_naive_h{h}"]114res = {k: (round(v, 10) if isinstance(v, float) else v) for k, v in res.items()}115116os.makedirs(OUT, exist_ok=True)117out = fc.drop(columns="error").copy()118for c in ["origin_date", "target_date", "last_target_date_used", "last_predictor_date_used"]:119 out[c] = out[c].dt.strftime("%Y-%m-%d")120out["last_predictor_date_used"] = out["last_predictor_date_used"].fillna("")121out["forecast"] = out["forecast"].round(8)122out.to_csv(f"{OUT}/forecasts.csv", index=False)123with open(f"{OUT}/results.json", "w") as fh:124 json.dump(res, fh, indent=2)125126# Exclusion ledger for the last estimation sample (last origin, h=1)127t = origins[-1]; h = 1128n1 = int((d.index <= M(t, -2)).sum())129gdp = d["gdp_yoy"][d.index <= M(t, -2)]; emp = d["emp_yoy"][d.index <= M(t, -1)]130n2 = len(design(gdp, emp, h, False)); n3 = n_last[1]131led = pd.DataFrame([132 ("1", "GDP month <= t-2 at last origin 2026-04 (publication lag); drops 2026-03..2026-09", n_raw, n1, "information-availability rule"),133 ("2", f"target month needs gdp_yoy(m-{h+2}) and emp_yoy(m-{h+1}) inside the file", n1, n2, "lags unavailable for first months of file"),134 ("3", "drop pairs with target/regressor month in 2020-03..2021-12", n2, n3, "pre-specified COVID outlier exclusion (model only)")],135 columns=["step", "rule", "n_before", "n_after", "rationale"])136led.insert(4, "n_removed", led.n_before - led.n_after)137led = led[["step", "rule", "n_before", "n_after", "n_removed", "rationale"]]138led.to_csv(f"{OUT}/exclusions.csv", index=False)139assert n_raw - led.n_removed.sum() == res["n_analysed"]140print(json.dumps({k: res[k] for k in res if k.startswith(("rmse", "mae", "rel", "dm", "sens", "n_an"))}, indent=1))141