1"""T1-forecast: real-time forecasting of US retail-sales YoY growth.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, out/exclusions.csv.4No randomness is used (random_seed = null)."""5import json, os6import numpy as np, pandas as pd7from scipy.stats import norm, binomtest8import matplotlib9matplotlib.use("Agg")10import matplotlib.pyplot as plt1112D, O = "data", "out"13ORIGINS = ["2025-10-31", "2025-11-30", "2025-12-31", "2026-01-31", "2026-02-28",14 "2026-03-31", "2026-04-30", "2026-05-31", "2026-06-30", "2026-07-31"]15HORIZONS = [1, 3]16MODELS = ["naive", "snaive", "ar1", "model"]17DRIFT_WINDOW = 36 # months, fixed a priori18n_input_rows = {}1920def read_series(path):21 df = pd.read_csv(path)22 n_input_rows[os.path.relpath(path, D)] = int(len(df))23 s = df.set_index(pd.to_datetime(df.iloc[:, 0])).iloc[:, 1]24 s.index = s.index.to_period("M")25 return s.astype(float)2627def yoy(level):28 l = np.log(level)29 return (100 * (l - l.shift(12))).dropna()3031fr = pd.read_csv(os.path.join(D, "first_release.csv"))32n_input_rows["first_release.csv"] = int(len(fr))33fr["tm"] = pd.to_datetime(fr.target_month).dt.to_period("M")34actual = fr.set_index("tm")["yoy_first_release"]35cur = read_series(os.path.join(D, "RSAFS_current.csv")) # used ONLY for estimand (d)36# Predictor files read for audit row counts only (the model does not use them)37for o in ORIGINS:38 read_series(os.path.join(D, "vintages", f"PAYEMS_vintage_{o}.csv"))39read_series(os.path.join(D, "UMCSENT.csv"))4041rows, data_issues, design = [], [], []42for o in ORIGINS:43 L = read_series(os.path.join(D, "vintages", f"RSAFS_vintage_{o}.csv"))44 assert L.index.is_monotonic_increasing and not L.index.duplicated().any() and L.notna().all()45 last = L.index[-1]46 vmonth = pd.Period(o, "M")47 g = yoy(L)48 # AR(1) on g over 1993-02 <= m <= last49 y = g[(g.index >= pd.Period("1993-02", "M"))]50 x = g.shift(1).reindex(y.index)51 X = np.column_stack([np.ones(len(y)), x.values]); 52 b = np.linalg.lstsq(X, y.values, rcond=None)[0]53 # model: random walk with drift on log level, base effect from vintage54 dl = np.log(L).diff().dropna()55 mu = dl.iloc[-DRIFT_WINDOW:].mean()56 design.append({"origin": o, "last_month": str(last), "delay_months": int((vmonth - last).n),57 "ar1_c": b[0], "ar1_phi": b[1], "drift_pct_per_month": 100 * mu})58 for h in HORIZONS:59 T = last + h60 assert T > last and T not in L.index # target not yet published in vintage V61 f_naive = g[last]62 f_snaive = g[T - 12]63 f = g[last]64 for _ in range(h):65 f = b[0] + b[1] * f66 f_ar1 = f67 f_model = 100 * (np.log(L[last]) + h * mu - np.log(L[T - 12]))68 a = float(actual[T])69 for m, fc in zip(MODELS, [f_naive, f_snaive, f_ar1, f_model]):70 rows.append({"origin_date": o, "target_date": T.to_timestamp().strftime("%Y-%m-%d"),71 "horizon": h, "forecast": round(float(fc), 10), "actual_first_release": a,72 "model": m, "vintage_used": o})73fc = pd.DataFrame(rows)74fc.to_csv(os.path.join(O, "forecasts.csv"), index=False)75fc["e"] = fc.forecast - fc.actual_first_release7677res = {"task_id": "T1-forecast", "model_name": "rw_drift_yoy",78 "model_description": ("Base-effect-aware random walk with drift on log levels. At origin V, using only the RSAFS vintage "79 "dated V: mu = mean of the last 36 monthly log changes of the level ending at last(V); the level forecast is "80 "ln L_{last+h} = ln L_last + h*mu, and the YoY forecast is 100*(ln L_{last+h} - ln L_{T-12}), where L_{T-12} "81 "(the base month, already published) is taken from vintage V. The 36-month window and the specification were fixed "82 "before any forecast error was computed; no parameter is tuned on the evaluation sample. PAYEMS and UMCSENT are not "83 "used, so their availability alignment is not an issue; RSAFS_current and first_release are never used to forecast."),84 "n_origins": len(ORIGINS), "first_origin": ORIGINS[0], "last_origin": ORIGINS[-1]}85table = []86for m in MODELS:87 for h in HORIZONS:88 e = fc[(fc.model == m) & (fc.horizon == h)].e.values89 rmse, mae = float(np.sqrt(np.mean(e ** 2))), float(np.mean(np.abs(e)))90 res[f"rmse_{m}_h{h}"], res[f"mae_{m}_h{h}"] = rmse, mae91 table.append({"model": m, "horizon": h, "rmse": rmse, "mae": mae, "n": len(e)})92pd.DataFrame(table).to_csv(os.path.join(O, "table1.csv"), index=False)9394def dm(d, h):95 n = len(d); dc = d - d.mean()96 lrv = np.sum(dc * dc) / n97 for k in range(1, h):98 lrv += 2 * (1 - k / h) * np.sum(dc[k:] * dc[:-k]) / n99 stat = d.mean() / np.sqrt(lrv / n)100 return float(stat), float(2 * (1 - norm.cdf(abs(stat)))), float(lrv)101102dm_extra = {}103for h in HORIZONS:104 en = fc[(fc.model == "naive") & (fc.horizon == h)].e.values105 em = fc[(fc.model == "model") & (fc.horizon == h)].e.values106 d = en ** 2 - em ** 2107 s, p, lrv = dm(d, h)108 res[f"rel_improvement_rmse_h{h}"] = 1 - res[f"rmse_model_h{h}"] / res[f"rmse_naive_h{h}"]109 res[f"dm_stat_h{h}"], res[f"dm_pvalue_h{h}"] = s, p110 better = "model_better" if res[f"rmse_model_h{h}"] < res[f"rmse_naive_h{h}"] else "naive_better"111 res[f"conclusion_code_h{h}"] = better + ("_significant" if p < 0.05 else "_not_significant")112 dm_extra[f"h{h}"] = {"mean_loss_diff": float(d.mean()), "lrv": lrv,113 "n_origins_model_lower_sq_error": int((d > 0).sum()),114 "sign_test_pvalue_two_sided": float(binomtest(int((d > 0).sum()), len(d), 0.5).pvalue)}115116# (d) revisions: h=1 targets117gcur = yoy(cur)118t1 = pd.to_datetime(fc[(fc.horizon == 1) & (fc.model == "naive")].target_date).dt.to_period("M")119rev = [abs(gcur[t] - actual[t]) for t in t1]120res["mean_abs_revision_h1_targets"] = float(np.mean(rev))121122delays = {x["origin"]: x["delay_months"] for x in design}123data_issues = [124 {"issue": "Irregular RSAFS publication timing: at the six origins 2025-10-31 ... 2026-03-31 the last published month is two "125 "months before the origin month (e.g. last(2025-10-31) = 2025-08) instead of the usual one; the 2026-04-30 vintage "126 "adds two months at once (2026-02 and 2026-03) and the last four origins have the usual one-month delay.",127 "handling": "Targets defined as last(V)+h per TASK.md, not by calendar. Consequently month 2026-03 is never an h=1 target "128 "(h=1 targets: 2026-02 from origin 2026-03-31, then 2026-04 from origin 2026-04-30), and for the first six origins the h=1 target month lies entirely before "129 "the origin date (a 'nowcast/backcast' of an unpublished month). Per-origin last(V) and delay are in real_time_design."},130 {"issue": "PAYEMS vintage 2025-10-31 ends in 2025-08 and vintage 2025-12-31 jumps to 2025-11 (missing employment releases); "131 "PAYEMS vintage 2026-02-28 contains a benchmark revision (level about 1 million lower).",132 "handling": "PAYEMS is not used by the model; noted for completeness."},133 {"issue": "UMCSENT has 210 empty values, all before 1978 (survey not monthly).",134 "handling": "UMCSENT not used by the model; no imputation."},135 {"issue": "RSAFS_current reflects a large benchmark revision (levels changed for 138 months since 2015-01, up to about 4.6%) relative "136 "to the 2026-07-31 vintage; first_release.csv levels for 2026-08 and 2026-09 come from post-benchmark vintages while "137 "earlier target months come from pre-benchmark vintages.",138 "handling": "Forecasts never use RSAFS_current; actuals are taken as given from first_release.csv (YoY is a within-vintage ratio, "139 "so level shifts largely cancel); the revision estimand (d) reflects these benchmark changes."},140 {"issue": "Monthly vintages revise past levels (up to 0.4% between consecutive vintages in the evaluation period).",141 "handling": "Each forecast uses only its own vintage; actual is the first release, so errors include revision noise in the base month."},142]143res.update({"feasibility": "complete", "feasibility_reason": None, "data_issues": data_issues,144 "n_input_rows": dict(sorted(n_input_rows.items())),145 "n_analysed": int(((fc.model == "model")).sum()), "random_seed": None,146 "real_time_design": design, "dm_details": dm_extra,147 "abs_revision_h1_by_target": {str(t): float(r) for t, r in zip(t1, rev)}})148with open(os.path.join(O, "results.json"), "w") as f:149 json.dump(res, f, indent=2, default=float)150151# exclusion ledger152n_fr = len(fr); used = sorted(set(pd.to_datetime(fc.target_date).dt.to_period("M")))153pd.DataFrame([154 {"step": 1, "rule": "candidate forecast-actual pairs = 10 origins x 2 horizons x 4 models", "n_before": 80, "n_after": 80,155 "n_removed": 0, "rationale": "all targets last(V)+h have a first-release actual; nothing excluded"},156 {"step": 2, "rule": "pairs for own model (evaluated for RMSE/MAE/DM)", "n_before": 80, "n_after": 20, "n_removed": 60,157 "rationale": "baseline rows (60) are evaluated separately in table1; n_analysed counts model pairs only"},158]).to_csv(os.path.join(O, "exclusions.csv"), index=False)159160# figure161fig, axes = plt.subplots(1, 2, figsize=(11, 4), sharey=True)162for ax, h in zip(axes, HORIZONS):163 sub = fc[fc.horizon == h]164 a = fr[(fr.tm >= pd.Period("2025-06", "M"))]165 ax.plot(a.tm.dt.to_timestamp(), a.yoy_first_release, "k-o", ms=3, label="first-release actual")166 for m, st in [("naive", "s--"), ("model", "^-"), ("ar1", "x:")]:167 s = sub[sub.model == m]168 ax.plot(pd.to_datetime(s.target_date), s.forecast, st, label=m)169 ax.set_title(f"h = {h}"); ax.set_ylabel("YoY growth, %"); ax.tick_params(axis="x", rotation=45)170axes[0].legend(fontsize=8)171fig.suptitle("RSAFS YoY growth: real-time forecasts vs first release")172fig.tight_layout(); fig.savefig(os.path.join(O, "figure1.png"), dpi=120)173174for k in ["rmse_naive_h1", "rmse_model_h1", "rmse_naive_h3", "rmse_model_h3", "dm_stat_h1", "dm_pvalue_h1",175 "dm_stat_h3", "dm_pvalue_h3", "conclusion_code_h1", "conclusion_code_h3", "mean_abs_revision_h1_targets"]:176 print(k, res[k])177