# Forecasting excess mortality with ARIMA model

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?

If auto.arima handles excess rates $$\Delta x_t$$ and produces interval forecasts, it will trivially handle cumulative excess rates $$x_t$$ (by first-differencing them), too. Supply the cumulatively summed series $$x_t$$ to auto.arima, and you will get the forecast interval for it.