"""T1-forecast: real-time forecasting of US retail-sales YoY growth from ALFRED vintages. Run from workspace root: ~/wp29/venv/bin/python out/analysis.py Reads only data/, writes out/results.json, out/forecasts.csv, out/table1.csv, out/figure1.png, out/exclusions.csv. Deterministic (no randomness; random_seed recorded as null).""" import json, os, glob import numpy as np import pandas as pd from scipy import stats import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt D = "data"; O = "out" os.makedirs(O, exist_ok=True) ORIGINS = ["2025-10-31", "2025-11-30", "2025-12-31", "2026-01-31", "2026-02-28", "2026-03-31", "2026-04-30", "2026-05-31", "2026-06-30", "2026-07-31"] HORIZONS = [1, 3] DRIFT_WINDOW = 12 # months of log-level growth used for the model's drift (pre-specified) def read_series(path): d = pd.read_csv(path) d["observation_date"] = pd.to_datetime(d["observation_date"]) s = d.set_index("observation_date").iloc[:, 0].astype(float) return s def yoy(level): # 100*ln(L_m/L_{m-12}); requires a gap-free monthly index (checked below) lv = np.log(level) return 100.0 * (lv - lv.shift(12)) n_input_rows = {} def count_rows(path): with open(path) as f: n = sum(1 for _ in f) - 1 n_input_rows[os.path.relpath(path, D)] = n # ---- read inputs ----------------------------------------------------------- fr = pd.read_csv(os.path.join(D, "first_release.csv"), parse_dates=["target_month"]) count_rows(os.path.join(D, "first_release.csv")) actual = fr.set_index("target_month")["yoy_first_release"] cur = read_series(os.path.join(D, "RSAFS_current.csv")); count_rows(os.path.join(D, "RSAFS_current.csv")) umc = read_series(os.path.join(D, "UMCSENT.csv")); count_rows(os.path.join(D, "UMCSENT.csv")) vint_r, vint_p = {}, {} for V in ORIGINS: pr = os.path.join(D, "vintages", f"RSAFS_vintage_{V}.csv") pp = os.path.join(D, "vintages", f"PAYEMS_vintage_{V}.csv") vint_r[V] = read_series(pr); count_rows(pr) vint_p[V] = read_series(pp); count_rows(pp) for V, s in vint_r.items(): idx = s.index assert (idx == pd.date_range(idx[0], idx[-1], freq="MS")).all(), f"gap in RSAFS {V}" assert s.notna().all() # ---- forecasts --------------------------------------------------------------- rows = [] delays = [] for V in ORIGINS: L = vint_r[V] g = yoy(L) last = L.index[-1] Vd = pd.Timestamp(V) lag_months = (Vd.year - last.year) * 12 + Vd.month - last.month delays.append({"origin": V, "last_month": last.strftime("%Y-%m"), "months_behind_origin": lag_months}) # AR(1) on g over 1993-02 <= m <= last gs = g.loc["1993-01-01":last] y = gs.loc["1993-02-01":].values x = gs.shift(1).loc["1993-02-01":].values X = np.column_stack([np.ones_like(x), x]) a, b = np.linalg.lstsq(X, y, rcond=None)[0] lnL = np.log(L) drift = (lnL.iloc[-1] - lnL.iloc[-1 - DRIFT_WINDOW]) / DRIFT_WINDOW for h in HORIZONS: T = last + pd.DateOffset(months=h) if T not in actual.index: raise RuntimeError(f"no first-release actual for {T}") f_naive = g.loc[last] f_snaive = g.loc[T - pd.DateOffset(months=12)] f_ar = g.loc[last] for _ in range(h): f_ar = a + b * f_ar # model: random walk with drift in log level, drift = mean monthly log growth over # last 12 months of the vintage; YoY forecast uses the vintage's level at T-12 lnT = lnL.iloc[-1] + h * drift f_model = 100.0 * (lnT - lnL.loc[T - pd.DateOffset(months=12)]) for name, f in [("naive", f_naive), ("snaive", f_snaive), ("ar1", f_ar), ("model", f_model)]: rows.append({"origin_date": V, "target_date": T.strftime("%Y-%m-%d"), "horizon": h, "forecast": float(f), "actual_first_release": float(actual.loc[T]), "model": name, "vintage_used": V}) fc = pd.DataFrame(rows) fc.to_csv(os.path.join(O, "forecasts.csv"), index=False, float_format="%.6f") # ---- evaluation --------------------------------------------------------------- fc["err"] = fc["actual_first_release"] - fc["forecast"] res = {} tab = [] for m in ["naive", "snaive", "ar1", "model"]: for h in HORIZONS: e = fc[(fc.model == m) & (fc.horizon == h)]["err"].values rmse = float(np.sqrt(np.mean(e ** 2))); mae = float(np.mean(np.abs(e))) res[f"rmse_{m}_h{h}"] = rmse; res[f"mae_{m}_h{h}"] = mae tab.append({"model": m, "horizon": h, "rmse": rmse, "mae": mae, "n": len(e)}) pd.DataFrame(tab).to_csv(os.path.join(O, "table1.csv"), index=False, float_format="%.6f") def dm_test(d, h): d = np.asarray(d, float); n = len(d); dbar = d.mean(); dc = d - dbar lrv = np.sum(dc ** 2) / n for k in range(1, h): w = 1 - k / h lrv += 2 * w * np.sum(dc[k:] * dc[:-k]) / n stat = dbar / np.sqrt(lrv / n) p = 2 * (1 - stats.norm.cdf(abs(stat))) return float(stat), float(p), float(lrv) extra = {} for h in HORIZONS: en = fc[(fc.model == "naive") & (fc.horizon == h)].sort_values("origin_date")["err"].values em = fc[(fc.model == "model") & (fc.horizon == h)].sort_values("origin_date")["err"].values d = en ** 2 - em ** 2 s, p, lrv = dm_test(d, h) res[f"rel_improvement_rmse_h{h}"] = 1 - res[f"rmse_model_h{h}"] / res[f"rmse_naive_h{h}"] res[f"dm_stat_h{h}"] = s; res[f"dm_pvalue_h{h}"] = p extra[f"dm_lrv_h{h}"] = lrv extra[f"n_model_better_sq_h{h}"] = int(np.sum(d > 0)) better = res[f"rmse_model_h{h}"] < res[f"rmse_naive_h{h}"] sig = p < 0.05 res[f"conclusion_code_h{h}"] = ("model_better_" if better else "naive_better_") + ("significant" if sig else "not_significant") # ---- (d) revisions ------------------------------------------------------------- gcur = yoy(cur) t1 = fc[(fc.model == "naive") & (fc.horizon == 1)] revs = [abs(gcur.loc[pd.Timestamp(t)] - a_) for t, a_ in zip(t1.target_date, t1.actual_first_release)] res["mean_abs_revision_h1_targets"] = float(np.mean(revs)) extra["revision_h1_targets"] = {t: float(gcur.loc[pd.Timestamp(t)] - a_) for t, a_ in zip(t1.target_date, t1.actual_first_release)} # ---- exclusion ledger ------------------------------------------------------------ n_pairs = len(ORIGINS) * len(HORIZONS) n_eval = int(fc[fc.model == "model"]["actual_first_release"].notna().sum()) pd.DataFrame([{"step": 1, "rule": "keep (origin,horizon) pairs whose target month has a first-release actual", "n_before": n_pairs, "n_after": n_eval, "n_removed": n_pairs - n_eval, "rationale": "all 20 targets (2025-09..2026-09) are in first_release.csv; nothing removed"}] ).to_csv(os.path.join(O, "exclusions.csv"), index=False) # ---- figure ------------------------------------------------------------------------ fig, axes = plt.subplots(1, 2, figsize=(11, 4), sharey=True) for ax, h in zip(axes, HORIZONS): sub = fc[fc.horizon == h] t = pd.to_datetime(sub[sub.model == "naive"].target_date) ax.plot(t, sub[sub.model == "naive"].actual_first_release, "k-o", label="first release (actual)") ax.plot(t, sub[sub.model == "naive"].forecast, "b--s", label="naive") ax.plot(t, sub[sub.model == "model"].forecast, "r--^", label="model (RW+drift, level)") ax.set_title(f"h = {h}"); ax.set_xlabel("target month"); ax.grid(alpha=.3) ax.tick_params(axis="x", rotation=45) axes[0].set_ylabel("YoY growth, % (100·log diff)"); axes[0].legend(fontsize=8) fig.suptitle("US retail & food services sales: real-time forecasts vs first release") fig.tight_layout(); fig.savefig(os.path.join(O, "figure1.png"), dpi=120); plt.close(fig) # ---- results.json ------------------------------------------------------------------ data_issues = [ {"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", "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); " "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).")}, {"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).", "handling": "Used as given in first_release.csv (first month-end vintage containing the month)."}, {"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.", "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)."}, {"issue": "UMCSENT has 210 empty fields (pre-1978 months without survey).", "handling": "UMCSENT not used by the model; read only for row counts/audit."}, {"issue": "PAYEMS vintages also show irregular timing (2025-12-31 vintage adds two months; ends 2025-11).", "handling": "PAYEMS not used by the model; read only for audit."}, ] out = { "task_id": "T1-forecast", "model_name": "rw_drift_loglevel_12m", "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, " "mu = mean monthly log growth over the last 12 months in vintage V (= g_last/12). The YoY forecast is " "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) " "enter the forecast, unlike the naive. No estimated parameters beyond mu; no PAYEMS/UMCSENT used. " "Specification fixed a priori (not tuned on evaluation actuals); every input is in the vintage dated V."), "n_origins": len(ORIGINS), "first_origin": ORIGINS[0], "last_origin": ORIGINS[-1], } out.update(res) out.update({ "feasibility": "complete", "feasibility_reason": None, "data_issues": data_issues, "n_input_rows": n_input_rows, "n_analysed": n_eval, "random_seed": None, "publication_delays": delays, **extra, }) with open(os.path.join(O, "results.json"), "w") as f: json.dump(out, f, indent=2) for 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"]: print(k, round(out[k], 4))