"""Real-time forecasts of US real GDP growth (F3-conf-02). Run from workspace root: ~/wp29/venv/bin/python out/analysis.py Reads only data/; writes out/results.json, out/forecasts.csv, out/exclusions.csv. Deterministic: no randomness used.""" import json, os, glob, hashlib import numpy as np, pandas as pd from scipy import stats DATA, OUT = "data", "out" ORIGINS = ["2024-01-31", "2024-04-30", "2024-07-31", "2024-10-31", "2025-01-31", "2025-04-30", "2025-07-31", "2025-12-31"] MODEL = "ar1_post1985_exCovid" SAMPLE_START = pd.Timestamp("1985-01-01") COVID_ROWS = [pd.Timestamp(x) for x in ("2020-04-01", "2020-07-01", "2020-10-01")] def read_series(path): d = pd.read_csv(path) d.columns = ["date", "value"] d["date"] = pd.to_datetime(d["date"]) return d def ar1_design(s): """s: DataFrame(date,value) from one vintage. Returns y, x, and ledger counts.""" n0 = len(s) w = s[s.date >= SAMPLE_START].reset_index(drop=True) n1 = len(w) rows = pd.DataFrame({"date": w.date.iloc[1:].values, "y": w.value.iloc[1:].values, "ylag": w.value.iloc[:-1].values}) # consecutive-quarter check assert (pd.to_datetime(rows.date).values - w.date.iloc[:-1].values).astype("timedelta64[D]").astype(int).max() <= 92 n2 = len(rows) rows = rows[~rows.date.isin(COVID_ROWS)].reset_index(drop=True) n3 = len(rows) return rows, (n0, n1, n2, n3) def fit_ar1(rows): X = np.column_stack([np.ones(len(rows)), rows.ylag.values]) beta, *_ = np.linalg.lstsq(X, rows.y.values, rcond=None) return beta # c, phi def dm_hln(e_model, e_naive, h): """Diebold-Mariano test (squared loss) with Harvey-Leybourne-Newbold small-sample correction; long-run variance uses autocovariances up to lag h-1 (rectangular).""" d = e_model ** 2 - e_naive ** 2 n = len(d) dbar = d.mean() dc = d - dbar lrv = np.sum(dc * dc) / n for k in range(1, h): lrv += 2 * np.sum(dc[k:] * dc[:-k]) / n if lrv <= 0: return None, None stat = dbar / np.sqrt(lrv / n) stat *= np.sqrt((n + 1 - 2 * h + h * (h - 1) / n) / n) p = 2 * stats.t.sf(abs(stat), df=n - 1) return float(stat), float(p) def main(): n_input_rows, data_issues = {}, [] fr = pd.read_csv(os.path.join(DATA, "first_release.csv"), parse_dates=["quarter", "first_release_vintage_date"]) n_input_rows["first_release.csv"] = int(len(fr)) cur = read_series(os.path.join(DATA, "current_vintage.csv")) n_input_rows["current_vintage.csv"] = int(len(cur)) frmap = dict(zip(fr.quarter, fr.first_release_value)) frdate = dict(zip(fr.quarter, fr.first_release_vintage_date)) recs, ledger = [], None for o in ORIGINS: rel = f"vintages/A191RL1Q225SBEA_vintage_{o}.csv" s = read_series(os.path.join(DATA, rel)) n_input_rows[rel] = int(len(s)) assert s.value.notna().all() and not s.date.duplicated().any() last_d, last_v = s.date.iloc[-1], float(s.value.iloc[-1]) # vintage's last value should equal that quarter's first release (sanity, not used) if last_d in frmap and abs(frmap[last_d] - last_v) > 1e-9: data_issues.append({"issue": f"vintage {o}: last value {last_v} != first release {frmap[last_d]}", "handling": "vintage value used as forecast input (real-time rule)"}) rows, counts = ar1_design(s) c, phi = fit_ar1(rows) f = last_v for h in (1, 2): tgt = last_d + pd.DateOffset(months=3 * h) actual = frmap.get(tgt, np.nan) assert frdate[tgt] > pd.Timestamp(o), "actual must be published after origin" f = c + phi * f base = dict(origin_date=o, target_date=tgt.strftime("%Y-%m-%d"), horizon=h, actual=actual, vintage_date=o, train_last_date=last_d.strftime("%Y-%m-%d"), train_last_value=last_v) recs.append({**base, "forecast": last_v, "model": "naive", "n_train": int(len(s))}) recs.append({**base, "forecast": float(f), "model": MODEL, "n_train": int(counts[3]), "_c": c, "_phi": phi}) if o == ORIGINS[-1]: ledger = counts data_issues.append({"issue": "Origins are not equally spaced: the 8th origin is 2025-12-31 (5 months after 2025-07-31) because the 2025Q3 advance estimate was delayed by the Oct-Nov 2025 federal shutdown (data/README.md).", "handling": "Used as given; horizons are defined relative to the last quarter in each vintage, so targets remain consecutive quarters."}) data_issues.append({"issue": "All vintages contain extreme COVID-19 observations (2020Q2 about -28, 2020Q3 about +35 pp SAAR).", "handling": "AR(1) regression rows with dependent quarter 2020Q2-2020Q4 excluded (pre-specified; see out/exclusions.csv). Naive forecast unaffected."}) data_issues.append({"issue": "Target quarter 2025Q4 first release became available only at the 2026-02-28 vintage (shutdown delay).", "handling": "first_release.csv used as given; it is after the origin, so no look-ahead."}) fc = pd.DataFrame(recs) fc["error"] = fc.actual - fc.forecast cols = ["origin_date", "target_date", "horizon", "forecast", "actual", "model", "vintage_date", "train_last_date", "train_last_value", "n_train"] os.makedirs(OUT, exist_ok=True) fc[cols].to_csv(os.path.join(OUT, "forecasts.csv"), index=False, float_format="%.6f") # exclusion ledger (last origin's model estimation sample) n0, n1, n2, n3 = ledger ex = pd.DataFrame([ dict(step=1, rule="keep quarters >= 1985Q1 (last-origin vintage 2025-12-31)", n_before=n0, n_after=n1, n_removed=n0 - n1, rationale="pre-specified post-Great-Moderation estimation window"), dict(step=2, rule="drop first quarter of window (no lag available for AR(1))", n_before=n1, n_after=n2, n_removed=n1 - n2, rationale="AR(1) regression requires y_{t-1}"), dict(step=3, rule="drop regression rows with t in 2020Q2,2020Q3,2020Q4", n_before=n2, n_after=n3, n_removed=n2 - n3, rationale="COVID outliers (-28/+35 pp) as dependent or lagged value would dominate OLS"), ]) ex.to_csv(os.path.join(OUT, "exclusions.csv"), index=False) res = {"task_id": "F3-conf-02", "model_name": MODEL, "model_description": ("AR(1) with intercept estimated by OLS, separately at each origin, on that origin's vintage " "from 1985Q1 to its last quarter, excluding regression rows whose dependent quarter is " "2020Q2-2020Q4; h=2 forecast by iterating the AR(1)."), "n_origins": len(ORIGINS), "first_origin": ORIGINS[0], "last_origin": ORIGINS[-1]} for m, tag in (("naive", "naive"), (MODEL, "model")): for h in (1, 2): e = fc[(fc.model == m) & (fc.horizon == h)].error.values res[f"rmse_{tag}_h{h}"] = float(np.sqrt(np.mean(e ** 2))) res[f"mae_{tag}_h{h}"] = float(np.mean(np.abs(e))) for h in (1, 2): res[f"rel_improvement_rmse_h{h}"] = 1 - res[f"rmse_model_h{h}"] / res[f"rmse_naive_h{h}"] t1 = fc[(fc.model == "naive") & (fc.horizon == 1)] curmap = dict(zip(cur.date, cur.value)) revs = [curmap[pd.Timestamp(t)] - a for t, a in zip(t1.target_date, t1.actual)] res["mean_abs_revision_h1_targets"] = float(np.mean(np.abs(revs))) res["feasibility"] = "complete" res["feasibility_reason"] = None res["data_issues"] = data_issues res["n_input_rows"] = n_input_rows res["n_analysed"] = int(n3) res["random_seed"] = None # supplementary for h in (1, 2): em = fc[(fc.model == MODEL) & (fc.horizon == h)].error.values en = fc[(fc.model == "naive") & (fc.horizon == h)].error.values st, p = dm_hln(em, en, h) res[f"supp_dm_hln_stat_h{h}"] = st res[f"supp_dm_hln_pvalue_h{h}"] = p res[f"supp_n_origins_model_abs_err_lt_naive_h{h}"] = int(np.sum(np.abs(em) < np.abs(en))) res[f"supp_mean_error_naive_h{h}"] = float(en.mean()) res[f"supp_mean_error_model_h{h}"] = float(em.mean()) mod = fc[(fc.model == MODEL) & (fc.horizon == 1)] res["supp_ar1_phi_range"] = [float(mod._phi.min()), float(mod._phi.max())] res["supp_ar1_implied_mean_range"] = [float((mod._c / (1 - mod._phi)).min()), float((mod._c / (1 - mod._phi)).max())] res["supp_revisions_h1_targets"] = {t: float(r) for t, r in zip(t1.target_date, revs)} with open(os.path.join(OUT, "results.json"), "w") as fh: json.dump(res, fh, indent=2) print(json.dumps({k: v for k, v in res.items() if k.startswith(("rmse", "mae", "rel", "mean_abs"))}, indent=1)) if __name__ == "__main__": main()