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.4Deterministic (no randomness)."""5import json, glob, os6import numpy as np7import pandas as pd89DATA, OUT = "data", "out"10MODEL = "ar1_post1985_excovid"11SAMPLE_START = pd.Timestamp("1985-01-01") # first dependent-variable quarter12COVID_EXCL = pd.to_datetime(["2020-04-01", "2020-07-01", "2020-10-01"]) # rows whose y_t or y_{t-1} is 2020Q2/Q31314def read_series(path):15 d = pd.read_csv(path, parse_dates=["observation_date"])16 return pd.Series(d.iloc[:, 1].astype(float).values, index=d["observation_date"])1718def qshift(ts, k):19 return (pd.Period(ts, freq="Q") + k).start_time2021def ar1_fit(s, ledger=None):22 """OLS AR(1) with intercept on the vintage series s."""23 n0 = len(s)24 df = pd.DataFrame({"y": s, "ylag": s.shift(1)})25 df1 = df[(df.index >= SAMPLE_START) & df["ylag"].notna()]26 df2 = df1[~df1.index.isin(COVID_EXCL)]27 if ledger is not None:28 ledger.append(("1", f"keep regression rows with dependent quarter >= {SAMPLE_START.date()} (lag needs prior quarter)",29 n0, len(df1), n0 - len(df1), "post-1984 (Great Moderation) sample; pre-1985 volatility regime not representative of 2024-26"))30 ledger.append(("2", "drop rows t in {2020Q2,2020Q3,2020Q4} (y_t or y_{t-1} is a COVID quarter)",31 len(df1), len(df2), len(df1) - len(df2), "pandemic outliers (~-28, ~+35 pp) would dominate OLS"))32 X = np.column_stack([np.ones(len(df2)), df2["ylag"].values])33 beta, *_ = np.linalg.lstsq(X, df2["y"].values, rcond=None)34 return beta, len(df2)3536def main():37 vfiles = sorted(glob.glob(os.path.join(DATA, "vintages", "A191RL1Q225SBEA_vintage_*.csv")))38 cur = read_series(os.path.join(DATA, "current_vintage.csv"))39 fr = pd.read_csv(os.path.join(DATA, "first_release.csv"),40 parse_dates=["quarter", "first_release_vintage_date"]).set_index("quarter")41 n_input_rows = {"current_vintage.csv": int(len(cur)), "first_release.csv": int(len(fr))}42 rows, h1_targets, coefs, ledger = [], [], {}, []43 for i, f in enumerate(vfiles):44 vdate = os.path.basename(f).split("_")[-1].replace(".csv", "")45 s = read_series(f)46 n_input_rows[os.path.relpath(f, DATA)] = int(len(s))47 assert s.index.is_monotonic_increasing and not s.index.duplicated().any() and s.notna().all()48 last_d, last_v = s.index[-1], float(s.iloc[-1])49 is_last = i == len(vfiles) - 150 beta, n_reg = ar1_fit(s, ledger if is_last else None)51 coefs[vdate] = [float(beta[0]), float(beta[1]), n_reg]52 f1 = beta[0] + beta[1] * last_v53 f2 = beta[0] + beta[1] * f154 for h, fm in ((1, f1), (2, f2)):55 tgt = qshift(last_d, h)56 actual = float(fr.loc[tgt, "first_release_value"])57 assert fr.loc[tgt, "first_release_vintage_date"] > pd.Timestamp(vdate) # chronology58 base = dict(origin_date=vdate, target_date=tgt.strftime("%Y-%m-%d"), horizon=h, actual=actual,59 vintage_date=vdate, train_last_date=last_d.strftime("%Y-%m-%d"), train_last_value=last_v)60 rows.append({**base, "forecast": last_v, "model": "naive", "n_train": int(len(s))})61 rows.append({**base, "forecast": float(fm), "model": MODEL, "n_train": int(n_reg)})62 if h == 1:63 h1_targets.append(tgt)64 if is_last:65 n_analysed = n_reg66 fc = pd.DataFrame(rows)[["origin_date", "target_date", "horizon", "forecast", "actual", "model",67 "vintage_date", "train_last_date", "train_last_value", "n_train"]]68 fc = fc.sort_values(["origin_date", "horizon", "model"]).reset_index(drop=True)69 fc.to_csv(os.path.join(OUT, "forecasts.csv"), index=False, float_format="%.10g")70 pd.DataFrame(ledger, columns=["step", "rule", "n_before", "n_after", "n_removed", "rationale"]).to_csv(71 os.path.join(OUT, "exclusions.csv"), index=False)7273 def err(model, h):74 d = fc[(fc.model == model) & (fc.horizon == h)]75 e = d["forecast"] - d["actual"]76 return float(np.sqrt(np.mean(e ** 2))), float(np.mean(np.abs(e)))7778 res = {"task_id": "F3-conf-02", "model_name": MODEL,79 "model_description": ("AR(1) with intercept estimated by OLS at each origin on that origin's vintage only, "80 "dependent quarters 1985Q1 to the last vintage quarter, excluding regression rows t=2020Q2-2020Q4 "81 "(COVID outliers in y_t or y_{t-1}); h=2 forecast iterated from the h=1 forecast."),82 "n_origins": len(vfiles),83 "first_origin": fc.origin_date.min(), "last_origin": fc.origin_date.max()}84 for h in (1, 2):85 res[f"rmse_naive_h{h}"], res[f"mae_naive_h{h}"] = err("naive", h)86 res[f"rmse_model_h{h}"], res[f"mae_model_h{h}"] = err(MODEL, h)87 for h in (1, 2):88 res[f"rel_improvement_rmse_h{h}"] = 1 - res[f"rmse_model_h{h}"] / res[f"rmse_naive_h{h}"]89 for h in (1, 2):90 a = fc[(fc.model == "naive") & (fc.horizon == h)].set_index("origin_date")91 b = fc[(fc.model == MODEL) & (fc.horizon == h)].set_index("origin_date")92 res[f"n_origins_model_abs_error_smaller_h{h}"] = int(((b.forecast - b.actual).abs() < (a.forecast - a.actual).abs()).sum())93 rev = [abs(float(cur.loc[t]) - float(fr.loc[t, "first_release_value"])) for t in h1_targets]94 res["mean_abs_revision_h1_targets"] = float(np.mean(rev))95 res.update({96 "feasibility": "complete", "feasibility_reason": None,97 "data_issues": [98 {"issue": "COVID-19 outliers 2020Q2 (~-28) and 2020Q3 (~+35) in every vintage",99 "handling": "AR regression rows t=2020Q2-2020Q4 excluded (dependent or lag is a COVID quarter); naive unaffected"},100 {"issue": "Irregular origin spacing: 8th origin 2025-12-31 instead of 2025-10-31 (2025Q3 advance estimate delayed by federal shutdown)",101 "handling": "Used as provided; horizons defined relative to last quarter in each vintage, so targets are 2025Q4 and 2026Q1"},102 {"issue": "Values of earlier quarters differ across vintages (revisions)",103 "handling": "Each origin uses only its own vintage; current vintage used only for revision statistic"}],104 "n_input_rows": n_input_rows, "n_analysed": int(n_analysed), "random_seed": None,105 "ar1_coefficients_by_origin": {k: {"intercept": v[0], "phi": v[1], "n_reg": v[2]} for k, v in coefs.items()},106 "abs_revision_h1_targets": {t.strftime("%Y-%m-%d"): r for t, r in zip(h1_targets, rev)},107 })108 with open(os.path.join(OUT, "results.json"), "w") as fh:109 json.dump(res, fh, indent=2)110 print(json.dumps({k: res[k] for k in res if k.startswith(("rmse", "mae", "rel", "mean"))}, indent=2))111112if __name__ == "__main__":113 main()114