1"""F3-conf-02: real-time forecasts of US real GDP growth from ALFRED vintages.2Run from workspace root: ~/wp29/venv/bin/python out/analysis.py3Reads only data/; writes out/results.json, out/forecasts.csv, out/exclusions.csv.4No randomness is used."""5import json, os6import numpy as np7import pandas as pd8from scipy import stats910DATA, OUT = "data", "out"11ORIGINS = ["2024-01-31", "2024-04-30", "2024-07-31", "2024-10-31",12 "2025-01-31", "2025-04-30", "2025-07-31", "2025-12-31"]13MODEL = "ar1_1985"14EST_START = pd.Timestamp("1985-01-01")15OUTLIERS = np.array([pd.Timestamp("2020-04-01"), pd.Timestamp("2020-07-01")], dtype="datetime64[ns]")161718def read_series(path):19 d = pd.read_csv(path)20 d.columns = ["date", "value"]21 d["date"] = pd.to_datetime(d["date"])22 return d.sort_values("date").reset_index(drop=True)232425def ar1_fit(s):26 """s: DataFrame(date,value) from the origin vintage only. Returns (c, phi, n_pairs_used, ledger)."""27 led = []28 w = s[s.date >= EST_START].reset_index(drop=True)29 led.append(("window", f"keep quarters >= {EST_START.date()}", len(s), len(w)))30 y, x = w.value.values[1:], w.value.values[:-1]31 dy, dx = w.date.values[1:], w.date.values[:-1]32 keep = ~(np.isin(dy, OUTLIERS) | np.isin(dx, OUTLIERS))33 led.append(("lag_pairs", "form (y_t, y_t-1) pairs; first window quarter used only as a lag", len(w), len(y)))34 led.append(("pandemic", "drop pairs with y_t or y_t-1 in 2020Q2-2020Q3", len(y), int(keep.sum())))35 X = np.column_stack([np.ones(keep.sum()), x[keep]])36 beta, *_ = np.linalg.lstsq(X, y[keep], rcond=None)37 return float(beta[0]), float(beta[1]), int(keep.sum()), led383940RATIONALE = {"window": "post-1984 (Great Moderation) sample, chosen a priori",41 "lag_pairs": "AR(1) needs one lag",42 "pandemic": "extreme COVID quarters would dominate OLS"}434445def main():46 n_rows = {}47 fr = pd.read_csv(os.path.join(DATA, "first_release.csv"), parse_dates=["quarter", "first_release_vintage_date"])48 n_rows["first_release.csv"] = len(fr)49 cur = read_series(os.path.join(DATA, "current_vintage.csv")) # revision statistic ONLY50 n_rows["current_vintage.csv"] = len(cur)51 frmap = dict(zip(fr.quarter, fr.first_release_value))52 frvint = dict(zip(fr.quarter, fr.first_release_vintage_date))5354 rows, excl = [], []55 for o in ORIGINS:56 rel = f"vintages/A191RL1Q225SBEA_vintage_{o}.csv"57 s = read_series(os.path.join(DATA, rel)) # the ONLY data used at this origin58 n_rows[rel] = len(s)59 assert s.value.notna().all() and not s.date.duplicated().any()60 last_d, last_v = s.date.iloc[-1], float(s.value.iloc[-1])61 c, phi, n_used, led = ar1_fit(s)62 for step, rule, nb, na in led:63 excl.append(dict(step=f"{o}:{step}", rule=rule, n_before=nb, n_after=na,64 n_removed=nb - na, rationale=RATIONALE[step]))65 f = last_v66 for h in (1, 2):67 tgt = last_d + pd.DateOffset(months=3 * h)68 assert tgt in frmap, f"no first release for {tgt}"69 assert frvint[tgt] > pd.Timestamp(o), "target first release not after origin"70 f = c + phi * f # iterated AR(1) forecast71 for model, fc, n in (("naive", last_v, len(s)), (MODEL, f, n_used)):72 rows.append(dict(origin_date=o, target_date=tgt.strftime("%Y-%m-%d"), horizon=h,73 forecast=round(fc, 6), actual=float(frmap[tgt]), model=model,74 vintage_date=o, train_last_date=last_d.strftime("%Y-%m-%d"),75 train_last_value=last_v, n_train=n))76 print(f"{o}: last {last_d.date()}={last_v}, AR(1) c={c:.3f} phi={phi:.3f} n_pairs={n_used}")7778 fc = pd.DataFrame(rows)79 fc.to_csv(os.path.join(OUT, "forecasts.csv"), index=False)80 pd.DataFrame(excl).to_csv(os.path.join(OUT, "exclusions.csv"), index=False)8182 res = dict(task_id="F3-conf-02", model_name=MODEL,83 model_description=("AR(1) with intercept estimated by OLS at each origin on that origin's vintage only, "84 "sample 1985Q1 to last published quarter, excluding regression pairs involving "85 "2020Q2-2020Q3; h=2 forecast iterated from the h=1 forecast."),86 n_origins=len(ORIGINS), first_origin=ORIGINS[0], last_origin=ORIGINS[-1])87 for h in (1, 2):88 for m, tag in (("naive", "naive"), (MODEL, "model")):89 e = fc[(fc.horizon == h) & (fc.model == m)].eval("forecast - actual").values90 assert len(e) == 891 res[f"rmse_{tag}_h{h}"] = round(float(np.sqrt(np.mean(e ** 2))), 6)92 res[f"mae_{tag}_h{h}"] = round(float(np.mean(np.abs(e))), 6)93 for h in (1, 2):94 res[f"rel_improvement_rmse_h{h}"] = round(1 - res[f"rmse_model_h{h}"] / res[f"rmse_naive_h{h}"], 6)95 for h in (1, 2): # descriptive loss-differential summary (n=8: indicative only)96 en = fc[(fc.horizon == h) & (fc.model == "naive")].eval("forecast - actual").values97 em = fc[(fc.horizon == h) & (fc.model == MODEL)].eval("forecast - actual").values98 d = en ** 2 - em ** 299 res[f"n_origins_model_better_sqerr_h{h}"] = int((d > 0).sum())100 res[f"dm_tstat_sqerr_h{h}"] = round(float(d.mean() / (d.std(ddof=1) / np.sqrt(len(d)))), 6)101 res[f"dm_pvalue_two_sided_h{h}"] = round(float(2 * stats.t.sf(abs(res[f"dm_tstat_sqerr_h{h}"]), len(d) - 1)), 6)102 t1 = pd.to_datetime(fc[(fc.horizon == 1) & (fc.model == "naive")].target_date)103 cmap = dict(zip(cur.date, cur.value))104 res["mean_abs_revision_h1_targets"] = round(float(np.mean([abs(cmap[t] - frmap[t]) for t in t1])), 6)105 res["feasibility"] = "complete"106 res["feasibility_reason"] = None107 res["data_issues"] = [108 {"issue": "Pandemic outliers 2020Q2 (about -28) and 2020Q3 (about +35) present in every vintage.",109 "handling": "Regression pairs involving these quarters excluded from AR(1) estimation; naive unaffected."},110 {"issue": "Origins unevenly spaced: 2025Q3 advance estimate delayed by the 2025 federal shutdown, so origin 8 is 2025-12-31.",111 "handling": "Used as specified; target 2025Q4 is the h=2 target of origin 7 and the h=1 target of origin 8."},112 {"issue": "current_vintage.csv contains revisions unknown at the origins.",113 "handling": "Used only for mean_abs_revision_h1_targets, never for estimation or forecast input."},114 ]115 res["n_input_rows"] = n_rows116 res["n_analysed"] = n_rows[f"vintages/A191RL1Q225SBEA_vintage_{ORIGINS[-1]}.csv"]117 res["n_model_estimation_pairs_last_origin"] = int(fc[fc.model == MODEL].n_train.iloc[-1])118 res["random_seed"] = None119 with open(os.path.join(OUT, "results.json"), "w") as fh:120 json.dump(res, fh, indent=2)121 print(json.dumps({k: v for k, v in res.items() if k.startswith(("rmse", "mae", "rel", "mean"))}, indent=1))122123124if __name__ == "__main__":125 main()126