1"""T1-forecast: real-time forecasting of US retail sales YoY growth (RSAFS).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.png4(and the auxiliary diagnostics file out/state/aux_diagnostics.json).5Deterministic: no random numbers are used (random_seed = null)."""6import json, os7import numpy as np8import pandas as pd9from scipy import stats10import matplotlib11matplotlib.use("Agg")12import matplotlib.pyplot as plt1314DATA, OUT = "data", "out"15ORIGINS = ["2025-10-31", "2025-11-30", "2025-12-31", "2026-01-31", "2026-02-28",16 "2026-03-31", "2026-04-30", "2026-05-31", "2026-06-30", "2026-07-31"]17HORIZONS = [1, 3]18MODELS = ["naive", "snaive", "ar1", "model"]19MODEL_NAME = "rw_drift_base_effect"20MODEL_DESCRIPTION = (21 "Random walk with drift for the log level of RSAFS, translated into YoY growth so that the known "22 "'base effects' (the monthly changes that drop out of the 12-month window) are used. At origin V, "23 "with x_m = ln L_m from the RSAFS vintage dated V and l = last(V): drift d = (x_l - x_{l-12})/12 "24 "(mean monthly log change over the last 12 published months, i.e. a rolling 12-month estimation window); "25 "forecast of the log level x_{l+h} = x_l + h*d; forecast YoY g_{l+h} = 100*(x_l + h*d - x_{l+h-12}), "26 "where x_{l+h-12} (h<=3) is already published in vintage V. Uses only the RSAFS vintage dated V; "27 "PAYEMS and UMCSENT are not used. The specification (including the 12-month drift window) was fixed "28 "before any forecast error was computed; no parameter was tuned on the evaluation targets.")293031def read_series(path):32 d = pd.read_csv(path)33 s = pd.Series(d.iloc[:, 1].values.astype(float), index=pd.to_datetime(d.iloc[:, 0]))34 return s, len(d)353637def yoy(level):38 x = np.log(level)39 return 100.0 * (x - x.shift(12, freq="MS")).reindex(level.index)404142n_input_rows = {}43fr = pd.read_csv(os.path.join(DATA, "first_release.csv"), parse_dates=["target_month"])44n_input_rows["first_release.csv"] = len(fr)45actual = fr.set_index("target_month")["yoy_first_release"]4647rows = []48last_months = {}49for V in ORIGINS:50 lv, n = read_series(os.path.join(DATA, "vintages", f"RSAFS_vintage_{V}.csv"))51 n_input_rows[f"vintages/RSAFS_vintage_{V}.csv"] = n52 # PAYEMS vintage is read only for the row count/audit; not used in any forecast.53 _, npay = read_series(os.path.join(DATA, "vintages", f"PAYEMS_vintage_{V}.csv"))54 n_input_rows[f"vintages/PAYEMS_vintage_{V}.csv"] = npay55 assert lv.index.max() <= pd.Timestamp(V), "vintage contains months after origin"56 assert not lv.isna().any()57 last = lv.index.max()58 last_months[V] = last59 g = yoy(lv)60 x = np.log(lv)61 # AR(1): OLS of g_m on const, g_{m-1}, 1993-02 <= m <= last62 gs = g.loc["1993-01-01":last]63 y = gs.loc["1993-02-01":].values64 ylag = gs.shift(1).loc["1993-02-01":].values65 X = np.column_stack([np.ones_like(ylag), ylag])66 a, b = np.linalg.lstsq(X, y, rcond=None)[0]67 d = (x[last] - x[last - pd.DateOffset(months=12)]) / 12.068 for h in HORIZONS:69 T = last + pd.DateOffset(months=h)70 T12 = T - pd.DateOffset(months=12)71 f = {}72 f["naive"] = g[last]73 f["snaive"] = g[T12]74 fa = g[last]75 for _ in range(h):76 fa = a + b * fa77 f["ar1"] = fa78 f["model"] = 100.0 * (x[last] + h * d - x[T12])79 for m in MODELS:80 rows.append(dict(origin_date=V, target_date=T.strftime("%Y-%m-%d"), horizon=h,81 forecast=float(f[m]), actual_first_release=float(actual[T]),82 model=m, vintage_used=V))8384fc = pd.DataFrame(rows)85fc = fc.sort_values(["horizon", "model", "origin_date"]).reset_index(drop=True)86fc["error"] = fc["actual_first_release"] - fc["forecast"]87fc.drop(columns="error").to_csv(os.path.join(OUT, "forecasts.csv"), index=False, float_format="%.6f")8889res = {"task_id": "T1-forecast", "model_name": MODEL_NAME, "model_description": MODEL_DESCRIPTION,90 "n_origins": len(ORIGINS), "first_origin": ORIGINS[0], "last_origin": ORIGINS[-1]}91tab = []92for m in MODELS:93 for h in HORIZONS:94 e = fc[(fc.model == m) & (fc.horizon == h)]["error"].values95 rmse, mae = float(np.sqrt(np.mean(e ** 2))), float(np.mean(np.abs(e)))96 res[f"rmse_{m}_h{h}"], res[f"mae_{m}_h{h}"] = rmse, mae97 tab.append(dict(model=m, horizon=h, rmse=rmse, mae=mae, n=len(e)))98pd.DataFrame(tab).to_csv(os.path.join(OUT, "table1.csv"), index=False, float_format="%.6f")99100101def dm_test(e_naive, e_model, h):102 dd = e_naive ** 2 - e_model ** 2103 n = len(dd)104 dbar = dd.mean()105 dc = dd - dbar106 lrv = np.sum(dc * dc) / n107 for j in range(1, h): # Bartlett weights 1 - j/h, truncation lag h-1108 w = 1.0 - j / h109 lrv += 2.0 * w * np.sum(dc[j:] * dc[:-j]) / n110 stat = dbar / np.sqrt(lrv / n)111 p = 2.0 * (1.0 - stats.norm.cdf(abs(stat)))112 return float(stat), float(p), float(lrv), float(dbar)113114115dm_extra = {}116for h in HORIZONS:117 en = fc[(fc.model == "naive") & (fc.horizon == h)].sort_values("origin_date")["error"].values118 em = fc[(fc.model == "model") & (fc.horizon == h)].sort_values("origin_date")["error"].values119 res[f"rel_improvement_rmse_h{h}"] = 1.0 - res[f"rmse_model_h{h}"] / res[f"rmse_naive_h{h}"]120 s, p, lrv, dbar = dm_test(en, em, h)121 assert lrv > 0122 res[f"dm_stat_h{h}"], res[f"dm_pvalue_h{h}"] = s, p123 nn = len(en)124 hln = np.sqrt((nn + 1 - 2 * h + h * (h - 1) / nn) / nn)125 dm_extra[h] = (lrv, dbar, s * hln, float(2 * (1 - stats.t.cdf(abs(s * hln), df=nn - 1))),126 [float(v) for v in (en ** 2 - em ** 2)], int(np.sum(em ** 2 < en ** 2)))127 better = "model" if res[f"rmse_model_h{h}"] < res[f"rmse_naive_h{h}"] else "naive"128 sig = "significant" if p < 0.05 else "not_significant"129 res[f"conclusion_code_h{h}"] = f"{better}_better_{sig}"130131# (d) revision: current-vintage growth vs first release, h=1 targets (RSAFS_current used ONLY here)132cur, ncur = read_series(os.path.join(DATA, "RSAFS_current.csv"))133n_input_rows["RSAFS_current.csv"] = ncur134gcur = yoy(cur)135t1 = fc[(fc.model == "naive") & (fc.horizon == 1)]["target_date"].map(pd.Timestamp)136rev = pd.DataFrame({"target": t1.values, "g_current": gcur.reindex(t1.values).values,137 "first_release": actual.reindex(t1.values).values})138rev["abs_rev"] = (rev.g_current - rev.first_release).abs()139res["mean_abs_revision_h1_targets"] = float(rev["abs_rev"].mean())140141res["feasibility"] = "complete"142res["feasibility_reason"] = None143lag = {V: (pd.Timestamp(V).to_period("M") - last_months[V].to_period("M")).n for V in ORIGINS}144res["data_issues"] = [145 {"issue": "Irregular RSAFS publication timing: last(V) is 2 months before the origin month for origins "146 + ", ".join(V for V in ORIGINS if lag[V] == 2)147 + " and 1 month before it for the others; the 2026-04-30 vintage added two months at once.",148 "handling": "Targets defined as last(V)+h exactly as TASK.md requires; each forecast uses only the vintage dated V (vintage_used = origin_date)."},149 {"issue": "First-release dating is at month-end granularity: 2025-09 first appears in the 2025-11-30 vintage, and 2026-02 and 2026-03 both first appear in the 2026-04-30 vintage (so the 2026-02 'first release' need not be the advance estimate).",150 "handling": "first_release.csv used as provided (task definition of the actual); documented in the report."},151 {"issue": "Large level revision: RSAFS levels in the current vintage are about 4.7-4.9% below those of the 2026-07-31 vintage for 2025-2026 months (first_release rows for 2026-08/09 already use the revised level base).",152 "handling": "YoY growth is always computed within a single vintage (levels never mixed across vintages); the revision enters only estimand (d)."},153 {"issue": "UMCSENT has 210 empty values, all before 1978 (quarterly survey era).",154 "handling": "UMCSENT is not used by the model; no imputation needed."},155 {"issue": "RSAFS_current.csv contains information not available at any forecast origin.",156 "handling": "Read only in the revision block for estimand (d), after all forecasts are produced."},157]158n_input_rows["UMCSENT.csv"] = len(pd.read_csv(os.path.join(DATA, "UMCSENT.csv"))) # audited, not used in forecasts159res["n_input_rows"] = dict(sorted(n_input_rows.items()))160res["n_analysed"] = int(len(fc[fc.model == "model"]))161res["random_seed"] = None162163with open(os.path.join(OUT, "results.json"), "w") as fh:164 json.dump(res, fh, indent=2)165166# auxiliary (non-required) diagnostics cited in the report167aux = {"dm_long_run_variance": {f"h{h}": dm_extra[h][0] for h in HORIZONS},168 "dm_mean_loss_diff": {f"h{h}": dm_extra[h][1] for h in HORIZONS},169 "sensitivity_dm_hln_stat": {f"h{h}": dm_extra[h][2] for h in HORIZONS},170 "sensitivity_dm_hln_pvalue_t_nminus1": {f"h{h}": dm_extra[h][3] for h in HORIZONS},171 "loss_diff_by_origin": {f"h{h}": dm_extra[h][4] for h in HORIZONS},172 "n_origins_model_sq_error_below_naive": {f"h{h}": dm_extra[h][5] for h in HORIZONS},173 "revision_by_target": rev.assign(target=rev.target.dt.strftime("%Y-%m")).round(6).to_dict("records"),174 "last_month_by_origin": {V: last_months[V].strftime("%Y-%m") for V in ORIGINS}}175os.makedirs(os.path.join(OUT, "state"), exist_ok=True)176with open(os.path.join(OUT, "state", "aux_diagnostics.json"), "w") as fh:177 json.dump(aux, fh, indent=2)178179# Figure180fig, axes = plt.subplots(2, 1, figsize=(8, 7), sharex=True)181for ax, h in zip(axes, HORIZONS):182 sub = fc[fc.horizon == h]183 a_ = sub[sub.model == "naive"].sort_values("target_date")184 ax.plot(pd.to_datetime(a_.target_date), a_.actual_first_release, "k-o", label="first release (actual)")185 for m, st in [("naive", "C0--s"), ("model", "C3-^"), ("ar1", "C2:x")]:186 s_ = sub[sub.model == m].sort_values("target_date")187 ax.plot(pd.to_datetime(s_.target_date), s_.forecast, st, label=m)188 ax.set_title(f"h = {h}: RSAFS YoY growth (%), by target month")189 ax.set_ylabel("percent"); ax.grid(alpha=.3)190axes[0].legend(fontsize=8)191fig.tight_layout()192fig.savefig(os.path.join(OUT, "figure1.png"), dpi=120)193194print(pd.DataFrame(tab).to_string(index=False))195for h in HORIZONS:196 print(f"h={h}: rel_impr={res[f'rel_improvement_rmse_h{h}']:.4f} DM={res[f'dm_stat_h{h}']:.3f} "197 f"p={res[f'dm_pvalue_h{h}']:.4f} -> {res[f'conclusion_code_h{h}']}")198print("mean_abs_revision_h1_targets", round(res["mean_abs_revision_h1_targets"], 6))199