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

# 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.

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