1"""F3-conf-01: pseudo-real-time forecasts of Canadian monthly GDP y/y growth.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"9MODEL = "ardl_emp"10COVID = (pd.Timestamp("2020-03-01"), pd.Timestamp("2021-06-01"))11HORIZONS = [1, 3]1213d = pd.read_csv(DATA, parse_dates=["month"]).sort_values("month").reset_index(drop=True)14n_raw = len(d)15assert not d.month.duplicated().any()16months = list(d.month)17idx = {m: i for i, m in enumerate(months)}18y = d.gdp_yoy.to_numpy(float); e = d.emp_yoy.to_numpy(float)1920def in_covid(i):21 return COVID[0] <= months[i] <= COVID[1]2223origins = pd.date_range("2022-05-01", "2026-04-01", freq="MS")24assert len(origins) == 4825fmt = lambda m: pd.Timestamp(m).strftime("%Y-%m-%d")2627def train_rows(t, h):28 """Training pairs for origin index t, horizon h: target y[s+h], regressors y[s-2], e[s-1].29 Constraint: s+h <= t-2 (target observed), s-1 <= t-1, s-2 >= 0; all values non-missing."""30 rows, n_covid = [], 031 for s in range(2, t - 2 - h + 1):32 if np.isnan(y[s + h]) or np.isnan(y[s - 2]) or np.isnan(e[s - 1]):33 continue34 if in_covid(s + h) or in_covid(s - 2) or in_covid(s - 1):35 n_covid += 1; continue36 rows.append(s)37 return rows, n_covid3839recs, last_sample = [], {}40for o in origins:41 t = idx[o]42 gdp_last, emp_last = t - 2, t - 143 for h in HORIZONS:44 actual = y[t + h]45 assert not np.isnan(actual)46 # naive47 recs.append(dict(origin_date=fmt(o), target_date=fmt(months[t + h]), horizon=h,48 forecast=y[gdp_last], actual=actual, model="naive",49 last_target_date_used=fmt(months[gdp_last]), last_predictor_date_used=""))50 # direct ARDL: y_{s+h} = a + b*y_{s-2} + c*e_{s-1}51 rows, n_cov = train_rows(t, h)52 S = np.array(rows)53 X = np.column_stack([np.ones(len(S)), y[S - 2], e[S - 1]])54 Y = y[S + h]55 beta, *_ = np.linalg.lstsq(X, Y, rcond=None)56 fc = float(beta @ np.array([1.0, y[gdp_last], e[emp_last]]))57 lt = max(S.max() + h, gdp_last); lp = max(S.max() - 1, emp_last)58 assert lt <= t - 2 and lp <= t - 159 recs.append(dict(origin_date=fmt(o), target_date=fmt(months[t + h]), horizon=h,60 forecast=fc, actual=actual, model=MODEL,61 last_target_date_used=fmt(months[lt]), last_predictor_date_used=fmt(months[lp])))62 if o == origins[-1]:63 last_sample[h] = dict(n=len(S), n_covid=n_cov, first=fmt(months[S.min() + h]), last=fmt(months[S.max() + h]),64 beta=[float(b) for b in beta])6566f = pd.DataFrame(recs)67f["forecast"] = f.forecast.round(6); f["error"] = f.forecast - f.actual68f.drop(columns="error").to_csv(os.path.join(OUT, "forecasts.csv"), index=False)6970def m(model, h):71 err = f[(f.model == model) & (f.horizon == h)].error.to_numpy()72 assert len(err) == 4873 return float(np.sqrt(np.mean(err ** 2))), float(np.mean(np.abs(err)))7475res = {"task_id": "F3-conf-01", "model_name": MODEL,76 "model_description": ("Direct h-step OLS regression of gdp_yoy(t+h) on a constant, the latest available gdp_yoy (t-2) "77 "and the latest available LFS emp_yoy (t-1), re-estimated at each origin on an expanding window "78 "starting 1998 using only pairs whose target is observed at the origin, excluding pairs touching 2020-03..2021-06 (COVID)."),79 "n_origins": int(len(origins)), "first_origin": fmt(origins[0]), "last_origin": fmt(origins[-1]),80 "last_target_date_available": fmt(months[idx[origins[-1]] - 2]),81 "last_predictor_date_available": fmt(months[idx[origins[-1]] - 1])}82for h in HORIZONS:83 rn, an = m("naive", h); rm, am = m(MODEL, h)84 res.update({f"rmse_naive_h{h}": rn, f"mae_naive_h{h}": an, f"rmse_model_h{h}": rm, f"mae_model_h{h}": am})85for h in HORIZONS:86 res[f"rel_improvement_rmse_h{h}"] = 1 - res[f"rmse_model_h{h}"] / res[f"rmse_naive_h{h}"]87from scipy import stats88for h in HORIZONS:89 en = f[(f.model == "naive") & (f.horizon == h)].error.to_numpy(); em = f[(f.model == MODEL) & (f.horizon == h)].error.to_numpy()90 dl = en ** 2 - em ** 2; T = len(dl); L = h + 1 # forecast gap is h+2 months -> MA(h+1) overlap91 u = dl - dl.mean(); lrv = u @ u / T92 for k in range(1, L + 1):93 lrv += 2 * (1 - k / (L + 1)) * (u[k:] @ u[:-k]) / T94 dm = dl.mean() / np.sqrt(lrv / T)95 hln = dm * np.sqrt((T + 1 - 2 * (h + 2) + (h + 2) * (h + 1) / T) / T) # Harvey-Leybourne-Newbold correction96 res[f"dm_stat_h{h}"] = float(hln); res[f"dm_pvalue_h{h}"] = float(2 * stats.t.sf(abs(hln), T - 1))97res.update({"feasibility": "complete", "feasibility_reason": None,98 "data_issues": [99 {"issue": "gdp_yoy/emp_yoy populated for 1998-01..1998-12 although README says empty for first 12 months (1997 levels not in file; cannot be verified)",100 "handling": "used as given; only affects the earliest training observations"},101 {"issue": "GDP columns missing for 2026-08 and 2026-09 (not yet published)", "handling": "expected; never needed (latest target 2026-07)"},102 {"issue": "COVID-19 extreme y/y swings 2020-03..2021-06 (gdp_yoy -16.2 to +18.5)",103 "handling": "training pairs whose target or regressor month is in 2020-03..2021-06 excluded from model estimation; naive unaffected"},104 {"issue": "data are current vintage, not first releases", "handling": "exercise simulates publication timing only; no real-time validity claimed"}],105 "n_input_rows": {"canada_gdp_employment.csv": int(n_raw)},106 "n_analysed": int(last_sample[1]["n"]),107 "n_analysed_h3": int(last_sample[3]["n"]),108 "last_estimation_sample": {f"h{h}": last_sample[h] for h in HORIZONS},109 "random_seed": None})110with open(os.path.join(OUT, "results.json"), "w") as fh:111 json.dump(res, fh, indent=2)112113# exclusion ledger for the last estimation sample (origin 2026-04-01, h=1)114t = idx[origins[-1]]; h = 1115n_avail = t - 2 + 1 # rows with month <= t-2 (GDP available)116n_pairs = len(range(2, t - 2 - h + 1))117ls = last_sample[1]118ex = pd.DataFrame([119 dict(step=1, rule="keep months <= last_origin-2 (2026-02) with GDP published at last origin", n_before=n_raw, n_after=n_avail,120 n_removed=n_raw - n_avail, rationale="information-availability rule (GDP <= t-2)"),121 dict(step=2, rule="form direct h=1 training pairs (need y(s-2), target y(s+1) <= t-2)", n_before=n_avail, n_after=n_pairs,122 n_removed=n_avail - n_pairs, rationale="lag/lead construction loses rows at sample edges"),123 dict(step=3, rule="drop pairs whose target or regressor month is in 2020-03..2021-06", n_before=n_pairs, n_after=ls["n"],124 n_removed=n_pairs - ls["n"], rationale="COVID outliers would dominate OLS; pre-specified")])125assert ex.n_removed.sum() == n_raw - ls["n"]126ex.to_csv(os.path.join(OUT, "exclusions.csv"), index=False)127print(json.dumps({k: v for k, v in res.items() if k.startswith(("rmse", "mae", "rel", "n_"))}, indent=1))128