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.4Deterministic (no randomness except the seeded bootstrap below)."""5import json, os6import numpy as np7import pandas as pd8from scipy import stats9import matplotlib10matplotlib.use("Agg")11import matplotlib.pyplot as plt1213DATA, OUT = "data", "out"14SEED = 2026103015ORIGINS = ["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(path, key):23 df = pd.read_csv(path)24 n_input_rows[key] = int(len(df))25 s = df.set_index(pd.PeriodIndex(pd.to_datetime(df.iloc[:, 0]), freq="M")).iloc[:, 1].astype(float)26 return s272829def yoy(level):30 lg = np.log(level)31 return 100.0 * (lg - lg.shift(12, freq="M").reindex(lg.index))323334def nw_dm(d, h):35 """DM stat with Bartlett long-run variance, truncation lag h-1, no HLN factor."""36 d = np.asarray(d, float); n = len(d); dbar = d.mean(); u = d - dbar37 lrv = np.sum(u * u) / n38 L = h - 139 for k in range(1, L + 1):40 w = 1.0 - k / (L + 1.0)41 lrv += 2.0 * w * np.sum(u[k:] * u[:-k]) / n42 stat = dbar / np.sqrt(lrv / n)43 p = 2.0 * (1.0 - stats.norm.cdf(abs(stat)))44 return float(stat), float(p), float(lrv)454647# ---------------------------------------------------------------- inputs48fr = pd.read_csv(os.path.join(DATA, "first_release.csv"))49n_input_rows["first_release.csv"] = int(len(fr))50actual = pd.Series(fr["yoy_first_release"].values,51 index=pd.PeriodIndex(pd.to_datetime(fr["target_month"]), freq="M"))52rs_vint = {V: read_series(os.path.join(DATA, "vintages", f"RSAFS_vintage_{V}.csv"),53 f"vintages/RSAFS_vintage_{V}.csv") for V in ORIGINS}54# PAYEMS / UMCSENT are read only for the audit row counts; they are not inputs to any forecast.55for V in ORIGINS:56 read_series(os.path.join(DATA, "vintages", f"PAYEMS_vintage_{V}.csv"), f"vintages/PAYEMS_vintage_{V}.csv")57read_series(os.path.join(DATA, "UMCSENT.csv"), "UMCSENT.csv")58rs_current = read_series(os.path.join(DATA, "RSAFS_current.csv"), "RSAFS_current.csv") # estimand (d) only5960# ---------------------------------------------------------------- forecasts61rows, ar_info, lastmonths = [], {}, {}62for V in ORIGINS:63 L = rs_vint[V]64 assert L.index.max() < pd.Period(V, "M") # last published month precedes the origin month65 g = yoy(L).dropna()66 last = L.index.max(); lastmonths[V] = str(last)67 lnL = np.log(L)68 # AR(1): g_m on const, g_{m-1}, 1993-02 <= m <= last69 y = g.loc[pd.Period("1993-02", "M"):last]70 x = g.shift(1, freq="M").reindex(y.index)71 X = np.column_stack([np.ones(len(y)), x.values])72 c, phi = np.linalg.lstsq(X, y.values, rcond=None)[0]73 ar_info[V] = {"const": float(c), "phi": float(phi), "n": int(len(y))}74 for h in HORIZONS:75 T = last + h76 fc = {}77 fc["naive"] = g[last]78 fc["snaive"] = g[T - 12]79 f = g[last]80 for _ in range(h):81 f = c + phi * f82 fc["ar1"] = f83 # Own model: base-effect random walk with drift in log level.84 # ln L_T forecast = ln L_last + h * drift, drift = g_last/12 (avg monthly log growth over last 12 months);85 # YoY forecast = 100*(ln L_T_hat - ln L_{T-12}), with L_{T-12} known in the same vintage.86 drift = g[last] / 100.0 / 12.087 fc["model"] = 100.0 * (lnL[last] + h * drift - lnL[T - 12])88 for m in MODELS:89 rows.append({"origin_date": V, "target_date": T.to_timestamp().strftime("%Y-%m-%d"),90 "horizon": h, "forecast": float(fc[m]),91 "actual_first_release": float(actual[T]) if T in actual.index else np.nan,92 "model": m, "vintage_used": V})93fc_df = pd.DataFrame(rows)94fc_df["error"] = fc_df["actual_first_release"] - fc_df["forecast"]95fc_df[["origin_date", "target_date", "horizon", "forecast", "actual_first_release", "model", "vintage_used"]] \96 .to_csv(os.path.join(OUT, "forecasts.csv"), index=False, float_format="%.6f")9798# ---------------------------------------------------------------- evaluation99res = {"task_id": "T1-forecast",100 "model_name": "BE-RWD (base-effect random walk with drift in log level)",101 "model_description": (102 "Univariate real-time model on the RSAFS vintage dated V only. The log level is forecast as a random walk "103 "with drift: ln L_hat(last+h) = ln L(last) + h*d, with d = g_last/1200, i.e. the average monthly log growth "104 "over the 12 months ending at last(V) in that vintage. The YoY forecast is 100*(ln L_hat(T) - ln L(T-12)), "105 "where L(T-12) is read from the same vintage, so the forecast equals the naive forecast plus a known "106 "base-effect correction (the change in the year-earlier comparison level) and drift. No parameters are "107 "estimated; PAYEMS and UMCSENT are not used. Specification fixed before looking at any evaluation outcome; "108 "no alternative models were tried on the evaluation sample."),109 "n_origins": len(ORIGINS), "first_origin": ORIGINS[0], "last_origin": ORIGINS[-1]}110table = []111for m in MODELS:112 for h in HORIZONS:113 e = fc_df.query("model==@m and horizon==@h")["error"].values114 rmse, mae = float(np.sqrt(np.mean(e ** 2))), float(np.mean(np.abs(e)))115 res[f"rmse_{m}_h{h}"], res[f"mae_{m}_h{h}"] = rmse, mae116 table.append({"model": m, "horizon": h, "rmse": rmse, "mae": mae, "n": int(len(e))})117pd.DataFrame(table).to_csv(os.path.join(OUT, "table1.csv"), index=False, float_format="%.6f")118119extra = {}120rng = np.random.default_rng(SEED)121for h in HORIZONS:122 en = fc_df.query("model=='naive' and horizon==@h").sort_values("origin_date")["error"].values123 em = fc_df.query("model=='model' and horizon==@h").sort_values("origin_date")["error"].values124 res[f"rel_improvement_rmse_h{h}"] = 1.0 - res[f"rmse_model_h{h}"] / res[f"rmse_naive_h{h}"]125 d = en ** 2 - em ** 2126 s, p, lrv = nw_dm(d, h)127 res[f"dm_stat_h{h}"], res[f"dm_pvalue_h{h}"] = s, p128 better = "model" if res[f"rmse_model_h{h}"] < res[f"rmse_naive_h{h}"] else "naive"129 res[f"conclusion_code_h{h}"] = f"{better}_better_{'significant' if p < 0.05 else 'not_significant'}"130 # supplementary: moving-block bootstrap (block length max(h,2)) of RMSE_naive - RMSE_model131 B, bl, n = 5000, max(h, 2), len(d)132 diffs = []133 for _ in range(B):134 starts = rng.integers(0, n - bl + 1, size=int(np.ceil(n / bl)))135 idx = np.concatenate([np.arange(s0, s0 + bl) for s0 in starts])[:n]136 diffs.append(np.sqrt(np.mean(en[idx] ** 2)) - np.sqrt(np.mean(em[idx] ** 2)))137 lo, hi = np.percentile(diffs, [2.5, 97.5])138 res[f"boot_rmse_diff_ci95_lo_h{h}"], res[f"boot_rmse_diff_ci95_hi_h{h}"] = float(lo), float(hi)139 res[f"share_origins_model_smaller_sq_error_h{h}"] = float(np.mean(d > 0))140 extra[f"h{h}"] = {"dm_lrv": lrv, "mean_loss_diff": float(d.mean()),141 "rmse_diff_naive_minus_model": res[f"rmse_naive_h{h}"] - res[f"rmse_model_h{h}"],142 "boot_rmse_diff_ci95": [float(lo), float(hi)], "boot_B": B, "boot_block": bl,143 "share_origins_model_smaller_sq_error": float(np.mean(d > 0))}144145# (d) revisions: h=1 targets, current-vintage YoY vs first release146g_cur = yoy(rs_current)147t1 = fc_df.query("model=='naive' and horizon==1")148tm = pd.PeriodIndex(pd.to_datetime(t1["target_date"]), freq="M")149rev = np.abs(g_cur.reindex(tm).values - actual.reindex(tm).values)150res["mean_abs_revision_h1_targets"] = float(np.mean(rev))151extra["revision_by_target"] = {str(k): float(v) for k, v in zip(tm, g_cur.reindex(tm).values - actual.reindex(tm).values)}152extra["ar1_params"] = ar_info153extra["last_month_by_origin"] = lastmonths154155res["feasibility"] = "complete"156res["feasibility_reason"] = None157res["data_issues"] = [158 {"issue": "Irregular RSAFS publication timing: last(V)=2025-08 at V=2025-10-31 and 2026-01 at V=2026-03-31 (two-month lag); V=2026-04-30 adds two months at once.",159 "handling": "Targets defined as last(V)+h per task rule; the 2025-10-31 and 2026-03-31 origins therefore forecast months that are 2 and 4 months past the last calendar month."},160 {"issue": "h=3 forecast errors of consecutive origins overlap (and some h=1 and h=3 target months coincide), inducing serial correlation in the loss differentials; only 10 origins are available.",161 "handling": "DM long-run variance uses Bartlett kernel with truncation lag h-1 as specified."},162 {"issue": "Large level revision of RSAFS after the 2026-07-31 vintage (current levels ~4-5% below vintage levels); first_release.csv shows the break between targets 2026-07 and 2026-08.",163 "handling": "YoY always computed within one vintage; h=3 actuals for targets 2026-08/2026-09 are post-revision first releases, used as given; estimand (d) reflects this revision."},164 {"issue": "PAYEMS benchmark revision visible in the 2026-02-28 vintage; PAYEMS 2025-10 first appears in the 2025-12-31 vintage.",165 "handling": "PAYEMS not used in any forecast; documented only."},166 {"issue": "UMCSENT has 210 empty values before 1978-01.",167 "handling": "UMCSENT not used in any forecast; documented only."},168]169res["n_input_rows"] = dict(sorted(n_input_rows.items()))170res["n_analysed"] = int(fc_df.query("model=='model'")["error"].notna().sum())171res["random_seed"] = SEED172173with open(os.path.join(OUT, "results.json"), "w") as fh:174 json.dump(res, fh, indent=1)175with open(os.path.join(OUT, "supplementary.json"), "w") as fh:176 json.dump(extra, fh, indent=1)177178# ---------------------------------------------------------------- figure179fig, axes = plt.subplots(1, 2, figsize=(12, 4.5), sharey=True)180for ax, h in zip(axes, HORIZONS):181 sub = fc_df.query("horizon==@h")182 for m, sty in [("naive", "s--"), ("model", "o-"), ("ar1", "^:")]:183 s = sub.query("model==@m").sort_values("target_date")184 ax.plot(pd.to_datetime(s["target_date"]), s["forecast"], sty, label=m, alpha=.8)185 s = sub.query("model=='naive'").sort_values("target_date")186 ax.plot(pd.to_datetime(s["target_date"]), s["actual_first_release"], "k-", lw=2.2, label="first release")187 ax.set_title(f"h = {h}: RSAFS YoY growth (%), real-time forecasts")188 ax.set_xlabel("target month"); ax.grid(alpha=.3); ax.tick_params(axis="x", rotation=45)189axes[0].set_ylabel("percent"); axes[0].legend()190fig.tight_layout(); fig.savefig(os.path.join(OUT, "figure1.png"), dpi=120)191192print(pd.DataFrame(table).to_string(index=False))193for h in HORIZONS:194 print(f"h={h}: rel_impr={res[f'rel_improvement_rmse_h{h}']:.4f} DM={res[f'dm_stat_h{h}']:.4f} p={res[f'dm_pvalue_h{h}']:.4f} {res[f'conclusion_code_h{h}']}")195print("mean_abs_revision_h1_targets", round(res["mean_abs_revision_h1_targets"], 4))196