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,4out/exclusions.csv. Deterministic (no randomness; random_seed recorded as null)."""5import json, os, glob6import numpy as np7import pandas as pd8from scipy import stats9import matplotlib10matplotlib.use("Agg")11import matplotlib.pyplot as plt1213D = "data"; O = "out"14os.makedirs(O, exist_ok=True)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]18DRIFT_WINDOW = 12 # months of log-level growth used for the model's drift (pre-specified)192021def read_series(path):22 d = pd.read_csv(path)23 d["observation_date"] = pd.to_datetime(d["observation_date"])24 s = d.set_index("observation_date").iloc[:, 0].astype(float)25 return s262728def yoy(level):29 # 100*ln(L_m/L_{m-12}); requires a gap-free monthly index (checked below)30 lv = np.log(level)31 return 100.0 * (lv - lv.shift(12))323334n_input_rows = {}35def count_rows(path):36 with open(path) as f:37 n = sum(1 for _ in f) - 138 n_input_rows[os.path.relpath(path, D)] = n3940# ---- read inputs -----------------------------------------------------------41fr = pd.read_csv(os.path.join(D, "first_release.csv"), parse_dates=["target_month"])42count_rows(os.path.join(D, "first_release.csv"))43actual = fr.set_index("target_month")["yoy_first_release"]44cur = read_series(os.path.join(D, "RSAFS_current.csv")); count_rows(os.path.join(D, "RSAFS_current.csv"))45umc = read_series(os.path.join(D, "UMCSENT.csv")); count_rows(os.path.join(D, "UMCSENT.csv"))46vint_r, vint_p = {}, {}47for V in ORIGINS:48 pr = os.path.join(D, "vintages", f"RSAFS_vintage_{V}.csv")49 pp = os.path.join(D, "vintages", f"PAYEMS_vintage_{V}.csv")50 vint_r[V] = read_series(pr); count_rows(pr)51 vint_p[V] = read_series(pp); count_rows(pp)5253for V, s in vint_r.items():54 idx = s.index55 assert (idx == pd.date_range(idx[0], idx[-1], freq="MS")).all(), f"gap in RSAFS {V}"56 assert s.notna().all()5758# ---- forecasts ---------------------------------------------------------------59rows = []60delays = []61for V in ORIGINS:62 L = vint_r[V]63 g = yoy(L)64 last = L.index[-1]65 Vd = pd.Timestamp(V)66 lag_months = (Vd.year - last.year) * 12 + Vd.month - last.month67 delays.append({"origin": V, "last_month": last.strftime("%Y-%m"), "months_behind_origin": lag_months})68 # AR(1) on g over 1993-02 <= m <= last69 gs = g.loc["1993-01-01":last]70 y = gs.loc["1993-02-01":].values71 x = gs.shift(1).loc["1993-02-01":].values72 X = np.column_stack([np.ones_like(x), x])73 a, b = np.linalg.lstsq(X, y, rcond=None)[0]74 lnL = np.log(L)75 drift = (lnL.iloc[-1] - lnL.iloc[-1 - DRIFT_WINDOW]) / DRIFT_WINDOW76 for h in HORIZONS:77 T = last + pd.DateOffset(months=h)78 if T not in actual.index:79 raise RuntimeError(f"no first-release actual for {T}")80 f_naive = g.loc[last]81 f_snaive = g.loc[T - pd.DateOffset(months=12)]82 f_ar = g.loc[last]83 for _ in range(h):84 f_ar = a + b * f_ar85 # model: random walk with drift in log level, drift = mean monthly log growth over86 # last 12 months of the vintage; YoY forecast uses the vintage's level at T-1287 lnT = lnL.iloc[-1] + h * drift88 f_model = 100.0 * (lnT - lnL.loc[T - pd.DateOffset(months=12)])89 for name, f in [("naive", f_naive), ("snaive", f_snaive), ("ar1", f_ar), ("model", f_model)]:90 rows.append({"origin_date": V, "target_date": T.strftime("%Y-%m-%d"), "horizon": h,91 "forecast": float(f), "actual_first_release": float(actual.loc[T]),92 "model": name, "vintage_used": V})93fc = pd.DataFrame(rows)94fc.to_csv(os.path.join(O, "forecasts.csv"), index=False, float_format="%.6f")9596# ---- evaluation ---------------------------------------------------------------97fc["err"] = fc["actual_first_release"] - fc["forecast"]98res = {}99tab = []100for m in ["naive", "snaive", "ar1", "model"]:101 for h in HORIZONS:102 e = fc[(fc.model == m) & (fc.horizon == h)]["err"].values103 rmse = float(np.sqrt(np.mean(e ** 2))); mae = float(np.mean(np.abs(e)))104 res[f"rmse_{m}_h{h}"] = rmse; res[f"mae_{m}_h{h}"] = mae105 tab.append({"model": m, "horizon": h, "rmse": rmse, "mae": mae, "n": len(e)})106pd.DataFrame(tab).to_csv(os.path.join(O, "table1.csv"), index=False, float_format="%.6f")107108109def dm_test(d, h):110 d = np.asarray(d, float); n = len(d); dbar = d.mean(); dc = d - dbar111 lrv = np.sum(dc ** 2) / n112 for k in range(1, h):113 w = 1 - k / h114 lrv += 2 * w * np.sum(dc[k:] * dc[:-k]) / n115 stat = dbar / np.sqrt(lrv / n)116 p = 2 * (1 - stats.norm.cdf(abs(stat)))117 return float(stat), float(p), float(lrv)118119extra = {}120for h in HORIZONS:121 en = fc[(fc.model == "naive") & (fc.horizon == h)].sort_values("origin_date")["err"].values122 em = fc[(fc.model == "model") & (fc.horizon == h)].sort_values("origin_date")["err"].values123 d = en ** 2 - em ** 2124 s, p, lrv = dm_test(d, h)125 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}"] = s; res[f"dm_pvalue_h{h}"] = p127 extra[f"dm_lrv_h{h}"] = lrv128 extra[f"n_model_better_sq_h{h}"] = int(np.sum(d > 0))129 better = res[f"rmse_model_h{h}"] < res[f"rmse_naive_h{h}"]130 sig = p < 0.05131 res[f"conclusion_code_h{h}"] = ("model_better_" if better else "naive_better_") + ("significant" if sig else "not_significant")132133# ---- (d) revisions -------------------------------------------------------------134gcur = yoy(cur)135t1 = fc[(fc.model == "naive") & (fc.horizon == 1)]136revs = [abs(gcur.loc[pd.Timestamp(t)] - a_) for t, a_ in zip(t1.target_date, t1.actual_first_release)]137res["mean_abs_revision_h1_targets"] = float(np.mean(revs))138extra["revision_h1_targets"] = {t: float(gcur.loc[pd.Timestamp(t)] - a_) for t, a_ in zip(t1.target_date, t1.actual_first_release)}139140# ---- exclusion ledger ------------------------------------------------------------141n_pairs = len(ORIGINS) * len(HORIZONS)142n_eval = int(fc[fc.model == "model"]["actual_first_release"].notna().sum())143pd.DataFrame([{"step": 1, "rule": "keep (origin,horizon) pairs whose target month has a first-release actual",144 "n_before": n_pairs, "n_after": n_eval, "n_removed": n_pairs - n_eval,145 "rationale": "all 20 targets (2025-09..2026-09) are in first_release.csv; nothing removed"}]146 ).to_csv(os.path.join(O, "exclusions.csv"), index=False)147148# ---- figure ------------------------------------------------------------------------149fig, axes = plt.subplots(1, 2, figsize=(11, 4), sharey=True)150for ax, h in zip(axes, HORIZONS):151 sub = fc[fc.horizon == h]152 t = pd.to_datetime(sub[sub.model == "naive"].target_date)153 ax.plot(t, sub[sub.model == "naive"].actual_first_release, "k-o", label="first release (actual)")154 ax.plot(t, sub[sub.model == "naive"].forecast, "b--s", label="naive")155 ax.plot(t, sub[sub.model == "model"].forecast, "r--^", label="model (RW+drift, level)")156 ax.set_title(f"h = {h}"); ax.set_xlabel("target month"); ax.grid(alpha=.3)157 ax.tick_params(axis="x", rotation=45)158axes[0].set_ylabel("YoY growth, % (100·log diff)"); axes[0].legend(fontsize=8)159fig.suptitle("US retail & food services sales: real-time forecasts vs first release")160fig.tight_layout(); fig.savefig(os.path.join(O, "figure1.png"), dpi=120); plt.close(fig)161162# ---- results.json ------------------------------------------------------------------163data_issues = [164 {"issue": "Irregular RSAFS publication timing: in the six vintages 2025-10-31..2026-03-31 the last published month is two months before the origin month (e.g. 2025-10-31 ends 2025-08), versus one month in 2026-04-30..2026-07-31; the 2026-04-30 vintage adds two months (2026-02 and 2026-03) at once",165 "handling": ("Targets defined as last(V)+h per TASK.md: h=1 targets are 2025-09..2026-02 and 2026-04..2026-07 (2026-03 is never an h=1 target); "166 "h=3 targets 2025-11..2026-04 and 2026-06..2026-09 (2026-05 skipped). All baselines and the model use the same last(V).")},167 {"issue": "First releases of 2025-09 (vintage 2025-11-30) and 2026-02/2026-03 (both first in 2026-04-30) were delayed relative to the normal calendar (government shutdown / disrupted schedule).",168 "handling": "Used as given in first_release.csv (first month-end vintage containing the month)."},169 {"issue": "Large level revision between first release and current vintage (e.g. 2026-07 first released at 763602, current 729915; levels of 2025 also lowered), indicating a benchmark/annual revision around 2026-08/09.",170 "handling": "YoY growth is computed within a single vintage for both numerator and denominator, so level shifts largely cancel; reported via estimand (d). RSAFS_current used only for (d)."},171 {"issue": "UMCSENT has 210 empty fields (pre-1978 months without survey).",172 "handling": "UMCSENT not used by the model; read only for row counts/audit."},173 {"issue": "PAYEMS vintages also show irregular timing (2025-12-31 vintage adds two months; ends 2025-11).",174 "handling": "PAYEMS not used by the model; read only for audit."},175]176out = {177 "task_id": "T1-forecast",178 "model_name": "rw_drift_loglevel_12m",179 "model_description": ("Random walk with drift on the log level of RSAFS from vintage V only: ln L_{last+h} = ln L_last + h*mu, "180 "mu = mean monthly log growth over the last 12 months in vintage V (= g_last/12). The YoY forecast is "181 "100*(ln L_{last+h} - ln L_{T-12}) with L_{T-12} taken from vintage V, so known base effects (the year-ago months) "182 "enter the forecast, unlike the naive. No estimated parameters beyond mu; no PAYEMS/UMCSENT used. "183 "Specification fixed a priori (not tuned on evaluation actuals); every input is in the vintage dated V."),184 "n_origins": len(ORIGINS), "first_origin": ORIGINS[0], "last_origin": ORIGINS[-1],185}186out.update(res)187out.update({188 "feasibility": "complete", "feasibility_reason": None,189 "data_issues": data_issues, "n_input_rows": n_input_rows,190 "n_analysed": n_eval, "random_seed": None,191 "publication_delays": delays, **extra,192})193with open(os.path.join(O, "results.json"), "w") as f:194 json.dump(out, f, indent=2)195for k in ["rmse_naive_h1", "rmse_model_h1", "rmse_naive_h3", "rmse_model_h3", "dm_stat_h1", "dm_pvalue_h1", "dm_stat_h3", "dm_pvalue_h3"]:196 print(k, round(out[k], 4))197