"""T1-forecast: real-time forecasting of US retail-sales YoY growth. 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. No randomness is used (random_seed = null).""" import json, os import numpy as np, pandas as pd from scipy.stats import norm, binomtest import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt D, O = "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"] DRIFT_WINDOW = 36 # months, fixed a priori n_input_rows = {} def read_series(path): df = pd.read_csv(path) n_input_rows[os.path.relpath(path, D)] = int(len(df)) s = df.set_index(pd.to_datetime(df.iloc[:, 0])).iloc[:, 1] s.index = s.index.to_period("M") return s.astype(float) def yoy(level): l = np.log(level) return (100 * (l - l.shift(12))).dropna() fr = pd.read_csv(os.path.join(D, "first_release.csv")) n_input_rows["first_release.csv"] = int(len(fr)) fr["tm"] = pd.to_datetime(fr.target_month).dt.to_period("M") actual = fr.set_index("tm")["yoy_first_release"] cur = read_series(os.path.join(D, "RSAFS_current.csv")) # used ONLY for estimand (d) # Predictor files read for audit row counts only (the model does not use them) for o in ORIGINS: read_series(os.path.join(D, "vintages", f"PAYEMS_vintage_{o}.csv")) read_series(os.path.join(D, "UMCSENT.csv")) rows, data_issues, design = [], [], [] for o in ORIGINS: L = read_series(os.path.join(D, "vintages", f"RSAFS_vintage_{o}.csv")) assert L.index.is_monotonic_increasing and not L.index.duplicated().any() and L.notna().all() last = L.index[-1] vmonth = pd.Period(o, "M") g = yoy(L) # AR(1) on g over 1993-02 <= m <= last y = g[(g.index >= pd.Period("1993-02", "M"))] x = g.shift(1).reindex(y.index) X = np.column_stack([np.ones(len(y)), x.values]); b = np.linalg.lstsq(X, y.values, rcond=None)[0] # model: random walk with drift on log level, base effect from vintage dl = np.log(L).diff().dropna() mu = dl.iloc[-DRIFT_WINDOW:].mean() design.append({"origin": o, "last_month": str(last), "delay_months": int((vmonth - last).n), "ar1_c": b[0], "ar1_phi": b[1], "drift_pct_per_month": 100 * mu}) for h in HORIZONS: T = last + h assert T > last and T not in L.index # target not yet published in vintage V f_naive = g[last] f_snaive = g[T - 12] f = g[last] for _ in range(h): f = b[0] + b[1] * f f_ar1 = f f_model = 100 * (np.log(L[last]) + h * mu - np.log(L[T - 12])) a = float(actual[T]) for m, fc in zip(MODELS, [f_naive, f_snaive, f_ar1, f_model]): rows.append({"origin_date": o, "target_date": T.to_timestamp().strftime("%Y-%m-%d"), "horizon": h, "forecast": round(float(fc), 10), "actual_first_release": a, "model": m, "vintage_used": o}) fc = pd.DataFrame(rows) fc.to_csv(os.path.join(O, "forecasts.csv"), index=False) fc["e"] = fc.forecast - fc.actual_first_release res = {"task_id": "T1-forecast", "model_name": "rw_drift_yoy", "model_description": ("Base-effect-aware random walk with drift on log levels. At origin V, using only the RSAFS vintage " "dated V: mu = mean of the last 36 monthly log changes of the level ending at last(V); the level forecast is " "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} " "(the base month, already published) is taken from vintage V. The 36-month window and the specification were fixed " "before any forecast error was computed; no parameter is tuned on the evaluation sample. PAYEMS and UMCSENT are not " "used, so their availability alignment is not an issue; RSAFS_current and first_release are never used to forecast."), "n_origins": len(ORIGINS), "first_origin": ORIGINS[0], "last_origin": ORIGINS[-1]} table = [] for m in MODELS: for h in HORIZONS: e = fc[(fc.model == m) & (fc.horizon == h)].e.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": len(e)}) pd.DataFrame(table).to_csv(os.path.join(O, "table1.csv"), index=False) def dm(d, h): n = len(d); dc = d - d.mean() lrv = np.sum(dc * dc) / n for k in range(1, h): lrv += 2 * (1 - k / h) * np.sum(dc[k:] * dc[:-k]) / n stat = d.mean() / np.sqrt(lrv / n) return float(stat), float(2 * (1 - norm.cdf(abs(stat)))), float(lrv) dm_extra = {} for h in HORIZONS: en = fc[(fc.model == "naive") & (fc.horizon == h)].e.values em = fc[(fc.model == "model") & (fc.horizon == h)].e.values d = en ** 2 - em ** 2 s, p, lrv = dm(d, h) res[f"rel_improvement_rmse_h{h}"] = 1 - res[f"rmse_model_h{h}"] / res[f"rmse_naive_h{h}"] res[f"dm_stat_h{h}"], res[f"dm_pvalue_h{h}"] = s, p better = "model_better" if res[f"rmse_model_h{h}"] < res[f"rmse_naive_h{h}"] else "naive_better" res[f"conclusion_code_h{h}"] = better + ("_significant" if p < 0.05 else "_not_significant") dm_extra[f"h{h}"] = {"mean_loss_diff": float(d.mean()), "lrv": lrv, "n_origins_model_lower_sq_error": int((d > 0).sum()), "sign_test_pvalue_two_sided": float(binomtest(int((d > 0).sum()), len(d), 0.5).pvalue)} # (d) revisions: h=1 targets gcur = yoy(cur) t1 = pd.to_datetime(fc[(fc.horizon == 1) & (fc.model == "naive")].target_date).dt.to_period("M") rev = [abs(gcur[t] - actual[t]) for t in t1] res["mean_abs_revision_h1_targets"] = float(np.mean(rev)) delays = {x["origin"]: x["delay_months"] for x in design} data_issues = [ {"issue": "Irregular RSAFS publication timing: at the six origins 2025-10-31 ... 2026-03-31 the last published month is two " "months before the origin month (e.g. last(2025-10-31) = 2025-08) instead of the usual one; the 2026-04-30 vintage " "adds two months at once (2026-02 and 2026-03) and the last four origins have the usual one-month delay.", "handling": "Targets defined as last(V)+h per TASK.md, not by calendar. Consequently month 2026-03 is never an h=1 target " "(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 " "the origin date (a 'nowcast/backcast' of an unpublished month). Per-origin last(V) and delay are in real_time_design."}, {"issue": "PAYEMS vintage 2025-10-31 ends in 2025-08 and vintage 2025-12-31 jumps to 2025-11 (missing employment releases); " "PAYEMS vintage 2026-02-28 contains a benchmark revision (level about 1 million lower).", "handling": "PAYEMS is not used by the model; noted for completeness."}, {"issue": "UMCSENT has 210 empty values, all before 1978 (survey not monthly).", "handling": "UMCSENT not used by the model; no imputation."}, {"issue": "RSAFS_current reflects a large benchmark revision (levels changed for 138 months since 2015-01, up to about 4.6%) relative " "to the 2026-07-31 vintage; first_release.csv levels for 2026-08 and 2026-09 come from post-benchmark vintages while " "earlier target months come from pre-benchmark vintages.", "handling": "Forecasts never use RSAFS_current; actuals are taken as given from first_release.csv (YoY is a within-vintage ratio, " "so level shifts largely cancel); the revision estimand (d) reflects these benchmark changes."}, {"issue": "Monthly vintages revise past levels (up to 0.4% between consecutive vintages in the evaluation period).", "handling": "Each forecast uses only its own vintage; actual is the first release, so errors include revision noise in the base month."}, ] res.update({"feasibility": "complete", "feasibility_reason": None, "data_issues": data_issues, "n_input_rows": dict(sorted(n_input_rows.items())), "n_analysed": int(((fc.model == "model")).sum()), "random_seed": None, "real_time_design": design, "dm_details": dm_extra, "abs_revision_h1_by_target": {str(t): float(r) for t, r in zip(t1, rev)}}) with open(os.path.join(O, "results.json"), "w") as f: json.dump(res, f, indent=2, default=float) # exclusion ledger n_fr = len(fr); used = sorted(set(pd.to_datetime(fc.target_date).dt.to_period("M"))) pd.DataFrame([ {"step": 1, "rule": "candidate forecast-actual pairs = 10 origins x 2 horizons x 4 models", "n_before": 80, "n_after": 80, "n_removed": 0, "rationale": "all targets last(V)+h have a first-release actual; nothing excluded"}, {"step": 2, "rule": "pairs for own model (evaluated for RMSE/MAE/DM)", "n_before": 80, "n_after": 20, "n_removed": 60, "rationale": "baseline rows (60) are evaluated separately in table1; n_analysed counts model pairs only"}, ]).to_csv(os.path.join(O, "exclusions.csv"), index=False) # figure fig, axes = plt.subplots(1, 2, figsize=(11, 4), sharey=True) for ax, h in zip(axes, HORIZONS): sub = fc[fc.horizon == h] a = fr[(fr.tm >= pd.Period("2025-06", "M"))] ax.plot(a.tm.dt.to_timestamp(), a.yoy_first_release, "k-o", ms=3, label="first-release actual") for m, st in [("naive", "s--"), ("model", "^-"), ("ar1", "x:")]: s = sub[sub.model == m] ax.plot(pd.to_datetime(s.target_date), s.forecast, st, label=m) ax.set_title(f"h = {h}"); ax.set_ylabel("YoY growth, %"); ax.tick_params(axis="x", rotation=45) axes[0].legend(fontsize=8) fig.suptitle("RSAFS YoY growth: real-time forecasts vs first release") fig.tight_layout(); fig.savefig(os.path.join(O, "figure1.png"), dpi=120) for k in ["rmse_naive_h1", "rmse_model_h1", "rmse_naive_h3", "rmse_model_h3", "dm_stat_h1", "dm_pvalue_h1", "dm_stat_h3", "dm_pvalue_h3", "conclusion_code_h1", "conclusion_code_h3", "mean_abs_revision_h1_targets"]: print(k, res[k])