1"""T1-forecast: real-time forecasting of US retail-sales YoY growth from ALFRED vintages.2Run from workspace root: ~/wp29/venv/bin/python out/analysis.py3Reads only data/, writes out/results.json, out/forecasts.csv, out/table1.csv,4out/figure1.png, out/exclusions.csv. Deterministic (bootstrap seed fixed)."""5import json, os6import numpy as np7import pandas as pd8from scipy.stats import norm9import matplotlib10matplotlib.use("Agg")11import matplotlib.pyplot as plt1213DATA, OUT = "data", "out"14SEED = 2026103015ORIGINS = ["2025-10-31", "2025-11-30", "2025-12-31", "2026-01-31", "2026-02-28",16 "2026-03-31", "2026-04-30", "2026-05-31", "2026-06-30", "2026-07-31"]17HORIZONS = [1, 3]18MODELS = ["naive", "snaive", "ar1", "model"]19n_input_rows = {}202122def read_series(rel):23 d = pd.read_csv(os.path.join(DATA, rel), parse_dates=["observation_date"])24 n_input_rows[rel] = int(len(d))25 return d.set_index("observation_date").iloc[:, 0].astype(float)262728def yoy(level):29 """g_m = 100 ln(L_m / L_{m-12}) on a gap-free monthly series."""30 s = level.asfreq("MS")31 assert not s.isna().any(), "gap in level series"32 return (100 * np.log(s / s.shift(12))).dropna()333435def mshift(t, k):36 return t + pd.DateOffset(months=k)373839# ---------------------------------------------------------------- inputs40fr = pd.read_csv(os.path.join(DATA, "first_release.csv"), parse_dates=["target_month"])41n_input_rows["first_release.csv"] = int(len(fr))42fr = fr.set_index("target_month")4344# PAYEMS vintages and UMCSENT are candidate predictors; they are read for the45# audit row counts only and do NOT enter the chosen model.46for V in ORIGINS:47 read_series(f"vintages/PAYEMS_vintage_{V}.csv")48read_series("UMCSENT.csv")4950# ---------------------------------------------------------------- forecasts51rows, origin_info = [], []52for V in ORIGINS:53 L = read_series(f"vintages/RSAFS_vintage_{V}.csv") # only vintage V54 g = yoy(L)55 last = L.index[-1]56 assert last < pd.Timestamp(V), "vintage contains month at/after origin"57 origin_info.append({"origin": V, "last": last.strftime("%Y-%m"),58 "lag_months": (pd.Timestamp(V).year - last.year) * 12 + pd.Timestamp(V).month - last.month})59 g_last = g[last]60 # AR(1): g_m on const + g_{m-1}, 1993-02 <= m <= last61 gs = g[g.index >= "1993-01-01"]62 y = gs[gs.index >= "1993-02-01"].values63 x = gs.shift(1)[gs.index >= "1993-02-01"].values64 X = np.column_stack([np.ones_like(x), x])65 (c, phi), *_ = np.linalg.lstsq(X, y, rcond=None)66 for h in HORIZONS:67 T = mshift(last, h)68 actual = float(fr.loc[T, "yoy_first_release"])69 # naive70 f_naive = g_last71 # seasonal naive72 f_snaive = g[mshift(T, -12)]73 # AR(1) iterated74 f = g_last75 for _ in range(h):76 f = c + phi * f77 f_ar1 = f78 # own model: base-effect drift79 f_model = 100 * np.log(L[last] / L[mshift(T, -12)]) + h * g_last / 12.080 for m, fc in zip(MODELS, [f_naive, f_snaive, f_ar1, f_model]):81 rows.append({"origin_date": V, "target_date": T.strftime("%Y-%m-%d"), "horizon": h,82 "forecast": round(float(fc), 6), "actual_first_release": actual,83 "model": m, "vintage_used": V})8485fc = pd.DataFrame(rows)86assert len(fc) == 8087fc.to_csv(os.path.join(OUT, "forecasts.csv"), index=False)88fc["err"] = fc["actual_first_release"] - fc["forecast"]8990# ---------------------------------------------------------------- evaluation91res = {}92tab = []93for m in MODELS:94 for h in HORIZONS:95 e = fc[(fc.model == m) & (fc.horizon == h)]["err"].values96 rmse, mae = float(np.sqrt(np.mean(e ** 2))), float(np.mean(np.abs(e)))97 res[f"rmse_{m}_h{h}"], res[f"mae_{m}_h{h}"] = rmse, mae98 tab.append({"model": m, "horizon": h, "rmse": rmse, "mae": mae, "n": len(e)})99pd.DataFrame(tab).to_csv(os.path.join(OUT, "table1.csv"), index=False)100101102def dm_test(e_naive, e_model, h):103 d = e_naive ** 2 - e_model ** 2104 n = len(d)105 dbar = d.mean()106 dc = d - dbar107 lrv = np.sum(dc * dc) / n108 for k in range(1, h): # Bartlett, truncation lag h-1109 w = 1 - k / h110 lrv += 2 * w * np.sum(dc[k:] * dc[:-k]) / n111 stat = dbar / np.sqrt(lrv / n)112 return float(stat), float(2 * (1 - norm.cdf(abs(stat)))), float(lrv)113114115boot = {}116rng = np.random.default_rng(SEED)117for h in HORIZONS:118 sub = fc[fc.horizon == h]119 en = sub[sub.model == "naive"].sort_values("origin_date")["err"].values120 em = sub[sub.model == "model"].sort_values("origin_date")["err"].values121 stat, p, lrv = dm_test(en, em, h)122 res[f"dm_stat_h{h}"], res[f"dm_pvalue_h{h}"] = stat, p123 res[f"rel_improvement_rmse_h{h}"] = 1 - res[f"rmse_model_h{h}"] / res[f"rmse_naive_h{h}"]124 better = res[f"rmse_model_h{h}"] < res[f"rmse_naive_h{h}"]125 sig = p < 0.05126 res[f"conclusion_code_h{h}"] = ("model_better" if better else "naive_better") + \127 ("_significant" if sig else "_not_significant")128 # supplementary: moving-block bootstrap of RMSE_model - RMSE_naive (block length max(h,2))129 n, b = len(en), max(h, 2)130 starts = np.arange(n - b + 1)131 draws = []132 for _ in range(9999):133 idx = np.concatenate([np.arange(s, s + b) for s in rng.choice(starts, size=int(np.ceil(n / b)))])[:n]134 draws.append(np.sqrt(np.mean(em[idx] ** 2)) - np.sqrt(np.mean(en[idx] ** 2)))135 draws = np.array(draws)136 boot[f"h{h}"] = {"rmse_diff_model_minus_naive": float(np.sqrt(np.mean(em ** 2)) - np.sqrt(np.mean(en ** 2))),137 "ci95": [float(np.percentile(draws, 2.5)), float(np.percentile(draws, 97.5))],138 "share_draws_model_better": float(np.mean(draws < 0)),139 "block_length": b, "n_boot": 9999, "dm_lrv": lrv,140 "n_origins_model_smaller_sq_error": int(np.sum(em ** 2 < en ** 2))}141142# ---------------------------------------------------------------- revision (estimand d)143# RSAFS_current is used ONLY here; it never enters a forecast.144cur = read_series("RSAFS_current.csv")145gcur = yoy(cur)146t1 = fc[(fc.horizon == 1) & (fc.model == "naive")]147targets = pd.to_datetime(t1["target_date"]).tolist()148revs = [float(gcur[T] - fr.loc[T, "yoy_first_release"]) for T in targets]149res["mean_abs_revision_h1_targets"] = float(np.mean(np.abs(revs)))150151# ---------------------------------------------------------------- exclusions ledger152excl = pd.DataFrame([153 {"step": 1, "rule": "first_release.csv: keep only rows that are targets last(V)+h of some origin",154 "n_before": len(fr), "n_after": int(fc.target_date.nunique()),155 "n_removed": len(fr) - int(fc.target_date.nunique()),156 "rationale": "rows outside the evaluation design are not needed as actuals (not an error-driven exclusion)"},157 {"step": 2, "rule": "forecast-actual pairs for own model: drop pairs with missing forecast or actual",158 "n_before": 20, "n_after": int(fc[fc.model == "model"].dropna(subset=["forecast", "actual_first_release"]).shape[0]),159 "n_removed": 20 - int(fc[fc.model == "model"].dropna(subset=["forecast", "actual_first_release"]).shape[0]),160 "rationale": "no pair is missing; nothing removed"},161])162excl.to_csv(os.path.join(OUT, "exclusions.csv"), index=False)163164# ---------------------------------------------------------------- figure165fig, axes = plt.subplots(2, 1, figsize=(9, 7), sharex=True)166frp = fr.loc["2025-06-01":"2026-09-01", "yoy_first_release"]167for ax, h in zip(axes, HORIZONS):168 ax.plot(frp.index, frp.values, "k-o", ms=4, label="first-release actual")169 for m, sty, lab in [("naive", "s--", "naive"), ("model", "^-", "BED model"), ("ar1", "x:", "AR(1)")]:170 s = fc[(fc.horizon == h) & (fc.model == m)]171 ax.plot(pd.to_datetime(s.target_date), s.forecast, sty, ms=5, label=lab)172 ax.set_title(f"h = {h}: forecasts by target month (10 origins)")173 ax.set_ylabel("YoY growth, % (100·log)")174 ax.grid(alpha=.3)175 ax.legend(fontsize=8)176fig.tight_layout()177fig.savefig(os.path.join(OUT, "figure1.png"), dpi=130)178179# ---------------------------------------------------------------- results.json180lag_issue = [o for o in origin_info if o["lag_months"] != 1]181results = {182 "task_id": "T1-forecast",183 "model_name": "BED (base-effect drift)",184 "model_description": (185 "Base-effect drift model on the RSAFS vintage dated V only. Log level forecast "186 "ln L_{last+h} = ln L_last + h*mu, with mu = mean monthly log growth over the last 12 published months "187 "= g_last/12 (percent). Hence forecast g_T = 100*ln(L_last / L_{T-12}) + h*g_last/12 for T = last(V)+h, "188 "where L_last and L_{T-12} are read from vintage V (T-12 <= last(V) always). It exploits that 12-h of the "189 "12 monthly growth terms in the target YoY are already published at the origin (known base effects), "190 "replacing only the h unknown months by their trailing 12-month average. No parameters are estimated; "191 "estimation window = last 13 levels of the vintage. PAYEMS and UMCSENT were not used. The specification "192 "was fixed before any out-of-sample error was computed; no tuning on the evaluation origins."),193 "n_origins": len(ORIGINS),194 "first_origin": ORIGINS[0], "last_origin": ORIGINS[-1],195}196for k in ["rmse_naive_h1", "mae_naive_h1", "rmse_naive_h3", "mae_naive_h3",197 "rmse_snaive_h1", "mae_snaive_h1", "rmse_snaive_h3", "mae_snaive_h3",198 "rmse_ar1_h1", "mae_ar1_h1", "rmse_ar1_h3", "mae_ar1_h3",199 "rmse_model_h1", "mae_model_h1", "rmse_model_h3", "mae_model_h3",200 "rel_improvement_rmse_h1", "rel_improvement_rmse_h3",201 "dm_stat_h1", "dm_pvalue_h1", "dm_stat_h3", "dm_pvalue_h3",202 "conclusion_code_h1", "conclusion_code_h3", "mean_abs_revision_h1_targets"]:203 results[k] = res[k]204results.update({205 "feasibility": "complete",206 "feasibility_reason": None,207 "data_issues": [208 {"issue": "Irregular publication timing: last(V) is two months before the origin for "209 + ", ".join(f"{o['origin']} (last {o['last']})" for o in lag_issue)210 + "; vintage 2026-04-30 adds 2026-02 and 2026-03 at once, so 2026-03 is never an h=1 target (it is the h=1 target of no origin) and 2026-02 is reached from the 2026-03-31 vintage.",211 "handling": "Targets defined as last(V)+h from the vintage content, per the task rule, not from the calendar."},212 {"issue": "Benchmark/annual revision of RSAFS levels: the current vintage is ~4% lower in level and its YoY growth is "213 "systematically lower than first release; the base level used in first_release jumps between targets 2026-07 and 2026-08 "214 "(rebasing between the 2026-08-31 and 2026-09-30 vintages), affecting h=3 targets 2026-08 and 2026-09.",215 "handling": "No adjustment: actuals are first-release values by definition; the revision is quantified in mean_abs_revision_h1_targets and discussed."},216 {"issue": "UMCSENT has 210 empty values (all before 1978, quarterly survey era) and is a current (not real-time) vintage.",217 "handling": "Not used by the model; would have been usable under the stated no-revision rule."},218 {"issue": "PAYEMS vintages sometimes extend beyond last(V) of RSAFS (e.g. 2025-12-31 vintage ends 2025-11 while RSAFS ends 2025-10).",219 "handling": "Not used by the model."},220 {"issue": "first_release.csv dates the first release of 2025-08 (last month of the 2025-10-31 vintage) to the unshipped 2025-09-30 vintage; its level in the 2025-10-31 vintage still equals level_first_release.",221 "handling": "None needed; checked for consistency."},222 ],223 "n_input_rows": n_input_rows,224 "n_analysed": int(fc[fc.model == "model"].dropna(subset=["forecast", "actual_first_release"]).shape[0]),225 "random_seed": SEED,226 "supplementary": {"bootstrap_rmse_diff": boot, "origins": origin_info,227 "revisions_h1_targets": {T.strftime("%Y-%m"): r for T, r in zip(targets, revs)},228 "ar1_note": "AR(1) re-estimated by OLS at each origin on the vintage, 1993-02..last(V)"},229})230with open(os.path.join(OUT, "results.json"), "w") as f:231 json.dump(results, f, indent=2)232233print(pd.DataFrame(tab).to_string(index=False))234for h in HORIZONS:235 print(f"h={h}: rel_impr={res[f'rel_improvement_rmse_h{h}']:.4f} DM={res[f'dm_stat_h{h}']:.3f} "236 f"p={res[f'dm_pvalue_h{h}']:.4f} -> {res[f'conclusion_code_h{h}']}")237print("mean_abs_revision_h1_targets =", round(res["mean_abs_revision_h1_targets"], 4))238