"""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. Deterministic: no randomness is used (random_seed = null).""" import json, os import numpy as np import pandas as pd from scipy import stats import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt DATA, OUT = "data", "out" 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] MODELS = ["naive", "snaive", "ar1", "model"] n_input_rows = {} def read_series(path): d = pd.read_csv(path) n_input_rows[os.path.relpath(path, DATA)] = int(len(d)) s = pd.Series(pd.to_numeric(d.iloc[:, 1], errors="coerce").values, index=pd.to_datetime(d.iloc[:, 0]), name=d.columns[1]) return s def yoy(level): lv = level.asfreq("MS") return 100.0 * np.log(lv / lv.shift(12)) def add_months(ts, k): return ts + pd.DateOffset(months=k) # ---------------- forecasting functions (receive only vintage-V information) ---------------- def f_naive(L, last, h): return float(yoy(L)[last]) def f_snaive(L, last, h): T = add_months(last, h) return float(yoy(L)[add_months(T, -12)]) def f_ar1(L, last, h): g = yoy(L) y = g.loc["1993-02-01":last] x = g.shift(1).loc["1993-02-01":last] ok = y.notna() & x.notna() X = np.column_stack([np.ones(ok.sum()), x[ok].values]) beta = np.linalg.lstsq(X, y[ok].values, rcond=None)[0] f = float(g[last]) for _ in range(h): f = beta[0] + beta[1] * f return float(f) def f_model(L, last, h): """Base-effect-aware random walk with drift in log levels. drift = mean monthly log change over the 12 months ending at last(V) (= g_last/1200); ln L_T = ln L_last + h*drift; YoY forecast uses the known base L_{T-12} of the same vintage.""" lv = np.log(L.asfreq("MS")) drift = (lv[last] - lv[add_months(last, -12)]) / 12.0 T = add_months(last, h) return float(100.0 * (lv[last] + h * drift - lv[add_months(T, -12)])) FUNCS = {"naive": f_naive, "snaive": f_snaive, "ar1": f_ar1, "model": f_model} def dm_test(e_naive, e_model, h): d = np.asarray(e_naive) ** 2 - np.asarray(e_model) ** 2 n = len(d) dbar = d.mean() dc = d - dbar lrv = np.sum(dc * dc) / n for k in range(1, h): w = 1.0 - k / h lrv += 2.0 * w * np.sum(dc[k:] * dc[:-k]) / n stat = dbar / np.sqrt(lrv / n) p = 2.0 * (1.0 - stats.norm.cdf(abs(stat))) return float(stat), float(p), float(lrv) def main(): fr = pd.read_csv(os.path.join(DATA, "first_release.csv"), parse_dates=["target_month"]) n_input_rows["first_release.csv"] = int(len(fr)) actual = fr.set_index("target_month")["yoy_first_release"] # predictors read and audited (not used by the chosen model); availability rule respected umc = read_series(os.path.join(DATA, "UMCSENT.csv")) rows, design = [], [] for V in ORIGINS: L = read_series(os.path.join(DATA, "vintages", f"RSAFS_vintage_{V}.csv")) P = read_series(os.path.join(DATA, "vintages", f"PAYEMS_vintage_{V}.csv")) last = L.index.max() Vts = pd.Timestamp(V) assert last < Vts and P.index.max() < Vts design.append({"origin": V, "last_rsafs": str(last.date()), "last_payems": str(P.index.max().date()), "lag_months": (Vts.year - last.year) * 12 + Vts.month - last.month}) for h in HORIZONS: T = add_months(last, h) for m in MODELS: rows.append({"origin_date": V, "target_date": T.strftime("%Y-%m-%d"), "horizon": h, "forecast": FUNCS[m](L, last, h), "actual_first_release": float(actual[T]), "model": m, "vintage_used": V}) fc = pd.DataFrame(rows) fc["error"] = fc["actual_first_release"] - fc["forecast"] fc.drop(columns="error").to_csv(os.path.join(OUT, "forecasts.csv"), index=False, float_format="%.6f") res = {"task_id": "T1-forecast", "model_name": "RW-drift-base-effect", "model_description": ( "Base-effect-aware random walk with drift in log levels, applied to the RSAFS vintage dated V only. " "drift = mean monthly log change of RSAFS over the 12 months ending at last(V) (equal to g_last/1200); " "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}), " "where L_{T-12} is the (already published) base-month level in the same vintage. No estimated parameters, " "no estimation window beyond the last 13 levels of the vintage; PAYEMS and UMCSENT are not used, so the " "only information is the RSAFS vintage V (satisfies the availability rule). Specification fixed before any " "evaluation error was computed.")} res["n_origins"] = len(ORIGINS) res["first_origin"], res["last_origin"] = ORIGINS[0], ORIGINS[-1] table = [] for m in MODELS: for h in HORIZONS: e = fc[(fc.model == m) & (fc.horizon == h)]["error"].values rmse, mae = float(np.sqrt(np.mean(e ** 2))), float(np.mean(np.abs(e))) res[f"rmse_{m}_h{h}"], res[f"mae_{m}_h{h}"] = rmse, mae table.append({"model": m, "horizon": h, "rmse": rmse, "mae": mae, "n": int(len(e))}) pd.DataFrame(table).to_csv(os.path.join(OUT, "table1.csv"), index=False, float_format="%.6f") extra = {} for h in HORIZONS: res[f"rel_improvement_rmse_h{h}"] = 1.0 - res[f"rmse_model_h{h}"] / res[f"rmse_naive_h{h}"] en = fc[(fc.model == "naive") & (fc.horizon == h)].sort_values("origin_date")["error"].values em = fc[(fc.model == "model") & (fc.horizon == h)].sort_values("origin_date")["error"].values s, p, lrv = dm_test(en, em, h) res[f"dm_stat_h{h}"], res[f"dm_pvalue_h{h}"] = s, p extra[f"dm_lrv_h{h}"] = lrv extra[f"mean_loss_diff_h{h}"] = float(np.mean(en ** 2 - em ** 2)) n = len(en) hln = np.sqrt((n + 1 - 2 * h + h * (h - 1) / n) / n) # sensitivity only, not the required statistic extra[f"dm_stat_hln_h{h}"] = float(s * hln) extra[f"dm_pvalue_hln_t_h{h}"] = float(2 * (1 - stats.t.cdf(abs(s * hln), df=n - 1))) extra[f"n_origins_model_sqerr_lower_h{h}"] = int(np.sum(em ** 2 < en ** 2)) better = "model" if res[f"rmse_model_h{h}"] < res[f"rmse_naive_h{h}"] else "naive" sig = "significant" if p < 0.05 else "not_significant" res[f"conclusion_code_h{h}"] = f"{better}_better_{sig}" # ---------------- estimand (d): revision analysis (only place RSAFS_current is used) ---------------- cur = read_series(os.path.join(DATA, "RSAFS_current.csv")) gcur = yoy(cur) t1 = fc[(fc.model == "naive") & (fc.horizon == 1)] tgt = pd.to_datetime(t1["target_date"]).values rev = np.abs(gcur.loc[tgt].values - t1["actual_first_release"].values) res["mean_abs_revision_h1_targets"] = float(np.mean(rev)) extra["mean_revision_h1_targets_signed"] = float(np.mean(gcur.loc[tgt].values - t1["actual_first_release"].values)) res["feasibility"] = "complete" res["feasibility_reason"] = None res["data_issues"] = [ {"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.", "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)."}, {"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).", "handling": "Used first_release.csv as given (month-end granularity definition)."}, {"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.", "handling": "PAYEMS read and audited but not used by the model; no effect on forecasts."}, {"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.", "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."}, {"issue": "UMCSENT has 210 empty values before 1978 (quarterly survey era); current vintage used.", "handling": "Not used by the model; admissible only for months <= month(V) by rule."}, {"issue": "COVID-19 period (2020-2021) produces extreme YoY swings inside the AR(1) estimation sample.", "handling": "AR(1) baseline estimated on the full sample 1993-02..last(V) as the task prescribes; no trimming."}, ] res["n_input_rows"] = dict(sorted(n_input_rows.items())) res["n_analysed"] = int(len(fc[fc.model == "model"])) res["random_seed"] = None res["supplementary"] = {"real_time_design": design, **extra} with open(os.path.join(OUT, "results.json"), "w") as f: json.dump(res, f, indent=2) # ---------------- figure ---------------- fig, axes = plt.subplots(1, 2, figsize=(12, 4.5), sharey=True) for ax, h in zip(axes, HORIZONS): sub = fc[fc.horizon == h] for m, st in [("naive", "s--"), ("model", "o-")]: s = sub[sub.model == m] ax.plot(pd.to_datetime(s.target_date), s.forecast, st, label=f"{m} forecast") a = fr[(fr.target_month >= "2025-06-01")] ax.plot(a.target_month, a.yoy_first_release, "k-", lw=2, label="first-release actual") ax.set_title(f"h = {h}"); ax.set_ylabel("YoY growth, %"); ax.grid(alpha=.3) ax.tick_params(axis="x", rotation=45) axes[0].legend() fig.suptitle("RSAFS YoY growth: first-release actuals vs real-time forecasts (target months)") fig.tight_layout(); fig.savefig(os.path.join(OUT, "figure1.png"), dpi=120) print(pd.DataFrame(table).to_string(index=False)) print({k: res[k] for k in res if k.startswith(("rel_", "dm_", "conclusion", "mean_abs"))}) if __name__ == "__main__": main()