1"""Real-time forecasts of US real GDP growth (F3-conf-02).2Run from workspace root: ~/wp29/venv/bin/python out/analysis.py3Reads only data/; writes out/results.json, out/forecasts.csv, out/exclusions.csv.4Deterministic: no randomness used."""5import json, os, glob, hashlib6import numpy as np, pandas as pd7from scipy import stats89DATA, OUT = "data", "out"10ORIGINS = ["2024-01-31", "2024-04-30", "2024-07-31", "2024-10-31",11 "2025-01-31", "2025-04-30", "2025-07-31", "2025-12-31"]12MODEL = "ar1_post1985_exCovid"13SAMPLE_START = pd.Timestamp("1985-01-01")14COVID_ROWS = [pd.Timestamp(x) for x in ("2020-04-01", "2020-07-01", "2020-10-01")]151617def read_series(path):18 d = pd.read_csv(path)19 d.columns = ["date", "value"]20 d["date"] = pd.to_datetime(d["date"])21 return d222324def ar1_design(s):25 """s: DataFrame(date,value) from one vintage. Returns y, x, and ledger counts."""26 n0 = len(s)27 w = s[s.date >= SAMPLE_START].reset_index(drop=True)28 n1 = len(w)29 rows = pd.DataFrame({"date": w.date.iloc[1:].values,30 "y": w.value.iloc[1:].values,31 "ylag": w.value.iloc[:-1].values})32 # consecutive-quarter check33 assert (pd.to_datetime(rows.date).values - w.date.iloc[:-1].values).astype("timedelta64[D]").astype(int).max() <= 9234 n2 = len(rows)35 rows = rows[~rows.date.isin(COVID_ROWS)].reset_index(drop=True)36 n3 = len(rows)37 return rows, (n0, n1, n2, n3)383940def fit_ar1(rows):41 X = np.column_stack([np.ones(len(rows)), rows.ylag.values])42 beta, *_ = np.linalg.lstsq(X, rows.y.values, rcond=None)43 return beta # c, phi444546def dm_hln(e_model, e_naive, h):47 """Diebold-Mariano test (squared loss) with Harvey-Leybourne-Newbold small-sample correction;48 long-run variance uses autocovariances up to lag h-1 (rectangular)."""49 d = e_model ** 2 - e_naive ** 250 n = len(d)51 dbar = d.mean()52 dc = d - dbar53 lrv = np.sum(dc * dc) / n54 for k in range(1, h):55 lrv += 2 * np.sum(dc[k:] * dc[:-k]) / n56 if lrv <= 0:57 return None, None58 stat = dbar / np.sqrt(lrv / n)59 stat *= np.sqrt((n + 1 - 2 * h + h * (h - 1) / n) / n)60 p = 2 * stats.t.sf(abs(stat), df=n - 1)61 return float(stat), float(p)626364def main():65 n_input_rows, data_issues = {}, []66 fr = pd.read_csv(os.path.join(DATA, "first_release.csv"), parse_dates=["quarter", "first_release_vintage_date"])67 n_input_rows["first_release.csv"] = int(len(fr))68 cur = read_series(os.path.join(DATA, "current_vintage.csv"))69 n_input_rows["current_vintage.csv"] = int(len(cur))70 frmap = dict(zip(fr.quarter, fr.first_release_value))71 frdate = dict(zip(fr.quarter, fr.first_release_vintage_date))7273 recs, ledger = [], None74 for o in ORIGINS:75 rel = f"vintages/A191RL1Q225SBEA_vintage_{o}.csv"76 s = read_series(os.path.join(DATA, rel))77 n_input_rows[rel] = int(len(s))78 assert s.value.notna().all() and not s.date.duplicated().any()79 last_d, last_v = s.date.iloc[-1], float(s.value.iloc[-1])80 # vintage's last value should equal that quarter's first release (sanity, not used)81 if last_d in frmap and abs(frmap[last_d] - last_v) > 1e-9:82 data_issues.append({"issue": f"vintage {o}: last value {last_v} != first release {frmap[last_d]}",83 "handling": "vintage value used as forecast input (real-time rule)"})84 rows, counts = ar1_design(s)85 c, phi = fit_ar1(rows)86 f = last_v87 for h in (1, 2):88 tgt = last_d + pd.DateOffset(months=3 * h)89 actual = frmap.get(tgt, np.nan)90 assert frdate[tgt] > pd.Timestamp(o), "actual must be published after origin"91 f = c + phi * f92 base = dict(origin_date=o, target_date=tgt.strftime("%Y-%m-%d"), horizon=h, actual=actual,93 vintage_date=o, train_last_date=last_d.strftime("%Y-%m-%d"), train_last_value=last_v)94 recs.append({**base, "forecast": last_v, "model": "naive", "n_train": int(len(s))})95 recs.append({**base, "forecast": float(f), "model": MODEL, "n_train": int(counts[3]),96 "_c": c, "_phi": phi})97 if o == ORIGINS[-1]:98 ledger = counts99 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).",100 "handling": "Used as given; horizons are defined relative to the last quarter in each vintage, so targets remain consecutive quarters."})101 data_issues.append({"issue": "All vintages contain extreme COVID-19 observations (2020Q2 about -28, 2020Q3 about +35 pp SAAR).",102 "handling": "AR(1) regression rows with dependent quarter 2020Q2-2020Q4 excluded (pre-specified; see out/exclusions.csv). Naive forecast unaffected."})103 data_issues.append({"issue": "Target quarter 2025Q4 first release became available only at the 2026-02-28 vintage (shutdown delay).",104 "handling": "first_release.csv used as given; it is after the origin, so no look-ahead."})105 fc = pd.DataFrame(recs)106 fc["error"] = fc.actual - fc.forecast107 cols = ["origin_date", "target_date", "horizon", "forecast", "actual", "model", "vintage_date",108 "train_last_date", "train_last_value", "n_train"]109 os.makedirs(OUT, exist_ok=True)110 fc[cols].to_csv(os.path.join(OUT, "forecasts.csv"), index=False, float_format="%.6f")111112 # exclusion ledger (last origin's model estimation sample)113 n0, n1, n2, n3 = ledger114 ex = pd.DataFrame([115 dict(step=1, rule="keep quarters >= 1985Q1 (last-origin vintage 2025-12-31)", n_before=n0, n_after=n1,116 n_removed=n0 - n1, rationale="pre-specified post-Great-Moderation estimation window"),117 dict(step=2, rule="drop first quarter of window (no lag available for AR(1))", n_before=n1, n_after=n2,118 n_removed=n1 - n2, rationale="AR(1) regression requires y_{t-1}"),119 dict(step=3, rule="drop regression rows with t in 2020Q2,2020Q3,2020Q4", n_before=n2, n_after=n3,120 n_removed=n2 - n3, rationale="COVID outliers (-28/+35 pp) as dependent or lagged value would dominate OLS"),121 ])122 ex.to_csv(os.path.join(OUT, "exclusions.csv"), index=False)123124 res = {"task_id": "F3-conf-02", "model_name": MODEL,125 "model_description": ("AR(1) with intercept estimated by OLS, separately at each origin, on that origin's vintage "126 "from 1985Q1 to its last quarter, excluding regression rows whose dependent quarter is "127 "2020Q2-2020Q4; h=2 forecast by iterating the AR(1)."),128 "n_origins": len(ORIGINS), "first_origin": ORIGINS[0], "last_origin": ORIGINS[-1]}129 for m, tag in (("naive", "naive"), (MODEL, "model")):130 for h in (1, 2):131 e = fc[(fc.model == m) & (fc.horizon == h)].error.values132 res[f"rmse_{tag}_h{h}"] = float(np.sqrt(np.mean(e ** 2)))133 res[f"mae_{tag}_h{h}"] = float(np.mean(np.abs(e)))134 for h in (1, 2):135 res[f"rel_improvement_rmse_h{h}"] = 1 - res[f"rmse_model_h{h}"] / res[f"rmse_naive_h{h}"]136 t1 = fc[(fc.model == "naive") & (fc.horizon == 1)]137 curmap = dict(zip(cur.date, cur.value))138 revs = [curmap[pd.Timestamp(t)] - a for t, a in zip(t1.target_date, t1.actual)]139 res["mean_abs_revision_h1_targets"] = float(np.mean(np.abs(revs)))140 res["feasibility"] = "complete"141 res["feasibility_reason"] = None142 res["data_issues"] = data_issues143 res["n_input_rows"] = n_input_rows144 res["n_analysed"] = int(n3)145 res["random_seed"] = None146 # supplementary147 for h in (1, 2):148 em = fc[(fc.model == MODEL) & (fc.horizon == h)].error.values149 en = fc[(fc.model == "naive") & (fc.horizon == h)].error.values150 st, p = dm_hln(em, en, h)151 res[f"supp_dm_hln_stat_h{h}"] = st152 res[f"supp_dm_hln_pvalue_h{h}"] = p153 res[f"supp_n_origins_model_abs_err_lt_naive_h{h}"] = int(np.sum(np.abs(em) < np.abs(en)))154 res[f"supp_mean_error_naive_h{h}"] = float(en.mean())155 res[f"supp_mean_error_model_h{h}"] = float(em.mean())156 mod = fc[(fc.model == MODEL) & (fc.horizon == 1)]157 res["supp_ar1_phi_range"] = [float(mod._phi.min()), float(mod._phi.max())]158 res["supp_ar1_implied_mean_range"] = [float((mod._c / (1 - mod._phi)).min()), float((mod._c / (1 - mod._phi)).max())]159 res["supp_revisions_h1_targets"] = {t: float(r) for t, r in zip(t1.target_date, revs)}160 with open(os.path.join(OUT, "results.json"), "w") as fh:161 json.dump(res, fh, indent=2)162 print(json.dumps({k: v for k, v in res.items() if k.startswith(("rmse", "mae", "rel", "mean_abs"))}, indent=1))163164165if __name__ == "__main__":166 main()167