Why do the simulations of my SARIMA model not resemble my original data?

Question

I want to simulate a SARIMA model I obtained using the auto.arima function from the R package "forecast". My objective is to be able to do a lot of simulations in order to "predict" for example the centennial flood.

My model's name is "final.model", and here is the output of the auto.arima function:

> summary(final.model) # -> ARIMA(1,0,0)(4,1,0)[12] 
Series: mlog 
ARIMA(1,0,0)(4,1,0)[12] 

Coefficients:
         ar1     sar1     sar2     sar3     sar4
      0.3317  -0.7356  -0.5288  -0.3651  -0.2587
s.e.  0.0388   0.0402   0.0489   0.0489   0.0409

sigma^2 = 0.361:  log likelihood = -536.49
AIC=1084.98   AICc=1085.12   BIC=1111.23

Training set error measures:
                      ME      RMSE       MAE       MPE     MAPE      MASE        ACF1
Training set 0.006117115 0.5922533 0.4615495 -41.48468 111.4677 0.7431494 -0.01109284 

When I try to do simulations, I obtain simulations that are too large in comparison with my original time series, here is an example:

Here is the code I used:

n_simulations=10 # number of replication
q=nsim/12
colors=c("red","blue","green","orange","purple")

plot(exp(mlog),xlim=c(2015,2035),ylim=c(0,100),type="l",col="black",
main="Original Data and Simulations",xlab="Year",ylab="Values")

simulated_values=replicate(n=n_simulations,expr=simulate(final.model,nsim=nsim))

for (i in 1:n_simulations) {
    lines(time(monthly_ts)[length(monthly_ts)] + 1:(q*12)/12+1/12, exp(simulated_values[,i]), col=alpha("red", 0.2))
}

Three last remarks: I cannot provide the data (as it is not allowed by the source), in the auto.arima function I provided the logarithm of the original time series, and it is environmental data (i.e. stream flow of a river).

My questions are thus: is my way of doing correct? Why do I obtain such weird results?

Answer 1

Your model is $\operatorname{ARIMA}(p,d,q)(P,D,Q)_m$ where $p$ is the number of AR terms, $d$ is the number of differences, $q$ is the number of MA terms, $P$ is the number of seasonal AR terms, $D$ is the number of seasonal (i.e., year-over-year) differences, and $Q$ is the number of seasonal MA terms, and $m$ is the frequency of the data.

The model you estimate is: $\operatorname{ARIMA}(1,0,0)(4,1,0)_{12}$, so you are modeling the seasonally differenced data (since $D=1$). In ARIMA modeling, whenever you take a difference of your data (be it a 'regular' difference of a 'seasonal' difference), you are making an assumption that your data has a unit root. A key feature of unit roots is that the unconditional variance increases linearly over time. As you forecast further into the future, the variance will increase, and your prediction intervals will get steadily wider over time. The point-forecasts (the mean of the forecast distribution at each point in time) will still stay roughly constant over time (with seasonal variation in your case).

For many reasons, fitting an ARIMA model on the seasonal difference does not seem appropriate for your data. Instead, you may want to experiment with using a regression model with constant seasonal dummies to capture the seasonality, rather than differencing. You can specify a linear regression model with ARIMA errors by specify your dummies as part of the xreg option in the Arima() function from the fpp2 package.

It seems like you are particularly interested in modeling the probability of rare events ("centennial flood"). If the errors from the above approach has fatter tails than implied by a Normal distribution, you may want to fit a model that has student-t distributed errors. I'm not aware of a package that allows you to easily do this in the context of an Arima model with seasonal dummies. Perhaps rugarch, but I'm not sure.

Finally, I don't know much about the field, but you probably want to look into extreme value theory if your interest lies in predicting extreme events.