1"""F3-conf-01: GDP y/y forecasts under publication lags. Run from workspace root:2 ~/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 pd67DATA = "data/canada_gdp_employment.csv"8OUT = "out"9MODEL = "dir_ols_emp"10COVID = (pd.Timestamp("2020-03-01"), pd.Timestamp("2021-12-01"))1112d = pd.read_csv(DATA, parse_dates=["month"]).sort_values("month").reset_index(drop=True)13n_raw = len(d)14assert d.month.is_unique15d = d.set_index("month")16M = lambda t, k: t + pd.DateOffset(months=k)17def ds(t): return t.strftime("%Y-%m-%d")18def in_covid(m): return (m >= COVID[0]) & (m <= COVID[1])1920origins = pd.date_range("2022-05-01", "2026-04-01", freq="MS")21assert len(origins) == 482223def build_rows(t, h):24 """Direct regression rows y_s ~ 1 + y_{s-h-2} + x_{s-h-1} using only info at origin t."""25 gdp = d.loc[d.index <= M(t, -2), "gdp_yoy"].dropna() # GDP known up to t-226 emp = d.loc[d.index <= M(t, -1), "emp_yoy"].dropna() # emp known up to t-127 s = gdp.index28 df = pd.DataFrame({"y": gdp.values,29 "ylag": gdp.reindex(s - pd.DateOffset(months=h + 2)).values,30 "xlag": emp.reindex(s - pd.DateOffset(months=h + 1)).values}, index=s)31 n_avail = len(df)32 df = df.dropna()33 n_lagged = len(df)34 bad = in_covid(df.index) | in_covid(df.index - pd.DateOffset(months=h + 2)) | in_covid(df.index - pd.DateOffset(months=h + 1))35 df = df[~bad]36 return df, gdp, emp, n_avail, n_lagged3738rows = []39last_info = {}40for t in origins:41 for h in (1, 3):42 tgt = M(t, h)43 actual = float(d.loc[tgt, "gdp_yoy"])44 df, gdp, emp, n_avail, n_lagged = build_rows(t, h)45 y_last, x_last = gdp.index.max(), emp.index.max()46 assert y_last <= M(t, -2) and x_last <= M(t, -1)47 # naive48 rows.append(dict(origin_date=ds(t), target_date=ds(tgt), horizon=h, forecast=float(gdp.loc[M(t, -2)]),49 actual=actual, model="naive", last_target_date_used=ds(M(t, -2)), last_predictor_date_used=""))50 # model51 X = np.column_stack([np.ones(len(df)), df.ylag.values, df.xlag.values])52 beta, *_ = np.linalg.lstsq(X, df.y.values, rcond=None)53 xnew = np.array([1.0, gdp.loc[M(t, -2)], emp.loc[M(t, -1)]])54 fc = float(xnew @ beta)55 rows.append(dict(origin_date=ds(t), target_date=ds(tgt), horizon=h, forecast=fc, actual=actual, model=MODEL,56 last_target_date_used=ds(max(df.index.max(), M(t, -2))),57 last_predictor_date_used=ds(M(t, -1))))58 last_info[h] = dict(n=len(df), n_avail=n_avail, n_lagged=n_lagged, beta=beta.tolist())5960fc = pd.DataFrame(rows)61fc.to_csv(os.path.join(OUT, "forecasts.csv"), index=False, float_format="%.6f")6263def err(model, h):64 e = fc[(fc.model == model) & (fc.horizon == h)]65 e = e.forecast - e.actual66 assert len(e) == 4867 return float(np.sqrt(np.mean(e ** 2))), float(np.mean(np.abs(e)))6869res = {"task_id": "F3-conf-01", "model_name": MODEL,70 "model_description": ("Direct h-step OLS of gdp_yoy(s) on a constant, gdp_yoy(s-h-2) and emp_yoy(s-h-1), "71 "re-estimated at each origin on an expanding window from 1998 with targets <= t-2, "72 "excluding rows whose target or regressor month lies in 2020-03..2021-12; "73 "forecast plugs in gdp_yoy(t-2) and emp_yoy(t-1)."),74 "n_origins": int(len(origins)), "first_origin": ds(origins[0]), "last_origin": ds(origins[-1]),75 "last_target_date_available": ds(M(origins[-1], -2)),76 "last_predictor_date_available": ds(M(origins[-1], -1))}77for m, tag in (("naive", "naive"), (MODEL, "model")):78 for h in (1, 3):79 r, a = err(m, h)80 res[f"rmse_{tag}_h{h}"], res[f"mae_{tag}_h{h}"] = r, a81for h in (1, 3):82 res[f"rel_improvement_rmse_h{h}"] = 1 - res[f"rmse_model_h{h}"] / res[f"rmse_naive_h{h}"]83res.update({"feasibility": "complete", "feasibility_reason": None,84 "data_issues": [85 {"issue": "gdp_level/gdp_yoy missing for 2026-08 and 2026-09 (not yet published)",86 "handling": "Not needed: latest target is 2026-07; rows never used as GDP information."},87 {"issue": "gdp_yoy and emp_yoy for 1998 depend on 1997 levels absent from the file and cannot be re-verified (all later yoy values match levels exactly)",88 "handling": "Accepted as given."},89 {"issue": "COVID-19 shock 2020-03..2021-12 with yoy values between -16 and +19 percent (extreme leverage)",90 "handling": "Regression rows touching this window excluded from model estimation (pre-specified); naive unaffected."},91 {"issue": "Current-vintage (revised) data, not historical first releases",92 "handling": "Only publication timing simulated; no real-time validity claimed."}],93 "n_input_rows": {"canada_gdp_employment.csv": int(n_raw)},94 "n_analysed": int(last_info[1]["n"]), "random_seed": None})95# Uncertainty: Diebold-Mariano (squared loss, HAC variance with h-1 Newey-West lags, Harvey et al. small-sample correction)96from scipy import stats97for h in (1, 3):98 a = fc[(fc.model == "naive") & (fc.horizon == h)].sort_values("origin_date")99 b = fc[(fc.model == MODEL) & (fc.horizon == h)].sort_values("origin_date")100 dd = (a.forecast.values - a.actual.values) ** 2 - (b.forecast.values - b.actual.values) ** 2101 n = len(dd); dm = dd - dd.mean()102 lrv = np.mean(dm ** 2) + 2 * sum((1 - k / h) * np.mean(dm[k:] * dm[:-k]) for k in range(1, h))103 stat = dd.mean() / np.sqrt(lrv / n) * np.sqrt((n + 1 - 2 * h + h * (h - 1) / n) / n)104 res[f"dm_stat_h{h}"] = float(stat)105 res[f"dm_pvalue_two_sided_h{h}"] = float(2 * stats.t.sf(abs(stat), df=n - 1))106# Ablation (diagnostic only): same direct regression without the employment regressor107for h in (1, 3):108 e = []109 for t in origins:110 df, gdp, emp, *_ = build_rows(t, h)111 X = np.column_stack([np.ones(len(df)), df.ylag.values])112 beta, *_ = np.linalg.lstsq(X, df.y.values, rcond=None)113 e.append(beta[0] + beta[1] * gdp.loc[M(t, -2)] - d.loc[M(t, h), "gdp_yoy"])114 res[f"rmse_ablation_no_emp_h{h}"] = float(np.sqrt(np.mean(np.square(e))))115with open(os.path.join(OUT, "results.json"), "w") as f:116 json.dump(res, f, indent=2)117118# exclusion ledger for the last estimation sample (origin 2026-04, h=1)119li = last_info[1]120n_info = int((d.index <= M(origins[-1], -2)).sum())121steps = [("info_filter", "keep months <= t-2 (GDP published) at last origin 2026-04-01", n_raw, n_info,122 "publication-lag rule"),123 ("lag_availability", "drop rows lacking gdp_yoy(s-3) or emp_yoy(s-2) (start of sample)", n_info, li["n_lagged"],124 "regressors not in file"),125 ("covid_window", "drop rows with target or regressor month in 2020-03..2021-12", li["n_lagged"], li["n"],126 "pandemic outliers; pre-specified")]127assert li["n_avail"] == n_info128pd.DataFrame([dict(step=s, rule=r, n_before=a, n_after=b, n_removed=a - b, rationale=q) for s, r, a, b, q in steps]129 ).to_csv(os.path.join(OUT, "exclusions.csv"), index=False)130print(json.dumps({k: v for k, v in res.items() if k.startswith(("rmse", "mae", "rel"))}, indent=1))131print("last-origin coefficients [const, gdp lag, emp lag]: h1", np.round(last_info[1]["beta"], 4).tolist(),132 "h3", np.round(last_info[3]["beta"], 4).tolist())133