Fix silent NaN forecast from zero-variance lag-subsample in NNS.ARMA

gists/timeseries/m4/.gitignore

@@ -0,0 +1,3 @@
+# Downloaded M4 data and transient benchmark outputs -- regenerated by the script.
+m4_data/
+results/

gists/timeseries/m4/_modonly_probe.py

@@ -0,0 +1,105 @@
+"""Local experimentation harness for NNS.ARMA on the M4 seasonal subsets.
+
+Self-contained: downloads the M4 data on first run (cached in ./m4_data), then
+scores NNS forecasts vs the Naive (last-value) baseline on a sample of series.
+
+Tweak the CONFIG block, run `python _modonly_probe.py`, compare OWA/sMAPE, then
+revert to m4_benchmark.py for the full committed setup.
+"""
+
+import warnings
+from urllib.request import urlretrieve
+
+import numpy as np
+import pandas as pd
+
+warnings.filterwarnings("ignore")
+import nns as NNS
+
+# ── CONFIG — edit and re-run ────────────────────────────────────────────────
+N_SERIES = 30          # series per subset (None = all)
+MODULO_FROM_M = True   # use M4 seasonality (Hourly=24) as nns_seas modulo; else MODULO below
+MODULO = {"Hourly": 24, "Daily": 7, "Weekly": 12}
+MOD_ONLY = True        # nns_seas mod_only
+USE_OPTIM = True       # True: nns_arma_optim (searches lin/nonlin/both); False: plain nns_arma
+TRAINING_SET = None    # int to pass an explicit training_set to the optimizer; None = default
+OBJ_FN = None          # e.g. lambda p, a: np.mean(np.abs(p - a)); None = optimizer default
+MAX_PERIODS = 25
+SUBSETS = {"Hourly": (24, 48), "Weekly": (1, 13), "Daily": (1, 14)}  # name -> (m, h)
+# ─────────────────────────────────────────────────────────────────────────────
+
+BASE = "https://raw.githubusercontent.com/Mcompetitions/M4-methods/master/Dataset"
+
+
+def fetch(freq: str) -> tuple[pd.DataFrame, pd.DataFrame]:
+    import os
+    os.makedirs("m4_data", exist_ok=True)
+    paths = {}
+    for split in ("Train", "Test"):
+        dest = f"m4_data/{freq}-{split.lower()}.csv"
+        if not os.path.exists(dest):
+            print(f"  downloading {freq}-{split.lower()}.csv ...")
+            urlretrieve(f"{BASE}/{split}/{freq}-{split.lower()}.csv", dest)
+        paths[split] = dest
+    return pd.read_csv(paths["Train"]), pd.read_csv(paths["Test"])
+
+
+def smape(a, f):
+    d = np.abs(a) + np.abs(f)
+    return 200 * np.mean(np.where(d == 0, 0, np.abs(a - f) / d))
+
+
+def mase(a, f, tr, m=1):
+    s = np.mean(np.abs(tr[m:] - tr[:-m])) or 1e-8
+    return np.mean(np.abs(a - f)) / s
+
+
+def periods_for(y, freq):
+    modulo = (SUBSETS[freq][0] if MODULO_FROM_M else MODULO[freq])
+    modulo = modulo if modulo and modulo > 1 else None
+    tn = int(0.8 * len(y))
+    lim = tn / min(4, max(3, round(tn / 100)))
+    p = np.asarray(NNS.nns_seas(y, modulo=modulo, mod_only=MOD_ONLY).get("periods", []), int)
+    p = np.unique(p[(p > 1) & (p < lim)])[:MAX_PERIODS]
+    return p if p.size else np.array([2])
+
+
+def forecast(y, h, freq):
+    sf = periods_for(y, freq)
+    kw = dict(variable=y, h=h, seasonal_factor=sf, negative_values=True, print_trace=False)
+    if OBJ_FN is not None:
+        kw.update(obj_fn=OBJ_FN, objective="min")
+    if USE_OPTIM:
+        if TRAINING_SET is not None:
+            kw["training_set"] = TRAINING_SET
+        return np.asarray(NNS.nns_arma_optim(**kw)["results"], float)
+    kw.pop("print_trace", None)
+    out = NNS.nns_arma(variable=y, h=h, seasonal_factor=sf, negative_values=True)
+    return np.asarray(out["results"] if isinstance(out, dict) else out, float)
+
+
+def run(freq):
+    m, h = SUBSETS[freq]
+    tr, te = fetch(freq)
+    n = len(tr) if N_SERIES is None else min(N_SERIES, len(tr))
+    ns, n2s, nm, n2m, fail = [], [], [], [], 0
+    for i in range(n):
+        y = tr.iloc[i, 1:].dropna().astype(float).values
+        test = te.iloc[i, 1:].dropna().astype(float).values[:h]
+        try:
+            fc = forecast(y, h, freq)
+        except Exception:
+            fc = np.repeat(y[-1], h); fail += 1
+        n2 = np.repeat(y[-1], h)
+        ns.append(smape(test, fc)); nm.append(mase(test, fc, y))
+        n2s.append(smape(test, n2)); n2m.append(mase(test, n2, y))
+    owa = 0.5 * (np.mean(ns) / np.mean(n2s) + np.mean(nm) / np.mean(n2m))
+    print(f"{freq:7s} n={n:<4} OWA={owa:.3f}  sMAPE={np.mean(ns):.2f}  MASE={np.mean(nm):.2f}  fail={fail}")
+
+
+if __name__ == "__main__":
+    print(f"config: optim={USE_OPTIM} mod_only={MOD_ONLY} modulo_from_m={MODULO_FROM_M} "
+          f"obj_fn={'custom' if OBJ_FN else 'default'} training_set={TRAINING_SET}")
+    for freq in SUBSETS:
+        run(freq)
+    print("\nOWA < 1.000 beats the Naive baseline; lower is better.")

gists/timeseries/m4/m4_benchmark.py

@@ -0,0 +1,234 @@
+"""NNS.ARMA on the M4 competition (Hourly + Weekly + Daily subsets).
+
+The M4 competition (Makridakis et al., 2018-2020) is the recognised benchmark for
+univariate forecasting: 100,000 series scored by OWA -- the average of sMAPE and
+MASE, each normalised by the competition's Naive2 baseline (so Naive2 == 1.000).
+
+This script runs the three *seasonal / longer-series* subsets, where a seasonal
+ARMA method is best suited:
+
+    Hourly   414 series   m=24   h=48
+    Weekly   359 series   m=1    h=13     (M4 treats Weekly as non-seasonal)
+    Daily   4227 series   m=1    h=14     (M4 treats Daily  as non-seasonal)
+
+For every series it discovers seasonal periods with nns_seas, forecasts h steps
+with NNS.ARMA.optim, and scores sMAPE / MASE / OWA against M4's own Naive2. No
+neural network, no training loop.
+
+Usage:
+    python m4_benchmark.py                 # all three subsets, all series
+    python m4_benchmark.py --freq Hourly   # one subset
+    python m4_benchmark.py --limit 50      # first 50 series per subset (smoke test)
+    python m4_benchmark.py --workers 4
+
+Results are checkpointed to results/m4_<freq>.csv (resumable) and summarised to
+results/m4_summary.csv.
+"""
+
+from __future__ import annotations
+
+import argparse
+import os
+import warnings
+from concurrent.futures import ProcessPoolExecutor, as_completed
+from urllib.request import urlretrieve
+
+import numpy as np
+import pandas as pd
+
+warnings.filterwarnings("ignore")
+
+import nns as NNS
+
+BASE = "https://raw.githubusercontent.com/Mcompetitions/M4-methods/master/Dataset"
+CACHE = "m4_data"
+RESULTS = "results"
+
+# M4 official seasonality (m) and forecast horizon (h) per frequency.
+SUBSETS = {
+    "Hourly": {"m": 24, "h": 48},
+    "Weekly": {"m": 1, "h": 13},
+    "Daily": {"m": 1, "h": 14},
+}
+
+
+# ── Data ──────────────────────────────────────────────────────────────────────
+
+def _download(rel: str) -> str:
+    os.makedirs(CACHE, exist_ok=True)
+    dest = os.path.join(CACHE, os.path.basename(rel))
+    if not os.path.exists(dest):
+        urlretrieve(f"{BASE}/{rel}", dest)
+    return dest
+
+
+def load_subset(freq: str) -> tuple[list[str], list[np.ndarray], np.ndarray]:
+    train = pd.read_csv(_download(f"Train/{freq}-train.csv"))
+    test = pd.read_csv(_download(f"Test/{freq}-test.csv"))
+    ids = train.iloc[:, 0].astype(str).to_numpy()
+    y_train = [row[~np.isnan(row)] for row in train.iloc[:, 1:].to_numpy(dtype=float)]
+    y_test = test.iloc[:, 1:].to_numpy(dtype=float)
+    return ids, y_train, y_test
+
+
+# ── Metrics ─────────────────────────────────────────────────────────────────
+
+def smape(actual: np.ndarray, forecast: np.ndarray) -> float:
+    denom = np.abs(actual) + np.abs(forecast)
+    diff = np.abs(actual - forecast)
+    return float(200.0 * np.mean(np.where(denom == 0, 0.0, diff / denom)))
+
+
+def mase(actual: np.ndarray, forecast: np.ndarray, train: np.ndarray, m: int) -> float:
+    m = m if len(train) > m else 1
+    scale = np.mean(np.abs(train[m:] - train[:-m]))
+    if scale == 0:
+        scale = 1e-8
+    return float(np.mean(np.abs(actual - forecast)) / scale)
+
+
+def _is_seasonal(train: np.ndarray, m: int) -> bool:
+    """M4's 90% autocorrelation seasonality test at lag m."""
+    n = len(train)
+    if m <= 1 or n < 3 * m:
+        return False
+    x = train - train.mean()
+    acf = np.array([np.sum(x[k:] * x[:-k]) / np.sum(x**2) for k in range(1, m + 1)])
+    limit = 1.645 * np.sqrt((1.0 + 2.0 * np.sum(acf[:-1] ** 2)) / n)
+    return bool(np.abs(acf[-1]) > limit)
+
+
+def naive2_forecast(train: np.ndarray, h: int, m: int) -> np.ndarray:
+    """M4 Naive2: seasonally-adjusted naive (==Naive1 when non-seasonal)."""
+    if not _is_seasonal(train, m):
+        return np.repeat(train[-1], h)
+    # Multiplicative seasonal indices via ratio-to-moving-average.
+    n = len(train)
+    ma = pd.Series(train).rolling(m, center=True).mean()
+    if m % 2 == 0:
+        ma = ma.rolling(2).mean().shift(-1)
+    ratio = train / ma.to_numpy()
+    idx = np.array([np.nanmean(ratio[i::m]) for i in range(m)])
+    idx *= m / idx.sum()
+    seas_train = np.tile(idx, n // m + 1)[:n]
+    deseason = train / seas_train
+    fut_seas = np.tile(idx, h // m + 1)[:h]
+    return deseason[-1] * fut_seas
+
+
+# ── Per-series NNS forecast ──────────────────────────────────────────────────
+
+def nns_forecast(train: np.ndarray, h: int, m: int) -> np.ndarray:
+    # Use M4's declared seasonality as the nns_seas modulo where it exists
+    # (Hourly=24); leave it to auto-detect on the non-seasonal subsets.
+    modulo = m if m > 1 else None
+    seas = NNS.nns_seas(train, modulo=modulo, plot=False)
+    periods = np.asarray(seas.get("periods", []), dtype=np.int64)
+    # Cap to the same limit nns_arma_optim enforces (period < train_n / denom).
+    train_n = int(0.8 * len(train))
+    limit = train_n / min(4, max(3, round(train_n / 100)))
+    periods = np.unique(periods[(periods > 1) & (periods < limit)])[:25]
+    if periods.size == 0:
+        periods = np.array([m if 1 < m < limit else 2], dtype=np.int64)
+    fit = NNS.nns_arma_optim(  # default objective; let it search lin/nonlin/both
+        variable=train,
+        h=h,
+        seasonal_factor=periods,
+        negative_values=True,
+        print_trace=False,
+        plot=False,
+    )
+    return np.asarray(fit["results"], dtype=np.float64)
+
+
+def _score_one(args: tuple) -> dict:
+    sid, train, test, m = args
+    h = len(test)
+    try:
+        fc = nns_forecast(train, h, m)
+        ok = True
+    except Exception:  # noqa: BLE001 -- a failed series falls back to Naive2
+        fc = naive2_forecast(train, h, m)
+        ok = False
+    n2 = naive2_forecast(train, h, m)
+    return {
+        "id": sid,
+        "nns_smape": smape(test, fc),
+        "nns_mase": mase(test, fc, train, m),
+        "naive2_smape": smape(test, n2),
+        "naive2_mase": mase(test, n2, train, m),
+        "nns_ok": ok,
+    }
+
+
+# ── Run one frequency ────────────────────────────────────────────────────────
+
+def run_freq(freq: str, limit: int | None, workers: int) -> pd.DataFrame:
+    m, h = SUBSETS[freq]["m"], SUBSETS[freq]["h"]
+    ids, y_train, y_test = load_subset(freq)
+    if limit:
+        ids, y_train, y_test = ids[:limit], y_train[:limit], y_test[:limit]
+
+    os.makedirs(RESULTS, exist_ok=True)
+    ckpt = os.path.join(RESULTS, f"m4_{freq}.csv")
+    done = set(pd.read_csv(ckpt)["id"]) if os.path.exists(ckpt) else set()
+
+    tasks = [
+        (sid, y_train[i], y_test[i][: h], m)
+        for i, sid in enumerate(ids)
+        if sid not in done
+    ]
+    print(f"[{freq}] {len(tasks)} series to run ({len(done)} cached), m={m}, h={h}")
+
+    rows: list[dict] = []
+    with ProcessPoolExecutor(max_workers=workers) as ex:
+        futures = {ex.submit(_score_one, t): t[0] for t in tasks}
+        for n_done, fut in enumerate(as_completed(futures), 1):
+            rows.append(fut.result())
+            if n_done % 50 == 0 or n_done == len(tasks):
+                pd.DataFrame(rows).to_csv(
+                    ckpt, mode="a", header=not os.path.exists(ckpt), index=False
+                )
+                print(f"  [{freq}] {n_done}/{len(tasks)}")
+                rows = []
+    if rows:
+        pd.DataFrame(rows).to_csv(ckpt, mode="a", header=not os.path.exists(ckpt), index=False)
+    return pd.read_csv(ckpt)
+
+
+def summarise(freq: str, df: pd.DataFrame) -> dict:
+    nns_s, nns_m = df["nns_smape"].mean(), df["nns_mase"].mean()
+    n2_s, n2_m = df["naive2_smape"].mean(), df["naive2_mase"].mean()
+    owa = 0.5 * (nns_s / n2_s + nns_m / n2_m)
+    return {
+        "freq": freq,
+        "series": len(df),
+        "nns_fail": int((~df["nns_ok"]).sum()),
+        "nns_sMAPE": round(nns_s, 3),
+        "nns_MASE": round(nns_m, 3),
+        "naive2_sMAPE": round(n2_s, 3),
+        "naive2_MASE": round(n2_m, 3),
+        "NNS_OWA": round(owa, 3),
+    }
+
+
+def main() -> None:
+    ap = argparse.ArgumentParser()
+    ap.add_argument("--freq", choices=list(SUBSETS), help="run a single subset")
+    ap.add_argument("--limit", type=int, default=None, help="first N series per subset")
+    ap.add_argument("--workers", type=int, default=max(1, (os.cpu_count() or 2)))
+    args = ap.parse_args()
+
+    freqs = [args.freq] if args.freq else list(SUBSETS)
+    summary = [summarise(f, run_freq(f, args.limit, args.workers)) for f in freqs]
+
+    os.makedirs(RESULTS, exist_ok=True)
+    out = pd.DataFrame(summary)
+    out.to_csv(os.path.join(RESULTS, "m4_summary.csv"), index=False)
+    print("\n=== M4 BENCHMARK (NNS.ARMA.optim) ===\n")
+    print(out.to_string(index=False))
+    print("\nOWA < 1.000 beats the Naive2 baseline; lower is better.")
+
+
+if __name__ == "__main__":
+    main()

src/nns/arma.py

@@ -792,7 +792,13 @@ def _numeric_seasonal_weights(
             np.float64(np.std(variable, ddof=1)) / np.float64(np.mean(variable))
         )
         relative = output / baseline_cv
-        seasonal_weighting = 1.0 / relative
+        # A perfectly stable lag-subsample in an otherwise-varying series has CV 0, so
+        # relative 0 -> 1/relative = Inf -> Inf/Inf = NaN in the normalisation below,
+        # poisoning the forecast. Floor relative so that (maximally seasonal) lag gets a
+        # large but finite weight. A fully constant series (baseline CV 0 -> relative
+        # NaN) is intentionally left to propagate NaN, matching the R reference.
+        floor = np.finfo(np.float64).eps
+        seasonal_weighting = 1.0 / np.maximum(relative, floor)
 
     observation_weighting = 1.0 / np.sqrt(lags.astype(np.float64))
     denom = float(np.sum(observation_weighting * seasonal_weighting))

tests/invariants/test_arma.py

@@ -33,6 +33,22 @@ def test_nns_arma_numeric_seasonal_dynamic_raises() -> None:
         nns_arma(variable, h=3, seasonal_factor=5, dynamic=True)
 
 
+def test_numeric_seasonal_weights_constant_subsample_is_finite() -> None:
+    # Regression: a lag-subsample with zero variance gives CV 0, so the seasonal
+    # weighting 1 / (CV / baseline) = 1 / 0 = Inf and the normalisation collapsed to
+    # Inf / Inf = NaN -- which produced an all-NaN forecast. Now guarded.
+    rng = np.random.default_rng(0)
+    variable = rng.uniform(50.0, 150.0, size=40)
+    variable[::-4] = 10.0  # every 4th point (from the end) identical -> std 0
+
+    weights = _numeric_seasonal_weights(variable, np.array([4], dtype=np.int64))
+    assert np.isfinite(weights).all()
+    np.testing.assert_allclose(weights.sum(), 1.0)
+
+    forecast = nns_arma(variable, h=8, seasonal_factor=4, method="lin")
+    assert np.isfinite(forecast).all()
+
+
 @pytest.mark.stochastic
 def test_nns_arma_pred_int_returns_interval_dict() -> None:
     variable = np.sin(np.arange(1, 40, dtype=np.float64) / 3.0) + 2.0