1"""T1-forecast: real-time forecasting of US retail-sales YoY growth from ALFRED vintages.23Run from the workspace root: ~/wp29/venv/bin/python out/analysis.py4Reads only data/, writes out/results.json, out/forecasts.csv, out/table1.csv, out/figure1.png.5"""6import json7import os8import glob910import numpy as np11import pandas as pd12from scipy import stats13import matplotlib14matplotlib.use("Agg")15import matplotlib.pyplot as plt1617DATA = "data"18OUT = "out"19SEED = 2026103020B_BOOT = 999921BLOCK = 322ORIGINS = ["2025-10-31", "2025-11-30", "2025-12-31", "2026-01-31", "2026-02-28",23 "2026-03-31", "2026-04-30", "2026-05-31", "2026-06-30", "2026-07-31"]24HORIZONS = [1, 3]25MODELS = ["naive", "snaive", "ar1", "model"]262728def read_series(path):29 d = pd.read_csv(path, parse_dates=[0], index_col=0)30 s = d.iloc[:, 0].astype(float)31 s.index = s.index.to_period("M")32 return s333435def n_rows(path):36 with open(path) as f:37 return sum(1 for _ in f) - 1383940def yoy(level):41 return 100.0 * (np.log(level) - np.log(level.shift(12)))424344def forecasts_for_origin(V):45 """All forecasts at origin V use ONLY the RSAFS vintage dated V."""46 L = read_series(os.path.join(DATA, "vintages", f"RSAFS_vintage_{V}.csv"))47 assert L.index.is_monotonic_increasing and L.index.is_unique48 last = L.index[-1]49 g = yoy(L)50 lnL = np.log(L)51 out = {}52 # AR(1) on g over 1993-02 <= m <= last (consecutive pairs; g starts 1993-01)53 gs = g.loc[pd.Period("1993-01", "M"):last]54 y = gs.loc[pd.Period("1993-02", "M"):].values55 x = gs.shift(1).loc[pd.Period("1993-02", "M"):].values56 X = np.column_stack([np.ones_like(x), x])57 c, phi = np.linalg.lstsq(X, y, rcond=None)[0]58 # model: base-effect drift. log level path extrapolated with trailing-12m mean monthly log change59 mu = (lnL.loc[last] - lnL.loc[last - 12]) / 12.060 for h in HORIZONS:61 T = last + h62 f_naive = g.loc[last]63 f_snaive = g.loc[T - 12]64 f = g.loc[last]65 for _ in range(h):66 f = c + phi * f67 f_ar1 = f68 f_model = 100.0 * (lnL.loc[last] + h * mu - lnL.loc[T - 12])69 out[h] = dict(target=T, naive=f_naive, snaive=f_snaive, ar1=f_ar1, model=f_model)70 return last, out, (c, phi)717273def dm_test(e_naive, e_model, h):74 d = e_naive ** 2 - e_model ** 275 n = len(d)76 dbar = d.mean()77 dc = d - dbar78 lrv = np.sum(dc * dc) / n79 for k in range(1, h): # truncation lag h-1, Bartlett weights 1 - k/h80 w = 1.0 - k / h81 lrv += 2.0 * w * np.sum(dc[k:] * dc[:-k]) / n82 stat = dbar / np.sqrt(lrv / n)83 p = 2.0 * (1.0 - stats.norm.cdf(abs(stat)))84 return float(stat), float(p), float(dbar), float(lrv)858687def block_boot_rmse_diff(e_naive, e_model, rng):88 n = len(e_naive)89 nb = int(np.ceil(n / BLOCK))90 diffs = np.empty(B_BOOT)91 for b in range(B_BOOT):92 starts = rng.integers(0, n, nb)93 idx = np.concatenate([(s + np.arange(BLOCK)) % n for s in starts])[:n]94 diffs[b] = np.sqrt(np.mean(e_model[idx] ** 2)) - np.sqrt(np.mean(e_naive[idx] ** 2))95 return float(np.percentile(diffs, 2.5)), float(np.percentile(diffs, 97.5)), float(np.mean(diffs < 0))969798def main():99 os.makedirs(OUT, exist_ok=True)100 fr = pd.read_csv(os.path.join(DATA, "first_release.csv"), parse_dates=["target_month"])101 fr["T"] = fr["target_month"].dt.to_period("M")102 actual = dict(zip(fr["T"], fr["yoy_first_release"]))103104 rows = []105 lasts = {}106 arcoef = {}107 for V in ORIGINS:108 last, fc, cp = forecasts_for_origin(V)109 lasts[V] = last110 arcoef[V] = cp111 for h in HORIZONS:112 T = fc[h]["target"]113 if T not in actual:114 raise RuntimeError(f"missing first-release actual for {T}")115 for m in MODELS:116 rows.append(dict(origin_date=V, target_date=T.to_timestamp().strftime("%Y-%m-%d"),117 horizon=h, forecast=round(float(fc[h][m]), 10),118 actual_first_release=float(actual[T]), model=m, vintage_used=V))119 F = pd.DataFrame(rows)120 F["error"] = F["actual_first_release"] - F["forecast"]121 F.drop(columns="error").to_csv(os.path.join(OUT, "forecasts.csv"), index=False)122123 res = {"task_id": "T1-forecast",124 "model_name": "base-effect drift (BED)",125 "model_description": (126 "Forecast of the log level of RSAFS at T = last(V)+h from the RSAFS vintage dated V only: "127 "ln L_hat_T = ln L_last + h*mu, where mu = (ln L_last - ln L_{last-12})/12 is the trailing-12-month "128 "mean monthly log change in the same vintage (i.e. the latest YoY growth spread evenly over months). "129 "The YoY forecast is g_hat_T = 100*(ln L_hat_T - ln L_{T-12}), where L_{T-12} (always <= last(V) since h<=3) "130 "is the base-period level as published in vintage V. Equivalently: naive forecast plus the known "131 "base effect (drop-out of the monthly log changes of months last-11..last-12+h) plus h months of drift. "132 "No parameters are estimated on the evaluation sample; the specification was fixed before computing any "133 "out-of-sample error and no variant was compared. PAYEMS and UMCSENT are not used; RSAFS_current and "134 "first_release.csv never enter any forecast."),135 "n_origins": len(ORIGINS), "first_origin": ORIGINS[0], "last_origin": ORIGINS[-1]}136 t1 = []137 for m in MODELS:138 for h in HORIZONS:139 e = F[(F.model == m) & (F.horizon == h)]["error"].values140 rmse = float(np.sqrt(np.mean(e ** 2)))141 mae = float(np.mean(np.abs(e)))142 res[f"rmse_{m}_h{h}"] = rmse143 res[f"mae_{m}_h{h}"] = mae144 t1.append(dict(model=m, horizon=h, rmse=rmse, mae=mae, n=len(e)))145 pd.DataFrame(t1).to_csv(os.path.join(OUT, "table1.csv"), index=False)146147 rng = np.random.default_rng(SEED)148 extra = {}149 for h in HORIZONS:150 en = F[(F.model == "naive") & (F.horizon == h)].sort_values("origin_date")["error"].values151 em = F[(F.model == "model") & (F.horizon == h)].sort_values("origin_date")["error"].values152 res[f"rel_improvement_rmse_h{h}"] = 1.0 - res[f"rmse_model_h{h}"] / res[f"rmse_naive_h{h}"]153 stat, p, dbar, lrv = dm_test(en, em, h)154 res[f"dm_stat_h{h}"] = stat155 res[f"dm_pvalue_h{h}"] = p156 better = "model" if res[f"rmse_model_h{h}"] < res[f"rmse_naive_h{h}"] else "naive"157 sig = "significant" if p < 0.05 else "not_significant"158 res[f"conclusion_code_h{h}"] = f"{better}_better_{sig}"159 lo, hi, share = block_boot_rmse_diff(en, em, rng)160 extra[f"dm_mean_loss_diff_h{h}"] = dbar161 extra[f"dm_lrv_h{h}"] = lrv162 extra[f"boot_rmse_diff_lo_h{h}"] = lo163 extra[f"boot_rmse_diff_hi_h{h}"] = hi164 extra[f"boot_share_model_better_h{h}"] = share165 extra[f"n_model_wins_h{h}"] = int(np.sum(em ** 2 < en ** 2))166167 # estimand (d): revision of h=1 targets, current vintage vs first release168 cur = read_series(os.path.join(DATA, "RSAFS_current.csv"))169 gcur = yoy(cur)170 t_h1 = [pd.Period(t, "M") for t in F[(F.model == "naive") & (F.horizon == 1)]["target_date"].str[:7]]171 revs = np.array([abs(gcur.loc[T] - actual[T]) for T in t_h1])172 res["mean_abs_revision_h1_targets"] = float(revs.mean())173 extra["mean_signed_revision_h1_targets"] = float(np.mean([gcur.loc[T] - actual[T] for T in t_h1]))174 extra["max_abs_revision_h1_targets"] = float(revs.max())175176 res["feasibility"] = "complete"177 res["feasibility_reason"] = None178 res["data_issues"] = [179 {"issue": "Irregular RSAFS publication timing: vintages 2025-10-31 to 2026-03-31 end two months before the vintage date; the 2026-04-30 vintage adds 2026-02 and 2026-03 at once.",180 "handling": "Targets defined as last(V)+h from the last month actually present in each vintage, as TASK.md prescribes; target months are therefore unevenly spaced (h=1 skips 2026-03, h=3 skips 2026-05)."},181 {"issue": "Recent months are revised from vintage to vintage (e.g. 2026-05 level 763705 in the 2026-06-30 vintage, 766876 in 2026-07-31).",182 "handling": "Each forecast uses the values as published in its own vintage; actuals are first releases from first_release.csv."},183 {"issue": "Large downward level revision (about 3.5-4%) between the 2026-08-31 and 2026-09-30 vintages, visible in first_release.csv (2026-08 base-year level 700218 vs 731700 in shipped vintages) and in RSAFS_current.",184 "handling": "YoY growth is always computed within a single vintage so the level shift mostly cancels; not corrected. It affects the h=3 actuals of the last two origins (2026-08, 2026-09) and estimand (d)."},185 {"issue": "PAYEMS publication also delayed (2025-10-31 vintage ends 2025-08; 2025-12-31 adds two months).",186 "handling": "PAYEMS not used by the model; recorded only."},187 {"issue": "UMCSENT has 210 empty fields, all before 1978 (quarterly survey era); file is current vintage.",188 "handling": "UMCSENT not used by the model; recorded only."},189 {"issue": "RSAFS_current.csv and first_release.csv contain post-origin information.",190 "handling": "RSAFS_current used only for estimand (d); first_release only as evaluation actual."},191 ]192 files = sorted(glob.glob(os.path.join(DATA, "vintages", "*.csv"))) + [193 os.path.join(DATA, f) for f in ("RSAFS_current.csv", "UMCSENT.csv", "first_release.csv")]194 # PAYEMS and UMCSENT files are read only for the audit row counts (they do not enter forecasts)195 res["n_input_rows"] = {os.path.relpath(f, DATA): n_rows(f) for f in files}196 res["n_analysed"] = int(len(F[F.model == "model"]))197 res["random_seed"] = SEED198 res.update(extra)199 res["block_bootstrap"] = {"B": B_BOOT, "block_length": BLOCK, "type": "circular moving-block, paired errors"}200 res["last_month_by_origin"] = {V: str(lasts[V]) for V in ORIGINS}201 res["ar1_coef_by_origin"] = {V: {"const": float(arcoef[V][0]), "phi": float(arcoef[V][1])} for V in ORIGINS}202 res["abs_revision_by_h1_target"] = {str(T): float(r) for T, r in zip(t_h1, revs)}203204 with open(os.path.join(OUT, "results.json"), "w") as f:205 json.dump(res, f, indent=2)206207 # figure208 fig, axes = plt.subplots(1, 2, figsize=(12, 4.5), sharey=True)209 for ax, h in zip(axes, HORIZONS):210 sub = F[F.horizon == h]211 a = sub[sub.model == "naive"].sort_values("target_date")212 x = pd.to_datetime(a["target_date"])213 ax.plot(x, a["actual_first_release"], "k-o", label="first-release actual")214 for m, st in (("naive", "b--s"), ("model", "r-.^")):215 s = sub[sub.model == m].sort_values("target_date")216 ax.plot(pd.to_datetime(s["target_date"]), s["forecast"], st,217 label=f"{m} (RMSE {res[f'rmse_{m}_h{h}']:.2f})" if m == "naive" else218 f"BED model (RMSE {res[f'rmse_{m}_h{h}']:.2f})")219 ax.set_title(f"h = {h}")220 ax.set_ylabel("YoY growth, % (100 ln)")221 ax.grid(alpha=0.3)222 ax.legend(fontsize=8)223 ax.tick_params(axis="x", rotation=45)224 fig.suptitle("RSAFS YoY growth: real-time forecasts vs first release (10 origins, 2025-10-31 to 2026-07-31)")225 fig.tight_layout()226 fig.savefig(os.path.join(OUT, "figure1.png"), dpi=120)227 plt.close(fig)228229 for h in HORIZONS:230 print(f"h={h}: " + ", ".join(f"{m} RMSE={res[f'rmse_{m}_h{h}']:.4f}" for m in MODELS)231 + f"; DM={res[f'dm_stat_h{h}']:.3f} p={res[f'dm_pvalue_h{h}']:.3f} -> {res[f'conclusion_code_h{h}']}")232 print(f"mean abs revision (h=1 targets) = {res['mean_abs_revision_h1_targets']:.4f}")233234235if __name__ == "__main__":236 main()237