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.py3"""4import json, os5import numpy as np, pandas as pd6from scipy import stats7import matplotlib8matplotlib.use("Agg")9import matplotlib.pyplot as plt1011DATA, OUT = "data", "out"12ORIGINS = ["2025-10-31", "2025-11-30", "2025-12-31", "2026-01-31", "2026-02-28",13 "2026-03-31", "2026-04-30", "2026-05-31", "2026-06-30", "2026-07-31"]14HORIZONS = [1, 3]15DRIFT_WINDOW = 36 # months of m/m log changes used for the drift (fixed a priori)16SEED = 2026103017n_rows = {}1819def read(path):20 d = pd.read_csv(os.path.join(DATA, path), parse_dates=["observation_date"])21 n_rows[path] = int(len(d))22 s = d.set_index("observation_date").iloc[:, 0].astype(float)23 s.index.freq = None24 return s2526def yoy(level):27 return 100 * np.log(level / level.shift(12))2829fr = pd.read_csv(os.path.join(DATA, "first_release.csv"), parse_dates=["target_month"])30n_rows["first_release.csv"] = int(len(fr))31actual = fr.set_index("target_month")["yoy_first_release"]3233def add_months(ts, k):34 return ts + pd.DateOffset(months=k)3536rows, audit = [], []37for V in ORIGINS:38 L = read(f"vintages/RSAFS_vintage_{V}.csv")39 P = read(f"vintages/PAYEMS_vintage_{V}.csv") # read for audit; not used by the model40 assert L.isna().sum() == 041 g = yoy(L)42 last = L.index[-1]43 vdate = pd.Timestamp(V)44 lag = (vdate.year - last.year) * 12 + vdate.month - last.month45 audit.append({"origin": V, "last_month": str(last.date()), "months_behind": int(lag),46 "payems_last": str(P.index[-1].date())})47 # AR(1) on g, pairs with 1993-02 <= m <= last48 gg = g.dropna()49 y = gg.loc["1993-02-01":last]50 x = gg.shift(1).loc[y.index]51 X = np.column_stack([np.ones(len(y)), x.values])52 beta = np.linalg.lstsq(X, y.values, rcond=None)[0]53 # model: drift-in-levels with exact base effects54 dl = 100 * np.log(L).diff().dropna()55 mu = dl.iloc[-DRIFT_WINDOW:].mean()56 for h in HORIZONS:57 T = add_months(last, h)58 f_naive = g.loc[last]59 f_snaive = g.loc[add_months(T, -12)]60 f = g.loc[last]61 for _ in range(h):62 f = beta[0] + beta[1] * f63 f_ar1 = f64 # g_T = 100 ln(L_last/L_{T-12}) + h*mu (L_{T-12} is in the vintage)65 f_model = 100 * np.log(L.loc[last] / L.loc[add_months(T, -12)]) + h * mu66 for name, val in [("naive", f_naive), ("snaive", f_snaive), ("ar1", f_ar1), ("model", f_model)]:67 rows.append({"origin_date": V, "target_date": T.strftime("%Y-%m-%d"), "horizon": h,68 "forecast": float(val), "actual_first_release": float(actual.loc[T]),69 "model": name, "vintage_used": V})7071fc = pd.DataFrame(rows)72fc["error"] = fc["actual_first_release"] - fc["forecast"]73fc.drop(columns="error").to_csv(os.path.join(OUT, "forecasts.csv"), index=False, float_format="%.6f")7475tab = []76for m in ["naive", "snaive", "ar1", "model"]:77 for h in HORIZONS:78 e = fc[(fc.model == m) & (fc.horizon == h)]["error"].values79 tab.append({"model": m, "horizon": h, "rmse": float(np.sqrt(np.mean(e**2))),80 "mae": float(np.mean(np.abs(e))), "n": int(len(e))})81tab = pd.DataFrame(tab)82tab.to_csv(os.path.join(OUT, "table1.csv"), index=False, float_format="%.6f")8384def dm(e_naive, e_model, h):85 d = e_naive**2 - e_model**286 n = len(d); dbar = d.mean(); dc = d - dbar87 lrv = np.sum(dc**2) / n88 for j in range(1, h):89 w = 1 - j / h90 lrv += 2 * w * np.sum(dc[j:] * dc[:-j]) / n91 stat = dbar / np.sqrt(lrv / n)92 return float(stat), float(2 * (1 - stats.norm.cdf(abs(stat))))9394res = {"task_id": "T1-forecast",95 "model_name": "Drift-in-levels with exact base effects (RW-drift36)",96 "model_description": (97 "The YoY growth of target month T = last(V)+h is written exactly as "98 "g_T = 100*ln(L_last/L_{T-12}) + sum_{k=1..h} 100*dln L_{last+k}. The first term uses only "99 "levels published in the RSAFS vintage dated V (L_{T-12} is always in the vintage), so the "100 "'base effect' of months dropping out of the 12-month window is known exactly at the origin. "101 "The unknown future monthly log changes are each forecast by their mean over the last 36 "102 "monthly log changes available in the same vintage (a random walk with drift in log levels; "103 "window fixed a priori, not tuned on evaluation targets). Re-estimated at every origin from the "104 "vintage dated V only; PAYEMS and UMCSENT are not used (no information leakage possible); "105 "RSAFS_current is not used."),106 "n_origins": len(ORIGINS), "first_origin": ORIGINS[0], "last_origin": ORIGINS[-1]}107for _, r in tab.iterrows():108 res[f"rmse_{r.model}_h{r.horizon}"] = r.rmse109 res[f"mae_{r.model}_h{r.horizon}"] = r.mae110111rng = np.random.default_rng(SEED)112boot = {}113for h in HORIZONS:114 en = fc[(fc.model == "naive") & (fc.horizon == h)].sort_values("origin_date")["error"].values115 em = fc[(fc.model == "model") & (fc.horizon == h)].sort_values("origin_date")["error"].values116 res[f"rel_improvement_rmse_h{h}"] = float(1 - res[f"rmse_model_h{h}"] / res[f"rmse_naive_h{h}"])117 s, p = dm(en, em, h)118 res[f"dm_stat_h{h}"], res[f"dm_pvalue_h{h}"] = s, p119 better = res[f"rmse_model_h{h}"] < res[f"rmse_naive_h{h}"]120 sig = p < 0.05121 res[f"conclusion_code_h{h}"] = ("model_better_" if better else "naive_better_") + ("significant" if sig else "not_significant")122 # circular moving-block bootstrap (block length max(h,2)) of RMSE_naive - RMSE_model123 n = len(en); b = max(h, 2); B = 9999; diffs = np.empty(B)124 for i in range(B):125 starts = rng.integers(0, n, int(np.ceil(n / b)))126 idx = np.concatenate([(np.arange(s0, s0 + b) % n) for s0 in starts])[:n]127 diffs[i] = np.sqrt(np.mean(en[idx]**2)) - np.sqrt(np.mean(em[idx]**2))128 boot[h] = {"rmse_diff_naive_minus_model": float(np.sqrt(np.mean(en**2)) - np.sqrt(np.mean(em**2))),129 "ci95": [float(np.percentile(diffs, 2.5)), float(np.percentile(diffs, 97.5))],130 "share_positive": float(np.mean(diffs > 0))}131132# (d) revisions133cur = read("RSAFS_current.csv")134gcur = yoy(cur)135t1 = fc[(fc.model == "naive") & (fc.horizon == 1)]["target_date"].map(pd.Timestamp)136revs = [(str(T.date()), float(gcur.loc[T]), float(actual.loc[T]), float(gcur.loc[T] - actual.loc[T])) for T in t1]137res["mean_abs_revision_h1_targets"] = float(np.mean([abs(r[3]) for r in revs]))138139um = read("UMCSENT.csv") # read for audit only; not used by the model140res["feasibility"] = "complete"141res["feasibility_reason"] = None142lag_txt = "; ".join(f"{a['origin']}->last {a['last_month'][:7]}" for a in audit if a['months_behind'] > 1)143res["data_issues"] = [144 {"issue": "Irregular RSAFS publication timing (normal: the month-end vintage ends one calendar month before V). Delayed origins (vintage ends two months before V): " + lag_txt + ". The 2026-04-30 vintage added two months (2026-02 and 2026-03) at once; from then on timing is normal.",145 "handling": "Targets defined relative to last(V) as required (T = last(V)+h), never relative to the calendar; for delayed origins the h=1 target is a month that would normally already be published. Consequently target month 2026-03 is never an h=1 target and 2026-05 is never an h=3 target."},146 {"issue": "Two target months (2026-02, 2026-03) were first released in the same month-end vintage (2026-04-30); 2025-09 first appeared only in the 2025-11-30 vintage (publication disruption).",147 "handling": "Actuals taken as given from first_release.csv (first month-end appearance); no adjustment."},148 {"issue": "Annual/benchmark revision of RSAFS between the 2026-08-31 and 2026-09-30 vintages lowered the level of the whole series by roughly 4% (e.g. 2025-08: 731700 in the 2026-07-31 vintage vs 700407 in the current file); first_release levels for 2026-08 onward are on the new basis.",149 "handling": "All comparisons use log YoY growth computed within a single vintage, which is largely invariant to proportional level rebasing; documented as affecting estimand (d)."},150 {"issue": "PAYEMS vintage 2026-02-28 contains the annual benchmark revision (level about 1 million lower than in the 2026-01-31 vintage).",151 "handling": "PAYEMS not used by the model; noted only."},152 {"issue": "UMCSENT has 210 empty values, all before 1978 (quarterly/irregular survey period).",153 "handling": "UMCSENT not used by the model; no imputation needed."},154 {"issue": "RSAFS 2026-07 first release (763602, old basis) vs current value 729915 (new basis): level comparisons across vintages are not meaningful.",155 "handling": "Revision measure (d) computed on growth rates only."},156]157res["n_input_rows"] = n_rows158res["n_analysed"] = int(len(fc[fc.model == "model"]))159res["random_seed"] = SEED160with open(os.path.join(OUT, "results.json"), "w") as f:161 json.dump(res, f, indent=2)162sens = {}163for W in [12, 24, 36, 60, 120]:164 for h in HORIZONS:165 errs = []166 for V in ORIGINS:167 L = pd.read_csv(os.path.join(DATA, f"vintages/RSAFS_vintage_{V}.csv"), parse_dates=[0]).set_index("observation_date").iloc[:, 0]168 last = L.index[-1]; T = add_months(last, h)169 mu = (100 * np.log(L).diff().dropna()).iloc[-W:].mean()170 errs.append(actual.loc[T] - (100 * np.log(L.loc[last] / L.loc[add_months(T, -12)]) + h * mu))171 sens[f"W{W}_h{h}_rmse"] = float(np.sqrt(np.mean(np.square(errs))))172with open(os.path.join(OUT, "supplement.json"), "w") as f:173 json.dump({"publication_audit": audit, "bootstrap": boot, "drift_window_sensitivity_posthoc": sens, "revisions_h1": revs,174 "ar1_note": "see forecasts.csv"}, f, indent=2)175176fig, axes = plt.subplots(2, 1, figsize=(8, 7), sharex=True)177for ax, h in zip(axes, HORIZONS):178 sub = fc[fc.horizon == h]179 a = sub[sub.model == "naive"].sort_values("target_date")180 t = pd.to_datetime(a.target_date)181 ax.plot(t, a.actual_first_release, "k-o", label="first-release actual")182 for m, st in [("naive", "C0--s"), ("ar1", "C2:^"), ("model", "C3-d")]:183 s2 = sub[sub.model == m].sort_values("target_date")184 ax.plot(pd.to_datetime(s2.target_date), s2.forecast, st, label=m)185 ax.set_title(f"h = {h}"); ax.set_ylabel("YoY growth, %"); ax.grid(alpha=.3)186axes[0].legend(fontsize=8)187fig.suptitle("US retail & food services sales YoY growth: real-time forecasts vs first release")188fig.tight_layout(); fig.savefig(os.path.join(OUT, "figure1.png"), dpi=130)189print(json.dumps({k: v for k, v in res.items() if k not in ("data_issues", "model_description")}, indent=1))190