Add effective_sample_size() and document that overlapping backtest windows are not independent

CHANGELOG.md

@@ -11,6 +11,7 @@ but cannot always guarantee backwards compatibility. Changes that may **break co
 
 **Improved**
 
+- Added `effective_sample_size()` to `darts.utils.statistics`, reporting how many independent observations an autocorrelated sample is worth. Motivated by `backtest()`: when `stride` is smaller than `forecast_horizon` the evaluation windows overlap, the per-window error scores are strongly autocorrelated, and the window count overstates how much the backtest pins down. The `reduction` parameter documentation now points at it. [#3176](https://github.com/unit8co/darts/pull/3176) by [ipezygj](https://github.com/ipezygj).
 - Added support for per-timestep (non-aggregated) encoder and decoder variable importances in `TFTExplainer`, exposed as `TimeSeries` via `TFTExplainabilityResult.get_encoder_importance_over_time()` and `get_decoder_importance_over_time()`. [#3170](https://github.com/unit8co/darts/pull/3170) by [exactml](https://github.com/exactml).
 - Calling `TFTModel.fit_from_dataset()` on a dataset that does not have future covariates now raises an informative exception. [#3149](https://github.com/unit8co/darts/pull/3149) by [YOON KIWOONG](https://github.com/kiwoongyoon).
 - 🔴 Percentage and range-based metrics (`ape`, `mape`, `sape`, `smape`, `wmape`, `ope`, `arre`, `marre`, `coefficient_of_variation`) no longer raise a hard `ValueError` when the denominator is exactly zero. A new `zero_division` parameter controls the behavior: [#3122](https://github.com/unit8co/darts/pull/3122) by [Mahimn](https://github.com/mahimn01).

darts/models/forecasting/forecasting_model.py

@@ -1412,6 +1412,14 @@ def backtest(
             value for each metric function.
             If explicitly set to `None`, the method will return a list of the individual error scores instead.
             Set to ``np.mean`` by default.
+
+            .. note::
+                When `stride` is smaller than `forecast_horizon` the evaluation windows overlap, and the
+                individual error scores are strongly autocorrelated as a result. The number of windows then
+                overstates how much the backtest actually pins down, so a standard error formed as
+                ``std / sqrt(number of windows)`` is far too small. Pass ``reduction=None`` and see
+                :func:`~darts.utils.statistics.effective_sample_size` for how many independent observations
+                the scores are worth.
         verbose
             Whether to print the progress.
         show_warnings

darts/tests/utils/test_effective_sample_size.py

@@ -0,0 +1,132 @@
+import numpy as np
+import pytest
+
+from darts import TimeSeries
+from darts.metrics import mae
+from darts.models import LinearRegressionModel
+from darts.utils.statistics import effective_sample_size
+
+
+def _ar1(n: int, rho: float, seed: int = 42) -> np.ndarray:
+    rng = np.random.default_rng(seed)
+    x = np.zeros(n)
+    for i in range(1, n):
+        x[i] = rho * x[i - 1] + rng.normal(0, 1)
+    return x
+
+
+class TestEffectiveSampleSize:
+    def test_independent_sample_is_worth_its_own_size(self):
+        """Known value: with no autocorrelation, n_eff must come back as n."""
+        n = 4000
+        x = np.random.default_rng(0).normal(0, 1, n)
+
+        assert effective_sample_size(x) == pytest.approx(n, rel=0.15)
+
+    @pytest.mark.parametrize("rho", [0.3, 0.6, 0.85])
+    def test_matches_the_analytic_ar1_value(self, rho):
+        """Known value from a different formula than the one implemented.
+
+        For an AR(1) process the effective size is ``n * (1 - rho) / (1 + rho)``.
+        The implementation instead sums the estimated autocorrelations, so
+        agreement here is two routes meeting, not one route restated.
+        """
+        n = 8000
+        x = _ar1(n, rho)
+        analytic = n * (1 - rho) / (1 + rho)
+
+        assert effective_sample_size(x) == pytest.approx(analytic, rel=0.25)
+
+    def test_never_exceeds_the_sample_size_and_never_drops_below_one(self):
+        rng = np.random.default_rng(3)
+        samples = (
+            rng.normal(0, 1, 500),
+            _ar1(500, 0.95),
+            # rho < 0 alternates sign, which drives the correction below 1 and
+            # would report *more* independent observations than were collected.
+            # Negating a positively correlated series would not do this: the
+            # autocorrelation of -x equals that of x.
+            _ar1(2000, -0.7),
+        )
+        for x in samples:
+            n_eff = effective_sample_size(x)
+            assert 1.0 <= n_eff <= len(x)
+
+    def test_anti_correlated_sample_is_not_credited_with_extra_observations(self):
+        x = _ar1(2000, -0.7)
+
+        assert effective_sample_size(x) == pytest.approx(len(x))
+
+    def test_result_does_not_drift_when_max_lag_is_raised(self):
+        """Long lags carry estimation noise, not information.
+
+        Summing them instead of stopping at the first non-positive
+        autocorrelation makes the answer depend on where the sum was cut off.
+        """
+        x = _ar1(3000, 0.6)
+        at_default = effective_sample_size(x)
+        at_high_lag = effective_sample_size(x, max_lag=400)
+
+        assert at_high_lag == pytest.approx(at_default, rel=0.05)
+
+    def test_more_correlation_means_fewer_independent_observations(self):
+        sizes = [effective_sample_size(_ar1(4000, rho)) for rho in (0.0, 0.5, 0.9)]
+
+        assert sizes[0] > sizes[1] > sizes[2]
+
+    def test_constant_input_is_rejected_rather_than_silently_scored(self):
+        with pytest.raises(ValueError):
+            effective_sample_size(np.full(100, 2.5))
+
+    @pytest.mark.parametrize("values", [[1.0], [], [1.0, np.nan, 2.0]])
+    def test_unusable_input_raises(self, values):
+        with pytest.raises(ValueError):
+            effective_sample_size(np.asarray(values, dtype=float))
+
+    def test_accepts_a_plain_sequence(self):
+        x = list(_ar1(600, 0.5))
+
+        assert effective_sample_size(x) == pytest.approx(
+            effective_sample_size(np.asarray(x))
+        )
+
+
+class TestEffectiveSampleSizeOnBacktestWindows:
+    """The reason this function is in darts: overlapping backtest windows."""
+
+    @staticmethod
+    def _window_metrics(stride: int) -> np.ndarray:
+        n = 400
+        t = np.arange(n)
+        rng = np.random.default_rng(7)
+        series = TimeSeries.from_values(
+            10 * np.sin(2 * np.pi * t / 24) + 0.02 * t + rng.normal(0, 1, n)
+        )
+        model = LinearRegressionModel(lags=24)
+        model.fit(series[: n // 2])
+        return np.asarray(
+            model.backtest(
+                series,
+                start=0.5,
+                forecast_horizon=12,
+                stride=stride,
+                metric=mae,
+                reduction=None,
+                retrain=False,
+                last_points_only=False,
+                verbose=False,
+            ),
+            dtype=float,
+        ).ravel()
+
+    def test_overlapping_windows_are_worth_far_fewer_observations(self):
+        values = self._window_metrics(stride=1)
+
+        assert effective_sample_size(values) < len(values) / 3
+
+    def test_non_overlapping_windows_keep_most_of_their_size(self):
+        """Control: the shrinkage above is caused by the overlap, not by the
+        metric or the series."""
+        values = self._window_metrics(stride=12)
+
+        assert effective_sample_size(values) > len(values) / 2

darts/utils/statistics.py

@@ -1230,3 +1230,86 @@ def plot_tolerance_curve(
     axis.set_title("Tolerance Curve")
     axis.grid(True, alpha=0.3)
     axis.legend()
+
+
+def effective_sample_size(
+    values: Sequence[float] | np.ndarray,
+    max_lag: int | None = None,
+) -> float:
+    r"""
+    Number of independent observations a correlated sample is worth.
+
+    A sample of `n` autocorrelated values carries less information than `n`
+    independent ones, so a standard error formed as ``std / sqrt(n)`` is too
+    small. This returns
+
+    .. math:: n_{eff} = \frac{n}{1 + 2 \sum_{k=1}^{K} \rho_k}
+
+    where :math:`\rho_k` is the sample autocorrelation at lag `k`, truncated at
+    the first non-positive :math:`\rho_k` (the initial positive sequence rule of
+    Geyer, 1992). The result is clipped to ``[1, n]``.
+
+    The motivating case in Darts is
+    :func:`~darts.models.forecasting.forecasting_model.ForecastingModel.backtest`
+    with ``reduction=None``: when ``stride < forecast_horizon`` the evaluation
+    windows overlap, successive window metrics are strongly autocorrelated, and
+    the number of windows badly overstates how much the backtest pins down.
+
+    Note that a small effective size is a warning, not something to divide by
+    and move on. A backtest of one series cannot, by any correction, deliver a
+    valid interval for expected performance on a *new* series: the variation
+    between series realisations is not visible inside one of them.
+
+    Parameters
+    ----------
+    values
+        The (possibly autocorrelated) sample, e.g. the per-window metrics from
+        ``backtest(..., reduction=None)``. Must be finite and not constant.
+    max_lag
+        Highest lag considered. Defaults to ``min(len(values) - 1, 10 * log10(n))``,
+        the statsmodels convention.
+
+    Returns
+    -------
+    float
+        The effective sample size, between 1 and ``len(values)``.
+
+    Examples
+    --------
+    >>> import numpy as np
+    >>> from darts.utils.statistics import effective_sample_size
+    >>> rng = np.random.default_rng(0)
+    >>> round(effective_sample_size(rng.normal(size=2000)))  # doctest: +SKIP
+    2013
+    """
+    array = np.asarray(values, dtype=float).ravel()
+    n = array.size
+
+    if n < 3:
+        raise_log(
+            ValueError(f"`values` must contain at least 3 observations, received {n}.")
+        )
+    if not np.all(np.isfinite(array)):
+        raise_log(ValueError("`values` must be finite; found NaN or infinity."))
+    if np.all(array == array[0]):
+        raise_log(
+            ValueError("`values` is constant, so its autocorrelation is undefined.")
+        )
+
+    if max_lag is None:
+        max_lag = int(min(n - 1, 10 * math.log10(n)))
+    max_lag = max(1, min(max_lag, n - 1))
+
+    rho = acf(array, nlags=max_lag, fft=True)[1:]
+    # Initial positive sequence: stop at the first lag whose correlation is not
+    # positive, rather than summing noise from long lags.
+    non_positive = np.flatnonzero(rho <= 0)
+    cutoff = non_positive[0] if non_positive.size else rho.size
+    inflation = 1.0 + 2.0 * float(np.sum(rho[:cutoff]))
+
+    # The truncated sum runs over strictly positive autocorrelations, so
+    # `inflation >= 1` by construction and the result can never exceed `n`.
+    # No upper clamp is applied: one would be unreachable, and an unreachable
+    # guard reads as if it were doing work. The lower floor is not provable in
+    # the same way for very small `n`, so it stays.
+    return float(max(1.0, n / inflation))