1"""F3-conf-01: GDP y/y growth forecasts under publication lags.2Run from workspace root: ~/wp29/venv/bin/python out/analysis.py3Naive: y(t-2). Model 'adl_direct_emp': direct OLS per horizon4 y(tau+h) = a + b*y(tau-2) + c*e(tau-1), training rows tau with tau+h <= t-2,5 excluding rows touching 2020-03..2021-12 (COVID); forecast plugs y(t-2), e(t-1).6No randomness is used.7"""8import json, os9import numpy as np, pandas as pd1011ROOT = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))12DATA = os.path.join(ROOT, "data", "canada_gdp_employment.csv")13OUT = os.path.join(ROOT, "out")14MODEL = "adl_direct_emp"15COVID = (pd.Timestamp("2020-03-01"), pd.Timestamp("2021-12-01"))1617d = pd.read_csv(DATA, parse_dates=["month"]).sort_values("month").reset_index(drop=True)18n_raw = len(d)19d = d.set_index("month")20months = d.index21origins = pd.date_range("2022-05-01", "2026-04-01", freq="MS")22H = [1, 3]23M = lambda m, k: m + pd.DateOffset(months=k)24fmt = lambda x: x.strftime("%Y-%m-%d")25in_covid = lambda m: (m >= COVID[0]) & (m <= COVID[1])262728def training_rows(t, h, return_ledger=False, exclude_covid=True):29 """Training rows available at origin t for horizon h (vectorised)."""30 gdp_last = M(t, -2) # latest GDP month observed at origin t31 emp_last = M(t, -1) # latest employment month observed at origin t32 y = d["gdp_yoy"].where(months <= gdp_last) # information truncation33 e = d["emp_yoy"].where(months <= emp_last)34 df = pd.DataFrame({"y_tgt": y.shift(-h), "y_lag": y.shift(2), "e_lag": e.shift(1)}, index=months)35 ledger = []36 n0 = len(df)37 df1 = df.dropna()38 ledger.append(("drop_unavailable", f"target month tau+{h} <= {fmt(gdp_last)}, y(tau-2) and e(tau-1) observed",39 n0, len(df1), "information-availability rule / series edges"))40 tau = df1.index41 bad = in_covid(tau) | in_covid(tau + pd.DateOffset(months=h)) | in_covid(tau - pd.DateOffset(months=2)) | in_covid(tau - pd.DateOffset(months=1))42 df2 = df1[~np.asarray(bad)] if exclude_covid else df143 ledger.append(("exclude_covid", "drop rows whose target/predictor months fall in 2020-03..2021-12",44 len(df1), len(df2), "COVID outliers dominate OLS; decided with pre-2022 info"))45 # sanity: every used month respects the rule46 assert (df2.index + pd.DateOffset(months=h)).max() <= gdp_last47 return (df2, ledger) if return_ledger else df2484950rows = []51n_analysed = None52for t in origins:53 y_t2 = d.at[M(t, -2), "gdp_yoy"]54 e_t1 = d.at[M(t, -1), "emp_yoy"]55 for h in H:56 tgt = M(t, h)57 actual = d.at[tgt, "gdp_yoy"]58 rows.append(dict(origin_date=fmt(t), target_date=fmt(tgt), horizon=h, forecast=y_t2, actual=actual,59 model="naive", last_target_date_used=fmt(M(t, -2)), last_predictor_date_used=""))60 tr = training_rows(t, h)61 X = np.column_stack([np.ones(len(tr)), tr["y_lag"], tr["e_lag"]])62 beta, *_ = np.linalg.lstsq(X, tr["y_tgt"].values, rcond=None)63 fc = float(beta @ np.array([1.0, y_t2, e_t1]))64 rows.append(dict(origin_date=fmt(t), target_date=fmt(tgt), horizon=h, forecast=fc, actual=actual,65 model=MODEL, last_target_date_used=fmt(M(t, -2)), last_predictor_date_used=fmt(M(t, -1))))6667fc = pd.DataFrame(rows)68assert len(fc) == 192 and fc["actual"].notna().all()69fc.to_csv(os.path.join(OUT, "forecasts.csv"), index=False, float_format="%.6f")7071# exclusion ledger for the last estimation sample (last origin, h=1); n_raw = all file rows72t_last = origins[-1]73tr_last, led = training_rows(t_last, 1, return_ledger=True)74n_analysed = len(tr_last)75assert led[0][2] == n_raw76ex = pd.DataFrame([dict(step=s, rule=r, n_before=a, n_after=b, n_removed=a - b, rationale=w) for s, r, a, b, w in led])77ex.to_csv(os.path.join(OUT, "exclusions.csv"), index=False)78n_last_h3 = len(training_rows(t_last, 3))7980res = {"task_id": "F3-conf-01", "model_name": MODEL,81 "model_description": ("Direct (per-horizon) OLS regression of gdp_yoy(t+h) on a constant, gdp_yoy(t-2) and emp_yoy(t-1), "82 "re-estimated at each origin on an expanding window of training rows (from 1998-03) whose target month is <= t-2, "83 "excluding rows touching 2020-03..2021-12."),84 "n_origins": int(len(origins)), "first_origin": fmt(origins[0]), "last_origin": fmt(t_last),85 "last_target_date_available": fmt(M(t_last, -2)), "last_predictor_date_available": fmt(M(t_last, -1))}86for h in H:87 for m, lab in [("naive", "naive"), (MODEL, "model")]:88 s = fc[(fc.horizon == h) & (fc.model == m)]89 err = s["forecast"] - s["actual"]90 res[f"rmse_{lab}_h{h}"] = float(np.sqrt(np.mean(err ** 2)))91 res[f"mae_{lab}_h{h}"] = float(np.mean(np.abs(err)))92for h in H:93 res[f"rel_improvement_rmse_h{h}"] = 1 - res[f"rmse_model_h{h}"] / res[f"rmse_naive_h{h}"]9495# Uncertainty: Diebold-Mariano test (squared loss), HAC variance with h-1 lags (Newey-West), HLN small-sample correction96from scipy import stats97dm = {}98for h in H:99 en = fc[(fc.horizon == h) & (fc.model == "naive")].sort_values("origin_date")100 em = fc[(fc.horizon == h) & (fc.model == MODEL)].sort_values("origin_date")101 dl = (em["forecast"].values - em["actual"].values) ** 2 - (en["forecast"].values - en["actual"].values) ** 2102 T = len(dl); dbar = dl.mean(); u = dl - dbar103 L = h + 1 # Newey-West lags h+1 since effective GDP horizon is h+2 months104 v = np.mean(u * u) + 2 * sum((1 - k / (L + 1)) * np.mean(u[k:] * u[:-k]) for k in range(1, L + 1))105 stat = dbar / np.sqrt(v / T)106 hk = h + 2107 stat *= np.sqrt((T + 1 - 2 * hk + hk * (hk - 1) / T) / T)108 p = 2 * stats.t.sf(abs(stat), df=T - 1)109 dm[f"h{h}"] = {"mean_loss_diff_model_minus_naive": float(dbar), "dm_stat_hln": float(stat), "p_value_two_sided": float(p), "nw_lags": L}110res["dm_test"] = dm111res["n_training_last_origin_h3"] = int(n_last_h3)112# Sensitivity: same model without the COVID exclusion (not the headline model)113sens = {}114for h in H:115 errs = []116 for t in origins:117 tr = training_rows(t, h, exclude_covid=False)118 X = np.column_stack([np.ones(len(tr)), tr["y_lag"], tr["e_lag"]])119 b, *_ = np.linalg.lstsq(X, tr["y_tgt"].values, rcond=None)120 f_ = float(b @ np.array([1.0, d.at[M(t, -2), "gdp_yoy"], d.at[M(t, -1), "emp_yoy"]]))121 errs.append(f_ - d.at[M(t, h), "gdp_yoy"])122 errs = np.array(errs)123 sens[f"rmse_h{h}"] = float(np.sqrt(np.mean(errs ** 2)))124 sens[f"rel_improvement_rmse_h{h}"] = 1 - sens[f"rmse_h{h}"] / res[f"rmse_naive_h{h}"]125res["sensitivity_no_covid_exclusion"] = sens126res.update({"feasibility": "complete", "feasibility_reason": None,127 "data_issues": [128 {"issue": "README states gdp_yoy is empty for the first 12 months, but the file contains yoy values for 1998 (both series).",129 "handling": "Used as given; yoy recomputed from levels agrees to rounding (<5e-5) where possible."},130 {"issue": "gdp columns missing for 2026-08 and 2026-09 (not yet published).",131 "handling": "Not needed: latest target (2026-07) and latest GDP input (2026-02) are present."},132 {"issue": "COVID-19 outliers (gdp_yoy from -16.2 to +18.5 in 2020-2021).",133 "handling": "Model training rows touching 2020-03..2021-12 excluded; naive unaffected; see out/exclusions.csv."},134 {"issue": "Current-vintage data (revised), not first releases.",135 "handling": "Exercise simulates publication timing only; no real-time validity claimed."}],136 "n_input_rows": {"canada_gdp_employment.csv": int(n_raw)},137 "n_analysed": int(n_analysed), "random_seed": None})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", "n_"))}, indent=2))141print(json.dumps(dm, indent=2))142