Probabilistic forecast of ARIMA peak value I am using a seasonal ARIMA model to model seasonal influenza. Specifically, I want a probabilistic forecast of the peak value for flu-like illness in the coming season. While it's easy to get forecasts for individual weeks ahead in the forecast package, it's not so straightforward to get a forecast of what the peak value will be.
My thought is to take advantage of the asymptotic multivariate normality of the parameter estimates and randomly generate new parameter estimates using the estimated covariance. From those, I can simulate the upcoming season and record the peak value. If I do this 10,000 times I can get an estimated distribution of the peak value. 
My question is: Is there some better way to generate a forecast for this value that I'm missing?
Caveats: I know a seasonal ARIMA is an oversimplification - I'm using this as a baseline to compare other models to. I've included my code below if it's helpful for anyone. These are not the full set of data, but enough to create a reproducible example. I used the auto.arima function previously to determine that the ARIMA (2, 0, 0)(1, 0, 0) is the best fit to the full data.
Thanks!
Code:
library(cdcfluview)
library(MASS)
library(dplyr)
library(forecast)

set.seed(500)

# Create ARIMA model
usflu_Arima <- get_flu_data("national", "ilinet", years=2013:2015) %>%
            select(
              observation=X..WEIGHTED.ILI) %>%
            tail(
              n=156) %>%
            ts(
              start = c(2013, 1), end = c(2015, 52), frequency = 52) %>%
            Arima(
              order = c(2, 0, 0), seasonal = c(1, 0, 0))

# Object to store simulated peak values
peak <- NULL

# Simulate peak value
for (i in 1:10000) {

  # Generate random values of parameter estimates from normal distribution
  rand.phis <- mvrnorm(1, usflu_Arima$coef, usflu_Arima$var.coef)

  # Create value for constant
  c <- rand.phis[4]*(1-rand.phis[1]-rand.phis[2]-rand.phis[3]+
                  rand.phis[1]*rand.phis[3]+
                  rand.phis[2]*rand.phis[3])

  # Extract previous data used to model new values
  y <- usflu_Arima$x[(length(usflu_Arima$x)-53):length(usflu_Arima$x)]

  # Model one year out
    for(k in 55:106){
      y[k] =  c + rand.phis[1]*y[k-1] +
        rand.phis[2]*y[k-2] +
        rand.phis[3]*y[k-52] - 
        rand.phis[1]*rand.phis[3]*y[k-53] -
        rand.phis[2]*rand.phis[3]*y[k-54]
    }

  # Keep predicted max flu value
  peak <- rbind(peak, max(y[55:106]))
}

summary(peak)
# V1       
# Min.   :2.055  
# 1st Qu.:2.669  
# Median :2.797  
# Mean   :2.786  
# 3rd Qu.:2.915  
# Max.   :3.325 

 A: Your idea sounds sensible to me. It reminds me the idea of resampling or bootstrapping, i.e., generating replicates of the data in order to approximate the distribution of a sample statistic, e.g., the mean, variance,... You can try it on order to resample the forecasts.
A common approach for bootstrapping is to generate data from the chosen model using the estimated coefficients (which are kept fixed instead of drawing from their distribution) and resampling the residuals. Assuming that the chosen model leads to white noise residuals, they could be resampled for example as sample(residuals, replace=TRUE). In your context, you could generate replicates of the data in this way and then fit an ARIMA model to it and get forecasts.
The code and the figure below sketches this idea.
# fit ARIMA(2,0,0)(1,0,0) model to the original data
# and store residuals
x <- usflu_Arima$x
n <- length(x)
fit <- arima(x, order=c(2,0,0), seasonal=list(order=c(1,0,0)))
resid <- residuals(fit)
# compute and store also forecasts for completeness
fcast <- forecast(fit, 12)
flow <- ts(fcast$lower[,2], frequency=52, start=2016)
fup <- ts(fcast$upper[,2], frequency=52, start=2016)

# storage matrix and main loop
peaks <- matrix(NA, nrow=1000, ncol=3)
set.seed(500)
for (i in seq_len(nrow(peaks)))
{
  # resample the residuals and generate bootstrap replicates 
  # upon the fitted model
  be <- sample(resid, size=n, replace=TRUE)
  bx <- stats::filter(be, method="recursive", 
    filter=-c(-coef(fit)[1:2], rep(0,49), -coef(fit)[3], coef(fit)[1]*coef(fit)[3], coef(fit)[2]*coef(fit)[3]))
  # add the estimated mean
  bx <- bx + coef(fit)[4]
  tsp(bx) <- tsp(x)

  # fit the same ARIMA model to the bootstrapped data, 
  # get one-year-ahead forecasts and record date and value of peaks
  bfit <- try(arima(bx, order=c(2,0,0), seasonal=list(order=c(1,0,0))))
  if (inherits(bfit, "Arima"))
  {
    fcast <- forecast(bfit, 12)
    peaks[i,1] <- which.max(fcast$mean)
    peaks[i,2] <- time(fcast$mean)[peaks[i,1]]
    peaks[i,3] <- fcast$mean[peaks[i,1]]

    # display results
    plot(cbind(x, flow, fup), plot.type="single", type="n", ylim=c(0,6), xlab="", ylab="flu data")
    mtext(text=paste("iteration", i), side=3)
    lines(x, col="gray")
    lines(fcast$mean, col="blue")
    lines(flow, col="red", lty=2)
    lines(fup, col="red", lty=2) 
    abline(v=c(2014, 2015, 2016), col="gray", lty=2)
    legend("topleft", bty="n", legend=c("original data", "bootstrapped data", "forecasts for bootstrapped data", "95% CI of forecasts for original data", "peaks"),
      col=c("gray", "black", "blue", "red", "red"), lty=c(1,1,1,2,NA), pch=c(NA,NA,NA,NA,16))
    lines(bx)
    points(peaks[,2], peaks[,3], col="red", pch=16)
  }
}


We can see that the bootstrapped series do not always keep faithfully the overall features of the series (sometimes there seem to cycles of shorter 
period than in the original data, sometimes negative values). The presence of outliers or the small size of the sample (3 observations for each week) may affect the quality of the results shown in this illustration.
I run the simulation for around 250 iterations. This is the summary of the times at which the peaks are forecasted.
summary(peaks[,1], na.action=na.remove)
#   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
#  1.000   1.000   7.000   6.458  11.000  12.000     736 

As shown in the plot below 
plot(density(na.omit(peaks[,1])), main="")

the distribution of these peaks is bimodal, with peaks at around weeks 2 and 11. This is not what we may have expected, since the annual paths of the original data suggest peaks at around weeks 12 and 23.


Alternative approaches
You may be interested in other possible approaches and methodologies. 
I'm not very familiar with them, so I won't give any examples.
Extreme value theory 
studies the probability of events that are more extreme.
Functional data analysis analyses curves  that are taken as a single observation. For example, in this case you would have three units of observations (three years or curves, each one consisting of 52 points). Nevertheless, three years of data may be a small sample.
Tagging your question with extreme-value and functional-data-analysis may result in some ideas in this sense.
A: We have implemented the boot-strapping procedure within AUTOBOX even allowing multivariable boot-strapping where predictor series are also bootstrapped providing more honest confidence limits reflecting the uncertainty in exogenous stochastic series. Additionally we have fully implemented incorporating pulse anomalies to be included in the forecast distribution providing safety levels for exceptional data to be addressed and peak levels to be targeted. The confidence limits generated from this approach are quite reasonable and do approximate the theoretical/standard limits exactly when the residuals are "nearly normal" . No longer do we have the naive symmetrical limit assumption.
