1"""F3-conf-02: real-time forecasts of US real GDP growth (A191RL1Q225SBEA).2Run from workspace root: ~/wp29/venv/bin/python out/analysis.py3"""4import json, glob, os5import numpy as np6import pandas as pd78DATA = "data"; OUT = "out"9MODEL = "ar1_1985_excovid"10EST_START = pd.Timestamp("1985-01-01")11COVID = (pd.Timestamp("2020-01-01"), pd.Timestamp("2021-01-01")) # quarters excluded from regression pairs12issues = []13n_rows = {}1415def read_series(path):16 df = pd.read_csv(path)17 n_rows[os.path.relpath(path, DATA)] = len(df)18 df.columns = ["date", "value"]19 df["date"] = pd.to_datetime(df["date"])20 df["value"] = pd.to_numeric(df["value"], errors="coerce")21 nbad = int(df["value"].isna().sum())22 if nbad:23 issues.append({"issue": f"{nbad} non-numeric/missing values in {path}", "handling": "dropped"})24 df = df.dropna()25 if df["date"].duplicated().any():26 issues.append({"issue": f"duplicate dates in {path}", "handling": "kept last"})27 df = df.drop_duplicates("date", keep="last")28 df = df.sort_values("date").set_index("date")["value"]29 # check regular quarterly spacing30 exp = pd.date_range(df.index[0], df.index[-1], freq="QS")31 if len(exp) != len(df):32 issues.append({"issue": f"gaps in quarterly index of {path}", "handling": "noted; AR pairs use consecutive quarters only"})33 return df3435fr = pd.read_csv(os.path.join(DATA, "first_release.csv"))36n_rows["first_release.csv"] = len(fr)37fr["quarter"] = pd.to_datetime(fr["quarter"])38fr_map = fr.set_index("quarter")["first_release_value"]39cur = read_series(os.path.join(DATA, "current_vintage.csv"))4041files = sorted(glob.glob(os.path.join(DATA, "vintages", "A191RL1Q225SBEA_vintage_*.csv")))42rows = []43last_n = None44for f in files:45 vd = os.path.basename(f).split("_vintage_")[1][:-4]46 y = read_series(f)47 last_date, last_val = y.index[-1], float(y.iloc[-1])48 targets = [last_date + pd.DateOffset(months=3 * h) for h in (1, 2)]49 # AR(1) by OLS: y_t = c + phi*y_{t-1} + e, t >= 1985Q1, consecutive quarters, excluding50 # pairs where y_t or y_{t-1} falls in 2020Q1..2021Q1 (pandemic outliers)51 ys = y[y.index >= EST_START - pd.DateOffset(months=3)]52 prev = ys.shift(1); prev_d = pd.Series(ys.index, index=ys.index).shift(1)53 ok = (ys.index >= EST_START) & prev.notna().values54 ok &= ((ys.index - pd.DatetimeIndex(prev_d.fillna(ys.index[0]))).days.values > 80) & \55 ((ys.index - pd.DatetimeIndex(prev_d.fillna(ys.index[0]))).days.values < 100)56 inC = lambda d: (d >= COVID[0]) & (d <= COVID[1])57 ok &= ~inC(ys.index) & ~inC(pd.DatetimeIndex(prev_d.fillna(ys.index[0])))58 Y = ys.values[ok]; X = np.column_stack([np.ones(ok.sum()), prev.values[ok]])59 c, phi = np.linalg.lstsq(X, Y, rcond=None)[0]60 f1 = c + phi * last_val; f2 = c + phi * f161 for h, t, fm in ((1, targets[0], f1), (2, targets[1], f2)):62 act = float(fr_map[t]) if t in fr_map.index else np.nan63 base = dict(origin_date=vd, target_date=t.strftime("%Y-%m-%d"), horizon=h, actual=act,64 vintage_date=vd, train_last_date=last_date.strftime("%Y-%m-%d"),65 train_last_value=last_val)66 rows.append({**base, "forecast": last_val, "model": "naive", "n_train": len(y)})67 rows.append({**base, "forecast": round(float(fm), 6), "model": MODEL, "n_train": int(ok.sum())})68 last_n = len(y)69 print(vd, "c=%.3f phi=%.3f n=%d" % (c, phi, ok.sum()))7071fc = pd.DataFrame(rows)[["origin_date", "target_date", "horizon", "forecast", "actual", "model",72 "vintage_date", "train_last_date", "train_last_value", "n_train"]]73fc.to_csv(os.path.join(OUT, "forecasts.csv"), index=False)74if fc["actual"].isna().any():75 issues.append({"issue": "missing first-release actuals for some targets", "handling": "excluded from metrics"})7677def metrics(m, h):78 d = fc[(fc.model == m) & (fc.horizon == h)].dropna(subset=["actual"])79 e = d.forecast - d.actual80 return float(np.sqrt((e ** 2).mean())), float(e.abs().mean())8182res = {"task_id": "F3-conf-02", "model_name": MODEL,83 "model_description": ("AR(1) with intercept estimated by OLS on each origin's vintage, sample 1985Q1 to the last "84 "published quarter, excluding regression pairs touching 2020Q1-2021Q1; h=2 forecast by iteration."),85 "n_origins": len(files),86 "first_origin": fc.origin_date.min(), "last_origin": fc.origin_date.max()}87for m, tag in (("naive", "naive"), (MODEL, "model")):88 for h in (1, 2):89 r, a = metrics(m, h)90 res[f"rmse_{tag}_h{h}"] = round(r, 6); res[f"mae_{tag}_h{h}"] = round(a, 6)91for h in (1, 2):92 res[f"rel_improvement_rmse_h{h}"] = round(1 - res[f"rmse_model_h{h}"] / res[f"rmse_naive_h{h}"], 6)93t1 = pd.to_datetime(fc[(fc.horizon == 1) & (fc.model == "naive")].target_date)94rev = [abs(cur[t] - fr_map[t]) for t in t1]95res["mean_abs_revision_h1_targets"] = round(float(np.mean(rev)), 6)96res["feasibility"] = "complete"; res["feasibility_reason"] = None97issues.append({"issue": "2025Q3 advance estimate delayed by the Oct-Nov 2025 federal shutdown; 8th origin is 2025-12-31 and the 2025Q4 first release is dated 2026-02-28",98 "handling": "used origins as supplied; spacing between last two origins is 5 months, targets still defined relative to last quarter in the vintage"})99issues.append({"issue": "values of earlier quarters differ across vintages (revisions, incl. annual updates)",100 "handling": "each origin uses only its own vintage file; current vintage used only for the revision statistic"})101res["data_issues"] = issues102res["n_input_rows"] = n_rows103res["n_analysed"] = last_n104res["random_seed"] = None105with open(os.path.join(OUT, "results.json"), "w") as fh:106 json.dump(res, fh, indent=2)107print(json.dumps(res, indent=2))108print(fc.to_string())109print("h1 revisions:", dict(zip(t1.dt.strftime("%Y-%m-%d"), rev)))110111# ---- Uncertainty: loss differentials, sign counts, simple DM-type t-stat (HLN small-sample correction) ----112from scipy import stats113unc = {}114for h in (1, 2):115 n_ = fc[(fc.model == "naive") & (fc.horizon == h)].reset_index(drop=True)116 m_ = fc[(fc.model == MODEL) & (fc.horizon == h)].reset_index(drop=True)117 d = (n_.forecast - n_.actual) ** 2 - (m_.forecast - m_.actual) ** 2118 T = len(d); dbar = d.mean()119 g0 = ((d - dbar) ** 2).mean()120 lr = g0 + (2 * ((d - dbar).iloc[1:].values * (d - dbar).iloc[:-1].values).sum() / T if h == 2 else 0)121 lr = lr if lr > 0 else g0122 dm = dbar / np.sqrt(lr / T)123 hln = dm * np.sqrt((T + 1 - 2 * h + h * (h - 1) / T) / T)124 unc[f"h{h}"] = {"mean_sq_loss_diff_naive_minus_model": round(float(dbar), 4),125 "n_origins_model_better": int((d > 0).sum()), "T": T,126 "hln_dm_stat": round(float(hln), 3),127 "p_two_sided_t(T-1)": round(float(2 * stats.t.sf(abs(hln), T - 1)), 3),128 "loss_diffs": [round(float(x), 3) for x in d]}129with open(os.path.join(OUT, "uncertainty.json"), "w") as fh:130 json.dump(unc, fh, indent=2)131print(json.dumps(unc, indent=2))132