# Simulating from TBATS models

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
```