Sample path generation for arima with R language

Assume a given time series and let's denote it by "series". Also assuming h as a given nonnegative integer, the following code should give you an object called forecast:

ARIMA_MODEL = arima(series, order=c(1,1,1)) 
forecast = predict(ARIMA_MODEL, h)

such that:

forecast$pred: mean value of the prediction for next h steps.
forecast$se: standard error of the prediction at each time step.

Now, not only obtaining the mean results, but I also would like to generate a sample path (stochastic forecast), using the information of ARIMA_MODEL (do some sort of simulation).

How to do so? Note that the problem is different to using arima.sim function naively with the same (p,d,q) parameters, since I would like to use the information of calculated mean, variance, etc which are available in ARIMA_MODEL object.

Solution

This can be done using the forecast package.

library(forecast)
library(ggplot2)

model <- Arima(WWWusage, order=c(1,1,1))
forecast <- forecast(model, h=20)
sample_path <- simulate(model, nsim=20)
autoplot(forecast) +
  autolayer(sample_path)

^{Created on 2020-07-28 by the reprex package (v0.3.0)}