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."""4import json, os, glob5import numpy as np6import pandas as pd78DATA, OUT = "data", "out"9ORIGINS = ["2024-01-31", "2024-04-30", "2024-07-31", "2024-10-31",10 "2025-01-31", "2025-04-30", "2025-07-31", "2025-12-31"]11SAMPLE_START = pd.Timestamp("1985-01-01")12EXCLUDED_Q = [pd.Timestamp("2020-04-01"), pd.Timestamp("2020-07-01")]13MODEL = "ar1_post1985"14DESC = ("AR(1) with intercept estimated by OLS on the origin vintage, sample 1985Q1 to the last quarter "15 "of that vintage, with 2020Q2 and 2020Q3 set to missing (no regression pair uses them); "16 "h=2 forecast obtained by iterating the AR(1) recursion.")1718def 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.set_index("date")["value"].astype(float)2324def qshift(ts, k):25 return ts + pd.DateOffset(months=3 * k)2627n_input_rows = {}28vint = {}29for o in ORIGINS:30 rel = f"vintages/A191RL1Q225SBEA_vintage_{o}.csv"31 s = read_series(os.path.join(DATA, rel))32 vint[o] = s33 n_input_rows[rel] = int(len(s))34cur = read_series(os.path.join(DATA, "current_vintage.csv")) # used ONLY for revision statistic35n_input_rows["current_vintage.csv"] = int(len(cur))36fr = pd.read_csv(os.path.join(DATA, "first_release.csv"), parse_dates=["quarter", "first_release_vintage_date"])37n_input_rows["first_release.csv"] = int(len(fr))38fr_val = fr.set_index("quarter")["first_release_value"].astype(float)39fr_vd = fr.set_index("quarter")["first_release_vintage_date"]4041data_issues = []42for o, s in vint.items():43 assert s.index.is_monotonic_increasing and not s.index.duplicated().any() and s.notna().all()44 assert (s.index.to_series().diff().dropna().dt.days <= 92).all()45data_issues.append({"issue": "Pandemic outliers 2020Q2/2020Q3 (about -28% and +35% SAAR) in every vintage.",46 "handling": "Set to missing in the AR(1) estimation sample (documented in out/exclusions.csv); naive forecasts unaffected."})47data_issues.append({"issue": "Historical values differ across vintages (annual revisions, e.g. 2024Q4 2.3 -> 2.4 -> 1.9 across vintages).",48 "handling": "Each origin uses only its own vintage file; current vintage used only for the revision statistic."})49data_issues.append({"issue": "Eighth origin is 2025-12-31 rather than 2025-10-31 (2025Q3 advance estimate delayed by the federal shutdown); its h=1 target 2025Q4 first release dated 2026-02-28.",50 "handling": "Used as given in TASK.md; origins are not equally spaced in calendar time but each still forecasts 1 and 2 quarters past its last published quarter."})5152def fit_ar1(s):53 y = s[s.index >= SAMPLE_START].copy()54 n_after_start = len(y)55 y.loc[y.index.isin(EXCLUDED_Q)] = np.nan56 lag = y.shift(1)57 # require the lag to be the immediately preceding quarter (index is contiguous quarterly)58 m = y.notna() & lag.notna()59 X = np.column_stack([np.ones(m.sum()), lag[m].values])60 beta, *_ = np.linalg.lstsq(X, y[m].values, rcond=None)61 return beta, int(y.notna().sum()), n_after_start, int(m.sum())6263rows, excl = [], []64fits = {}65for o in ORIGINS:66 s = vint[o]67 T, yT = s.index[-1], float(s.iloc[-1])68 beta, n_used, n_after_start, n_pairs = fit_ar1(s)69 fits[o] = {"c": float(beta[0]), "phi": float(beta[1]), "n_pairs": n_pairs}70 excl.append({"step": f"{o}:sample_start", "rule": "drop quarters before 1985-01-01 (model only)",71 "n_before": len(s), "n_after": n_after_start, "n_removed": len(s) - n_after_start,72 "rationale": "post-Great-Moderation estimation sample for AR(1)"})73 excl.append({"step": f"{o}:pandemic", "rule": "set 2020-04-01 and 2020-07-01 to missing (model only)",74 "n_before": n_after_start, "n_after": n_used, "n_removed": n_after_start - n_used,75 "rationale": "extreme pandemic swings would dominate OLS AR coefficient"})76 f_ar = {}77 prev = yT78 for h in (1, 2):79 prev = beta[0] + beta[1] * prev80 f_ar[h] = float(prev)81 for h in (1, 2):82 tgt = qshift(T, h)83 assert tgt in fr_val.index, f"missing first release for {tgt}"84 assert fr_vd[tgt] >= pd.Timestamp(o), "actual published before origin"85 for model, fc, ntr in (("naive", yT, len(s)), (MODEL, f_ar[h], n_used)):86 rows.append({"origin_date": o, "target_date": tgt.strftime("%Y-%m-%d"), "horizon": h,87 "forecast": round(float(fc), 10), "actual": float(fr_val[tgt]), "model": model,88 "vintage_date": o, "train_last_date": T.strftime("%Y-%m-%d"),89 "train_last_value": yT, "n_train": int(ntr)})9091fc = pd.DataFrame(rows)92fc["error"] = fc["actual"] - fc["forecast"]93fc.drop(columns="error").to_csv(os.path.join(OUT, "forecasts.csv"), index=False)94# Main ledger: single chain for the last origin (whose model training set defines n_analysed);95# the same two rules applied at every origin are listed in exclusions_by_origin.csv.96pd.DataFrame([e for e in excl if e["step"].startswith(ORIGINS[-1])]).to_csv(os.path.join(OUT, "exclusions.csv"), index=False)97pd.DataFrame(excl).to_csv(os.path.join(OUT, "exclusions_by_origin.csv"), index=False)9899def metrics(model, h):100 e = fc[(fc.model == model) & (fc.horizon == h)]["error"].values101 assert len(e) == len(ORIGINS)102 return float(np.sqrt(np.mean(e ** 2))), float(np.mean(np.abs(e)))103104res = {"task_id": "F3-conf-02", "model_name": MODEL, "model_description": DESC,105 "n_origins": len(ORIGINS), "first_origin": ORIGINS[0], "last_origin": ORIGINS[-1]}106for model, tag in (("naive", "naive"), (MODEL, "model")):107 for h in (1, 2):108 r, a = metrics(model, h)109 res[f"rmse_{tag}_h{h}"] = round(r, 10)110 res[f"mae_{tag}_h{h}"] = round(a, 10)111for h in (1, 2):112 res[f"rel_improvement_rmse_h{h}"] = round(1 - res[f"rmse_model_h{h}"] / res[f"rmse_naive_h{h}"], 10)113h1t = pd.to_datetime(fc[(fc.model == "naive") & (fc.horizon == 1)]["target_date"])114rev = [abs(float(cur[t]) - float(fr_val[t])) for t in h1t]115res["mean_abs_revision_h1_targets"] = round(float(np.mean(rev)), 10)116res["feasibility"] = "complete"117res["feasibility_reason"] = None118res["data_issues"] = data_issues119res["n_input_rows"] = n_input_rows120res["n_analysed"] = int(fc[(fc.origin_date == ORIGINS[-1]) & (fc.model == MODEL)]["n_train"].iloc[0])121res["random_seed"] = None122# supplementary diagnostics (not required keys)123res["supp_ar1_coefficients_by_origin"] = fits124res["supp_abs_revision_h1_targets"] = {t.strftime("%Y-%m-%d"): round(r, 10) for t, r in zip(h1t, rev)}125for model, tag in (("naive", "naive"), (MODEL, "model")):126 for h in (1, 2):127 res[f"supp_mean_error_{tag}_h{h}"] = round(float(fc[(fc.model == model) & (fc.horizon == h)]["error"].mean()), 10)128# Diebold-Mariano-type comparison (squared-error loss), simple t-stat on loss differential, small-sample caveat129from scipy import stats130for h in (1, 2):131 en = fc[(fc.model == "naive") & (fc.horizon == h)]["error"].values132 em = fc[(fc.model == MODEL) & (fc.horizon == h)]["error"].values133 d = en ** 2 - em ** 2134 n = len(d)135 # HAC variance with h-1 lags (rectangular), fallback to plain variance if non-positive136 g0 = np.mean((d - d.mean()) ** 2)137 v = g0 + sum(2 * np.mean((d[k:] - d.mean()) * (d[:-k] - d.mean())) for k in range(1, h))138 if v <= 0:139 v = g0140 dm = float(d.mean() / np.sqrt(v / n))141 hln = dm * np.sqrt((n + 1 - 2 * h + h * (h - 1) / n) / n)142 res[f"supp_dm_hln_stat_h{h}"] = round(float(hln), 10)143 res[f"supp_dm_hln_pvalue_two_sided_h{h}"] = round(float(2 * stats.t.sf(abs(hln), df=n - 1)), 10)144145with open(os.path.join(OUT, "results.json"), "w") as f:146 json.dump(res, f, indent=2)147print(json.dumps({k: v for k, v in res.items() if k.startswith(("rmse", "mae", "rel", "mean_abs", "supp_dm"))}, indent=2))148