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.py3Regenerates out/results.json, out/forecasts.csv, out/table1.csv, out/figure1.png from data/ only."""4import json, os5import numpy as np6import pandas as pd7from scipy.stats import norm8import matplotlib9matplotlib.use("Agg")10import matplotlib.pyplot as plt1112D, O = "data", "out"13os.makedirs(O, exist_ok=True)14ORIGINS = ["2025-10-31", "2025-11-30", "2025-12-31", "2026-01-31", "2026-02-28",15 "2026-03-31", "2026-04-30", "2026-05-31", "2026-06-30", "2026-07-31"]16HORIZONS = [1, 3]17SEED = 2026103018n_rows = {}192021def read(rel):22 df = pd.read_csv(os.path.join(D, rel), parse_dates=["observation_date"]).set_index("observation_date")23 n_rows[rel] = len(df)24 s = df.iloc[:, 0].astype(float)25 assert s.index.is_monotonic_increasing and not s.index.duplicated().any()26 return s272829def yoy(level):30 return 100 * np.log(level / level.shift(12))313233fr = pd.read_csv(os.path.join(D, "first_release.csv"), parse_dates=["target_month", "first_release_vintage"])34n_rows["first_release.csv"] = len(fr)35actual = fr.set_index("target_month")["yoy_first_release"]3637data_issues = []38rows = []39last_months = {}40for V in ORIGINS:41 L = read(f"vintages/RSAFS_vintage_{V}.csv")42 P = read(f"vintages/PAYEMS_vintage_{V}.csv") # read for audit only; not used by the chosen model43 assert L.isna().sum() == 0 and (L.index.to_series().diff().dt.days.dropna() <= 31).all()44 g = yoy(L)45 last = L.index[-1]46 last_months[V] = last47 # AR(1) on g_m, 1993-02 <= m <= last48 y = g.loc["1993-02-01":last]49 x = g.shift(1).loc["1993-02-01":last]50 X = np.column_stack([np.ones(len(y)), x.values])51 a, b = np.linalg.lstsq(X, y.values, rcond=None)[0]52 lnL = np.log(L)53 mu = (lnL.loc[last] - lnL.loc[last - pd.DateOffset(months=12)]) / 12 # mean monthly log growth, last 12 months54 for h in HORIZONS:55 T = last + pd.DateOffset(months=h)56 f_naive = g.loc[last]57 f_snaive = g.loc[T - pd.DateOffset(months=12)]58 f = g.loc[last]59 for _ in range(h):60 f = a + b * f61 f_ar1 = f62 # Model: base-effect-aware drift. Forecast the level by extrapolating the last published level with the63 # trailing-12-month mean monthly log growth, then divide by the (already published) base-period level.64 f_model = 100 * (lnL.loc[last] + h * mu - lnL.loc[T - pd.DateOffset(months=12)])65 for m, fc in [("naive", f_naive), ("snaive", f_snaive), ("ar1", f_ar1), ("model", f_model)]:66 rows.append(dict(origin_date=V, target_date=T.strftime("%Y-%m-%d"), horizon=h, forecast=round(float(fc), 6),67 actual_first_release=float(actual.loc[T]), model=m, vintage_used=V))6869fc = pd.DataFrame(rows)70fc.to_csv(os.path.join(O, "forecasts.csv"), index=False)71fc["err"] = fc.actual_first_release - fc.forecast7273# publication-timing audit74for V in ORIGINS:75 lag = (pd.Timestamp(V).year - last_months[V].year) * 12 + pd.Timestamp(V).month - last_months[V].month76 if lag != 1:77 data_issues.append({"issue": f"Irregular publication timing: RSAFS vintage {V} ends at {last_months[V]:%Y-%m} ({lag} months before the origin month).",78 "handling": "Per task rules the target is defined relative to last(V), the last month present in the vintage, not to the calendar; forecasts use only that vintage."})79lm = pd.Series(last_months)80added = (lm.dt.year * 12 + lm.dt.month).diff()81for V, k in added.items():82 if k == 2:83 data_issues.append({"issue": f"Vintage {V} adds two new months at once (catch-up after delayed release); the month {lm[V] - pd.DateOffset(months=1):%Y-%m} is never an h=1 target and its first release is in the same vintage as the following month.",84 "handling": "No adjustment; targets follow last(V)+h. Horizon-1 targets therefore skip one month and h=3 targets skip one month, so the 10 origins are not equally spaced in target time."})8586table = []87for m in ["naive", "snaive", "ar1", "model"]:88 for h in HORIZONS:89 e = fc[(fc.model == m) & (fc.horizon == h)].err.values90 table.append(dict(model=m, horizon=h, rmse=float(np.sqrt(np.mean(e ** 2))), mae=float(np.mean(np.abs(e))), n=len(e)))91tab = pd.DataFrame(table)92tab.to_csv(os.path.join(O, "table1.csv"), index=False)939495def nw_lrv(d, L):96 d = d - d.mean()97 n = len(d)98 v = np.sum(d * d) / n99 for l in range(1, L + 1):100 v += 2 * (1 - l / (L + 1)) * np.sum(d[l:] * d[:-l]) / n101 return v102103104res = dict(task_id="T1-forecast",105 model_name="Base-effect drift (trailing-12m mean growth level extrapolation)",106 model_description=(107 "At origin V, using only the RSAFS vintage dated V: let last = last month in the vintage and "108 "mu = (ln L_last - ln L_{last-12})/12, the mean monthly log growth over the last 12 published months. "109 "The level of target month T = last+h is forecast as ln L_T = ln L_last + h*mu, and the YoY forecast is "110 "100*(ln L_last + h*mu - ln L_{T-12}), where L_{T-12} is the base-period level already published in the same vintage. "111 "Equivalently: naive forecast + h*g_last/12 - 100*(ln L_{T-12} - ln L_{last-12}), i.e. it corrects the naive forecast "112 "for known base effects (the months dropping out of the 12-month window). No parameters are estimated; "113 "PAYEMS and UMCSENT are not used (specification fixed a priori, before any evaluation against actuals). "114 "All inputs are from vintage V, so the availability rule holds trivially."),115 n_origins=len(ORIGINS), first_origin=ORIGINS[0], last_origin=ORIGINS[-1])116for r in table:117 res[f"rmse_{r['model']}_h{r['horizon']}"] = r["rmse"]118 res[f"mae_{r['model']}_h{r['horizon']}"] = r["mae"]119120rng = np.random.default_rng(SEED)121boot = {}122for h in HORIZONS:123 en = fc[(fc.model == "naive") & (fc.horizon == h)].sort_values("origin_date").err.values124 em = fc[(fc.model == "model") & (fc.horizon == h)].sort_values("origin_date").err.values125 res[f"rel_improvement_rmse_h{h}"] = 1 - res[f"rmse_model_h{h}"] / res[f"rmse_naive_h{h}"]126 d = en ** 2 - em ** 2127 n = len(d)128 lrv = nw_lrv(d, h - 1)129 dm = d.mean() / np.sqrt(lrv / n)130 p = 2 * (1 - norm.cdf(abs(dm)))131 res[f"dm_stat_h{h}"] = float(dm)132 res[f"dm_pvalue_h{h}"] = float(p)133 better = "model" if res[f"rmse_model_h{h}"] < res[f"rmse_naive_h{h}"] else "naive"134 res[f"conclusion_code_h{h}"] = f"{better}_better_{'significant' if p < 0.05 else 'not_significant'}"135 # supplementary: HLN-adjusted DM with t(n-1), and moving-block bootstrap (block 3) of RMSE_naive - RMSE_model136 k = np.sqrt((n + 1 - 2 * h + h * (h - 1) / n) / n)137 hln = dm * k138 B, bl = 10000, 3139 diffs = []140 for _ in range(B):141 starts = rng.integers(0, n - bl + 1, size=int(np.ceil(n / bl)))142 idx = np.concatenate([np.arange(s, s + bl) for s in starts])[:n]143 diffs.append(np.sqrt(np.mean(en[idx] ** 2)) - np.sqrt(np.mean(em[idx] ** 2)))144 diffs = np.array(diffs)145 from scipy.stats import t as tdist146 boot[h] = dict(hln_stat=float(hln), hln_p=float(2 * (1 - tdist.cdf(abs(hln), n - 1))),147 rmse_diff=float(res[f"rmse_naive_h{h}"] - res[f"rmse_model_h{h}"]),148 boot_ci95=[float(np.percentile(diffs, 2.5)), float(np.percentile(diffs, 97.5))],149 wins_model=int(np.sum(np.abs(em) < np.abs(en))), n=n)150151# (d) revision of h=1 targets152cur = read("RSAFS_current.csv")153um = read("UMCSENT.csv") # audited only (not used by the model)154n_um_missing = int(um.isna().sum())155gcur = yoy(cur)156t1 = pd.to_datetime(fc[(fc.horizon == 1) & (fc.model == "naive")].target_date)157revs = np.abs(gcur.loc[t1].values - actual.loc[t1].values)158res["mean_abs_revision_h1_targets"] = float(np.mean(revs))159lvl_rev = 100 * (cur.loc[fr.target_month].values / fr.level_first_release.values - 1)160data_issues.append({"issue": f"Large level revision: current RSAFS levels for 2024-06..2026-07 are {-lvl_rev[:-2].max():.1f}% to {-lvl_rev[:-2].min():.1f}% below their first-release levels (a benchmark/annual revision that first appears in the 2026-09-30 month-end vintage shifted the whole level; the 2026-08 and 2026-09 first releases are already on the new level); YoY growth is revised by up to {np.max(np.abs(gcur.loc[fr.target_month].values - fr.yoy_first_release.values)):.2f} pp.",161 "handling": "Forecasts are evaluated against first-release YoY (computed within one vintage, so level shifts cancel to first order). The current file is used only for estimand (d)."})162data_issues.append({"issue": "Real-time vintages differ from each other through routine revisions of the last 2-3 months and through annual revisions; values of a given month differ across origins.",163 "handling": "Each origin uses only its own vintage file; no splicing of vintages."})164data_issues.append({"issue": f"UMCSENT has empty values for {n_um_missing} months before 1978 (quarterly survey); PAYEMS vintages also show irregular timing (2025-10-31 ends at 2025-08; 2025-12-31 adds two months at once) due to delayed BLS releases.",165 "handling": "Neither predictor enters the chosen model; files were read, audited and row-counted only (PAYEMS, UMCSENT)."})166res["feasibility"] = "complete"167res["feasibility_reason"] = None168res["data_issues"] = data_issues169res["n_input_rows"] = n_rows170res["n_analysed"] = int(len(fc[fc.model == "model"]))171res["random_seed"] = SEED172res["supplementary"] = {"per_horizon": {f"h{h}": boot[h] for h in HORIZONS},173 "abs_revision_h1_targets": {t.strftime("%Y-%m"): float(r) for t, r in zip(t1, revs)},174 "ar1_note": "AR(1) re-estimated at every origin on that vintage only."}175with open(os.path.join(O, "results.json"), "w") as f:176 json.dump(res, f, indent=2)177178# figure179fig, axes = plt.subplots(1, 2, figsize=(12, 4.5), sharey=True)180for ax, h in zip(axes, HORIZONS):181 sub = fc[fc.horizon == h]182 a = actual.loc["2025-06-01":"2026-09-01"]183 ax.plot(a.index, a.values, "k-o", ms=4, label="first-release actual")184 for m, st in [("naive", "s--"), ("model", "^-"), ("ar1", "x:")]:185 s = sub[sub.model == m]186 ax.plot(pd.to_datetime(s.target_date), s.forecast, st, label=m)187 ax.set_title(f"h = {h}")188 ax.set_ylabel("YoY growth of RSAFS, %")189 ax.grid(alpha=.3)190 ax.tick_params(axis="x", rotation=45)191axes[0].legend()192fig.suptitle("Real-time forecasts of US retail sales YoY growth vs first release")193fig.tight_layout()194fig.savefig(os.path.join(O, "figure1.png"), dpi=130)195print(json.dumps({k: v for k, v in res.items() if k not in ("data_issues", "model_description")}, indent=1))196