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, out/figure1.png.4"""5import json, os6import numpy as np7import pandas as pd8from scipy.stats import norm, t as student_t9import matplotlib10matplotlib.use("Agg")11import matplotlib.pyplot as plt1213SEED = 2026103014DATA, OUT = "data", "out"15ORIGINS = ["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 df = pd.read_csv(os.path.join(DATA, rel))24 n_input_rows[rel] = int(len(df))25 s = df.set_index(pd.to_datetime(df.iloc[:, 0])).iloc[:, 1].astype(float)26 s.index = s.index.to_period("M")27 return s282930def yoy(level):31 return (100 * np.log(level / level.shift(12))).dropna()323334def fc_naive(L, g, last, h):35 return g[last]363738def fc_snaive(L, g, last, h):39 return g[last + h - 12]404142def fc_ar1(L, g, last, h):43 gg = g[(g.index >= pd.Period("1993-01", "M")) & (g.index <= last)]44 y = gg.values[1:] # m from 1993-02 to last(V)45 x = gg.values[:-1] # m-146 assert gg.index[1] == pd.Period("1993-02", "M")47 X = np.column_stack([np.ones_like(x), x])48 a, b = np.linalg.lstsq(X, y, rcond=None)[0]49 f = g[last]50 for _ in range(h):51 f = a + b * f52 return f535455def fc_model(L, g, last, h):56 """Base-effect drift model: future monthly log changes = mean of the last 12 monthly57 log changes in vintage V (= g_last/12); the base-period levels L_{T-12} are taken58 from vintage V (they are already published)."""59 drift = g[last] / 12.0 # percent per month60 lnL_T = 100 * np.log(L[last]) + h * drift61 return lnL_T - 100 * np.log(L[last + h - 12])626364FUNCS = {"naive": fc_naive, "snaive": fc_snaive, "ar1": fc_ar1, "model": fc_model}6566# ---------------- actuals ----------------67fr_df = pd.read_csv(os.path.join(DATA, "first_release.csv"))68n_input_rows["first_release.csv"] = int(len(fr_df))69fr = pd.Series(fr_df.yoy_first_release.values,70 index=pd.PeriodIndex(pd.to_datetime(fr_df.target_month), freq="M"))7172# ---------------- forecasts ----------------73rows, lasts = [], {}74for V in ORIGINS:75 L = read_series(f"vintages/RSAFS_vintage_{V}.csv") # only the vintage dated V76 read_series(f"vintages/PAYEMS_vintage_{V}.csv") # audited (row count) but not used by any model77 g = yoy(L)78 last = L.index.max()79 assert last.to_timestamp() < pd.Timestamp(V) # information rule sanity check80 lasts[V] = str(last)81 for h in HORIZONS:82 T = last + h83 for m in MODELS:84 rows.append(dict(origin_date=V, target_date=T.to_timestamp().strftime("%Y-%m-%d"),85 horizon=h, forecast=round(float(FUNCS[m](L, g, last, h)), 6),86 actual_first_release=float(fr[T]), model=m, vintage_used=V))87fc = pd.DataFrame(rows)88fc.to_csv(os.path.join(OUT, "forecasts.csv"), index=False)89fc["err"] = fc.actual_first_release - fc.forecast9091# ---------------- accuracy ----------------92res = {}93tab = []94for m in MODELS:95 for h in HORIZONS:96 e = fc[(fc.model == m) & (fc.horizon == h)].err.values97 rmse, mae = float(np.sqrt(np.mean(e ** 2))), float(np.mean(np.abs(e)))98 res[f"rmse_{m}_h{h}"], res[f"mae_{m}_h{h}"] = rmse, mae99 tab.append(dict(model=m, horizon=h, rmse=round(rmse, 6), mae=round(mae, 6), n=len(e)))100pd.DataFrame(tab).to_csv(os.path.join(OUT, "table1.csv"), index=False)101102103def nw_lrv(d, lag):104 d = np.asarray(d, float); n = len(d); dc = d - d.mean()105 v = np.sum(dc * dc) / n106 for k in range(1, lag + 1):107 v += 2 * (1 - k / (lag + 1)) * np.sum(dc[k:] * dc[:-k]) / n108 return v109110111rng = np.random.default_rng(SEED)112B = 10000113extra = {}114for h in HORIZONS:115 sub = fc[fc.horizon == h]116 en = sub[sub.model == "naive"].sort_values("origin_date").err.values117 em = sub[sub.model == "model"].sort_values("origin_date").err.values118 d = en ** 2 - em ** 2119 n = len(d)120 lrv = nw_lrv(d, h - 1)121 if lrv > 0:122 dm = float(d.mean() / np.sqrt(lrv / n)); p = float(2 * (1 - norm.cdf(abs(dm))))123 else:124 dm, p = None, None125 res[f"rel_improvement_rmse_h{h}"] = 1 - res[f"rmse_model_h{h}"] / res[f"rmse_naive_h{h}"]126 res[f"dm_stat_h{h}"], res[f"dm_pvalue_h{h}"] = dm, p127 better = "model" if res[f"rmse_model_h{h}"] < res[f"rmse_naive_h{h}"] else "naive"128 sig = "significant" if (p is not None and p < 0.05) else "not_significant"129 res[f"conclusion_code_h{h}"] = f"{better}_better_{sig}"130 extra[f"dm_lrv_h{h}"] = float(lrv)131 if dm is not None: # sensitivity only (not the required statistic): HLN small-sample correction, t(n-1)132 hln = np.sqrt((n + 1 - 2 * h + h * (h - 1) / n) / n)133 extra[f"sens_dm_hln_stat_h{h}"] = float(dm * hln)134 extra[f"sens_dm_hln_pvalue_t_h{h}"] = float(2 * student_t.sf(abs(dm * hln), n - 1))135 extra[f"dm_mean_loss_diff_h{h}"] = float(d.mean())136 # circular block bootstrap of RMSE_model - RMSE_naive (block length h; iid at h=1)137 bl = h138 nb = int(np.ceil(n / bl))139 diffs = np.empty(B)140 for b in range(B):141 starts = rng.integers(0, n, nb)142 idx = np.concatenate([(s + np.arange(bl)) % n for s in starts])[:n]143 diffs[b] = np.sqrt(np.mean(em[idx] ** 2)) - np.sqrt(np.mean(en[idx] ** 2))144 extra[f"boot_rmse_diff_model_minus_naive_h{h}"] = float(np.sqrt(np.mean(em ** 2)) - np.sqrt(np.mean(en ** 2)))145 extra[f"boot_ci95_lo_h{h}"] = float(np.percentile(diffs, 2.5))146 extra[f"boot_ci95_hi_h{h}"] = float(np.percentile(diffs, 97.5))147 extra[f"boot_share_model_better_h{h}"] = float(np.mean(diffs < 0))148149# ---------------- revision analysis (d) ----------------150cur = read_series("RSAFS_current.csv") # used ONLY here151gcur = yoy(cur)152h1 = fc[(fc.horizon == 1) & (fc.model == "naive")]153T1 = pd.PeriodIndex(pd.to_datetime(h1.target_date), freq="M")154rev = gcur[T1].values - fr[T1].values155res["mean_abs_revision_h1_targets"] = float(np.mean(np.abs(rev)))156extra["mean_revision_h1_targets_current_minus_first"] = float(np.mean(rev))157extra["revisions_h1_targets"] = {str(t): round(float(r), 6) for t, r in zip(T1, rev)}158159ucs = pd.read_csv(os.path.join(DATA, "UMCSENT.csv")) # read only for the row count; not used160n_input_rows["UMCSENT.csv"] = int(len(ucs))161162# ---------------- figure ----------------163fig, axes = plt.subplots(1, 2, figsize=(12, 4.5), sharey=True)164for ax, h in zip(axes, HORIZONS):165 sub = fc[fc.horizon == h].copy(); sub["t"] = pd.to_datetime(sub.target_date)166 a = sub[sub.model == "naive"].sort_values("t")167 ax.plot(a.t, a.actual_first_release, "k-o", label="first-release actual")168 for m, sty, lab in [("naive", "C0--s", "naive"), ("model", "C3-^", "base-effect drift model"),169 ("ar1", "C2:", "AR(1)")]:170 s = sub[sub.model == m].sort_values("t")171 ax.plot(s.t, s.forecast, sty, label=lab)172 ax.set_title(f"h = {h}: target month = last(V) + {h}")173 ax.set_xlabel("target month"); ax.grid(alpha=.3); ax.tick_params(axis="x", rotation=45)174axes[0].set_ylabel("YoY growth, % (100·log ratio)")175axes[0].legend(fontsize=8)176fig.suptitle("US retail & food services sales: real-time forecasts vs first release (10 origins)")177fig.tight_layout(); fig.savefig(os.path.join(OUT, "figure1.png"), dpi=130); plt.close(fig)178179# ---------------- results.json ----------------180issues = [181 {"issue": "Irregular RSAFS publication timing: vintages 2025-10-31 and 2026-03-31 end two months before the origin; vintage 2026-04-30 adds two months (2026-02, 2026-03) at once.",182 "handling": "Targets defined as last(V)+h from the vintage itself, per task rule; 2026-03 is never an h=1 target, h=1 target of 2025-10-31 is 2025-09."},183 {"issue": "PAYEMS vintage 2025-12-31 adds 2025-10 and 2025-11 at once (delayed release) and vintage 2026-02-28 contains a benchmark level shift (~-1 million).",184 "handling": "Documented; PAYEMS is not used by any forecast."},185 {"issue": "UMCSENT has 210 empty fields (pre-1978 non-survey months).",186 "handling": "Documented; UMCSENT is not used by any forecast."},187 {"issue": "Large annual benchmark revision of RSAFS between the last origin vintage and the current vintage (138 months changed back to 2015, up to 4.6% in level).",188 "handling": "Affects only estimand (d), which by design uses the current vintage; forecasts use vintage V only and are evaluated against first release."},189 {"issue": "first_release.csv dates first release at month-end granularity; 2026-02 and 2026-03 share the first-release vintage 2026-04-30.",190 "handling": "Used as given; yoy recomputed from levels (max abs diff < 5e-7) and levels verified against the shipped vintages."},191]192out = {193 "task_id": "T1-forecast",194 "model_name": "base-effect drift (BED)",195 "model_description": ("Univariate real-time model using only the RSAFS vintage dated V. The h-step-ahead log level is "196 "100*ln L_{last(V)} + h*mu, where the monthly drift mu = mean of the last 12 monthly log changes in vintage V "197 "(= g_{last(V)}/12); the forecast of YoY growth for T = last(V)+h is that log level minus 100*ln L_{T-12}, with L_{T-12} "198 "taken from vintage V (already published). It therefore keeps the recent trend but accounts for the known base "199 "effects (the monthly changes that drop out of the 12-month window). No parameters are estimated (estimation window: "200 "the 13 most recent levels of vintage V); PAYEMS and UMCSENT are not used. Specified before looking at evaluation errors; no tuning on the 10 origins."),201 "n_origins": len(ORIGINS), "first_origin": ORIGINS[0], "last_origin": ORIGINS[-1],202}203for m in MODELS:204 for h in HORIZONS:205 out[f"rmse_{m}_h{h}"] = res[f"rmse_{m}_h{h}"]; out[f"mae_{m}_h{h}"] = res[f"mae_{m}_h{h}"]206for k in ["rel_improvement_rmse_h1", "rel_improvement_rmse_h3", "dm_stat_h1", "dm_pvalue_h1",207 "dm_stat_h3", "dm_pvalue_h3", "conclusion_code_h1", "conclusion_code_h3", "mean_abs_revision_h1_targets"]:208 out[k] = res[k]209out.update({210 "feasibility": "complete", "feasibility_reason": None, "data_issues": issues,211 "n_input_rows": dict(sorted(n_input_rows.items())),212 "n_analysed": int(len(fc[fc.model == "model"])), "random_seed": SEED,213 "last_month_by_origin": lasts,214})215out.update(extra)216217218def rnd(o):219 if isinstance(o, float): return round(o, 6)220 if isinstance(o, dict): return {k: rnd(v) for k, v in o.items()}221 if isinstance(o, list): return [rnd(v) for v in o]222 return o223224225out = rnd(out)226with open(os.path.join(OUT, "results.json"), "w") as f:227 json.dump(out, f, indent=1)228print(pd.DataFrame(tab).to_string(index=False))229for h in HORIZONS:230 print(f"h={h}: rel_impr={out[f'rel_improvement_rmse_h{h}']}, DM={out[f'dm_stat_h{h}']}, p={out[f'dm_pvalue_h{h}']}, {out[f'conclusion_code_h{h}']}")231print("mean_abs_revision_h1_targets =", out["mean_abs_revision_h1_targets"])232