1"""F3-conf-01: forecasting Canadian monthly GDP y/y growth under publication lags.2Run from workspace root: ~/wp29/venv/bin/python out/analysis.py3Reads data/ only; writes out/results.json, out/forecasts.csv, out/exclusions.csv. Deterministic (no randomness)."""4import json, os5import numpy as np, pandas as pd67DATA = "data/canada_gdp_employment.csv"8OUT = "out"9H = [1, 3]10MODEL = "ardl_direct_emp"11COVID = (pd.Timestamp("2020-03-01"), pd.Timestamp("2021-12-01")) # pre-specified exclusion window1213raw = pd.read_csv(DATA, parse_dates=["month"])14n_raw = len(raw)15d = raw.set_index("month").sort_index()16months = d.index17y = d["gdp_yoy"]; e = d["emp_yoy"]1819def shift(m, k):20 return m + pd.DateOffset(months=k)2122def in_covid(m):23 return (m >= COVID[0]) & (m <= COVID[1])2425def design(t, h):26 """Training rows for origin t, horizon h (direct regression).27 Pseudo-origin s: target y[s+h], regressors const, y[s-2], e[s-1].28 Availability at origin t: y months <= t-2, e months <= t-1."""29 s = months30 tgt_m = s + pd.DateOffset(months=h)31 ym = s - pd.DateOffset(months=2); em = s - pd.DateOffset(months=1)32 steps = []33 n0 = len(s)34 keep = np.ones(n0, bool)35 k1 = keep & (tgt_m <= shift(t, -2)) & (ym <= shift(t, -2)) & (em <= shift(t, -1))36 steps.append(("availability", f"target month s+{h} <= t-2 and regressors observed at origin t", keep.sum(), k1.sum()))37 yt = y.reindex(tgt_m).values; yl = y.reindex(ym).values; el = e.reindex(em).values38 k2 = k1 & ~np.isnan(yt) & ~np.isnan(yl) & ~np.isnan(el)39 steps.append(("lags_in_sample", "target and lagged regressors present in file (drops first months lacking y[s-2]/e[s-1] and GDP-empty months)", k1.sum(), k2.sum()))40 cov = in_covid(tgt_m) | in_covid(ym) | in_covid(em)41 k3 = k2 & ~np.asarray(cov)42 steps.append(("covid_window", "drop rows whose target or regressor month lies in 2020-03..2021-12", k2.sum(), k3.sum()))43 X = np.column_stack([np.ones(k3.sum()), yl[k3], el[k3]])44 return X, yt[k3], tgt_m[k3], steps4546rows = []47origins = pd.date_range("2022-05-01", "2026-04-01", freq="MS")48assert len(origins) == 4849last_info = {}50ar_rows = []51for t in origins:52 yl_t = y.loc[shift(t, -2)]; el_t = e.loc[shift(t, -1)]53 for h in H:54 tm = shift(t, h)55 actual = y.loc[tm]56 rows.append(dict(origin_date=t, target_date=tm, horizon=h, forecast=yl_t, actual=actual, model="naive",57 last_target_date_used=shift(t, -2), last_predictor_date_used=None))58 X, Y, tms, steps = design(t, h)59 beta, *_ = np.linalg.lstsq(X, Y, rcond=None)60 f = float(beta @ np.array([1.0, yl_t, el_t]))61 rows.append(dict(origin_date=t, target_date=tm, horizon=h, forecast=f, actual=actual, model=MODEL,62 last_target_date_used=max(tms.max(), shift(t, -2)), last_predictor_date_used=shift(t, -1)))63 b_ar, *_ = np.linalg.lstsq(X[:, :2], Y, rcond=None) # diagnostic: same sample, no employment64 ar_rows.append(dict(origin_date=t, horizon=h, forecast=float(b_ar @ np.array([1.0, yl_t])), actual=actual))65 last_info[h] = dict(n=len(Y), steps=steps, beta=beta.tolist(), last_target=str(tms.max().date()))6667fc = pd.DataFrame(rows)68# information-rule checks69assert (fc.last_target_date_used <= fc.origin_date - pd.DateOffset(months=2)).all()70m = fc.model == MODEL71assert (fc.loc[m, "last_predictor_date_used"] <= fc.loc[m, "origin_date"] - pd.DateOffset(months=1)).all()72assert fc.actual.notna().all() and len(fc) == 1927374out = fc.copy()75for c in ["origin_date", "target_date", "last_target_date_used", "last_predictor_date_used"]:76 out[c] = out[c].apply(lambda v: "" if v is None or pd.isna(v) else pd.Timestamp(v).strftime("%Y-%m-%d"))77out["forecast"] = out["forecast"].round(6)78out.to_csv(os.path.join(OUT, "forecasts.csv"), index=False)7980def err(model, h):81 s = fc[(fc.model == model) & (fc.horizon == h)]82 u = (s.actual - s.forecast).values83 return float(np.sqrt(np.mean(u ** 2))), float(np.mean(np.abs(u)))8485res = {"task_id": "F3-conf-01", "model_name": MODEL,86 "model_description": ("Direct OLS regression per horizon h of gdp_yoy[s+h] on a constant, gdp_yoy[s-2] and emp_yoy[s-1], "87 "estimated at each origin on an expanding window from 1998 using only pairs whose target is published (<= t-2), "88 "excluding rows touching the pre-specified COVID window 2020-03..2021-12; forecast plugs in gdp_yoy[t-2] and emp_yoy[t-1]."),89 "n_origins": 48, "first_origin": "2022-05-01", "last_origin": "2026-04-01",90 "last_target_date_available": shift(origins[-1], -2).strftime("%Y-%m-%d"),91 "last_predictor_date_available": shift(origins[-1], -1).strftime("%Y-%m-%d")}92for h in H:93 rn, mn = err("naive", h); rm, mm = err(MODEL, h)94 res[f"rmse_naive_h{h}"] = rn; res[f"mae_naive_h{h}"] = mn95 res[f"rmse_model_h{h}"] = rm; res[f"mae_model_h{h}"] = mm96order = ["rmse_naive_h1", "mae_naive_h1", "rmse_naive_h3", "mae_naive_h3", "rmse_model_h1", "mae_model_h1", "rmse_model_h3", "mae_model_h3"]97res = {**res, **{k: res[k] for k in order}}98for h in H:99 res[f"rel_improvement_rmse_h{h}"] = 1 - res[f"rmse_model_h{h}"] / res[f"rmse_naive_h{h}"]100101# Diebold-Mariano (squared loss, HAC with h-1 lags, Harvey-Leybourne-Newbold correction) for uncertainty102from scipy import stats103for h in H:104 a = fc[(fc.model == "naive") & (fc.horizon == h)].sort_values("origin_date")105 b = fc[(fc.model == MODEL) & (fc.horizon == h)].sort_values("origin_date")106 dlt = (a.actual.values - a.forecast.values) ** 2 - (b.actual.values - b.forecast.values) ** 2107 T = len(dlt); dm = dlt - dlt.mean()108 L = max(h - 1, 0) + 2 # a few extra lags given overlapping y/y targets109 v = dm @ dm / T + 2 * sum((1 - l / (L + 1)) * (dm[l:] @ dm[:-l]) / T for l in range(1, L + 1))110 stat = dlt.mean() / np.sqrt(v / T)111 stat *= np.sqrt((T + 1 - 2 * h + h * (h - 1) / T) / T)112 res[f"dm_stat_h{h}"] = float(stat); res[f"dm_pvalue_two_sided_h{h}"] = float(2 * stats.t.sf(abs(stat), T - 1))113 res[f"dm_hac_lags_h{h}"] = L114115ar = pd.DataFrame(ar_rows)116for h in H:117 a = ar[ar.horizon == h]; u = (a.actual - a.forecast).values118 res[f"diag_rmse_ar_only_h{h}"] = float(np.sqrt(np.mean(u ** 2)))119 res[f"diag_mae_ar_only_h{h}"] = float(np.mean(np.abs(u)))120res["feasibility"] = "complete"; res["feasibility_reason"] = None121res["data_issues"] = [122 {"issue": "gdp_yoy is populated for 1998-01..1998-12 although data/README.md says it is empty for the first 12 months and the file has no 1997 levels to verify it.",123 "handling": "Values kept as given (they are the official target column); all gdp_yoy values computable from gdp_level in the file match to 4 decimals."},124 {"issue": "gdp_level/gdp_yoy missing for 2026-08 and 2026-09 (not yet published).", "handling": "Expected; never needed (latest target 2026-07)."},125 {"issue": "COVID-19 swings (gdp_yoy -16.2% in 2020-04, +18.5% in 2021-04) would dominate OLS estimates.",126 "handling": "Pre-specified exclusion of training rows whose target or regressor month lies in 2020-03..2021-12; naive forecast unaffected."},127 {"issue": "Current-vintage data, not real-time vintages.", "handling": "Exercise simulates publication timing only; no real-time validity claimed."}]128res["n_input_rows"] = {"canada_gdp_employment.csv": int(n_raw)}129res["n_analysed"] = int(last_info[1]["n"])130res["n_analysed_last_h3"] = int(last_info[3]["n"])131res["last_coefficients_h1"] = dict(zip(["const", "gdp_yoy_lag2", "emp_yoy_lag1"], last_info[1]["beta"]))132res["last_coefficients_h3"] = dict(zip(["const", "gdp_yoy_lag2", "emp_yoy_lag1"], last_info[3]["beta"]))133res["random_seed"] = None134135pd.DataFrame([dict(step=s, rule=r, n_before=int(a), n_after=int(b), n_removed=int(a - b),136 rationale="last estimation sample (origin 2026-04-01, h=1); rows = candidate pseudo-origins s, one per input month")137 for s, r, a, b in last_info[1]["steps"]]).to_csv(os.path.join(OUT, "exclusions.csv"), index=False)138with open(os.path.join(OUT, "results.json"), "w") as f:139 json.dump(res, f, indent=2)140print(json.dumps({k: v for k, v in res.items() if k.startswith(("rmse", "mae", "rel", "dm_", "n_"))}, indent=2))141