I am using the forecast package by Prof Hyndman, and have had success fitting ARIMA models to excess mortality (from the COVID-19 pandemic) data. I am currently trying to produce plots for cumulative excess deaths, but am unsure how to proceed about producing prediction intervals for these plots.
In more detail, I am fitting an ARIMA model using monthly mortalaity (rate) data from January 2000 to Feburary 2020 via the auto.arima()
function. I then use forecast()
to work out the counterfactual deaths for March 2020- March 2021 (deaths that would have occurred in the absence of the pandemic). This also gives me prediction intervals for the counterfactual deaths.
However, I want to produce a plot for cumulative excess death rates starting from March 2020. If this were a normal linear model, then to work out the confidence intervals I would simply use the following:
data$cumse[1]<-sqrt(data$se[1]^2)
for(i in 2:nrow(data)){
data$cumse[i]<-sqrt(data$cumse[i-1]^2 + data$se[i]^2)
}
where data$se
is the standard error for my counterfactual deaths for March 2020-March 2021, and data$cumse
would be the standard error for the cumulative excess deaths. To then produce the prediction intervals, I would simply take the bounds to be data$cumexcess
$\pm$1.96*data$cumse
.
However, my ARIMA model contains AR terms so I don't think the above would be appropriate as it assumes independence. Please could someone help me work out how to produce prediction intervals in this scenario?