1"""F3-conf-01: pseudo-real-time forecasts of Canadian monthly GDP YoY growth (h=1,3).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, hashlib5from pathlib import Path6import numpy as np, pandas as pd78DATA = Path("data/canada_gdp_employment.csv"); OUT = Path("out")9MODEL = "ardl_emp_direct"10COVID = (pd.Timestamp("2020-03-01"), pd.Timestamp("2021-12-01")) # pre-specified exclusion window11H = (1, 3)1213raw = pd.read_csv(DATA, parse_dates=["month"])14n_raw = len(raw)15d = raw.set_index("month").sort_index()16assert d.index.is_unique17idx = d.index18def shift(m, k): return m + pd.DateOffset(months=k)19def ymd(m): return m.strftime("%Y-%m-%d")2021origins = pd.date_range("2022-05-01", "2026-04-01", freq="MS")22assert len(origins) == 482324def design(h):25 """Direct-regression rows indexed by target month m: y(m) on [1, y(m-h-2), e(m-h-1), e(m-h-2)].26 At origin t=m-h these regressors are gdp at t-2 and emp at t-1, t-2."""27 rows = pd.DataFrame(index=idx)28 rows["y"] = d["gdp_yoy"]29 rows["y_lag"] = d["gdp_yoy"].shift(h + 2)30 rows["e_lag1"] = d["emp_yoy"].shift(h + 1)31 rows["e_lag2"] = d["emp_yoy"].shift(h + 2)32 # months involved for COVID exclusion: target month m and regressor months m-h-2 .. m-h-133 m = rows.index34 first_month = m - pd.DateOffset(months=h + 2)35 covid = np.array([(a <= COVID[1]) and (b >= COVID[0]) for a, b in zip(first_month, m)])36 rows["covid"] = covid37 return rows3839ledger, excl = [], []40last_n = {}41for h in H:42 X = design(h)43 for t in origins:44 tgt = shift(t, h)45 actual = float(d.loc[tgt, "gdp_yoy"])46 g_av, e_av = shift(t, -2), shift(t, -1) # information boundaries at origin t47 # naive48 naive = float(d.loc[g_av, "gdp_yoy"])49 ledger.append(dict(origin_date=ymd(t), target_date=ymd(tgt), horizon=h, forecast=naive,50 actual=actual, model="naive", last_target_date_used=ymd(g_av),51 last_predictor_date_used=""))52 # model: training targets must be observed at origin -> target month <= t-253 s0 = X[X.index <= g_av]54 s1 = s0.dropna(subset=["y", "y_lag", "e_lag1", "e_lag2"])55 s2 = s1[~s1["covid"]]56 # guard: regressor months of all training rows lie before the target month <= t-257 A = np.column_stack([np.ones(len(s2)), s2[["y_lag", "e_lag1", "e_lag2"]].values])58 beta, *_ = np.linalg.lstsq(A, s2["y"].values, rcond=None)59 xnew = np.array([1.0, d.loc[g_av, "gdp_yoy"], d.loc[e_av, "emp_yoy"], d.loc[shift(t, -2), "emp_yoy"]])60 assert not np.isnan(xnew).any()61 fc = float(xnew @ beta)62 ledger.append(dict(origin_date=ymd(t), target_date=ymd(tgt), horizon=h, forecast=fc,63 actual=actual, model=MODEL, last_target_date_used=ymd(g_av),64 last_predictor_date_used=ymd(e_av)))65 if t == origins[-1]:66 last_n[h] = len(s2)67 excl += [68 dict(step=f"h{h}_1", rule=f"origin {ymd(t)}, h={h}: keep rows with target month <= t-2 ({ymd(g_av)}) (information rule)",69 n_before=n_raw, n_after=len(s0), n_removed=n_raw - len(s0),70 rationale="targets after t-2 not yet published at origin (includes 2 rows with empty GDP)"),71 dict(step=f"h{h}_2", rule="drop rows whose lagged regressors fall before 1998-01 (not in file)",72 n_before=len(s0), n_after=len(s1), n_removed=len(s0) - len(s1),73 rationale="lagged regressors unavailable"),74 dict(step=f"h{h}_3", rule="drop rows whose target or regressor months fall in 2020-03..2021-12",75 n_before=len(s1), n_after=len(s2), n_removed=len(s1) - len(s2),76 rationale="pre-specified COVID outlier exclusion; extreme swings would dominate OLS"),77 ]7879fc = pd.DataFrame(ledger)80fc["forecast"] = fc["forecast"].round(6)81fc = fc.sort_values(["model", "horizon", "origin_date"]).reset_index(drop=True)82fc.to_csv(OUT / "forecasts.csv", index=False)83pd.DataFrame(excl).to_csv(OUT / "exclusions.csv", index=False)8485def err(model, h):86 s = fc[(fc.model == model) & (fc.horizon == h)]87 e = s.forecast - s.actual88 return float(np.sqrt(np.mean(e ** 2))), float(np.mean(np.abs(e)))8990res = dict(task_id="F3-conf-01", model_name=MODEL,91 model_description=("Direct h-step OLS, separately per horizon: gdp_yoy(t+h) on a constant, gdp_yoy(t-2), "92 "emp_yoy(t-1) and emp_yoy(t-2), re-estimated at each origin on an expanding window of all target months "93 "from 1998 up to t-2 (training rows with target or regressor months in 2020-03..2021-12 excluded)."),94 n_origins=len(origins), first_origin=ymd(origins[0]), last_origin=ymd(origins[-1]),95 last_target_date_available=ymd(shift(origins[-1], -2)),96 last_predictor_date_available=ymd(shift(origins[-1], -1)))97for h in H:98 rn, mn = err("naive", h); rm, mm = err(MODEL, h)99 res.update({f"rmse_naive_h{h}": round(rn, 6), f"mae_naive_h{h}": round(mn, 6),100 f"rmse_model_h{h}": round(rm, 6), f"mae_model_h{h}": round(mm, 6)})101for h in H:102 res[f"rel_improvement_rmse_h{h}"] = round(1 - res[f"rmse_model_h{h}"] / res[f"rmse_naive_h{h}"], 6)103res.update(feasibility="complete", feasibility_reason=None,104 data_issues=[105 {"issue": "gdp_level/gdp_yoy empty for 2026-08 and 2026-09 (not yet published)",106 "handling": "never needed: latest target is 2026-07 and latest GDP used is 2026-02"},107 {"issue": "README says gdp_yoy/emp_yoy are empty for the first 12 months, but the file has values from 1998-01 (computed from pre-1998 levels not in the file)",108 "handling": "values used as given; recomputation from levels matches exactly wherever possible (from 1999-01)"},109 {"issue": "COVID-19 swings 2020-03..2021-12 (gdp_yoy from -16.2 to +18.5 %)",110 "handling": "training rows touching this window excluded from model estimation (pre-specified); evaluation period 2022-06..2026-07 unaffected"},111 {"issue": "data are current vintage, not first releases",112 "handling": "exercise simulates publication timing only; no real-time validity claimed"}],113 n_input_rows={"canada_gdp_employment.csv": n_raw},114 n_analysed=int(last_n[1]),115 n_analysed_h3=int(last_n[3]),116 random_seed=None,117 input_sha256=hashlib.sha256(DATA.read_bytes()).hexdigest())118# --- Uncertainty: Diebold-Mariano (squared loss), Newey-West variance with h+1 lags119# (errors overlap: the naive information gap is h+2 months), Harvey-Leybourne-Newbold small-sample correction.120from scipy import stats121for h in H:122 a = fc[(fc.model == "naive") & (fc.horizon == h)].set_index("origin_date")123 b = fc[(fc.model == MODEL) & (fc.horizon == h)].set_index("origin_date")124 dl = ((b.forecast - b.actual) ** 2 - (a.forecast - a.actual) ** 2).values125 T, L = len(dl), h + 1126 u = dl - dl.mean()127 lrv = u @ u / T + 2 * sum((1 - k / (L + 1)) * (u[k:] @ u[:-k]) / T for k in range(1, L + 1))128 dm = dl.mean() / np.sqrt(lrv / T)129 hln = dm * np.sqrt((T + 1 - 2 * h + h * (h - 1) / T) / T)130 res[f"dm_stat_h{h}"] = round(float(hln), 4)131 res[f"dm_pvalue_two_sided_h{h}"] = round(float(2 * stats.t.cdf(-abs(hln), df=T - 1)), 4)132res["dm_note"] = "negative stat = model has lower squared loss than naive; HLN-corrected, NW lags h+1, t(T-1) reference"133(OUT / "results.json").write_text(json.dumps(res, indent=2))134print(json.dumps({k: v for k, v in res.items() if k.startswith(("rmse", "mae", "rel", "n_"))}, indent=1))135