1"""Real-time forecasts of US real GDP growth (F3-conf-02).2Reads only data/; writes out/results.json, out/forecasts.csv, out/exclusions.csv. Deterministic (no randomness)."""3import glob, json, os4import numpy as np, pandas as pd5from scipy import stats67ROOT = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))8D = os.path.join(ROOT, "data"); O = os.path.join(ROOT, "out")9MODEL = "ar1_1985_excovid"10EST_START = pd.Timestamp("1985-01-01") # first quarter allowed in the estimation sample (as lag)11COVID = {pd.Timestamp("2020-04-01"), pd.Timestamp("2020-07-01")} # pandemic outlier quarters1213def read_series(path):14 df = pd.read_csv(path)15 df.columns = ["date", "value"]16 df["date"] = pd.to_datetime(df["date"])17 return df.sort_values("date").reset_index(drop=True)1819def fit_ar1(df, ledger=None):20 """OLS AR(1) with intercept on y_t ~ y_{t-1}, t-1 >= EST_START, dropping pairs touching COVID quarters."""21 y = df.set_index("date")["value"]22 pairs = pd.DataFrame({"y": y, "ylag": y.shift(1), "lagdate": pd.Series(y.index, index=y.index).shift(1)}).dropna()23 n0 = len(df)24 p1 = pairs[pairs["lagdate"] >= EST_START]25 p2 = p1[~(p1.index.isin(list(COVID)) | p1["lagdate"].isin(list(COVID)))]26 if ledger is not None:27 ledger += [("1", "regression pairs (y_t, y_{t-1}) with y_{t-1} dated >= 1985-01-01 (pre-1985 quarters and the first, lag-less row dropped)", n0, len(p1), n0 - len(p1),28 "post-Great-Moderation estimation sample; avoids high-volatility early postwar regime"),29 ("2", "drop pairs in which y_t or y_{t-1} is 2020Q2 or 2020Q3", len(p1), len(p2), len(p1) - len(p2),30 "pandemic collapse/rebound (about -28/+35 % SAAR) would dominate OLS and bias the AR coefficient")]31 X = np.column_stack([np.ones(len(p2)), p2["ylag"].values])32 beta, *_ = np.linalg.lstsq(X, p2["y"].values, rcond=None)33 return beta, len(p2)3435def main():36 fr = pd.read_csv(os.path.join(D, "first_release.csv"), parse_dates=["quarter", "first_release_vintage_date"])37 frmap = dict(zip(fr["quarter"], fr["first_release_value"]))38 cur = read_series(os.path.join(D, "current_vintage.csv"))39 curmap = dict(zip(cur["date"], cur["value"]))40 files = sorted(glob.glob(os.path.join(D, "vintages", "*.csv")))41 n_input_rows = {"current_vintage.csv": int(len(cur)), "first_release.csv": int(len(fr))}42 rows, issues, ledger, coefs = [], [], [], {}43 for i, f in enumerate(files):44 origin = os.path.basename(f).split("_vintage_")[1].replace(".csv", "")45 v = read_series(f)46 n_input_rows[os.path.relpath(f, D)] = int(len(v))47 last_date, last_val = v["date"].iloc[-1], float(v["value"].iloc[-1])48 assert v["value"].notna().all() and not v["date"].duplicated().any()49 is_last = i == len(files) - 150 beta, n_ar = fit_ar1(v, ledger if is_last else None)51 coefs[origin] = [round(float(b), 6) for b in beta]52 if is_last:53 n_analysed = n_ar54 fc_ar, prev = {}, last_val55 for h in (1, 2):56 prev = beta[0] + beta[1] * prev57 fc_ar[h] = prev58 for h in (1, 2):59 tgt = last_date + pd.DateOffset(months=3 * h)60 act = frmap.get(tgt, np.nan)61 if np.isnan(act):62 issues.append({"issue": f"no first release for {tgt.date()}", "handling": "target excluded"})63 for m, fc, ntr in (("naive", last_val, len(v)), (MODEL, fc_ar[h], n_ar)):64 rows.append(dict(origin_date=origin, target_date=str(tgt.date()), horizon=h, forecast=round(float(fc), 10),65 actual=act, model=m, vintage_date=origin, train_last_date=str(last_date.date()),66 train_last_value=last_val, n_train=int(ntr)))67 fc = pd.DataFrame(rows)68 fc.to_csv(os.path.join(O, "forecasts.csv"), index=False)69 pd.DataFrame(ledger, columns=["step", "rule", "n_before", "n_after", "n_removed", "rationale"]).to_csv(70 os.path.join(O, "exclusions.csv"), index=False)7172 res = {"task_id": "F3-conf-02", "model_name": MODEL,73 "model_description": ("AR(1) with intercept estimated by OLS separately at each origin on that origin's vintage, "74 "sample of pairs (y_t, y_{t-1}) with y_{t-1} from 1985Q1 to the vintage's last quarter, excluding pairs "75 "involving 2020Q2 or 2020Q3; h=2 forecast by iterating the h=1 forecast."),76 "n_origins": len(files), "first_origin": fc["origin_date"].min(), "last_origin": fc["origin_date"].max()}77 for m, tag in (("naive", "naive"), (MODEL, "model")):78 for h in (1, 2):79 s = fc[(fc.model == m) & (fc.horizon == h)].dropna(subset=["actual"])80 e = s["forecast"] - s["actual"]81 res[f"rmse_{tag}_h{h}"] = float(np.sqrt(np.mean(e ** 2)))82 res[f"mae_{tag}_h{h}"] = float(np.mean(np.abs(e)))83 for h in (1, 2):84 res[f"rel_improvement_rmse_h{h}"] = 1 - res[f"rmse_model_h{h}"] / res[f"rmse_naive_h{h}"]85 for h in (1, 2): # Diebold-Mariano (squared loss) with Harvey-Leybourne-Newbold small-sample correction86 a = fc[(fc.model == "naive") & (fc.horizon == h)].reset_index(drop=True)87 b = fc[(fc.model == MODEL) & (fc.horizon == h)].reset_index(drop=True)88 d = ((a.forecast - a.actual) ** 2 - (b.forecast - b.actual) ** 2).values89 n = len(d); dc = d - d.mean()90 lrv = np.mean(dc ** 2) + 2 * sum(np.sum(dc[k:] * dc[:-k]) / n for k in range(1, h))91 dm = d.mean() / np.sqrt(lrv / n)92 hln = dm * np.sqrt((n + 1 - 2 * h + h * (h - 1) / n) / n)93 res[f"dm_hln_stat_h{h}"] = float(hln)94 res[f"dm_hln_pvalue_two_sided_h{h}"] = float(2 * stats.t.sf(abs(hln), df=n - 1))95 t1 = fc[(fc.model == "naive") & (fc.horizon == 1)]["target_date"]96 revs = [abs(curmap[pd.Timestamp(t)] - frmap[pd.Timestamp(t)]) for t in t1]97 res["mean_abs_revision_h1_targets"] = float(np.mean(revs))98 issues.append({"issue": "2025Q3 advance estimate delayed by Oct-Nov 2025 federal shutdown; 8th origin is 2025-12-31, not 2025-10-31",99 "handling": "used the provided 2025-12-31 vintage as the 8th origin; targets remain 2025Q4 (h=1) and 2026Q1 (h=2)"})100 issues.append({"issue": "pandemic outliers 2020Q2/2020Q3 in every vintage",101 "handling": "excluded from AR(1) estimation pairs (see exclusions.csv); not relevant to the naive forecast"})102 res.update({"feasibility": "complete" if fc["actual"].notna().all() else "partial",103 "feasibility_reason": None if fc["actual"].notna().all() else "missing first releases for some targets",104 "data_issues": issues, "n_input_rows": n_input_rows, "n_analysed": int(n_analysed), "random_seed": None,105 "ar1_coefficients_by_origin": coefs})106 with open(os.path.join(O, "results.json"), "w") as fh:107 json.dump(res, fh, indent=2)108 print(json.dumps({k: v for k, v in res.items() if k.startswith(("rmse", "mae", "rel", "mean"))}, indent=2))109110if __name__ == "__main__":111 main()112