"""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 (bootstrap seed fixed).""" import json, os import numpy as np import pandas as pd from scipy.stats import norm import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt DATA, OUT = "data", "out" SEED = 20261030 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(rel): d = pd.read_csv(os.path.join(DATA, rel), parse_dates=["observation_date"]) n_input_rows[rel] = int(len(d)) return d.set_index("observation_date").iloc[:, 0].astype(float) def yoy(level): """g_m = 100 ln(L_m / L_{m-12}) on a gap-free monthly series.""" s = level.asfreq("MS") assert not s.isna().any(), "gap in level series" return (100 * np.log(s / s.shift(12))).dropna() def mshift(t, k): return t + pd.DateOffset(months=k) # ---------------------------------------------------------------- inputs fr = pd.read_csv(os.path.join(DATA, "first_release.csv"), parse_dates=["target_month"]) n_input_rows["first_release.csv"] = int(len(fr)) fr = fr.set_index("target_month") # PAYEMS vintages and UMCSENT are candidate predictors; they are read for the # audit row counts only and do NOT enter the chosen model. for V in ORIGINS: read_series(f"vintages/PAYEMS_vintage_{V}.csv") read_series("UMCSENT.csv") # ---------------------------------------------------------------- forecasts rows, origin_info = [], [] for V in ORIGINS: L = read_series(f"vintages/RSAFS_vintage_{V}.csv") # only vintage V g = yoy(L) last = L.index[-1] assert last < pd.Timestamp(V), "vintage contains month at/after origin" origin_info.append({"origin": V, "last": last.strftime("%Y-%m"), "lag_months": (pd.Timestamp(V).year - last.year) * 12 + pd.Timestamp(V).month - last.month}) g_last = g[last] # AR(1): g_m on const + g_{m-1}, 1993-02 <= m <= last gs = g[g.index >= "1993-01-01"] y = gs[gs.index >= "1993-02-01"].values x = gs.shift(1)[gs.index >= "1993-02-01"].values X = np.column_stack([np.ones_like(x), x]) (c, phi), *_ = np.linalg.lstsq(X, y, rcond=None) for h in HORIZONS: T = mshift(last, h) actual = float(fr.loc[T, "yoy_first_release"]) # naive f_naive = g_last # seasonal naive f_snaive = g[mshift(T, -12)] # AR(1) iterated f = g_last for _ in range(h): f = c + phi * f f_ar1 = f # own model: base-effect drift f_model = 100 * np.log(L[last] / L[mshift(T, -12)]) + h * g_last / 12.0 for m, fc in zip(MODELS, [f_naive, f_snaive, f_ar1, f_model]): rows.append({"origin_date": V, "target_date": T.strftime("%Y-%m-%d"), "horizon": h, "forecast": round(float(fc), 6), "actual_first_release": actual, "model": m, "vintage_used": V}) fc = pd.DataFrame(rows) assert len(fc) == 80 fc.to_csv(os.path.join(OUT, "forecasts.csv"), index=False) fc["err"] = fc["actual_first_release"] - fc["forecast"] # ---------------------------------------------------------------- evaluation res = {} tab = [] for m in MODELS: for h in HORIZONS: e = fc[(fc.model == m) & (fc.horizon == h)]["err"].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 tab.append({"model": m, "horizon": h, "rmse": rmse, "mae": mae, "n": len(e)}) pd.DataFrame(tab).to_csv(os.path.join(OUT, "table1.csv"), index=False) def dm_test(e_naive, e_model, h): d = e_naive ** 2 - e_model ** 2 n = len(d) dbar = d.mean() dc = d - dbar lrv = np.sum(dc * dc) / n for k in range(1, h): # Bartlett, truncation lag h-1 w = 1 - k / h lrv += 2 * w * np.sum(dc[k:] * dc[:-k]) / n stat = dbar / np.sqrt(lrv / n) return float(stat), float(2 * (1 - norm.cdf(abs(stat)))), float(lrv) boot = {} rng = np.random.default_rng(SEED) for h in HORIZONS: sub = fc[fc.horizon == h] en = sub[sub.model == "naive"].sort_values("origin_date")["err"].values em = sub[sub.model == "model"].sort_values("origin_date")["err"].values stat, p, lrv = dm_test(en, em, h) res[f"dm_stat_h{h}"], res[f"dm_pvalue_h{h}"] = stat, p res[f"rel_improvement_rmse_h{h}"] = 1 - res[f"rmse_model_h{h}"] / res[f"rmse_naive_h{h}"] 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") # supplementary: moving-block bootstrap of RMSE_model - RMSE_naive (block length max(h,2)) n, b = len(en), max(h, 2) starts = np.arange(n - b + 1) draws = [] for _ in range(9999): idx = np.concatenate([np.arange(s, s + b) for s in rng.choice(starts, size=int(np.ceil(n / b)))])[:n] draws.append(np.sqrt(np.mean(em[idx] ** 2)) - np.sqrt(np.mean(en[idx] ** 2))) draws = np.array(draws) boot[f"h{h}"] = {"rmse_diff_model_minus_naive": float(np.sqrt(np.mean(em ** 2)) - np.sqrt(np.mean(en ** 2))), "ci95": [float(np.percentile(draws, 2.5)), float(np.percentile(draws, 97.5))], "share_draws_model_better": float(np.mean(draws < 0)), "block_length": b, "n_boot": 9999, "dm_lrv": lrv, "n_origins_model_smaller_sq_error": int(np.sum(em ** 2 < en ** 2))} # ---------------------------------------------------------------- revision (estimand d) # RSAFS_current is used ONLY here; it never enters a forecast. cur = read_series("RSAFS_current.csv") gcur = yoy(cur) t1 = fc[(fc.horizon == 1) & (fc.model == "naive")] targets = pd.to_datetime(t1["target_date"]).tolist() revs = [float(gcur[T] - fr.loc[T, "yoy_first_release"]) for T in targets] res["mean_abs_revision_h1_targets"] = float(np.mean(np.abs(revs))) # ---------------------------------------------------------------- exclusions ledger excl = pd.DataFrame([ {"step": 1, "rule": "first_release.csv: keep only rows that are targets last(V)+h of some origin", "n_before": len(fr), "n_after": int(fc.target_date.nunique()), "n_removed": len(fr) - int(fc.target_date.nunique()), "rationale": "rows outside the evaluation design are not needed as actuals (not an error-driven exclusion)"}, {"step": 2, "rule": "forecast-actual pairs for own model: drop pairs with missing forecast or actual", "n_before": 20, "n_after": int(fc[fc.model == "model"].dropna(subset=["forecast", "actual_first_release"]).shape[0]), "n_removed": 20 - int(fc[fc.model == "model"].dropna(subset=["forecast", "actual_first_release"]).shape[0]), "rationale": "no pair is missing; nothing removed"}, ]) excl.to_csv(os.path.join(OUT, "exclusions.csv"), index=False) # ---------------------------------------------------------------- figure fig, axes = plt.subplots(2, 1, figsize=(9, 7), sharex=True) frp = fr.loc["2025-06-01":"2026-09-01", "yoy_first_release"] for ax, h in zip(axes, HORIZONS): ax.plot(frp.index, frp.values, "k-o", ms=4, label="first-release actual") for m, sty, lab in [("naive", "s--", "naive"), ("model", "^-", "BED model"), ("ar1", "x:", "AR(1)")]: s = fc[(fc.horizon == h) & (fc.model == m)] ax.plot(pd.to_datetime(s.target_date), s.forecast, sty, ms=5, label=lab) ax.set_title(f"h = {h}: forecasts by target month (10 origins)") ax.set_ylabel("YoY growth, % (100·log)") ax.grid(alpha=.3) ax.legend(fontsize=8) fig.tight_layout() fig.savefig(os.path.join(OUT, "figure1.png"), dpi=130) # ---------------------------------------------------------------- results.json lag_issue = [o for o in origin_info if o["lag_months"] != 1] results = { "task_id": "T1-forecast", "model_name": "BED (base-effect drift)", "model_description": ( "Base-effect drift model on the RSAFS vintage dated V only. Log level forecast " "ln L_{last+h} = ln L_last + h*mu, with mu = mean monthly log growth over the last 12 published months " "= g_last/12 (percent). Hence forecast g_T = 100*ln(L_last / L_{T-12}) + h*g_last/12 for T = last(V)+h, " "where L_last and L_{T-12} are read from vintage V (T-12 <= last(V) always). It exploits that 12-h of the " "12 monthly growth terms in the target YoY are already published at the origin (known base effects), " "replacing only the h unknown months by their trailing 12-month average. No parameters are estimated; " "estimation window = last 13 levels of the vintage. PAYEMS and UMCSENT were not used. The specification " "was fixed before any out-of-sample error was computed; no tuning on the evaluation origins."), "n_origins": len(ORIGINS), "first_origin": ORIGINS[0], "last_origin": ORIGINS[-1], } for k in ["rmse_naive_h1", "mae_naive_h1", "rmse_naive_h3", "mae_naive_h3", "rmse_snaive_h1", "mae_snaive_h1", "rmse_snaive_h3", "mae_snaive_h3", "rmse_ar1_h1", "mae_ar1_h1", "rmse_ar1_h3", "mae_ar1_h3", "rmse_model_h1", "mae_model_h1", "rmse_model_h3", "mae_model_h3", "rel_improvement_rmse_h1", "rel_improvement_rmse_h3", "dm_stat_h1", "dm_pvalue_h1", "dm_stat_h3", "dm_pvalue_h3", "conclusion_code_h1", "conclusion_code_h3", "mean_abs_revision_h1_targets"]: results[k] = res[k] results.update({ "feasibility": "complete", "feasibility_reason": None, "data_issues": [ {"issue": "Irregular publication timing: last(V) is two months before the origin for " + ", ".join(f"{o['origin']} (last {o['last']})" for o in lag_issue) + "; vintage 2026-04-30 adds 2026-02 and 2026-03 at once, so 2026-03 is never an h=1 target (it is the h=1 target of no origin) and 2026-02 is reached from the 2026-03-31 vintage.", "handling": "Targets defined as last(V)+h from the vintage content, per the task rule, not from the calendar."}, {"issue": "Benchmark/annual revision of RSAFS levels: the current vintage is ~4% lower in level and its YoY growth is " "systematically lower than first release; the base level used in first_release jumps between targets 2026-07 and 2026-08 " "(rebasing between the 2026-08-31 and 2026-09-30 vintages), affecting h=3 targets 2026-08 and 2026-09.", "handling": "No adjustment: actuals are first-release values by definition; the revision is quantified in mean_abs_revision_h1_targets and discussed."}, {"issue": "UMCSENT has 210 empty values (all before 1978, quarterly survey era) and is a current (not real-time) vintage.", "handling": "Not used by the model; would have been usable under the stated no-revision rule."}, {"issue": "PAYEMS vintages sometimes extend beyond last(V) of RSAFS (e.g. 2025-12-31 vintage ends 2025-11 while RSAFS ends 2025-10).", "handling": "Not used by the model."}, {"issue": "first_release.csv dates the first release of 2025-08 (last month of the 2025-10-31 vintage) to the unshipped 2025-09-30 vintage; its level in the 2025-10-31 vintage still equals level_first_release.", "handling": "None needed; checked for consistency."}, ], "n_input_rows": n_input_rows, "n_analysed": int(fc[fc.model == "model"].dropna(subset=["forecast", "actual_first_release"]).shape[0]), "random_seed": SEED, "supplementary": {"bootstrap_rmse_diff": boot, "origins": origin_info, "revisions_h1_targets": {T.strftime("%Y-%m"): r for T, r in zip(targets, revs)}, "ar1_note": "AR(1) re-estimated by OLS at each origin on the vintage, 1993-02..last(V)"}, }) with open(os.path.join(OUT, "results.json"), "w") as f: json.dump(results, f, indent=2) print(pd.DataFrame(tab).to_string(index=False)) for h in HORIZONS: print(f"h={h}: rel_impr={res[f'rel_improvement_rmse_h{h}']:.4f} DM={res[f'dm_stat_h{h}']:.3f} " f"p={res[f'dm_pvalue_h{h}']:.4f} -> {res[f'conclusion_code_h{h}']}") print("mean_abs_revision_h1_targets =", round(res["mean_abs_revision_h1_targets"], 4))