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. Deterministic (no randomness)."""4import json, os, glob5import numpy as np, pandas as pd67DATA, OUT = "data", "out"8ORIGINS = ["2024-01-31", "2024-04-30", "2024-07-31", "2024-10-31",9 "2025-01-31", "2025-04-30", "2025-07-31", "2025-12-31"]10MODEL = "ar1_post1985_excovid"11SAMPLE_START = pd.Timestamp("1985-01-01")12COVID = [pd.Timestamp(x) for x in ("2020-04-01", "2020-07-01")]1314def vpath(d): return f"vintages/A191RL1Q225SBEA_vintage_{d}.csv"1516def read_series(rel):17 df = pd.read_csv(os.path.join(DATA, rel))18 df.columns = ["date", "value"]19 df["date"] = pd.to_datetime(df["date"])20 return df.sort_values("date").reset_index(drop=True)2122n_input_rows = {}23cur = read_series("current_vintage.csv"); n_input_rows["current_vintage.csv"] = len(cur)24fr = pd.read_csv(os.path.join(DATA, "first_release.csv")); n_input_rows["first_release.csv"] = len(fr)25fr["quarter"] = pd.to_datetime(fr["quarter"]); fr["first_release_vintage_date"] = pd.to_datetime(fr["first_release_vintage_date"])26fr_map = fr.set_index("quarter")27cur_map = cur.set_index("date")["value"]2829def ar1_design(v):30 """Regression obs t with y_t, y_{t-1}; keep t >= 1985Q1, drop if t or t-1 is a pandemic outlier quarter."""31 y = v["value"].to_numpy(); d = v["date"].to_numpy()32 df = pd.DataFrame({"date": d[1:], "lagdate": d[:-1], "y": y[1:], "ylag": y[:-1]})33 steps = [("raw_regression_pairs", len(df))]34 df = df[df["date"] >= SAMPLE_START]; steps.append(("t>=1985Q1", len(df)))35 df = df[~df["date"].isin(COVID) & ~df["lagdate"].isin(COVID)]; steps.append(("drop_pandemic_2020Q2-Q4", len(df)))36 return df, steps3738rows, excl, data_issues = [], [], []39for o in ORIGINS:40 rel = vpath(o); v = read_series(rel); n_input_rows[rel] = len(v)41 assert v["value"].notna().all() and not v["date"].duplicated().any()42 last_d, last_y = v["date"].iloc[-1], float(v["value"].iloc[-1])43 # model estimation (vintage o only)44 df, steps = ar1_design(v)45 X = np.column_stack([np.ones(len(df)), df["ylag"].to_numpy()])46 beta, *_ = np.linalg.lstsq(X, df["y"].to_numpy(), rcond=None)47 c, phi = float(beta[0]), float(beta[1])48 f_model, prev = {}, last_y49 for h in (1, 2):50 prev = c + phi * prev; f_model[h] = prev51 if o == ORIGINS[-1]:52 excl.append(dict(step=0, rule="rows in last-origin vintage file", n_before=len(v), n_after=len(v), n_removed=0,53 rationale="raw input (vintage 2025-12-31)"))54 excl.append(dict(step=1, rule="form AR(1) pairs (first quarter has no lag)", n_before=len(v), n_after=steps[0][1],55 n_removed=len(v) - steps[0][1], rationale="1947Q2 has no lagged value in file"))56 excl.append(dict(step=2, rule="keep t >= 1985Q1", n_before=steps[0][1], n_after=steps[1][1],57 n_removed=steps[0][1] - steps[1][1], rationale="post-Great-Moderation estimation sample (a priori)"))58 excl.append(dict(step=3, rule="drop t in {2020Q2,2020Q3,2020Q4} (pandemic outlier as y_t or y_{t-1})", n_before=steps[1][1],59 n_after=steps[2][1], n_removed=steps[1][1] - steps[2][1], rationale="extreme pandemic outliers (-28%, +35%) would dominate OLS"))60 n_analysed = steps[2][1]61 for h in (1, 2):62 tgt = last_d + pd.DateOffset(months=3 * h)63 if tgt not in fr_map.index:64 raise RuntimeError(f"no first release for {tgt}")65 actual = float(fr_map.loc[tgt, "first_release_value"])66 assert fr_map.loc[tgt, "first_release_vintage_date"] > pd.Timestamp(o), "actual published before origin"67 for m, f, n in (("naive", last_y, len(v)), (MODEL, f_model[h], len(df))):68 rows.append(dict(origin_date=o, target_date=tgt.strftime("%Y-%m-%d"), horizon=h, forecast=round(f, 10),69 actual=actual, model=m, vintage_date=o, train_last_date=last_d.strftime("%Y-%m-%d"),70 train_last_value=last_y, n_train=n, ar_const=c if m != "naive" else np.nan,71 ar_phi=phi if m != "naive" else np.nan))7273fc = pd.DataFrame(rows)74os.makedirs(OUT, exist_ok=True)75fc.drop(columns=["ar_const", "ar_phi"]).to_csv(os.path.join(OUT, "forecasts.csv"), index=False)76pd.DataFrame(excl).to_csv(os.path.join(OUT, "exclusions.csv"), index=False)7778def err(m, h):79 s = fc[(fc.model == m) & (fc.horizon == h)]; e = s["forecast"] - s["actual"]80 return float(np.sqrt(np.mean(e ** 2))), float(np.mean(np.abs(e)))8182res = {"task_id": "F3-conf-02", "model_name": MODEL,83 "model_description": ("AR(1) with intercept estimated by OLS separately on each origin's vintage, regression sample "84 "t = 1985Q1 to the last quarter in the vintage excluding t in 2020Q2-2020Q4 (pandemic outliers as "85 "y_t or y_{t-1}); h=2 forecast by iterating the h=1 forecast."),86 "n_origins": len(ORIGINS), "first_origin": ORIGINS[0], "last_origin": ORIGINS[-1]}87for m, k in (("naive", "naive"), (MODEL, "model")):88 for h in (1, 2):89 r, a = err(m, h); res[f"rmse_{k}_h{h}"] = r; res[f"mae_{k}_h{h}"] = a90for h in (1, 2):91 res[f"rel_improvement_rmse_h{h}"] = 1 - res[f"rmse_model_h{h}"] / res[f"rmse_naive_h{h}"]92t1 = fc[(fc.model == "naive") & (fc.horizon == 1)]93revs = [abs(float(cur_map.loc[pd.Timestamp(t)]) - a) for t, a in zip(t1.target_date, t1.actual)]94# Diebold-Mariano (squared loss) with Harvey-Leybourne-Newbold small-sample correction; supplementary95from scipy import stats96for h in (1, 2):97 en = fc[(fc.model == "naive") & (fc.horizon == h)].sort_values("origin_date")98 em = fc[(fc.model == MODEL) & (fc.horizon == h)].sort_values("origin_date")99 dl = (en["forecast"].to_numpy() - en["actual"].to_numpy()) ** 2 - (em["forecast"].to_numpy() - em["actual"].to_numpy()) ** 2100 T = len(dl); dbar = dl.mean(); g = [np.mean((dl[k:] - dbar) * (dl[:T - k] - dbar)) for k in range(h)]101 var = (g[0] + 2 * sum(g[1:])) / T102 dm = dbar / np.sqrt(var) if var > 0 else np.nan103 hln = dm * np.sqrt((T + 1 - 2 * h + h * (h - 1) / T) / T)104 res[f"supp_dm_hln_stat_h{h}"] = float(hln)105 res[f"supp_dm_hln_pvalue_two_sided_h{h}"] = float(2 * stats.t.sf(abs(hln), df=T - 1))106 res[f"supp_n_model_wins_sqloss_h{h}"] = int(np.sum(dl > 0))107res["mean_abs_revision_h1_targets"] = float(np.mean(revs))108res["feasibility"] = "complete"; res["feasibility_reason"] = None109res["data_issues"] = [110 {"issue": "Pandemic outlier quarters 2020Q2 (about -28%) and 2020Q3 (about +35%) present in every vintage",111 "handling": "AR(1) regression observations with 2020Q2 or 2020Q3 as dependent or lagged value are excluded (decided before evaluation); naive unaffected"},112 {"issue": "Eighth origin is 2025-12-31 instead of 2025-10-31 because the 2025Q3 advance estimate was delayed by the federal shutdown",113 "handling": "Used as given; targets defined relative to the last quarter in that vintage (2025Q4, 2026Q1)"},114 {"issue": "current_vintage.csv contains revisions unknown at the origins",115 "handling": "Used only for the revision statistic mean_abs_revision_h1_targets, never for estimation or forecast input"}]116res["n_input_rows"] = n_input_rows117res["n_analysed"] = int(n_analysed)118res["random_seed"] = None119with open(os.path.join(OUT, "results.json"), "w") as f:120 json.dump(res, f, indent=2)121# per-origin diagnostics (printed only; consistent with written results)122print(fc[["origin_date", "target_date", "horizon", "model", "forecast", "actual", "n_train", "ar_const", "ar_phi"]].to_string(index=False))123print("revisions h1 targets:", [round(x, 2) for x in revs])124print(json.dumps({k: v for k, v in res.items() if k.startswith(("rmse", "mae", "rel", "mean"))}, indent=1))125