# Simulating from TBATS models

**R on Rob J Hyndman**, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

I’ve had several requests for an R function to simulate future values from a TBATS model. We will eventually include TBATS in the `fable`

package, and the facilities will be added there. But in the meantime, if you are using the `forecast`

package and want to simulate from a fitted TBATS model, here is how do it.

## Simulating via one-step forecasts

Doing it efficiently would require a more complicated approach, but this is super easy if you are willing to sacrifice some speed. The trick is to realise that a simulation can be handled easily for almost any time series model using residuals and one-step forecasts. Note that a residual is given by $e_t = y_t – \hat{y}_{t|t-1}$ so we can write

$$y_t = \hat{y}_{t|t-1} + e_t.$$

Therefore, given data to time $T$, we can simulate iteratively using $$y_{T+i}^* = \hat{y}_{T+i|T+i-1} + \varepsilon_{T+i}, \qquad i=1,\dots,h,$$ where $\varepsilon_{T+i}$ is randomly generated from the error distribution, or bootstrapped by randomly sampling from past residuals. The value of $\hat{y}_{T+i|T+i-1}$ can be obtained by applying the model to the series $\{y_1,\dots,y_T,y^*_{T+1},\dots,y^*_{T+i-1}\}$ (without re-estimating the parameters) and forecasting one-step ahead. This is the same trick we use to get prediction intervals for neural network models.

## Implementation

Because `simulate()`

is an S3 method in R, we have to make sure the corresponding `simulate.tbats()`

function has all the necessary arguments to match other `simulate`

functions. It’s also good to make it as close as possible to other `simulate`

functions in the `forecast`

package, to make it easier for users when switching between them. The real work is done in the last few lines.

simulate.tbats <- function(object, nsim=length(object$y), seed = NULL, future=TRUE, bootstrap=FALSE, innov = NULL, ...) { if (is.null(innov)) { if (!exists(".Random.seed", envir = .GlobalEnv)) { runif(1) } if (is.null(seed)) { RNGstate <- .Random.seed } else { R.seed <- .Random.seed set.seed(seed) RNGstate <- structure(seed, kind = as.list(RNGkind())) on.exit(assign(".Random.seed", R.seed, envir = .GlobalEnv)) } } else { nsim <- length(innov) } if (bootstrap) { res <- residuals(object) res <- na.omit(res - mean(res, na.rm = TRUE)) e <- sample(res, nsim, replace = TRUE) } else if (is.null(innov)) { e <- rnorm(nsim, 0, sqrt(object$variance)) } else { e <- innov } x <- getResponse(object) y <- numeric(nsim) if(future) { dataplusy <- x } else { # Start somewhere in the original series dataplusy <- ts(sample(x, 1), start=-1/frequency(x), frequency = frequency(x)) } fitplus <- object for(i in seq_along(y)) { y[i] <- forecast(fitplus, h=1)$mean + e[i] dataplusy <- ts(c(dataplusy, y[i]), start=start(dataplusy), frequency=frequency(dataplusy)) fitplus <- tbats(dataplusy, model=fitplus) } return(tail(dataplusy, nsim)) }

I’ve added this to the `forecast`

package for the next version.

Something similar could be written for any other forecasting function that doesn’t already have a `simulate`

method. Just swap the `tbats`

call to the relevant modelling function.

## Example

library(forecast) library(ggplot2) fit <- tbats(USAccDeaths) p <- USAccDeaths %>% autoplot() + labs(x = "Year", y = "US Accidental Deaths", title = "TBATS simulations") for (i in seq(9)) { p <- p + autolayer(simulate(fit, nsim = 36), series = paste("Sim", i)) } p

**leave a comment**for the author, please follow the link and comment on their blog:

**R on Rob J Hyndman**.

R-bloggers.com offers

**daily e-mail updates**about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.