Source code for statsmodels.tsa.stl.mstl

"""
Author: Kishan Manani
License: BSD-3 Clause

An implementation of MSTL [1], an algorithm for time series decomposition when
there are multiple seasonal components.

This implementation has the following differences with the original algorithm:
- Missing data must be handled outside of this class.
- The algorithm proposed in the paper handles a case when there is no
seasonality. This implementation assumes that there is at least one seasonal
component.

[1] K. Bandura, R.J. Hyndman, and C. Bergmeir (2021)
MSTL: A Seasonal-Trend Decomposition Algorithm for Time Series with Multiple
Seasonal Patterns
https://arxiv.org/pdf/2107.13462.pdf
"""
from typing import Optional, Union
from collections.abc import Sequence
import warnings

import numpy as np
import pandas as pd
from scipy.stats import boxcox

from statsmodels.tools.typing import ArrayLike1D
from statsmodels.tsa.stl._stl import STL
from statsmodels.tsa.tsatools import freq_to_period


[docs] class MSTL: """ MSTL(endog, periods=None, windows=None, lmbda=None, iterate=2, stl_kwargs=None) Season-Trend decomposition using LOESS for multiple seasonalities. .. versionadded:: 0.14.0 Parameters ---------- endog : array_like Data to be decomposed. Must be squeezable to 1-d. periods : {int, array_like, None}, optional Periodicity of the seasonal components. If None and endog is a pandas Series or DataFrame, attempts to determine from endog. If endog is a ndarray, periods must be provided. windows : {int, array_like, None}, optional Length of the seasonal smoothers for each corresponding period. Must be an odd integer, and should normally be >= 7 (default). If None then default values determined using 7 + 4 * np.arange(1, n + 1, 1) where n is number of seasonal components. lmbda : {float, str, None}, optional The lambda parameter for the Box-Cox transform to be applied to `endog` prior to decomposition. If None, no transform is applied. If "auto", a value will be estimated that maximizes the log-likelihood function. iterate : int, optional Number of iterations to use to refine the seasonal component. stl_kwargs: dict, optional Arguments to pass to STL. See Also -------- statsmodels.tsa.seasonal.STL References ---------- .. [1] K. Bandura, R.J. Hyndman, and C. Bergmeir (2021) MSTL: A Seasonal-Trend Decomposition Algorithm for Time Series with Multiple Seasonal Patterns. arXiv preprint arXiv:2107.13462. Examples -------- Start by creating a toy dataset with hourly frequency and multiple seasonal components. >>> import numpy as np >>> import matplotlib.pyplot as plt >>> import pandas as pd >>> pd.plotting.register_matplotlib_converters() >>> np.random.seed(0) >>> t = np.arange(1, 1000) >>> trend = 0.0001 * t ** 2 + 100 >>> daily_seasonality = 5 * np.sin(2 * np.pi * t / 24) >>> weekly_seasonality = 10 * np.sin(2 * np.pi * t / (24 * 7)) >>> noise = np.random.randn(len(t)) >>> y = trend + daily_seasonality + weekly_seasonality + noise >>> index = pd.date_range(start='2000-01-01', periods=len(t), freq='h') >>> data = pd.DataFrame(data=y, index=index) Use MSTL to decompose the time series into two seasonal components with periods 24 (daily seasonality) and 24*7 (weekly seasonality). >>> from statsmodels.tsa.seasonal import MSTL >>> res = MSTL(data, periods=(24, 24*7)).fit() >>> res.plot() >>> plt.tight_layout() >>> plt.show() .. plot:: plots/mstl_plot.py """ def __init__( self, endog: ArrayLike1D, *, periods: Optional[Union[int, Sequence[int]]] = None, windows: Optional[Union[int, Sequence[int]]] = None, lmbda: Optional[Union[float, str]] = None, iterate: int = 2, stl_kwargs: Optional[dict[str, Union[int, bool, None]]] = None, ): self.endog = endog self._y = self._to_1d_array(endog) self.nobs = self._y.shape[0] self.lmbda = lmbda self.periods, self.windows = self._process_periods_and_windows( periods, windows ) self.iterate = iterate self._stl_kwargs = self._remove_overloaded_stl_kwargs( stl_kwargs if stl_kwargs else {} )
[docs] def fit(self): """ Estimate a trend component, multiple seasonal components, and a residual component. Returns ------- DecomposeResult Estimation results. """ num_seasons = len(self.periods) iterate = 1 if num_seasons == 1 else self.iterate # Box Cox if self.lmbda == "auto": y, lmbda = boxcox(self._y, lmbda=None) self.est_lmbda = lmbda elif self.lmbda: y = boxcox(self._y, lmbda=self.lmbda) else: y = self._y # Get STL fit params stl_inner_iter = self._stl_kwargs.pop("inner_iter", None) stl_outer_iter = self._stl_kwargs.pop("outer_iter", None) # Iterate over each seasonal component to extract seasonalities seasonal = np.zeros(shape=(num_seasons, self.nobs)) deseas = y for _ in range(iterate): for i in range(num_seasons): deseas = deseas + seasonal[i] res = STL( endog=deseas, period=self.periods[i], seasonal=self.windows[i], **self._stl_kwargs, ).fit(inner_iter=stl_inner_iter, outer_iter=stl_outer_iter) seasonal[i] = res.seasonal deseas = deseas - seasonal[i] seasonal = np.squeeze(seasonal.T) trend = res.trend rw = res.weights resid = deseas - trend # Return pandas if endog is pandas if isinstance(self.endog, (pd.Series, pd.DataFrame)): index = self.endog.index y = pd.Series(y, index=index, name="observed") trend = pd.Series(trend, index=index, name="trend") resid = pd.Series(resid, index=index, name="resid") rw = pd.Series(rw, index=index, name="robust_weight") cols = [f"seasonal_{period}" for period in self.periods] if seasonal.ndim == 1: seasonal = pd.Series(seasonal, index=index, name="seasonal") else: seasonal = pd.DataFrame(seasonal, index=index, columns=cols) # Avoid circular imports from statsmodels.tsa.seasonal import DecomposeResult return DecomposeResult(y, seasonal, trend, resid, rw)
def __str__(self): return ( "MSTL(endog," f" periods={self.periods}," f" windows={self.windows}," f" lmbda={self.lmbda}," f" iterate={self.iterate})" ) def _process_periods_and_windows( self, periods: Union[int, Sequence[int], None], windows: Union[int, Sequence[int], None], ) -> tuple[Sequence[int], Sequence[int]]: periods = self._process_periods(periods) if windows: windows = self._process_windows(windows, num_seasons=len(periods)) periods, windows = self._sort_periods_and_windows(periods, windows) else: windows = self._process_windows(windows, num_seasons=len(periods)) periods = sorted(periods) if len(periods) != len(windows): raise ValueError("Periods and windows must have same length") # Remove long periods from decomposition if any(period >= self.nobs / 2 for period in periods): warnings.warn( "A period(s) is larger than half the length of time series." " Removing these period(s).", UserWarning ) periods = tuple( period for period in periods if period < self.nobs / 2 ) windows = windows[: len(periods)] return periods, windows def _process_periods( self, periods: Union[int, Sequence[int], None] ) -> Sequence[int]: if periods is None: periods = (self._infer_period(),) elif isinstance(periods, int): periods = (periods,) else: pass return periods def _process_windows( self, windows: Union[int, Sequence[int], None], num_seasons: int, ) -> Sequence[int]: if windows is None: windows = self._default_seasonal_windows(num_seasons) elif isinstance(windows, int): windows = (windows,) else: pass return windows def _infer_period(self) -> int: freq = None if isinstance(self.endog, (pd.Series, pd.DataFrame)): freq = getattr(self.endog.index, "inferred_freq", None) if freq is None: raise ValueError("Unable to determine period from endog") period = freq_to_period(freq) return period @staticmethod def _sort_periods_and_windows( periods, windows ) -> tuple[Sequence[int], Sequence[int]]: if len(periods) != len(windows): raise ValueError("Periods and windows must have same length") periods, windows = zip(*sorted(zip(periods, windows))) return periods, windows @staticmethod def _remove_overloaded_stl_kwargs(stl_kwargs: dict) -> dict: args = ["endog", "period", "seasonal"] for arg in args: stl_kwargs.pop(arg, None) return stl_kwargs @staticmethod def _default_seasonal_windows(n: int) -> Sequence[int]: return tuple(7 + 4 * i for i in range(1, n + 1)) # See [1] @staticmethod def _to_1d_array(x): y = np.ascontiguousarray(np.squeeze(np.asarray(x)), dtype=np.double) if y.ndim != 1: raise ValueError("y must be a 1d array") return y

Last update: Oct 29, 2024