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

[ipezygj]

What

backtest() reduces the per-window error scores with np.nanmean. When stride < forecast_horizon the evaluation windows overlap, so those scores are not independent — and the number of windows badly overstates how much the backtest actually pins down.

This adds darts.utils.statistics.effective_sample_size(), which says how many independent observations a correlated sample is worth, and a note on backtest()'s reduction parameter pointing at it.

values = model.backtest(series, forecast_horizon=12, stride=1, reduction=None, ...)
len(values)                     # 190 windows
effective_sample_size(values)   # ~10

Measured, on darts' own backtest()

Setup: a seasonal series of 400 points, LinearRegressionModel(lags=24), forecast_horizon=12, metric=mae, reduction=None — 190 windows. Ground truth for coverage is the expected per-window MAE under the data-generating process, pooled over 60 independent series (11 400 windows), so it does not come from any interval being tested.

The per-window scores are heavily autocorrelated:

lag	autocorrelation of per-window MAE
1	+0.85
2	+0.72
5	+0.38

A 95% interval built from them as if they were independent covers the truth 27.5% of the time:

interval over 200 independent series	coverage (nominal 95%)	mean width
naive, std / sqrt(190)	27.5%	0.19
moving-block bootstrap	49.5%	0.42
effective-sample-size corrected	79.0%	0.97

Two controls, to show the cause is the overlap and not the data:

run	windows	lag-1 ACF	naive coverage
overlapping (stride=1), autocorrelated noise	190	+0.91	27.3%
non-overlapping (stride=forecast_horizon)	16	+0.08	62.7%
overlapping, white-noise series	190	+0.89	31.3%

Setting stride = forecast_horizon removes the autocorrelation almost entirely, and a white-noise series still shows it — so it comes from the windows sharing data, not from the series being autocorrelated.

What this PR deliberately does not add

My first plan was a confidence-interval helper. The measurements above killed it: the best within-series correction I tried reaches 79% coverage against a nominal 95%. Shipping that as a confidence interval would be shipping a number that lies, in a library where people compare models on exactly these scores.

The reason is visible in the last table: even with the overlap removed, coverage only reaches 62.7%. The residual is variation between series realisations, and no amount of resampling inside a single series can see it. A backtest of one series cannot yield a valid interval for expected performance on a new series — that is a property of the setup, not of the estimator.

So this adds a diagnostic, not an interval, and the docstring says so in those words. effective_sample_size answers "how much is this backtest worth", which is answerable. It does not pretend to answer "what is the error bar on this model's future performance", which from one series is not.

Changes

• darts/utils/statistics.py: effective_sample_size(values, max_lag=None). Autocorrelation sum truncated at the first non-positive lag (Geyer's initial positive sequence rule). Uses statsmodels.tsa.stattools.acf, already imported in that module — no new dependency. It sits beside _bartlett_formula, which already corrects ACF confidence bands for the same reason, so the module is not learning a new idea here.
• darts/models/forecasting/forecasting_model.py: a .. note:: under reduction explaining the overlap and pointing at the new function.
• darts/tests/utils/test_effective_sample_size.py: 15 tests.

Testing

The two that carry the weight are known-value checks against formulas the implementation does not use:

• independent sample → n_eff ≈ n (n = 4000, within 15%);
• AR(1) with ρ ∈ {0.3, 0.6, 0.85} → n_eff ≈ n(1-ρ)/(1+ρ), the analytic value, within 25%. The implementation sums estimated autocorrelations instead, so this is two routes meeting rather than one restated.

Plus two integration tests on a real backtest: overlapping windows must come out worth less than a third of their count, and non-overlapping windows must keep more than half — the control from the table above, pinned.

I then ran the suite against seven deliberately broken versions to check the tests can fail:

mutation	result
correlation correction removed	5 failed
factor of 2 dropped from the sum	3 failed
no initial-positive truncation	3 failed
only lag 1 used	3 failed
truncation condition inverted	3 failed
constant-input guard removed	1 failed
lower floor removed	0 failed — see below

Two things that mutation run found, both fixed in this branch:

1. A test of mine meant to exercise negative autocorrelation used -ar1(rho=0.9). Negating a series does not change its autocorrelation, so that case was testing nothing. It now uses ar1(rho=-0.7).
2. I had a np.clip(..., 1.0, n) on the result. The upper bound is unreachable: the truncated sum runs over strictly positive autocorrelations, so the inflation factor is ≥ 1 by construction and the result can never exceed n. An unreachable guard reads as if it were doing work, so it is gone, replaced by a comment stating the invariant — and the test now asserts the invariant itself rather than the clamp.

The remaining floor (max(1.0, ...)) is not exercised by any test: I could not construct an input that reaches it (the most extreme small-n, strongly-correlated case I searched gave n_eff = 2.48 at n = 4). It is not provably unreachable the way the upper bound is, so I kept it, but I would rather say plainly that it is untested than leave you to find that out.

ruff check and ruff format --check are clean on all three files. pytest darts/tests/utils/test_effective_sample_size.py → 15 passed in 1.6s.

I will add the CHANGELOG.md entry once this has a PR number.

Notes

• Nothing existing changes behaviour; the function is additive and backtest() gains only documentation.
• If you would rather this lived somewhere other than utils/statistics.py, or under a different name, say the word.
• The scripts that produced the tables are throwaway but I am happy to attach them, or to fold the coverage experiment into the test suite as a slow-marked test if you would want it guarded against regression.

[dennisbader]

Hi @ipezygj and thanks for the PR. Usually, we prefer to first open an issue before a PR so we can discuss whether we will address the issue or not.

Having said that, I'm not sure I quite follow what this PR is trying to solve. The user is free to chose forecast_horizon, stride and also the reduction function (can even be a custom function). Every evaluated window comes from a new model forecast. Of course, if the windows are overlapping there might be auto-correlation, but it is not "wrong" per-se.

For example if our use case was to predict the next 7 days each day (with daily frequency), then forecast_horizon=7 and stride=1 would simulate this scenario.

[ipezygj]

Thanks for looking, and you're right — I made the claim too broadly. Overlapping windows are the correct simulation of "predict the next 7 days, each day", nothing about forecast_horizon=7, stride=1 is wrong, and this PR changes no behaviour. Also noted on issue-before-PR; I'll follow that next time.

The narrower claim I do want to defend is this: the number of windows backtest() returns is not the number of independent observations behind the number it returns, and nothing in darts currently says so. Here it is measured, with darts doing the work.

Setup: LinearRegressionModel(lags=24) fitted once on a held-out realisation and never refitted (retrain=False), stationary seasonal DGP, forecast_horizon=12. The estimand is that fixed model's expected MAE. For each of 800 independent realisations I call backtest(..., reduction=None) and take the mean of the window scores. The spread of that mean across the 800 realisations is the true standard error — no correction involved, just repetition.

	stride=1	stride=12 (control)
windows per backtest	590	50
naive SE = std / sqrt(windows)	0.0147	0.0506
true SE (across 800 realisations)	0.0535	0.0565
true / naive	×3.65	×1.12
actual coverage of a nominal 95% interval	45.0%	90.4%

Two things in that table:

1. 590 windows and 50 windows pin the answer down equally well (0.0535 vs 0.0565). The extra 540 windows are not extra evidence; they are the same evidence recounted. That is the part the window count hides.
2. Anyone forming an interval from the window count is off by a factor of 3.65, and their 95% interval covers 45% of the time.

Substituting effective_sample_size(scores) for the window count: SE 0.0500 against the true 0.0535, coverage 91.0%. Not exact — it is still a single-series correction — but the failure mode is gone. The control matters as much: at stride=12 the correction gives 42.3 against 50 windows, i.e. it is ≈1 when the windows are disjoint, so it isn't just shrinking everything it touches.

Where this touches the library rather than a user's own arithmetic: gridsearch() defaults to stride=1 (forecasting_model.py:1659), and in expanding-window mode it selects argmin over exactly these means and returns the winner's score at face value. The comparison it makes is judged against noise ~3.6× larger than the window count suggests. I have measured the standard error of one such mean, above; I have not measured the selection rate of a grid of equivalent candidates, so I'm stating that as the consequence, not as a result.

There's also precedent inside darts.utils.statistics: _bartlett_formula (statistics.py:117) exists precisely so that plot_acf's confidence band accounts for autocorrelation instead of assuming independence. This is the same correction in scalar form, for the sample the user gets back from backtest.

So the PR is not claiming backtest is wrong — it adds the one thing that lets a user tell 590 windows from 590 observations, plus a note on reduction pointing at it. If you'd prefer this discussed as an issue first, say the word and I'll open one and close the PR. And if a utility of this kind is simply out of scope for darts.utils.statistics, I'd rather hear that than leave it open.

Repro script (runs in ~70s)

import numpy as np
from darts import TimeSeries
from darts.metrics import mae
from darts.models import LinearRegressionModel
from darts.utils.statistics import effective_sample_size

N, LAGS, HORIZON, SEASON, R = 1000, 24, 12, 24, 800


def realisation(seed):
    rng = np.random.default_rng(seed)
    t = np.arange(N)
    noise = np.zeros(N)
    for i in range(1, N):
        noise[i] = 0.6 * noise[i - 1] + rng.normal(0, 1)
    return TimeSeries.from_values(10 * np.sin(2 * np.pi * t / SEASON) + noise)


def window_scores(model, series, stride):
    return np.asarray(
        model.backtest(
            series, start=0.4, forecast_horizon=HORIZON, stride=stride, metric=mae,
            reduction=None, retrain=False, last_points_only=False, verbose=False,
        ),
        dtype=float,
    ).ravel()


model = LinearRegressionModel(lags=LAGS)
model.fit(realisation(seed=0))

for stride in (1, HORIZON):
    means, naive, ess = [], [], []
    for r in range(R):
        s = window_scores(model, realisation(1000 + r), stride)
        sd = s.std(ddof=1)
        means.append(s.mean())
        naive.append(sd / np.sqrt(s.size))
        ess.append(sd / np.sqrt(effective_sample_size(s)))
    means, naive, ess = map(np.asarray, (means, naive, ess))
    truth, true_se = means.mean(), means.std(ddof=1)
    cover = lambda se: np.mean(np.abs(means - truth) <= 1.96 * se)
    print(f"stride={stride}: naive {naive.mean():.4f} | ESS {ess.mean():.4f} | "
          f"true {true_se:.4f} | coverage naive {cover(naive):.1%} ESS {cover(ess):.1%}")