Fix Bartlett significance threshold in seasonality detection

CHANGELOG.md

@@ -13,6 +13,7 @@ but cannot always guarantee backwards compatibility. Changes that may **break co
 
 **Fixed**
 
+- Fixed `check_seasonality()` using an incorrect significance threshold and omitting the preceding lag from Bartlett's standard error, which caused false seasonality detections. This can change automatically inferred periods in `Theta` and `FourTheta`; set `seasonality_period` when the period is known. [#3286](https://github.com/unit8co/darts/issues/3286) by [Aman Agnihotri](https://github.com/its-amann).
 - Fixed `CatBoostModel` with a `"gaussian"` likelihood reporting the variance instead of the standard deviation when predicting with `predict_likelihood_parameters=True`. The `*_sigma` component now holds the standard deviation, consistent with the torch `GaussianLikelihood`. Sampled forecasts are unaffected. [#3302](https://github.com/unit8co/darts/pull/3302) by [Miral Amin](https://github.com/aminmiral).
 - Fixed compatibility with pandas 3.1: fitting and predicting models on series with a daily frequency failed with a `TypeError`, and `MIDAS` accepted a series whose frequency was already the target low frequency. [#3289](https://github.com/unit8co/darts/pull/3289) by [Nicolò Angileri](https://github.com/nicoloangileri).
 

darts/tests/models/forecasting/test_4theta.py

@@ -4,6 +4,7 @@
 import numpy as np
 import pytest
 
+from darts import TimeSeries
 from darts.metrics import mape
 from darts.models import FourTheta, Theta
 from darts.utils.timeseries_generation import linear_timeseries as lt
@@ -15,6 +16,18 @@
 
 
 class TestFourTheta:
+    @pytest.mark.parametrize("model_cls", [Theta, FourTheta])
+    @pytest.mark.parametrize("period", [0, 12])
+    def test_inferred_seasonality(self, model_cls, period):
+        values = 10 + np.random.default_rng(0).normal(scale=0.1, size=120)
+        if period:
+            values += np.cos(2 * np.pi * np.arange(len(values)) / period)
+        model = model_cls().fit(TimeSeries.from_values(values))
+
+        assert model.is_seasonal == (period > 0)
+        assert model.season_period == period
+        assert np.isfinite(model.predict(12).values()).all()
+
     def test_input(self):
         with pytest.raises(ValueError):
             FourTheta(model_mode=SeasonalityMode.ADDITIVE)

darts/tests/models/forecasting/test_local_forecasting_models.py

@@ -70,14 +70,15 @@
     (ExponentialSmoothing(), 5.4),
     (ARIMA(12, 2, 1), 5.2),
     (ARIMA(1, 1, 1), 24),
-    (Theta(), 11),
-    (Theta(1), 17),
-    (Theta(-1), 12),
-    (FourTheta(1), 17),
-    (FourTheta(-1), 12),
-    (FourTheta(trend_mode=TrendMode.EXPONENTIAL), 6.0),
-    (FourTheta(model_mode=ModelMode.MULTIPLICATIVE), 10),
-    (FourTheta(season_mode=SeasonalityMode.ADDITIVE), 12),
+    # Use the known annual period for the AirPassengers performance benchmark.
+    (Theta(seasonality_period=12), 11),
+    (Theta(1, seasonality_period=12), 17),
+    (Theta(-1, seasonality_period=12), 12),
+    (FourTheta(1, seasonality_period=12), 17),
+    (FourTheta(-1, seasonality_period=12), 12),
+    (FourTheta(trend_mode=TrendMode.EXPONENTIAL, seasonality_period=12), 6.0),
+    (FourTheta(model_mode=ModelMode.MULTIPLICATIVE, seasonality_period=12), 10),
+    (FourTheta(season_mode=SeasonalityMode.ADDITIVE, seasonality_period=12), 12),
     (FFT(trend="poly"), 13),
     (KalmanForecaster(dim_x=3), 20),
     (LinearRegressionModel(lags=12), 13),

darts/tests/utils/test_statistics.py

@@ -4,6 +4,7 @@
 import numpy as np
 import pandas as pd
 import pytest
+from statsmodels.tsa.stattools import acf
 
 from darts import TimeSeries
 from darts.datasets import AirPassengersDataset
@@ -46,6 +47,42 @@ def test_check_seasonality(self):
         with pytest.raises(AssertionError):
             check_seasonality(series.stack(series))
 
+    @pytest.mark.parametrize("m", [None, 4])
+    def test_check_seasonality_white_noise(self, m):
+        series = TimeSeries.from_values(np.random.default_rng(0).normal(size=100))
+
+        assert check_seasonality(series, m=m) == (False, 0)
+
+    @pytest.mark.parametrize("m", [None, 2])
+    @pytest.mark.parametrize(
+        "alpha,expected", [(0.01, (False, 0)), (0.05, (False, 0)), (0.1, (True, 2))]
+    )
+    def test_check_seasonality_bartlett_confidence_interval(self, m, alpha, expected):
+        values = np.random.default_rng(13).normal(size=40)
+        for i in range(1, len(values)):
+            values[i] -= 0.4 * values[i - 1]
+
+        r, confint = acf(values, nlags=24, fft=False, alpha=alpha)
+        # Lag 2 is a local peak; its interval must include the strong lag-1 correlation.
+        assert r[2] > max(r[1], r[3])
+        assert (confint[2, 0] > 0) == expected[0]
+        assert (
+            check_seasonality(TimeSeries.from_values(values), m=m, alpha=alpha)
+            == expected
+        )
+
+    @pytest.mark.parametrize("period", [2, 6, 12])
+    @pytest.mark.parametrize("fixed_period", [False, True])
+    def test_check_seasonality_noisy_periodic_series(self, period, fixed_period):
+        values = np.cos(2 * np.pi * np.arange(240) / period)
+        values += np.random.default_rng(0).normal(scale=0.1, size=len(values))
+        series = TimeSeries.from_values(values)
+
+        assert check_seasonality(series, m=period if fixed_period else None) == (
+            True,
+            period,
+        )
+
     def test_granger_causality(self):
         series_cause_1 = constant_timeseries(start=0, end=9999).stack(
             constant_timeseries(start=0, end=9999)

darts/utils/statistics.py

@@ -56,7 +56,8 @@ def check_seasonality(
     max_lag
         The maximal lag allowed in the ACF.
     alpha
-        The desired confidence level (default 5%).
+        The significance level for each tested lag (default 5%). No multiple-testing correction is applied
+        when `m` is None.
 
     Returns
     -------
@@ -100,13 +101,12 @@ def check_seasonality(
     # Remove r[0], the auto-correlation at lag order 0, that introduces bias.
     r = r[1:]
 
-    # The non-adjusted upper limit of the significance interval.
-    band_upper = r.mean() + norm.ppf(1 - alpha / 2) * r.var()
+    # Normal critical value; Bartlett's formula supplies the standard error.
+    band_upper = norm.ppf(1 - alpha / 2)
 
-    # Significance test, stops at first admissible value. The two '-1' below
-    # compensate for the index change due to the restriction of the original r to r[1:].
+    # Test each lag using all preceding nonzero-lag autocorrelations.
     for candidate in candidates:
-        stat = _bartlett_formula(r, candidate - 1, len(ts))
+        stat = _bartlett_formula(r, candidate, len(ts))
         if r[candidate - 1] > stat * band_upper:
             return True, candidate
     return False, 0
@@ -119,7 +119,7 @@ def _bartlett_formula(r: np.ndarray, m: int, length: int) -> float:
     Parameters
     ----------
     r
-        The array whose standard error is to be computed.
+        The autocorrelations at lags 1, 2, ... (excluding lag 0).
     m
         The order of the standard error.
     length