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 pd67DATA = "data"; OUT = "out"8MODEL = "ar1_1985_xcovid"9EST_START = pd.Timestamp("1985-01-01")10COVID = {pd.Timestamp("2020-04-01"), pd.Timestamp("2020-07-01")}11issues = []12n_rows = {}1314def read_series(rel):15 df = pd.read_csv(os.path.join(DATA, rel))16 n_rows[rel] = int(len(df))17 df.columns = ["date", "value"]18 df["date"] = pd.to_datetime(df["date"])19 df["value"] = pd.to_numeric(df["value"], errors="coerce")20 nmiss = int(df["value"].isna().sum())21 if nmiss:22 issues.append({"issue": f"{rel}: {nmiss} missing/non-numeric values",23 "handling": "dropped"})24 df = df.dropna()25 s = df.set_index("date")["value"].sort_index()26 exp = pd.date_range(s.index[0], s.index[-1], freq="QS")27 if len(exp) != len(s) or not (exp == s.index).all():28 issues.append({"issue": f"{rel}: irregular/missing quarters", "handling": "noted; AR uses only consecutive pairs"})29 if s.index.duplicated().any():30 issues.append({"issue": f"{rel}: duplicate dates", "handling": "kept last"})31 s = s[~s.index.duplicated(keep="last")]32 return s3334fr = pd.read_csv(os.path.join(DATA, "first_release.csv"))35n_rows["first_release.csv"] = int(len(fr))36fr["quarter"] = pd.to_datetime(fr["quarter"])37fr_map = dict(zip(fr["quarter"], fr["first_release_value"]))38cur = read_series("current_vintage.csv")3940def fit_ar1(s):41 """OLS AR(1) with intercept on vintage s, from EST_START, dropping pairs42 where y_t or y_{t-1} is a COVID quarter (2020Q2, 2020Q3)."""43 s = s[s.index >= EST_START - pd.offsets.QuarterBegin(startingMonth=1)]44 y = s.values[1:]; x = s.values[:-1]45 dt = s.index[1:]; dl = s.index[:-1]46 consecutive = np.array([(a + pd.offsets.QuarterBegin(startingMonth=1)) == b for a, b in zip(dl, dt)])47 keep = consecutive & np.array([(a not in COVID) and (b not in COVID) for a, b in zip(dl, dt)]) & (dt >= EST_START)48 y, x = y[keep], x[keep]49 X = np.column_stack([np.ones_like(x), x])50 beta, *_ = np.linalg.lstsq(X, y, rcond=None)51 return beta, int(keep.sum())5253files = sorted(glob.glob(os.path.join(DATA, "vintages", "A191RL1Q225SBEA_vintage_*.csv")))54rows = []55for f in files:56 rel = os.path.relpath(f, DATA)57 vdate = os.path.basename(f).split("_")[-1].replace(".csv", "")58 s = read_series(rel)59 last_d, last_v = s.index[-1], float(s.iloc[-1])60 if last_d in fr_map and abs(fr_map[last_d] - last_v) > 1e-9:61 issues.append({"issue": f"{rel}: last value {last_v} differs from first release {fr_map[last_d]}",62 "handling": "vintage value used as training data (per rules)"})63 beta, n_est = fit_ar1(s)64 f_prev = last_v65 for h in (1, 2):66 tgt = last_d + pd.DateOffset(months=3 * h)67 f_prev = beta[0] + beta[1] * f_prev68 actual = fr_map.get(tgt, np.nan)69 base = dict(origin_date=vdate, target_date=tgt.strftime("%Y-%m-%d"), horizon=h,70 actual=actual, vintage_date=vdate,71 train_last_date=last_d.strftime("%Y-%m-%d"), train_last_value=last_v)72 rows.append({**base, "forecast": last_v, "model": "naive", "n_train": int(len(s))})73 rows.append({**base, "forecast": round(float(f_prev), 6), "model": MODEL, "n_train": n_est})74 print(vdate, "beta", beta.round(3), "n_est", n_est)7576fc = pd.DataFrame(rows)[["origin_date", "target_date", "horizon", "forecast", "actual", "model",77 "vintage_date", "train_last_date", "train_last_value", "n_train"]]78fc.to_csv(os.path.join(OUT, "forecasts.csv"), index=False)79if fc["actual"].isna().any():80 issues.append({"issue": "missing first-release actuals", "handling": "excluded from metrics"})8182def err(model, h):83 d = fc[(fc.model == model) & (fc.horizon == h)].dropna(subset=["actual"])84 e = d["forecast"] - d["actual"]85 return float(np.sqrt((e ** 2).mean())), float(e.abs().mean())8687res = {"task_id": "F3-conf-02", "model_name": MODEL,88 "model_description": ("AR(1) with intercept estimated by OLS on each origin's vintage from 1985Q1 to the "89 "last published quarter, excluding regression pairs involving the 2020Q2/2020Q3 "90 "COVID quarters; h=2 forecast by iterating the AR(1)."),91 "n_origins": len(files),92 "first_origin": fc.origin_date.min(), "last_origin": fc.origin_date.max()}93for m, tag in (("naive", "naive"), (MODEL, "model")):94 for h in (1, 2):95 r, a = err(m, h)96 res[f"rmse_{tag}_h{h}"] = round(r, 6); res[f"mae_{tag}_h{h}"] = round(a, 6)97for h in (1, 2):98 res[f"rel_improvement_rmse_h{h}"] = round(1 - res[f"rmse_model_h{h}"] / res[f"rmse_naive_h{h}"], 6)99t1 = pd.to_datetime(fc[(fc.horizon == 1) & (fc.model == "naive")]["target_date"])100revs = [abs(cur[t] - fr_map[t]) for t in t1]101res["mean_abs_revision_h1_targets"] = round(float(np.mean(revs)), 6)102# Origin-spacing check (shutdown-delayed 2025Q3 advance estimate)103od = pd.to_datetime(sorted(fc.origin_date.unique()))104gaps = [(b - a).days for a, b in zip(od[:-1], od[1:])]105for a, b, g in zip(od[:-1], od[1:], gaps):106 if g > 100:107 issues.append({"issue": f"irregular origin spacing: {a.date()} -> {b.date()} ({g} days), due to delayed BEA advance estimate (2025 shutdown)",108 "handling": "kept as given; h is defined relative to last quarter in the vintage, so targets stay consecutive quarters"})109# Small-sample Diebold-Mariano (Harvey-Leybourne-Newbold) test, squared-error loss110from scipy import stats111for h in (1, 2):112 a = fc[(fc.model == "naive") & (fc.horizon == h)].sort_values("origin_date")113 b = fc[(fc.model == MODEL) & (fc.horizon == h)].sort_values("origin_date")114 d = (a.forecast.values - a.actual.values) ** 2 - (b.forecast.values - b.actual.values) ** 2115 T = len(d); dbar = d.mean(); dc = d - dbar116 lrv = (dc @ dc) / T + 2 * sum((dc[k:] @ dc[:-k]) / T for k in range(1, h))117 dm = dbar / np.sqrt(lrv / T)118 hln = dm * np.sqrt((T + 1 - 2 * h + h * (h - 1) / T) / T)119 res[f"dm_hln_stat_h{h}"] = round(float(hln), 4)120 res[f"dm_hln_pvalue_h{h}"] = round(float(2 * stats.t.sf(abs(hln), T - 1)), 4)121res["feasibility"] = "complete"; res["feasibility_reason"] = None122res["data_issues"] = issues123res["n_input_rows"] = n_rows124res["n_analysed"] = int(len(read_series(os.path.relpath(files[-1], DATA))))125res["random_seed"] = None126with open(os.path.join(OUT, "results.json"), "w") as fh:127 json.dump(res, fh, indent=2)128print(json.dumps(res, indent=2))129