Aggregation of Interval Predictions Given a montly time series, my objective is to provide my client with the next 12 point forecasts along with a yearly forecast. To obtain the yearly forecast, I simply sum up the 12 points forecasts. However, I don't know how to estimate the prediction interval associated with this sum.
Any ideas on how to obtain a prediction interval for the yearly forecast based on the prediction intervals estimated for each month?
Thanks
 A: When creating the yearly forecasts by summing "up the 12 points forecasts", we actually sum the predicted trajectory across those 12 forward steps. It is important to emphasise this because it highlights that we do not have independent predictions. (i.e. our predictions for $\hat{y}_{t=5}$ is contingent to our $\hat{y}_{t=4}$ prediction, etc.) With that in mind, we want to create multiply trajectories (e.g. via bootstrapping) and then sum them up to get our final interval predictions for our yearly aggregates.
A quick R example can be found here:
library(forecast); 
ets_fit = ets(USAccDeaths, model="ZZZ"); 
horizon = 12; 
deaths.fcast = forecast(ets_fit, h=horizon);  
n_sim=5000;
set.seed(123);
sims = replicate(n_sim, sum(simulate(ets_fit, future=TRUE, nsim=horizon))); 
# Proportional difference between predicted mean and simulated mean
100*(mean(sims) - sum(deaths.fcast$mean)) / sum(deaths.fcast$mean) 

Python has a number of time-series libraries that are all great but none of them has emerged as the dominant one as forecast in R has. I have used PyFlux in the past. It has the ability to create posterior predictive checks for the mean prediction (e.g. look at plot_ppc() member methods which are available for most models) allowing you to specify nsims, etc. sktime also seems to have the ability to simulate a new time series following a state-space model but I don't see any worked examples in its documentation, seems possible though.
