API Reference / algorithms / timeseries /

tbats.py

TBATS: trigonometric seasonal exponential smoothing with Box-Cox and ARMA errors.

Classes

TBATS

class algorithms.timeseries.tbats.TBATS(Regressor)

T\ rigonometric seasonality, B\ ox-Cox transform, A\ RMA errors, T\ rend and S\ easonal components -- an exponential smoothing state-space model for series with complex, multiple, high frequency or non-integer seasonal periods.

Classical seasonal models carry one state per seasonal index, so a daily series with a yearly cycle needs 365 states and a period of 365.25 cannot be expressed at all. TBATS instead represents each seasonal pattern by a handful of trigonometric (Fourier) terms, so the state count depends on how many harmonics the pattern needs rather than on the length of the period. That single change is what makes seasonal_periods=[7, 365.25] both representable and cheap.
Constructor 
__init__(
    self,
    seasonal_periods: Optional[Sequence[float] | float] = None,
    n_harmonics: Optional[Sequence[int] | int] = None,
    use_trend: bool = True,
    damped_trend: bool = True,
    box_cox: bool = False,
    box_cox_lambda: Optional[float] = None,
    use_arma_errors: bool = True,
    arma_order: Tuple[int, int] = (),
    maxiter: int = 40,
    tol: float = 1e-06,
)

Overview

Fitting proceeds as follows:

  1. Optionally Box-Cox transform the series to stabilise the variance,
either at a fixed \lambda or by a small grid search.
  1. Seed the level, trend and seasonal states by ordinary least squares
on a design of [1, t, cos/sin harmonics], which gives the recursion a starting point already close to the data.
  1. Run the state-space smoothing recursion, updating level, damped trend
and every trigonometric seasonal pair from the one-step error.
  1. Minimise the sum of squared one-step errors over the smoothing
parameters (\alpha, \beta, \phi, \gamma_1, \gamma_2) with a bounded optimiser.
  1. Fit ARMA errors to the residual series and add their forecast.
  2. Forecast by iterating the recursion with zero errors and inverting
the Box-Cox transform.

Theory

With y_t^{(\lambda)} the Box-Cox transformed observation, the model is

\begin{aligned} y_t^{(\lambda)} &= \ell_{t-1} + \phi b_{t-1} + \sum_{i=1}^{M} s^{(i)}_{t-1} + d_t \\ \ell_t &= \ell_{t-1} + \phi b_{t-1} + \alpha d_t \\ b_t &= \phi b_{t-1} + \beta d_t \end{aligned}

where \phi is the damping parameter and d_t is an ARMA error process. Each seasonal component is the sum of k_i harmonic pairs that rotate at the seasonal frequencies \lambda^{(i)}_j = 2\pi j / m_i:

\begin{aligned} s^{(i)}_{j,t} &= s^{(i)}_{j,t-1}\cos\lambda^{(i)}_j + s^{*(i)}_{j,t-1}\sin\lambda^{(i)}_j + \gamma^{(i)}_1 d_t \\ s^{*(i)}_{j,t} &= -s^{(i)}_{j,t-1}\sin\lambda^{(i)}_j + s^{*(i)}_{j,t-1}\cos\lambda^{(i)}_j + \gamma^{(i)}_2 d_t \end{aligned}

with s^{(i)}_t = \sum_{j=1}^{k_i} s^{(i)}_{j,t}. Because m_i enters only through the angle 2\pi j/m_i, it need not be an integer -- 365.25 is as valid as 12. The Box-Cox transform is

y^{(\lambda)} = \begin{cases} (y^{\lambda} - 1)/\lambda, & \lambda \neq 0 \\ \log y, & \lambda = 0 . \end{cases}

Parameters

seasonal_periods
sequence of float or float or None = None
Seasonal period lengths. May be non-integer (365.25) and there may be several ([7, 365.25]). None disables seasonality.
n_harmonics
sequence of int or int or None = None
Number of harmonic pairs per seasonal period. None picks min(floor(m / 2), 5) for each period, which keeps the state small for very long periods.
use_trend
bool = True
Include a local linear trend component.
damped_trend
bool = True
Damp the trend with a fitted :math:`\phi \in [0.8, 1]`. Ignored when use_trend is False.
box_cox
bool = False
Apply a Box-Cox transform. Requires strictly positive data.
box_cox_lambda
float or None = None
Fixed :math:`\\lambda`. When None and box_cox is True, :math:`\\lambda` is chosen from a small grid by profile likelihood.
use_arma_errors
bool = True
Fit ARMA errors to the smoothing residuals and add their forecast.
arma_order
tuple of (int, int), 0) = (1
The (p, q) order of the residual ARMA model.
maxiter
int = 40
Maximum optimiser iterations for the smoothing parameters.
tol
float = 1e-6
Optimiser convergence tolerance.

Attributes

params_
dict
Fitted smoothing parameters: alpha, beta, phi and the per-period gamma1/gamma2.
lambda_
float or None
Box-Cox parameter actually used, None when box_cox is False.
harmonics_
list of int
Number of harmonic pairs used for each seasonal period.
level_
float
Final level state.
trend_
float
Final trend state (0.0 when use_trend is False).
seasonal_
np.ndarray
Final trigonometric seasonal states, stacked as :math:`(s_1, \dots, s_K, s^{}_1, \dots, s^{}_K)`.
ar_params_
np.ndarray
Fitted AR coefficients of the residual ARMA model.
ma_params_
np.ndarray
Fitted MA coefficients of the residual ARMA model.
fitted_values_
np.ndarray of shape (n_samples,)
One-step-ahead in-sample predictions on the original scale.
resid_
np.ndarray of shape (n_samples,)
One-step-ahead errors on the transformed scale.
sse_
float
Minimised sum of squared one-step errors.
aic_
float
Akaike information criterion computed from sse_.
n_obs_
int
Number of observations supplied to fit.

Notes

Complexity:

  • Training: O(\text{maxiter} \cdot k \cdot n K) where
K = \sum_i k_i is the total number of harmonic pairs and k the number of free smoothing parameters. Crucially it does not grow with the seasonal period lengths.
  • Prediction: O(h K) for h steps.
When to use TBATS:
  • Multiple simultaneous seasonalities -- daily plus weekly plus yearly.
  • Non-integer periods such as 365.25 (leap years) or 52.18 (weeks
per year), which seasonal ARIMA and Holt-Winters cannot represent.
  • High-frequency seasonality where one state per seasonal index
would be prohibitive.
  • Multiplicative-looking variance that a Box-Cox transform can tame.
Simplifications relative to De Livera et al. (2011). This implementation is deliberately a well-tested subset rather than a partial version of the whole paper:
  • The ARMA error stage is estimated in a second pass (Hannan-Rissanen
on the smoothing residuals) rather than jointly with the smoothing parameters inside a single likelihood. Point forecasts are almost unaffected; standard errors from a joint fit would be tighter.
  • Model selection over the discrete choices (trend on/off, damping
on/off, ARMA order, number of harmonics) is not automated by AIC as in the paper; those are constructor parameters.
  • Estimation minimises the sum of squared errors rather than the exact
Gaussian likelihood, and the Box-Cox \lambda is chosen over a small grid rather than jointly optimised.

Seasonality here is additive on the (optionally Box-Cox transformed) scale, which is the standard TBATS formulation -- multiplicative behaviour is obtained through the transform, not through a separate multiplicative seasonal form.

References

DeLivera2011
De Livera, A. M., Hyndman, R. J., & Snyder, R. D. (2011). Forecasting time series with complex seasonal patterns using exponential smoothing. Journal of the American Statistical Association, 106(496), 1513-1527. :doi:`10.1198/jasa.2011.tm09771`
Hyndman2008
Hyndman, R. J., Koehler, A. B., Ord, J. K., & Snyder, R. D. (2008). Forecasting with Exponential Smoothing: The State Space Approach. Springer. :doi:`10.1007/978-3-540-71918-2`
BoxCox1964
Box, G. E. P., & Cox, D. R. (1964). An analysis of transformations. Journal of the Royal Statistical Society: Series B, 26(2), 211-243. :doi:`10.1111/j.2517-6161.1964.tb00553.x`

See Also

ExponentialSmoothing

Holt-Winters smoothing, one state per seasonal index, integer periods only.
SARIMAX

Seasonal ARIMA with exogenous regressors, for a single integer seasonal period.
Prophet

Additive regression with Fourier seasonality and holiday effects.
python 
>>> import numpy as np
>>> from tuiml.algorithms.timeseries.tbats import TBATS
>>> t = np.arange(120)
>>> y = 10.0 + 0.05 * t + 3 * np.sin(2 * np.pi * t / 12)
>>> model = TBATS(seasonal_periods=[12]).fit(y)
>>> forecast = model.predict(steps=12)
>>> truth = 10.0 + 0.05 * np.arange(120, 132) + 3 * np.sin(
...     2 * np.pi * np.arange(120, 132) / 12)
>>> bool(np.mean(np.abs(forecast - truth)) < 0.2)
True

A non-integer period is handled exactly like an integer one:

python 
>>> t = np.arange(200)
>>> y = 5.0 + 2 * np.sin(2 * np.pi * t / 52.18)
>>> model = TBATS(seasonal_periods=52.18).fit(y)
>>> model.predict(steps=4).shape
(4,)
>>> model.harmonics_
[5]

Methods

get_parameter_schema (cls) -> Dict[str, Dict[str, Any]]

Return JSON Schema for constructor parameters.

get_capabilities (cls) -> List[str]

Return supported capabilities.

get_complexity (cls) -> str

Return complexity analysis.

get_references (cls) -> List[str]

Return academic citations.

fit (self, y: np.ndarray, X: Optional[np.ndarray]=None) -> 'TBATS'

Fit the TBATS model to a time series.

Parameters
y
np.ndarray of shape (n_samples,)
The time series to model.
X
np.ndarray = None
Ignored. Present for API consistency with other regressors; use SARIMAX when exogenous regressors are needed.
Returns
self
TBATS
Fitted estimator.
predict (self, steps: int=1, X: Optional[np.ndarray]=None) -> np.ndarray

Forecast future values.

Parameters
steps
int = 1
Number of future time steps to forecast.
X
np.ndarray = None
Ignored. Present for API consistency with other regressors.
Returns
forecast
np.ndarray of shape (steps,)
Point forecasts on the scale of the original series.
fit_predict (self, y: np.ndarray, steps: int=1) -> np.ndarray

Fit the model and forecast in one call.

Parameters
y
np.ndarray of shape (n_samples,)
Time series to fit.
steps
int = 1
Number of future time steps to forecast.
Returns
forecast
np.ndarray of shape (steps,)
Forecasted values.