How to get a density from a forecast with prediction interval

Question

Some reproducible code to have in your environment a time series and a possible forecast:

packages <- c('forecast',
              'robets',
              'quantmod')
lapply(X = packages,
       FUN = function(package){
         if (!require(package = package,
                      character.only = TRUE))
         {
           install.packages(pkgs = package,
                            repos = "https://cloud.r-project.org")
           library(package = package,
                   character.only = TRUE)
         } else {
           library(package = package,
                   character.only = TRUE)    
         }
       })
getSymbols('GLD')
adjustOHLC(GLD)
GLD %>%
  Cl() %>%
  log() %>%
  robets() %>%
  forecast(h = 250, level = seq(from = 51, to = 99, by = 1)) %>%
  autoplot()

I would like to draw a random sample from the distribution of forecasted values, in order to possibly fit a probability density function. The value returned by forecast has some slots which might be useful, like $mean and $residuals, but I don't know how to draw such a sample. Furthermore, aside from this specific example, I would like to understand the theoretical procedure to get a density from a forecast with expected value and prediction interval only.

Answer 1

If you can assume a normal distribution of errors, then you can construct the distribution of the forecasts from the mean and the prediction intervals.

See this blog post by Rob Hyndman for details.

Alternatively, you can generate sample paths for the forecasts, and then calculate the quantiles - in this case you don't even need the prediction intervals (but you still need to make assumptions on the distribution).