Skip to content

TBATS

Jump to bottom
Ivan Svetunkov edited this page Oct 8, 2026 · 1 revision

TBATS - Trigonometric Box-Cox ARMA Trend Seasonal model

TBATS (De Livera, Hyndman & Snyder, 2011) models seasonality with trigonometric terms (harmonics) instead of one state per season. This makes it suitable for long seasonal periods (e.g. 336 half-hours in a week), several periods at once (daily and weekly patterns in hourly data) and fractional periods (e.g. lags=c(1, 7, 365.25) for daily data).

smooth implements TBATS in ADAM's Single Source of Error framework: the level and trend of ETS, a trigonometric seasonality for each period, and ADAM's ARMA, all in the space of the Box-Cox transformed data. It is available in R (tbats()) and Python (TBATS) from smooth 4.6.0 (R) and 1.2.0 (Python). The two share the C++ core and fit the same models (see R and Python below).

Function signatures

R

tbats(y, lags = c(1, frequency(y)), harmonics = NULL,
      trend = c("auto", "none", "additive", "damped"),
      lambda = NULL, orders = list(ar = 3, ma = 3, select = TRUE),
      xreg = NULL, regressors = c("use", "select", "adapt"),
      occurrence = c("none", "auto", "fixed", "general", "odds-ratio",
                     "inverse-odds-ratio", "direct"),
      distribution = c("auto", "dnorm", "dlaplace", "ds", "dgnorm"),
      loss = c("likelihood", "MSE", "MAE", "HAM", "MSEh", "TMSE",
               "GTMSE", "MSCE", "GPL"),
      ic = c("AICc", "AIC", "BIC", "BICc"), h = 0, holdout = FALSE,
      initial = c("backcasting", "optimal", "two-stage", "complete", "gradient"),
      bounds = c("admissible", "usual", "none"), silent = TRUE, model = NULL, ...)

Python

class TBATS:
    def __init__(
        self,
        lags: list[float] | None = None,
        harmonics: list[int] | None = None,
        trend: str = "auto",
        lambda_bc: float | None = None,
        orders: dict | None = None,
        regressors: str = "use",
        occurrence: Any = "none",
        distribution: str = "auto",
        loss: str | Callable = "likelihood",
        ic: str = "AICc",
        h: int = 0,
        holdout: bool = False,
        initial: str = "backcasting",
        bounds: str = "admissible",
        verbose: int = 0,
        # ... B, lb, ub, NLopt settings (maxeval, maxtime, algorithm, tolerances),
        # n_iterations, head_length, fi, step_size, shape
    ) -> None: ...

    def fit(self, y, X=None) -> "TBATS": ...

The Box-Cox parameter is lambda in R and lambda_bc in Python (lambda is a reserved word there). Explanatory variables go to xreg in R and to fit(y, X) in Python.

Overview

The model works in the space of the Box-Cox transformed data, with lambda estimated in [0, 1] (with loss="likelihood") unless it is provided. Its components are:

  • the level and the trend of ETS, the trend being none, additive or damped;
  • for each seasonal period m, k harmonics: pairs of states that rotate at the frequencies 2πj/m, j = 1, ..., k, with the smoothing parameters γ1 and γ2 shared by the harmonics of the same period;
  • ARMA errors in ADAM's state-space form.

The model's name follows De Livera et al. (2011): TBATS(λ, {p,q}, φ, <m1,k1>, <m2,k2>, ...), where p and q are the AR and MA orders, φ the damping parameter (- without damping) and k the number of harmonics of the period m. Harmonics that coincide across periods are dropped: with periods 48 and 336, the weekly harmonics 7 and 14 have the frequencies of the daily harmonics 1 and 2.

The forecasts are produced in the transformed space and transformed back. The point forecasts are by default the medians (the inverse transform of the point forecasts of the transformed data); point="mean" gives the means.

Automatic selection

By default, tbats() / TBATS selects the structure itself:

  1. the number of harmonics for each period, by the information criterion of the global model (a regression on a trend and the Fourier terms);
  2. the trend, comparing no trend, an additive and a damped one by the information criterion;
  3. the ARMA orders (up to orders$ar and orders$ma), screened with Hannan-Rissanen on the residuals of the global model; the model with the ARMA is kept only if it improves the information criterion;
  4. the distribution of the error term (see below).

Specifying harmonics, trend, orders=list(ar=p, ma=q, select=FALSE), lambda and distribution skips the corresponding step, which is much faster on long series.

distribution="auto"

The Generalised Normal distribution nests the others: its shape β is 2 for the normal, 1 for the Laplace and 0.5 for the S distribution. The default distribution="auto" uses this to choose the distribution without fitting a model for each one:

  1. the structure is selected with the Generalised Normal distribution, its shape estimated;
  2. the named distribution closest to the estimated shape on the log scale is taken: S for β < 0.71, Laplace for 0.71 ≤ β < 1.41, normal otherwise;
  3. that distribution is fitted on the selected structure, from the default starting values and from the Generalised Normal estimates, and the fit with the higher likelihood is returned.

The chosen distribution is in model$distribution (R) / model.distribution_ (Python), and its information criterion is added to ICs / ics. With a loss other than the likelihood, "auto" is the distribution the loss implies, as in ADAM: "dlaplace" for "MAE", "ds" for "HAM", "dnorm" otherwise.

R Usage

library(smooth)

# Automatic model: harmonics, trend, ARMA orders and distribution selected
model <- tbats(y, lags=c(1, 48, 336), h=336, holdout=TRUE)
model
forecast(model, h=336, interval="prediction") |> plot()

# Known structure: 8 and 16 harmonics, no trend, no ARMA, no transformation
tbats(y, lags=c(1, 48, 336), harmonics=c(8, 16), trend="none",
      orders=list(ar=0, ma=0, select=FALSE), lambda=1,
      distribution="dnorm", h=336, holdout=TRUE)

# Fractional period, e.g. daily data with the annual seasonality
tbats(yDaily, lags=c(1, 7, 365.25), h=28)

Python Usage

from smooth import TBATS

# Automatic model
model = TBATS(lags=[1, 48, 336], h=336, holdout=True)
model.fit(y)
print(model)              # TBATS(lambda, {p,q}, phi, <48,k1>, <336,k2>), AICc ...
model.ics                 # the information criteria of the candidates
model.distribution_       # the distribution chosen by "auto"

fc = model.predict(h=336, interval="prediction", level=0.95)
fc.mean, fc.lower, fc.upper

# Known structure
model = TBATS(lags=[1, 48, 336], harmonics=[8, 16], trend="none",
              orders={"ar": 0, "ma": 0, "select": False}, lambda_bc=1,
              distribution="dnorm", h=336, holdout=True)
model.fit(y)

Initialisation

The initial states are backcast by default (initial="backcasting"), which keeps the number of estimated parameters low even with many harmonics. "optimal" estimates them as deviations from the global model, "two-stage" starts the optimisation of all parameters from the backcast fit, "complete" backcasts all the states including the coefficients of the explanatory variables, and "gradient" solves for the initial states by least squares for each set of the other parameters (slow with many harmonics). See Initialisation.

Explanatory variables and intermittent demand

Explanatory variables enter in the transformed space, as the regressors of ADAM: regressors="use" keeps their coefficients constant, "select" runs stepwise() on the errors of the model chosen without them (refitted with the selected ones and kept if it improves the information criterion), and "adapt" gives each coefficient a smoothing parameter. Their future values come from newdata (R) / predict(h, X) (Python), else the holdout, else their ADAM forecasts with a warning. See Explanatory-Variables.

tbats(y, xreg=X, regressors="select", h=12, holdout=TRUE)
TBATS(lags=[1, 12], regressors="select", h=12, holdout=True).fit(y, X)

occurrence= models intermittent demand: a fitted OM model, the probabilities (or 0/1) of occurrence, or a type of occurrence model to construct (a level-only om(model="ZXN") with the trend selected). The sizes are modelled on the non-zero values, the states evolving through the zeros, and the forecasts combine the probability with the sizes. Missing values are treated as gaps.

Parameters

Parameter Type (R) Type (Python) Default Description
y vector/ts NDArray/Series (in fit) - Time series
lags numeric vector list[float] c(1, frequency(y)) (R) / [1] (Python) 1 and the seasonal periods, possibly fractional
harmonics integer vector list[int]/None NULL Harmonics per period above 1; selected if NULL/None
trend character str "auto" "none" / "additive" / "damped" / "auto"
lambda (R) / lambda_bc (Python) numeric float/None NULL Box-Cox parameter; estimated in [0, 1] if NULL/None
orders list dict/None list(ar=3, ma=3, select=TRUE) Maximum ARMA orders and whether to select them
xreg (R) / X in fit (Python) matrix/data.frame NDArray/DataFrame NULL/None Explanatory variables
regressors character str "use" "use" / "select" / "adapt"
occurrence character/om/numeric str/OM/NDArray "none" Occurrence model for intermittent demand
distribution character str "auto" "auto" / "dnorm" / "dlaplace" / "ds" / "dgnorm"
loss character/function str/callable "likelihood" Loss function, see Loss-Functions
ic character str "AICc" Information criterion of the selection
h, holdout integer, logical int, bool 0, FALSE Forecast horizon and holdout
initial character str "backcasting" See Initialisation
bounds character str "admissible" "admissible" (stable model), "usual" or "none", see Bounds
model adamTBATS — NULL A previous model to reuse (R)

Fitted Attributes

Element (R) Element (Python) Description
model model_name Model name, e.g. "TBATS(0.13, {1,1}, -, <24,4>, <168,7>)"
lambda lambda_ Box-Cox parameter
harmonics, periods harmonics_, periods_ Harmonics and seasonal periods
trendType trend_type_ "none" / "additive" / "damped"
orders orders_ ARMA orders
distribution distribution_ Distribution of the error term (the one chosen with "auto")
phi phi_ Damping parameter
persistence persistence_vector Smoothing parameters
transition, measurement transition, measurement State-space matrices
initial initial_value Initial states
states states States over time
fitted, residuals fitted, residuals Fitted values and residuals (of the transformed data for the residuals)
forecast forecast_ / predict() Point forecasts
logLik loglik Log-likelihood, with the Jacobian of the transform
ICs ics Information criteria of the fitted candidates
B coef, coef_names Estimated parameters
scale scale Scale of the error term
timeElapsed time_elapsed Estimation time

R returns an object of class c("adamTBATS", "adam", "smooth"), so the methods that the forecast package registers for its own class "tbats" never dispatch on it.

Methods

The methods of ADAM apply, working in the space of the transformed data and transforming the results back:

R and Python

R and Python share the C++ core, the least squares of the global model (a Householder QR in src/headers/olsCore.h) and the eigenvalues of the admissible bounds (eigenModuliCore() in src/headers/eigenCalc.h), all written without BLAS or LAPACK so that both languages round identically. The likelihood surfaces of TBATS are flat enough for a last-bit difference to change the estimates, so the densities matter too: with R's greybox 2.0.10 and Python's greybox 1.0.9, which compute the log-densities of the Laplace, S and Generalised Normal distributions identically, the two fit the same model to all 414 hourly series of the M4 competition, with identical forecasts. With earlier greybox versions, distribution="auto" gave different models on a few series. See R-Python-differences.

Performance

On the 414 hourly series of M4 (h=48, periods 24 and 168), RMSSE and pinball scaled by the in-sample one-step differences:

smooth (auto) smooth (dnorm) forecast::tbats() statsforecast AutoTBATS
Mean RMSSE 0.829 0.845 0.909 1.041
Mean pinball 0.316 0.323 0.362 0.408
Calibration error 0.041 0.051 0.083 0.067
Seconds per series 20.0 9.7 40.6 34.2

The notebook, the runners and the per-series results are in python/tests/notebooks/ of the repository.

Relationship to ADAM

TBATS uses ADAM's C++ core for the fit, the forecasts and the simulations: the harmonics take the place of ADAM's ARIMA states, with a sparse transition matrix that makes the fit faster with many harmonics. Unlike ADAM with multiple seasonal lags, it needs a pair of states per harmonic instead of one state per season, so a weekly period of half-hourly data costs a few dozen states instead of 336.

References

  • 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
  • Svetunkov, I. (2023). Forecasting and Analytics with the Augmented Dynamic Adaptive Model (ADAM). Chapman and Hall/CRC. DOI: 10.1201/9781003452652

See Also

Related Functions

  • ADAM - Main unified framework
  • MSARIMA - Multiple seasonal ARIMA
  • msdecompose - Multiple seasonal decomposition
  • OM - Occurrence models for intermittent demand

Parameter Documentation

Clone this wiki locally