1"""F3-conf-02: real-time forecasts of US real GDP growth (A191RL1Q225SBEA).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 (random_seed = null)."""5import json, os, glob6import numpy as np7import pandas as pd89DATA = "data"; OUT = "out"10MODEL = "ar1"11EST_START = pd.Timestamp("1985-01-01")12COVID_DEP = pd.to_datetime(["2020-04-01", "2020-07-01", "2020-10-01"]) # dependent quarters dropped13DESC = ("AR(1) with intercept, y_t = c + phi*y_{t-1} + e_t, estimated by OLS on the origin vintage "14 "from 1985Q1 (first lag) to the vintage's last quarter, dropping observations whose dependent "15 "quarter is 2020Q2-2020Q4 (COVID outliers as dependent or lag); h=2 by iterating the 1-step forecast.")1617def read_series(path):18 d = pd.read_csv(path)19 d.columns = ["date", "value"]20 d["date"] = pd.to_datetime(d["date"])21 return d.sort_values("date").reset_index(drop=True)2223def add_q(ts, h):24 return ts + pd.DateOffset(months=3 * h)2526def ar1_fit(s):27 """s: vintage DataFrame. Returns (c, phi, n_obs, exclusion steps)."""28 steps = []29 n0 = len(s)30 w = s[s.date >= EST_START].reset_index(drop=True)31 steps.append(("window", "keep quarters >= 1985-01-01", n0, len(w), "Great-Moderation-era estimation window"))32 y = w.value.values[1:]; x = w.value.values[:-1]; dts = w.date.values[1:]33 steps.append(("lag", "first quarter of window used only as lag", len(w), len(y), "AR(1) loses one observation"))34 keep = ~np.isin(dts, COVID_DEP.values)35 steps.append(("covid", "drop dependent quarters 2020Q2-2020Q4", len(y), int(keep.sum()),36 "COVID outliers (-28/+35 % SAAR) would dominate OLS"))37 y, x = y[keep], x[keep]38 X = np.column_stack([np.ones_like(x), x])39 beta, *_ = np.linalg.lstsq(X, y, rcond=None)40 return float(beta[0]), float(beta[1]), int(len(y)), steps4142def main():43 vfiles = sorted(glob.glob(os.path.join(DATA, "vintages", "A191RL1Q225SBEA_vintage_*.csv")))44 fr = pd.read_csv(os.path.join(DATA, "first_release.csv"), parse_dates=["quarter", "first_release_vintage_date"])45 cur = read_series(os.path.join(DATA, "current_vintage.csv")) # used ONLY for revision statistic46 n_input_rows = {"current_vintage.csv": len(cur), "first_release.csv": len(fr)}47 issues = []48 rows = []; excl = None; h1_targets = []; coefs = {}49 for f in vfiles:50 vdate = pd.Timestamp(os.path.basename(f).split("_vintage_")[1][:-4])51 s = read_series(f)52 n_input_rows[os.path.relpath(f, DATA)] = len(s)53 assert s.value.notna().all() and not s.date.duplicated().any()54 last_d, last_v = s.date.iloc[-1], float(s.value.iloc[-1])55 # sanity: last vintage value equals first release of that quarter56 frl = fr.loc[fr.quarter == last_d, "first_release_value"]57 if len(frl) and abs(float(frl.iloc[0]) - last_v) > 1e-9:58 issues.append({"issue": f"vintage {vdate.date()} last value differs from first release",59 "handling": "vintage value used for forecasting, first release as actual"})60 c, phi, n_ar, steps = ar1_fit(s)61 excl = steps # keep last origin's ledger62 coefs[vdate.strftime('%Y-%m-%d')] = {'c': round(c, 6), 'phi': round(phi, 6)}63 f1 = c + phi * last_v; f2 = c + phi * f164 for h, fm in [(1, f1), (2, f2)]:65 tgt = add_q(last_d, h)66 a = fr.loc[fr.quarter == tgt]67 if len(a) == 0:68 actual = np.nan69 issues.append({"issue": f"no first release for {tgt.date()}", "handling": "excluded from errors"})70 else:71 actual = float(a.first_release_value.iloc[0])72 assert a.first_release_vintage_date.iloc[0] > vdate # actual published after origin73 if h == 1: h1_targets.append(tgt)74 base = dict(origin_date=vdate.strftime("%Y-%m-%d"), target_date=tgt.strftime("%Y-%m-%d"), horizon=h,75 actual=actual, vintage_date=vdate.strftime("%Y-%m-%d"),76 train_last_date=last_d.strftime("%Y-%m-%d"), train_last_value=last_v)77 rows.append({**base, "forecast": last_v, "model": "naive", "n_train": len(s)})78 rows.append({**base, "forecast": round(fm, 6), "model": MODEL, "n_train": n_ar})79 fc = pd.DataFrame(rows)[["origin_date", "target_date", "horizon", "forecast", "actual", "model",80 "vintage_date", "train_last_date", "train_last_value", "n_train"]]81 fc = fc.sort_values(["origin_date", "horizon", "model"]).reset_index(drop=True)82 fc.to_csv(os.path.join(OUT, "forecasts.csv"), index=False)8384 def err(model, h):85 d = fc[(fc.model == model) & (fc.horizon == h)].dropna(subset=["actual"])86 e = d.forecast - d.actual87 return float(np.sqrt(np.mean(e ** 2))), float(np.mean(np.abs(e)))88 res = {"task_id": "F3-conf-02", "model_name": MODEL, "model_description": DESC,89 "n_origins": len(vfiles),90 "first_origin": os.path.basename(vfiles[0]).split("_vintage_")[1][:-4],91 "last_origin": os.path.basename(vfiles[-1]).split("_vintage_")[1][:-4]}92 for m, tag in [("naive", "naive"), (MODEL, "model")]:93 for h in (1, 2):94 r, a = err(m, h)95 res[f"rmse_{tag}_h{h}"] = round(r, 6); res[f"mae_{tag}_h{h}"] = round(a, 6)96 for h in (1, 2):97 res[f"rel_improvement_rmse_h{h}"] = round(1 - res[f"rmse_model_h{h}"] / res[f"rmse_naive_h{h}"], 6)98 revs = []99 for t in h1_targets:100 cv = float(cur.loc[cur.date == t, "value"].iloc[0])101 fv = float(fr.loc[fr.quarter == t, "first_release_value"].iloc[0])102 revs.append(abs(cv - fv))103 res["mean_abs_revision_h1_targets"] = round(float(np.mean(revs)), 6)104 res["feasibility"] = "complete"; res["feasibility_reason"] = None105 issues.append({"issue": "COVID quarters 2020Q2-2020Q3 are extreme outliers in every vintage",106 "handling": "AR(1) estimation drops observations with dependent quarter 2020Q2-2020Q4; naive unaffected"})107 issues.append({"issue": "8th origin is 2025-12-31 (2025Q3 advance estimate delayed by the 2025 shutdown); origins are not equally spaced",108 "handling": "used as given in TASK.md; horizons defined relative to last quarter in each vintage"})109 res["data_issues"] = issues110 res["n_input_rows"] = n_input_rows111 res["n_analysed"] = int(fc[(fc.origin_date == res["last_origin"]) & (fc.model == MODEL)].n_train.iloc[0])112 res["random_seed"] = None113 res["ar1_coefficients_by_origin"] = coefs114 for h in (1, 2):115 a = fc[(fc.horizon == h)].pivot(index="origin_date", columns="model", values="forecast")116 act = fc[(fc.horizon == h) & (fc.model == "naive")].set_index("origin_date").actual117 res[f"n_origins_model_abs_error_lower_h{h}"] = int(((a[MODEL] - act).abs() < (a["naive"] - act).abs()).sum())118 with open(os.path.join(OUT, "results.json"), "w") as fh:119 json.dump(res, fh, indent=2)120 ex = pd.DataFrame(excl, columns=["step", "rule", "n_before", "n_after", "rationale"])121 ex["n_removed"] = ex.n_before - ex.n_after122 ex[["step", "rule", "n_before", "n_after", "n_removed", "rationale"]].to_csv(123 os.path.join(OUT, "exclusions.csv"), index=False)124 print(json.dumps({k: v for k, v in res.items() if k.startswith(("rmse", "mae", "rel", "mean"))}, indent=1))125126if __name__ == "__main__":127 main()128