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 is used (random_seed = null)."""5import json, os6import numpy as np7import pandas as pd8from scipy import stats9import matplotlib10matplotlib.use("Agg")11import matplotlib.pyplot as plt1213DATA, OUT = "data", "out"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]17MODELS = ["naive", "snaive", "ar1", "model"]18n_input_rows = {}192021def read_series(path):22 d = pd.read_csv(path)23 n_input_rows[os.path.relpath(path, DATA)] = int(len(d))24 s = pd.Series(pd.to_numeric(d.iloc[:, 1], errors="coerce").values,25 index=pd.to_datetime(d.iloc[:, 0]), name=d.columns[1])26 return s272829def yoy(level):30 lv = level.asfreq("MS")31 return 100.0 * np.log(lv / lv.shift(12))323334def add_months(ts, k):35 return ts + pd.DateOffset(months=k)363738# ---------------- forecasting functions (receive only vintage-V information) ----------------39def f_naive(L, last, h):40 return float(yoy(L)[last])414243def f_snaive(L, last, h):44 T = add_months(last, h)45 return float(yoy(L)[add_months(T, -12)])464748def f_ar1(L, last, h):49 g = yoy(L)50 y = g.loc["1993-02-01":last]51 x = g.shift(1).loc["1993-02-01":last]52 ok = y.notna() & x.notna()53 X = np.column_stack([np.ones(ok.sum()), x[ok].values])54 beta = np.linalg.lstsq(X, y[ok].values, rcond=None)[0]55 f = float(g[last])56 for _ in range(h):57 f = beta[0] + beta[1] * f58 return float(f)596061def f_model(L, last, h):62 """Base-effect-aware random walk with drift in log levels.63 drift = mean monthly log change over the 12 months ending at last(V) (= g_last/1200);64 ln L_T = ln L_last + h*drift; YoY forecast uses the known base L_{T-12} of the same vintage."""65 lv = np.log(L.asfreq("MS"))66 drift = (lv[last] - lv[add_months(last, -12)]) / 12.067 T = add_months(last, h)68 return float(100.0 * (lv[last] + h * drift - lv[add_months(T, -12)]))697071FUNCS = {"naive": f_naive, "snaive": f_snaive, "ar1": f_ar1, "model": f_model}727374def dm_test(e_naive, e_model, h):75 d = np.asarray(e_naive) ** 2 - np.asarray(e_model) ** 276 n = len(d)77 dbar = d.mean()78 dc = d - dbar79 lrv = np.sum(dc * dc) / n80 for k in range(1, h):81 w = 1.0 - k / h82 lrv += 2.0 * w * np.sum(dc[k:] * dc[:-k]) / n83 stat = dbar / np.sqrt(lrv / n)84 p = 2.0 * (1.0 - stats.norm.cdf(abs(stat)))85 return float(stat), float(p), float(lrv)868788def main():89 fr = pd.read_csv(os.path.join(DATA, "first_release.csv"), parse_dates=["target_month"])90 n_input_rows["first_release.csv"] = int(len(fr))91 actual = fr.set_index("target_month")["yoy_first_release"]9293 # predictors read and audited (not used by the chosen model); availability rule respected94 umc = read_series(os.path.join(DATA, "UMCSENT.csv"))9596 rows, design = [], []97 for V in ORIGINS:98 L = read_series(os.path.join(DATA, "vintages", f"RSAFS_vintage_{V}.csv"))99 P = read_series(os.path.join(DATA, "vintages", f"PAYEMS_vintage_{V}.csv"))100 last = L.index.max()101 Vts = pd.Timestamp(V)102 assert last < Vts and P.index.max() < Vts103 design.append({"origin": V, "last_rsafs": str(last.date()), "last_payems": str(P.index.max().date()),104 "lag_months": (Vts.year - last.year) * 12 + Vts.month - last.month})105 for h in HORIZONS:106 T = add_months(last, h)107 for m in MODELS:108 rows.append({"origin_date": V, "target_date": T.strftime("%Y-%m-%d"), "horizon": h,109 "forecast": FUNCS[m](L, last, h),110 "actual_first_release": float(actual[T]), "model": m, "vintage_used": V})111 fc = pd.DataFrame(rows)112 fc["error"] = fc["actual_first_release"] - fc["forecast"]113 fc.drop(columns="error").to_csv(os.path.join(OUT, "forecasts.csv"), index=False, float_format="%.6f")114115 res = {"task_id": "T1-forecast",116 "model_name": "RW-drift-base-effect",117 "model_description": (118 "Base-effect-aware random walk with drift in log levels, applied to the RSAFS vintage dated V only. "119 "drift = mean monthly log change of RSAFS over the 12 months ending at last(V) (equal to g_last/1200); "120 "forecast log level ln L_T = ln L_last(V) + h*drift for T = last(V)+h; YoY forecast = 100*ln(L_T / L_{T-12}), "121 "where L_{T-12} is the (already published) base-month level in the same vintage. No estimated parameters, "122 "no estimation window beyond the last 13 levels of the vintage; PAYEMS and UMCSENT are not used, so the "123 "only information is the RSAFS vintage V (satisfies the availability rule). Specification fixed before any "124 "evaluation error was computed.")}125 res["n_origins"] = len(ORIGINS)126 res["first_origin"], res["last_origin"] = ORIGINS[0], ORIGINS[-1]127 table = []128 for m in MODELS:129 for h in HORIZONS:130 e = fc[(fc.model == m) & (fc.horizon == h)]["error"].values131 rmse, mae = float(np.sqrt(np.mean(e ** 2))), float(np.mean(np.abs(e)))132 res[f"rmse_{m}_h{h}"], res[f"mae_{m}_h{h}"] = rmse, mae133 table.append({"model": m, "horizon": h, "rmse": rmse, "mae": mae, "n": int(len(e))})134 pd.DataFrame(table).to_csv(os.path.join(OUT, "table1.csv"), index=False, float_format="%.6f")135 extra = {}136 for h in HORIZONS:137 res[f"rel_improvement_rmse_h{h}"] = 1.0 - res[f"rmse_model_h{h}"] / res[f"rmse_naive_h{h}"]138 en = fc[(fc.model == "naive") & (fc.horizon == h)].sort_values("origin_date")["error"].values139 em = fc[(fc.model == "model") & (fc.horizon == h)].sort_values("origin_date")["error"].values140 s, p, lrv = dm_test(en, em, h)141 res[f"dm_stat_h{h}"], res[f"dm_pvalue_h{h}"] = s, p142 extra[f"dm_lrv_h{h}"] = lrv143 extra[f"mean_loss_diff_h{h}"] = float(np.mean(en ** 2 - em ** 2))144 n = len(en)145 hln = np.sqrt((n + 1 - 2 * h + h * (h - 1) / n) / n) # sensitivity only, not the required statistic146 extra[f"dm_stat_hln_h{h}"] = float(s * hln)147 extra[f"dm_pvalue_hln_t_h{h}"] = float(2 * (1 - stats.t.cdf(abs(s * hln), df=n - 1)))148 extra[f"n_origins_model_sqerr_lower_h{h}"] = int(np.sum(em ** 2 < en ** 2))149 better = "model" if res[f"rmse_model_h{h}"] < res[f"rmse_naive_h{h}"] else "naive"150 sig = "significant" if p < 0.05 else "not_significant"151 res[f"conclusion_code_h{h}"] = f"{better}_better_{sig}"152153 # ---------------- estimand (d): revision analysis (only place RSAFS_current is used) ----------------154 cur = read_series(os.path.join(DATA, "RSAFS_current.csv"))155 gcur = yoy(cur)156 t1 = fc[(fc.model == "naive") & (fc.horizon == 1)]157 tgt = pd.to_datetime(t1["target_date"]).values158 rev = np.abs(gcur.loc[tgt].values - t1["actual_first_release"].values)159 res["mean_abs_revision_h1_targets"] = float(np.mean(rev))160 extra["mean_revision_h1_targets_signed"] = float(np.mean(gcur.loc[tgt].values - t1["actual_first_release"].values))161162 res["feasibility"] = "complete"163 res["feasibility_reason"] = None164 res["data_issues"] = [165 {"issue": "Irregular publication timing: RSAFS 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.",166 "handling": "Targets defined as last(V)+h per the task rule; h=1 targets therefore skip 2026-03 and the 2026-03-31 origin forecasts 2026-02 (h=1) and 2026-04 (h=3)."},167 {"issue": "First releases of 2026-02 and 2026-03 both date from the 2026-04-30 vintage; 2025-09 first appears in the 2025-11-30 vintage (publication disruption).",168 "handling": "Used first_release.csv as given (month-end granularity definition)."},169 {"issue": "PAYEMS vintage 2025-12-31 adds two months at once; PAYEMS annual benchmark revision (~ -800 thousand) in vintage 2026-02-28; ragged edge between PAYEMS and RSAFS.",170 "handling": "PAYEMS read and audited but not used by the model; no effect on forecasts."},171 {"issue": "Large RSAFS benchmark level revision after the last origin (current levels ~4% below first-release levels; first-release levels drop between 2026-07 and 2026-08 targets); annual revision visible in vintage 2026-07-31.",172 "handling": "Growth rates are always computed within a single vintage, never across vintages; the revision is quantified by estimand (d) using RSAFS_current only there."},173 {"issue": "UMCSENT has 210 empty values before 1978 (quarterly survey era); current vintage used.",174 "handling": "Not used by the model; admissible only for months <= month(V) by rule."},175 {"issue": "COVID-19 period (2020-2021) produces extreme YoY swings inside the AR(1) estimation sample.",176 "handling": "AR(1) baseline estimated on the full sample 1993-02..last(V) as the task prescribes; no trimming."},177 ]178 res["n_input_rows"] = dict(sorted(n_input_rows.items()))179 res["n_analysed"] = int(len(fc[fc.model == "model"]))180 res["random_seed"] = None181 res["supplementary"] = {"real_time_design": design, **extra}182 with open(os.path.join(OUT, "results.json"), "w") as f:183 json.dump(res, f, indent=2)184185 # ---------------- figure ----------------186 fig, axes = plt.subplots(1, 2, figsize=(12, 4.5), sharey=True)187 for ax, h in zip(axes, HORIZONS):188 sub = fc[fc.horizon == h]189 for m, st in [("naive", "s--"), ("model", "o-")]:190 s = sub[sub.model == m]191 ax.plot(pd.to_datetime(s.target_date), s.forecast, st, label=f"{m} forecast")192 a = fr[(fr.target_month >= "2025-06-01")]193 ax.plot(a.target_month, a.yoy_first_release, "k-", lw=2, label="first-release actual")194 ax.set_title(f"h = {h}"); ax.set_ylabel("YoY growth, %"); ax.grid(alpha=.3)195 ax.tick_params(axis="x", rotation=45)196 axes[0].legend()197 fig.suptitle("RSAFS YoY growth: first-release actuals vs real-time forecasts (target months)")198 fig.tight_layout(); fig.savefig(os.path.join(OUT, "figure1.png"), dpi=120)199200 print(pd.DataFrame(table).to_string(index=False))201 print({k: res[k] for k in res if k.startswith(("rel_", "dm_", "conclusion", "mean_abs"))})202203204if __name__ == "__main__":205 main()206