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.py3"""4import json, glob, os5import numpy as np, pandas as pd6from scipy import stats78DATA = "data"; OUT = "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"]11MODEL = "ar1_covid_excl"12# Equations whose dependent variable or lag falls in 2020Q2/2020Q3 are dropped13COVID = pd.to_datetime(["2020-04-01", "2020-07-01"])1415n_rows, issues = {}, []1617def read_series(rel):18 d = pd.read_csv(os.path.join(DATA, rel))19 n_rows[rel] = len(d)20 d.columns = ["date", "value"]21 d["date"] = pd.to_datetime(d["date"])22 d["value"] = pd.to_numeric(d["value"], errors="coerce")23 exp = pd.date_range(d.date.iloc[0], d.date.iloc[-1], freq="QS")24 if d.value.isna().any() or d.date.duplicated().any() or len(exp) != len(d):25 issues.append({"issue": f"{rel}: missing/duplicate/gapped quarters",26 "handling": "rows with missing values dropped"})27 d = d.dropna()28 return d.set_index("date")["value"]2930fr = pd.read_csv(os.path.join(DATA, "first_release.csv"))31n_rows["first_release.csv"] = len(fr)32fr["quarter"] = pd.to_datetime(fr["quarter"])33first = fr.set_index("quarter")["first_release_value"]34cur = read_series("current_vintage.csv") # used ONLY for the revision statistic3536def fit_ar1(y):37 """OLS AR(1) with intercept; drop equations touching COVID quarters."""38 yt, yl = y.values[1:], y.values[:-1]39 dt, dl = y.index[1:], y.index[:-1]40 keep = ~(dt.isin(COVID) | dl.isin(COVID))41 X = np.column_stack([np.ones(keep.sum()), yl[keep]])42 b, *_ = np.linalg.lstsq(X, yt[keep], rcond=None)43 return b, int(keep.sum())4445rows, fit_info = [], []46for o in ORIGINS:47 rel = f"vintages/A191RL1Q225SBEA_vintage_{o}.csv"48 y = read_series(rel)49 last_d, last_v = y.index[-1], float(y.iloc[-1])50 (c, phi), neq = fit_ar1(y)51 fit_info.append({"origin": o, "const": c, "phi": phi, "n_eq": neq})52 f = last_v53 for h in (1, 2):54 tgt = last_d + pd.DateOffset(months=3 * h)55 f = c + phi * f56 act = float(first[tgt]) if tgt in first.index else np.nan57 base = dict(origin_date=o, target_date=tgt.strftime("%Y-%m-%d"), horizon=h,58 actual=act, vintage_date=o, train_last_date=last_d.strftime("%Y-%m-%d"),59 train_last_value=last_v, n_train=len(y))60 rows.append({**base, "forecast": last_v, "model": "naive"})61 rows.append({**base, "forecast": round(float(f), 6), "model": MODEL})62 n_last = len(y)6364fc = pd.DataFrame(rows)[["origin_date", "target_date", "horizon", "forecast", "actual", "model",65 "vintage_date", "train_last_date", "train_last_value", "n_train"]]66fc = fc.sort_values(["origin_date", "horizon", "model"]).reset_index(drop=True)67fc.to_csv(os.path.join(OUT, "forecasts.csv"), index=False)68fc["err"] = fc["actual"] - fc["forecast"]6970res = {}71for h in (1, 2):72 for m, lab in (("naive", "naive"), (MODEL, "model")):73 e = fc[(fc.horizon == h) & (fc.model == m)]["err"]74 res[f"rmse_{lab}_h{h}"] = float(np.sqrt((e ** 2).mean()))75 res[f"mae_{lab}_h{h}"] = float(e.abs().mean())76 res[f"rel_improvement_rmse_h{h}"] = 1 - res[f"rmse_model_h{h}"] / res[f"rmse_naive_h{h}"]7778t1 = pd.to_datetime(fc[(fc.horizon == 1) & (fc.model == "naive")]["target_date"])79rev = [abs(float(cur[t]) - float(first[t])) for t in t1]80mean_rev = float(np.mean(rev))8182# Diebold-Mariano with Harvey-Leybourne-Newbold small-sample correction (squared loss)83dm = {}84for h in (1, 2):85 en = fc[(fc.horizon == h) & (fc.model == "naive")]["err"].values86 em = fc[(fc.horizon == h) & (fc.model == MODEL)]["err"].values87 d = en ** 2 - em ** 2; T = len(d); dbar = d.mean()88 g = [np.mean((d[k:] - dbar) * (d[:T - k] - dbar)) for k in range(h)]89 v = (g[0] + 2 * sum(g[1:])) / T90 if v <= 0:91 v = g[0] / T92 stat = dbar / np.sqrt(v) * np.sqrt((T + 1 - 2 * h + h * (h - 1) / T) / T)93 dm[f"h{h}"] = {"stat": float(stat), "p_two_sided": float(2 * stats.t.sf(abs(stat), T - 1)),94 "mean_loss_diff_naive_minus_model": float(dbar)}9596# Descriptive robustness: per-origin wins and leave-one-origin-out relative improvement97robust = {}98for h in (1, 2):99 a = fc[(fc.horizon == h)].pivot(index="origin_date", columns="model", values="err")100 wins = int((a[MODEL].abs() < a["naive"].abs()).sum())101 loo = {o: float(1 - np.sqrt((a[MODEL].drop(o) ** 2).mean()) / np.sqrt((a["naive"].drop(o) ** 2).mean()))102 for o in a.index}103 robust[f"h{h}"] = {"n_origins_model_abs_error_smaller": wins, "loo_rel_improvement_rmse": loo}104105issues.append({"issue": "2020Q2 (-28) and 2020Q3 (+35) COVID outliers in every vintage would dominate OLS AR estimates",106 "handling": "AR(1) equations whose dependent value or lag is 2020Q2 or 2020Q3 are excluded from estimation (equivalent to impulse dummies for 2020Q2-Q4); naive unaffected"})107issues.append({"issue": "2025Q3 advance estimate delayed by the Oct-Nov 2025 federal shutdown; 8th origin is 2025-12-31, and the h=1 target 2025Q4 first release is dated 2026-02-28",108 "handling": "used as provided; origin spacing is irregular for the last origin"})109issues.append({"issue": "Past values differ across vintages (e.g. 2023Q4: 3.3 in 2024-01-31 vintage, 3.4 in 2024-04-30, 3.2 in 2024-10-31) and annual revisions alter whole history",110 "handling": "each origin uses only its own vintage file; actuals are first releases"})111112out = {"task_id": "F3-conf-02", "model_name": MODEL,113 "model_description": ("AR(1) with intercept estimated by OLS on the full origin vintage (1947Q2 to last published quarter), "114 "excluding equations whose dependent value or lag is 2020Q2 or 2020Q3, re-estimated at each origin; "115 "h=2 forecast iterated."),116 "n_origins": len(ORIGINS), "first_origin": ORIGINS[0], "last_origin": ORIGINS[-1],117 **{k: res[k] for k in ["rmse_naive_h1", "mae_naive_h1", "rmse_naive_h2", "mae_naive_h2",118 "rmse_model_h1", "mae_model_h1", "rmse_model_h2", "mae_model_h2",119 "rel_improvement_rmse_h1", "rel_improvement_rmse_h2"]},120 "mean_abs_revision_h1_targets": mean_rev,121 "feasibility": "complete", "feasibility_reason": None,122 "data_issues": issues, "n_input_rows": n_rows, "n_analysed": n_last, "random_seed": None,123 "diebold_mariano_hln": dm, "robustness": robust, "ar_fits": fit_info,124 "revisions_h1_targets": {t.strftime("%Y-%m-%d"): {"first": float(first[t]), "current": float(cur[t])} for t in t1}}125with open(os.path.join(OUT, "results.json"), "w") as fh:126 json.dump(out, fh, indent=2)127print(json.dumps(out, indent=2))128print(fc.to_string())129