tbats() - Trigonometric Box-Cox ARMA Trend Seasonal model

This vignette covers the tbats() function, which is part of smooth package. It implements the TBATS model of De Livera et al. (2011) in the Single Source of Error state space framework of ADAM: the level and trend of ETS, a trigonometric seasonality for each seasonal period, and an ARMA, all in the space of the Box-Cox transformed data. The trigonometric seasonality makes the model suitable for long and several seasonal periods (e.g. daily and weekly patterns in half-hourly data), and the periods can be fractional (e.g. lags=c(1, 7, 365.25) for daily data).

Let’s load the package:

require(smooth)

The data

We use the same artificial half-hourly series as in the adam() vignette, generated with sim.gum(), with the daily (48) and weekly (336) seasonalities:

set.seed(41)
ordersGUM <- c(1,1,1)
lagsGUM <- c(1,48,336)
initialGUM1 <- -25381.7
initialGUM2 <- c(23955.09, 24248.75, 24848.54, 25012.63, 24634.14, 24548.22, 24544.63, 24572.77,
                 24498.33, 24250.94, 24545.44, 25005.92, 26164.65, 27038.55, 28262.16, 28619.83,
                 28892.19, 28575.07, 28837.87, 28695.12, 28623.02, 28679.42, 28682.16, 28683.40,
                 28647.97, 28374.42, 28261.56, 28199.69, 28341.69, 28314.12, 28252.46, 28491.20,
                 28647.98, 28761.28, 28560.11, 28059.95, 27719.22, 27530.23, 27315.47, 27028.83,
                 26933.75, 26961.91, 27372.44, 27362.18, 27271.31, 26365.97, 25570.88, 25058.01)
initialGUM3 <- c(23920.16, 23026.43, 22812.23, 23169.52, 23332.56, 23129.27, 22941.20, 22692.40,
                 22607.53, 22427.79, 22227.64, 22580.72, 23871.99, 25758.34, 28092.21, 30220.46,
                 31786.51, 32699.80, 33225.72, 33788.82, 33892.25, 34112.97, 34231.06, 34449.53,
                 34423.61, 34333.93, 34085.28, 33948.46, 33791.81, 33736.17, 33536.61, 33633.48,
                 33798.09, 33918.13, 33871.41, 33403.75, 32706.46, 31929.96, 31400.48, 30798.24,
                 29958.04, 30020.36, 29822.62, 30414.88, 30100.74, 29833.49, 28302.29, 26906.72,
                 26378.64, 25382.11, 25108.30, 25407.07, 25469.06, 25291.89, 25054.11, 24802.21,
                 24681.89, 24366.97, 24134.74, 24304.08, 25253.99, 26950.23, 29080.48, 31076.33,
                 32453.20, 33232.81, 33661.61, 33991.21, 34017.02, 34164.47, 34398.01, 34655.21,
                 34746.83, 34596.60, 34396.54, 34236.31, 34153.32, 34102.62, 33970.92, 34016.13,
                 34237.27, 34430.08, 34379.39, 33944.06, 33154.67, 32418.62, 31781.90, 31208.69,
                 30662.59, 30230.67, 30062.80, 30421.11, 30710.54, 30239.27, 28949.56, 27506.96,
                 26891.75, 25946.24, 25599.88, 25921.47, 26023.51, 25826.29, 25548.72, 25405.78,
                 25210.45, 25046.38, 24759.76, 24957.54, 25815.10, 27568.98, 29765.24, 31728.25,
                 32987.51, 33633.74, 34021.09, 34407.19, 34464.65, 34540.67, 34644.56, 34756.59,
                 34743.81, 34630.05, 34506.39, 34319.61, 34110.96, 33961.19, 33876.04, 33969.95,
                 34220.96, 34444.66, 34474.57, 34018.83, 33307.40, 32718.90, 32115.27, 31663.53,
                 30903.82, 31013.83, 31025.04, 31106.81, 30681.74, 30245.70, 29055.49, 27582.68,
                 26974.67, 25993.83, 25701.93, 25940.87, 26098.63, 25771.85, 25468.41, 25315.74,
                 25131.87, 24913.15, 24641.53, 24807.15, 25760.85, 27386.39, 29570.03, 31634.00,
                 32911.26, 33603.94, 34020.90, 34297.65, 34308.37, 34504.71, 34586.78, 34725.81,
                 34765.47, 34619.92, 34478.54, 34285.00, 34071.90, 33986.48, 33756.85, 33799.37,
                 33987.95, 34047.32, 33924.48, 33580.82, 32905.87, 32293.86, 31670.02, 31092.57,
                 30639.73, 30245.42, 30281.61, 30484.33, 30349.51, 29889.23, 28570.31, 27185.55,
                 26521.85, 25543.84, 25187.82, 25371.59, 25410.07, 25077.67, 24741.93, 24554.62,
                 24427.19, 24127.21, 23887.55, 24028.40, 24981.34, 26652.32, 28808.00, 30847.09,
                 32304.13, 33059.02, 33562.51, 33878.96, 33976.68, 34172.61, 34274.50, 34328.71,
                 34370.12, 34095.69, 33797.46, 33522.96, 33169.94, 32883.32, 32586.24, 32380.84,
                 32425.30, 32532.69, 32444.24, 32132.49, 31582.39, 30926.58, 30347.73, 29518.04,
                 29070.95, 28586.20, 28416.94, 28598.76, 28529.75, 28424.68, 27588.76, 26604.13,
                 26101.63, 25003.82, 24576.66, 24634.66, 24586.21, 24224.92, 23858.42, 23577.32,
                 23272.28, 22772.00, 22215.13, 21987.29, 21948.95, 22310.79, 22853.79, 24226.06,
                 25772.55, 27266.27, 28045.65, 28606.14, 28793.51, 28755.83, 28613.74, 28376.47,
                 27900.76, 27682.75, 27089.10, 26481.80, 26062.94, 25717.46, 25500.27, 25171.05,
                 25223.12, 25634.63, 26306.31, 26822.46, 26787.57, 26571.18, 26405.21, 26148.41,
                 25704.47, 25473.10, 25265.97, 26006.94, 26408.68, 26592.04, 26224.64, 25407.27,
                 25090.35, 23930.21, 23534.13, 23585.75, 23556.93, 23230.25, 22880.24, 22525.52,
                 22236.71, 21715.08, 21051.17, 20689.40, 20099.18, 19939.71, 19722.69, 20421.58,
                 21542.03, 22962.69, 23848.69, 24958.84, 25938.72, 26316.56, 26742.61, 26990.79,
                 27116.94, 27168.78, 26464.41, 25703.23, 25103.56, 24891.27, 24715.27, 24436.51,
                 24327.31, 24473.02, 24893.89, 25304.13, 25591.77, 25653.00, 25897.55, 25859.32,
                 25918.32, 25984.63, 26232.01, 26810.86, 27209.70, 26863.50, 25734.54, 24456.96)
y <- sim.gum(orders=ordersGUM, lags=lagsGUM, nsim=1, frequency=336, obs=3360,
             measurement=rep(1,3), transition=diag(3), persistence=c(0.045,0.162,0.375),
             initial=rbind(initialGUM1,initialGUM2,initialGUM3))$data

Automatic model

By default, tbats() selects the structure of the model 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 the models with no trend, an additive and a damped one by the information criterion;
  3. the orders of the ARMA (up to orders$ar and orders$ma), screened 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, explained in Distribution.

The parameter lambda of the Box-Cox transform is estimated together with the other parameters (with loss="likelihood"), unless it is provided. We hold out the last week of data:

tbatsModel <- tbats(y, lags=c(1,48,336), h=336, holdout=TRUE)
tbatsModel
## Time elapsed: 81.56 seconds
## Model estimated using tbats() function: TBATS(0.899, {3,0}, 0.57, <48,9>, <336,36>)
## With backcasting initialisation
## Distribution assumed in the model: Normal
## Loss function type: likelihood; Loss function value: 21003.93
## Smoothing parameters:
##       alpha        beta  gamma1[48]  gamma2[48] gamma1[336] gamma2[336] 
##     -0.0710      0.0767      0.0038     -0.0033      0.0000      0.0000 
## Damping parameter: 0.5695 
## Box-Cox lambda: 0.8988
## ARMA parameters of the model:
##         Lag 1
## AR(1) -0.0984
## AR(2) -0.2464
## AR(3) -0.1020
## 
## Sample size: 3024
## Number of estimated parameters: 107
## Number of degrees of freedom: 2917
## Information criteria:
##      AIC     AICc      BIC     BICc 
## 42221.86 42229.78 42865.39 42897.15 
## 
## Forecast errors:
## ME: 45.4; MAE: 246.083; RMSE: 310.739
## sCE: 50.33%; Asymmetry: 17%; sMAE: 0.812%; sMSE: 0.011%
## MASE: 0.331; RMSSE: 0.302; rMAE: 0.038; rRMSE: 0.039

The name of the model follows De Livera et al. (2011): TBATS(\(\lambda\), {\(p\),\(q\)}, \(\phi\), <\(m_1\),\(k_1\)>, <\(m_2\),\(k_2\)>), where \(\lambda\) is the Box-Cox parameter, \(p\) and \(q\) are the AR and MA orders, \(\phi\) is the damping parameter (“-” for a model without damping) and \(k_i\) is the number of harmonics for the period \(m_i\). The smoothing parameters \(\gamma_1\) and \(\gamma_2\) are shared by the harmonics of the same period.

The forecasts are produced in the transformed space and then transformed back. The point forecasts are by default the medians of the forecast distribution (the inverse transform of the point forecasts of the transformed data):

forecast(tbatsModel, h=336, interval="prediction") |> plot()

Specifying the structure

The automatic selection takes time on long series, because several models are estimated. If the structure is known, it can be specified directly, which is much faster. Here is the model with 8 and 16 harmonics, no trend, no ARMA, no transformation (lambda=1) and the normal distribution:

tbatsModel <- 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)
tbatsModel
## Time elapsed: 2.92 seconds
## Model estimated using tbats() function: TBATS(1, {0,0}, -, <48,8>, <336,14>)
## With backcasting initialisation
## Distribution assumed in the model: Normal
## Loss function type: likelihood; Loss function value: 22623.76
## Smoothing parameters:
##       alpha  gamma1[48]  gamma2[48] gamma1[336] gamma2[336] 
##      0.0371      0.0024      0.0060      0.0508     -0.0064 
## Box-Cox lambda: 1
## Sample size: 3024
## Number of estimated parameters: 51
## Number of degrees of freedom: 2973
## Information criteria:
##      AIC     AICc      BIC     BICc 
## 45349.52 45351.30 45656.25 45663.40 
## 
## Forecast errors:
## ME: 69.14; MAE: 488.601; RMSE: 627.919
## sCE: 76.649%; Asymmetry: 16.9%; sMAE: 1.612%; sMSE: 0.043%
## MASE: 0.658; RMSSE: 0.611; rMAE: 0.076; rRMSE: 0.08

The weekly harmonics 7 and 14 have the same frequencies as the daily harmonics 1 and 2, so they are dropped, leaving 14 harmonics for the weekly period. The trend can be "none", "additive", "damped" or "auto", and orders=list(ar=1, ma=1, select=FALSE) fits the ARMA(1,1) instead of selecting the orders.

Distribution

The error term in the space of the transformed data follows the normal ("dnorm"), Laplace ("dlaplace"), S ("ds") or Generalised Normal ("dgnorm") distribution. The last one nests the others: its shape \(\beta\) 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 of the model is selected as above, with the Generalised Normal distribution and its shape estimated;
  2. the named distribution closest to the estimated shape on the log scale is taken: the S distribution for \(\beta < 0.71\), the Laplace for \(0.71 \leq \beta < 1.41\) and the normal otherwise;
  3. this distribution is fitted on the selected structure, from the default starting values and from the estimates of the Generalised Normal model, and the fit with the higher likelihood is returned.

The shape tells how heavy the tails of the error term are, so the data with outliers get the Laplace or the S distribution, which are less sensitive to them, and the data without them get the normal. The estimate of the shape is noisy, so the named distribution, which has one parameter less, is used even when the Generalised Normal one has the lower information criterion: on the hourly series of the M4 competition, its forecasts were more accurate than those of the Generalised Normal distribution and of the choice between the two by the information criterion. The chosen distribution and its information criterion are in the model:

tbatsModel$distribution
tbatsModel$ICs

The two steps after the selection add to the time of the estimation, so a known distribution is better specified directly, e.g. distribution="dnorm". With a loss other than the likelihood, "auto" is the distribution that the loss implies, as in adam(): "dlaplace" for loss="MAE", "ds" for loss="HAM" and "dnorm" otherwise.

Initialisation

The initial states are by default backcasted (initial="backcasting"), which keeps the number of estimated parameters low even with many harmonics. The other options are:

  • initial="optimal" - the initial states are estimated together with the other parameters, as deviations from the global model;
  • initial="two-stage" - backcasting first, then the optimisation of all parameters starting from its results;
  • initial="complete" - backcasting of all the states, including the coefficients of the explanatory variables;
  • initial="gradient" - the initial states are solved for exactly, by least squares, for each set of the other parameters.

The last option is slower than backcasting on long series with many harmonics (several minutes for the model above instead of a few seconds), because the solve involves all the states of the model:

tbats(y, lags=c(1,48,336), harmonics=c(8,16), trend="none",
      orders=list(ar=0, ma=0, select=FALSE), lambda=1,
      initial="gradient", h=336, holdout=TRUE)

Provided parameters

As in adam(), the smoothing parameters (persistence), the damping parameter (phi), the ARMA parameters (arma) and the initial states (initial) can be provided, and only the rest is estimated. The lists name the parts that are fixed; the seasonal elements have one vector per period of lags: the pair \(\gamma_1, \gamma_2\) of the smoothing parameters, and the sine and then the cosine coefficients of the harmonics of the initial seasonality. The initial states are in the space of the Box-Cox transformed data, so they are meaningful with a provided lambda:

tbats(AirPassengers, harmonics=2, trend="damped", lambda=0,
      orders=list(ar=1, ma=0, select=FALSE),
      persistence=list(level=0.3, seasonal=list(c(0.001, 0))), phi=0.98,
      arma=list(ar=0.5), initial=list(level=4.7), h=12, holdout=TRUE)

The vectors without names are matched with the parameters by their positions, so they need a structure that is not selected (the trend, the harmonics and the ARMA orders provided), and phi needs a provided damped trend: otherwise, they are estimated with a warning. Providing arma switches off the selection of the ARMA orders.

Explanatory variables and intermittent demand

Explanatory variables are included via xreg, and regressors defines what to do with them: "use" all of them, "select" the relevant ones based on the information criterion, or "adapt" their parameters over time, as in adam().

The function also handles intermittent data via the occurrence parameter, which is either a type of the occurrence model to construct (see om()), a fitted occurrence model, or the provided probabilities of occurrence. The sizes are then modelled on the non-zero values:

set.seed(41)
yIntermittent <- ts(rpois(120, 3) * rbinom(120, 1, 0.7), frequency=12)
tbats(yIntermittent, lags=c(1,12), occurrence="auto", h=12, holdout=TRUE)
## Time elapsed: 0.45 seconds
## Model estimated using tbats() function: TBATS(0.418, {0,0}, -)
## With backcasting initialisation
## Occurrence model type: Fixed probability
## Distribution assumed in the model: Mixture of Bernoulli and Normal
## Loss function type: likelihood; Loss function value: 115.7044
## Smoothing parameters:
## alpha 
## 0.062 
## Box-Cox lambda: 0.4177
## Sample size: 108
## Number of estimated parameters: 5
## Number of degrees of freedom: 103
## Information criteria:
##      AIC     AICc      BIC     BICc 
## 388.1145 388.5028 401.5251 397.7521 
## 
## Forecast errors:
## Asymmetry: -26.908%; sMSE: 34.656%; rRMSE: 0.761; sPIS: 653.675%; sCE: -44.272%

Methods

Since tbats() returns an object of the class adam, the methods of adam() work for it as well, including summary(), plot(), vcov(), confint(), simulate(), reapply() and reforecast(). For example, here are the diagnostic plots of the model: actuals vs fitted, the QQ-plot of the residuals, the fit over time and the ACF of the residuals:

par(mfcol=c(2,2))
plot(tbatsModel, which=c(1,6,7,10))

References

De Livera, Alysha M., Rob J. Hyndman, and Ralph D. Snyder. 2011. “Forecasting time series with complex seasonal patterns using exponential smoothing.” Journal of the American Statistical Association 106 (496): 1513–27. https://doi.org/10.1198/jasa.2011.tm09771.