1"""Real-time forecasts of US real GDP growth (F3-conf-02).2Run from workspace root: ~/wp29/venv/bin/python out/analysis.py3"""4import json, glob, os5import numpy as np, pandas as pd67DATA = "data"; OUT = "out"8EST_START = pd.Timestamp("1985-01-01")9COVID = [pd.Timestamp("2020-04-01"), pd.Timestamp("2020-07-01")]10MODEL = "ar1_1985_excovid"1112n_rows, issues = {}, []1314def read_series(path):15 df = pd.read_csv(path)16 n_rows[os.path.relpath(path, DATA)] = 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"{nmiss} missing values in {path}", "handling": "dropped"})23 df = df.dropna()24 if df["date"].duplicated().any():25 issues.append({"issue": f"duplicate dates in {path}", "handling": "kept last"})26 df = df.drop_duplicates("date", keep="last")27 df = df.sort_values("date").reset_index(drop=True)28 exp = pd.date_range(df["date"].iloc[0], df["date"].iloc[-1], freq="QS")29 if len(exp) != len(df):30 issues.append({"issue": f"quarterly gaps in {path}", "handling": "lag pairs require consecutive quarters"})31 return df3233files = sorted(glob.glob(os.path.join(DATA, "vintages", "*.csv")))34vint = {os.path.basename(f).split("_vintage_")[1][:-4]: read_series(f) for f in files}35cur = read_series(os.path.join(DATA, "current_vintage.csv"))36fr = pd.read_csv(os.path.join(DATA, "first_release.csv"))37n_rows["first_release.csv"] = len(fr)38fr["quarter"] = pd.to_datetime(fr["quarter"])39frmap = dict(zip(fr["quarter"], fr["first_release_value"]))4041def fit_ar1(df):42 s = df.set_index("date")["value"]43 prev = s.shift(1); prev_date = pd.Series(s.index, index=s.index).shift(1)44 d = pd.DataFrame({"y": s, "x": prev, "pd": prev_date})45 # consecutive quarter requirement46 d = d[d["pd"] == (d.index - pd.offsets.QuarterBegin(1, startingMonth=1))]47 d = d[d.index >= EST_START]48 d = d[~d.index.isin(COVID) & ~d["pd"].isin(COVID)]49 d = d.dropna()50 X = np.column_stack([np.ones(len(d)), d["x"].values])51 beta, *_ = np.linalg.lstsq(X, d["y"].values, rcond=None)52 return beta, len(d)5354rows = []55for origin, df in sorted(vint.items()):56 last_date, last_val = df["date"].iloc[-1], float(df["value"].iloc[-1])57 (c, phi), n_reg = fit_ar1(df)58 f = last_val59 for h in (1, 2):60 tgt = last_date + pd.DateOffset(months=3 * h)61 f = c + phi * f62 actual = frmap.get(tgt, np.nan)63 base = dict(origin_date=origin, target_date=tgt.strftime("%Y-%m-%d"), horizon=h,64 actual=actual, vintage_date=origin,65 train_last_date=last_date.strftime("%Y-%m-%d"), train_last_value=last_val)66 rows.append({**base, "forecast": last_val, "model": "naive", "n_train": len(df)})67 rows.append({**base, "forecast": round(float(f), 6), "model": MODEL, "n_train": n_reg,68 "ar_const": c, "ar_phi": phi})69fc = pd.DataFrame(rows)70coefs = fc[fc.model == MODEL].groupby("origin_date")[["ar_const", "ar_phi"]].first()71fc = fc.drop(columns=["ar_const", "ar_phi"])72cols = ["origin_date", "target_date", "horizon", "forecast", "actual", "model",73 "vintage_date", "train_last_date", "train_last_value", "n_train"]74fc = fc[cols]75fc.to_csv(os.path.join(OUT, "forecasts.csv"), index=False)76if fc["actual"].isna().any():77 issues.append({"issue": "missing first-release actuals", "handling": "excluded from metrics"})7879res = {"task_id": "F3-conf-02", "model_name": MODEL,80 "model_description": ("AR(1) with intercept estimated by OLS at each origin on the origin's vintage, "81 "sample 1985Q1 to last published quarter, excluding regression pairs involving "82 "2020Q2-2020Q3 (COVID outliers); h=2 forecast iterated."),83 "n_origins": len(vint), "first_origin": min(vint), "last_origin": max(vint)}84for m, tag in (("naive", "naive"), (MODEL, "model")):85 for h in (1, 2):86 e = fc[(fc.model == m) & (fc.horizon == h)].eval("forecast-actual").dropna()87 res[f"rmse_{tag}_h{h}"] = round(float(np.sqrt((e**2).mean())), 4)88 res[f"mae_{tag}_h{h}"] = round(float(e.abs().mean()), 4)89for h in (1, 2):90 res[f"rel_improvement_rmse_h{h}"] = round(1 - res[f"rmse_model_h{h}"] / res[f"rmse_naive_h{h}"], 4)91curmap = dict(zip(cur["date"], cur["value"]))92t1 = pd.to_datetime(fc[(fc.model == "naive") & (fc.horizon == 1)]["target_date"])93revs = [abs(curmap[t] - frmap[t]) for t in t1]94res["mean_abs_revision_h1_targets"] = round(float(np.mean(revs)), 4)95issues.append({"issue": "No 2025-10-31 vintage: 2025Q3 advance estimate delayed by federal shutdown; 8th origin is 2025-12-31 (gap of 5 months between last two origins)",96 "handling": "used origins as provided; horizons defined relative to last quarter in each file"})97issues.append({"issue": "Extreme COVID observations 2020Q2 (~-28%) and 2020Q3 (~+35%) dominate OLS AR estimates",98 "handling": "regression pairs with dependent or lagged value in 2020Q2-2020Q3 excluded from AR estimation (pre-specified, not tuned on outcomes)"})99res.update({"feasibility": "complete", "feasibility_reason": None, "data_issues": issues,100 "n_input_rows": n_rows,101 "n_analysed": int(fc[(fc.model == MODEL) & (fc.origin_date == max(vint))]["n_train"].iloc[0]),102 "random_seed": None})103with open(os.path.join(OUT, "results.json"), "w") as fh:104 json.dump(res, fh, indent=2)105# diagnostics (not a required deliverable)106coefs.to_csv(os.path.join(OUT, "ar_coefficients.csv"))107print(json.dumps(res, indent=2)); print(fc.to_string()); print(coefs)108print("h1 target revisions:", dict(zip([t.strftime('%Y-%m-%d') for t in t1], [round(r,2) for r in revs])))109110# ---- supplementary: uncertainty and sensitivity (written to out/supplementary.json) ----111from scipy import stats112sup = {}113for h in (1, 2):114 a = fc[(fc.model == "naive") & (fc.horizon == h)].reset_index(drop=True)115 b = fc[(fc.model == MODEL) & (fc.horizon == h)].reset_index(drop=True)116 d = (a.forecast - a.actual) ** 2 - (b.forecast - b.actual) ** 2117 n = len(d); dbar = d.mean()118 # HAC variance with h-1 lags (rectangular), Harvey-Leybourne-Newbold correction119 g = [((d - dbar) * (d - dbar).shift(k)).sum() / n for k in range(h)]120 v = (g[0] + 2 * sum(g[1:])) / n121 dm = dbar / np.sqrt(v) if v > 0 else np.nan122 k = np.sqrt((n + 1 - 2 * h + h * (h - 1) / n) / n)123 stat = dm * k124 sup[f"h{h}"] = {"mean_sq_loss_diff_naive_minus_model": round(float(dbar), 4),125 "dm_hln_stat": round(float(stat), 4),126 "dm_pvalue_two_sided_t(n-1)": round(float(2 * stats.t.sf(abs(stat), n - 1)), 4),127 "n_origins_model_abs_error_smaller": int(((b.forecast - b.actual).abs() < (a.forecast - a.actual).abs()).sum()),128 "n": n}129# sensitivity: AR(1) without excluding COVID quarters130COVID_SAVE = COVID[:]; COVID.clear()131sens = {1: [], 2: []}132for origin, df in sorted(vint.items()):133 (c, phi), _ = fit_ar1(df); f = float(df["value"].iloc[-1]); ld = df["date"].iloc[-1]134 for h in (1, 2):135 f = c + phi * f136 sens[h].append(f - frmap[ld + pd.DateOffset(months=3 * h)])137COVID.extend(COVID_SAVE)138for h in (1, 2):139 e = np.array(sens[h])140 sup[f"sensitivity_ar1_incl_covid_h{h}"] = {"rmse": round(float(np.sqrt((e**2).mean())), 4),141 "mae": round(float(np.abs(e).mean()), 4)}142sup["h1_target_abs_revisions"] = dict(zip([t.strftime('%Y-%m-%d') for t in t1], [round(float(r), 2) for r in revs]))143with open(os.path.join(OUT, "supplementary.json"), "w") as fh:144 json.dump(sup, fh, indent=2)145print(json.dumps(sup, indent=2))146