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. No randomness."""4import json, os, glob5import numpy as np, pandas as pd67D, O = "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_1985_excovid"11SAMPLE_START = pd.Timestamp("1985-04-01") # first dependent-variable quarter12COVID_EXCL = pd.to_datetime(["2020-04-01", "2020-07-01", "2020-10-01"]) # dependent quarters dropped1314def read_series(path):15 d = pd.read_csv(path)16 d.columns = ["date", "value"]17 d["date"] = pd.to_datetime(d["date"])18 return d1920n_input_rows, issues = {}, []21fr = pd.read_csv(os.path.join(D, "first_release.csv"), parse_dates=["quarter", "first_release_vintage_date"])22n_input_rows["first_release.csv"] = len(fr)23cur = read_series(os.path.join(D, "current_vintage.csv"))24n_input_rows["current_vintage.csv"] = len(cur)25frmap = dict(zip(fr["quarter"], fr["first_release_value"]))26frdate = dict(zip(fr["quarter"], fr["first_release_vintage_date"]))2728def ar1_fit(s):29 """s: vintage DataFrame. OLS of y_t on (1, y_{t-1}) over the restricted sample."""30 y = s.set_index("date")["value"]31 df = pd.DataFrame({"y": y, "ylag": y.shift(1)}).dropna()32 df = df[(df.index >= SAMPLE_START) & (~df.index.isin(COVID_EXCL))]33 X = np.column_stack([np.ones(len(df)), df["ylag"].values])34 beta, *_ = np.linalg.lstsq(X, df["y"].values, rcond=None)35 return beta, len(df), df3637rows, excl_rows = [], []38for o in ORIGINS:39 rel = f"vintages/A191RL1Q225SBEA_vintage_{o}.csv"40 s = read_series(os.path.join(D, rel))41 n_input_rows[rel] = len(s)42 assert s["value"].notna().all() and not s["date"].duplicated().any()43 last_date, last_val = s["date"].iloc[-1], float(s["value"].iloc[-1])44 beta, n_ar, df_ar = ar1_fit(s)45 if o == ORIGINS[-1]:46 n_raw = len(s)47 n1 = int((s["date"] >= SAMPLE_START).sum())48 n2 = n_ar49 excl_rows += [50 dict(step=1, rule="drop dependent quarters before 1985-04-01", n_before=n_raw, n_after=n1, n_removed=n_raw - n1,51 rationale="AR estimated on post-1985 (Great Moderation) sample, fixed a priori; 1947Q2 also lacks a lag"),52 dict(step=2, rule="drop dependent quarters 2020Q2-2020Q4", n_before=n1, n_after=n2, n_removed=n1 - n2,53 rationale="COVID outliers (-28%, +35%) would dominate OLS as y_t or y_{t-1}"),54 ]55 n_analysed = n256 path_model, prev = [], last_val57 for h in (1, 2):58 prev = beta[0] + beta[1] * prev59 path_model.append(prev)60 for h in (1, 2):61 tgt = last_date + pd.DateOffset(months=3 * h)62 assert tgt in frmap, f"no first release for {tgt}"63 assert frdate[tgt] > pd.Timestamp(o), "actual must be published after origin"64 for model, fc, ntr in (("naive", last_val, len(s)), (MODEL, path_model[h - 1], n_ar)):65 rows.append(dict(origin_date=o, target_date=tgt.strftime("%Y-%m-%d"), horizon=h,66 forecast=round(float(fc), 10), actual=float(frmap[tgt]), model=model,67 vintage_date=o, train_last_date=last_date.strftime("%Y-%m-%d"),68 train_last_value=last_val, n_train=int(ntr)))69 if abs(last_val - frmap.get(last_date, np.nan)) > 1e-9:70 issues.append({"issue": f"vintage {o} last value differs from first_release for {last_date.date()}",71 "handling": "vintage value used as input (information rule); first release used only as actual"})7273fc = pd.DataFrame(rows)74fc.to_csv(os.path.join(O, "forecasts.csv"), index=False)75pd.DataFrame(excl_rows).to_csv(os.path.join(O, "exclusions.csv"), index=False)7677def err(model, h):78 e = fc[(fc.model == model) & (fc.horizon == h)]79 e = e["forecast"] - e["actual"]80 assert len(e) == 881 return float(np.sqrt(np.mean(e ** 2))), float(np.mean(np.abs(e)))8283res = {"task_id": "F3-conf-02", "model_name": MODEL,84 "model_description": ("AR(1) with intercept estimated by OLS at each origin on that origin's vintage only, dependent sample "85 "1985Q2 to the last vintage quarter excluding dependent quarters 2020Q2-2020Q4 (COVID outliers); "86 "h=2 forecast by iterating the AR(1)."),87 "n_origins": len(ORIGINS), "first_origin": ORIGINS[0], "last_origin": ORIGINS[-1]}88for h in (1, 2):89 rn, mn = err("naive", h); rm, mm = err(MODEL, h)90 res[f"rmse_naive_h{h}"], res[f"mae_naive_h{h}"] = rn, mn91 res[f"rmse_model_h{h}"], res[f"mae_model_h{h}"] = rm, mm92for h in (1, 2):93 res[f"rel_improvement_rmse_h{h}"] = 1 - res[f"rmse_model_h{h}"] / res[f"rmse_naive_h{h}"]94tg = pd.to_datetime(fc[(fc.horizon == 1) & (fc.model == "naive")]["target_date"])95curmap = dict(zip(cur["date"], cur["value"]))96res["mean_abs_revision_h1_targets"] = float(np.mean([abs(curmap[t] - frmap[t]) for t in tg]))97res["feasibility"] = "complete"; res["feasibility_reason"] = None98issues.append({"issue": "COVID-19 outliers 2020Q2 (about -28%) and 2020Q3 (about +35%) in all vintages",99 "handling": "dependent quarters 2020Q2-2020Q4 excluded from AR(1) estimation (decided before evaluation)"})100issues.append({"issue": "8th origin is 2025-12-31 (2025Q3 advance estimate delayed by the 2025 federal shutdown); 2025Q4 first release dated 2026-02-28",101 "handling": "used as provided; origin dates and first-release dates checked to respect chronology"})102res["data_issues"] = issues103res["n_input_rows"] = n_input_rows104res["n_analysed"] = int(n_analysed)105res["random_seed"] = None106# supplementary uncertainty summaries (n=8 per horizon: descriptive only)107from scipy import stats108for h in (1, 2):109 a = fc[(fc.model == "naive") & (fc.horizon == h)].reset_index(drop=True)110 b = fc[(fc.model == MODEL) & (fc.horizon == h)].reset_index(drop=True)111 d = (a.forecast - a.actual) ** 2 - (b.forecast - b.actual) ** 2 # >0: model better112 n = len(d); dc = d - d.mean()113 lrv = np.sum(dc ** 2) / n + sum(2 * (1 - k / h) * np.sum(dc[k:].values * dc[:-k].values) / n for k in range(1, h))114 t = d.mean() / np.sqrt(lrv / n)115 res[f"supp_n_origins_model_lower_sqerr_h{h}"] = int((d > 0).sum())116 res[f"supp_dm_stat_h{h}"] = float(t)117 res[f"supp_dm_pvalue_t7_h{h}"] = float(2 * stats.t.sf(abs(t), df=n - 1))118res["ar1_coefficients_last_origin"] = {"intercept": float(beta[0]), "phi": float(beta[1])}119with open(os.path.join(O, "results.json"), "w") as f:120 json.dump(res, f, indent=2)121print(json.dumps({k: v for k, v in res.items() if k.startswith(("rmse", "mae", "rel", "mean"))}, indent=2))122