# How do I calculate estimated variance for an ensemble forecast?

I have several (n) different forecasts of comparable quality for a variable, based on the same data but using wildly different statistical models. For each, I have generated an estimate for m periods in the future, using the standard R tools for producing both a point estimate and a confidence interval for each forecast. I wish to combine these estimates into a single point estimate for m periods ahead, and a single confidence interval for the combined point estimate.

If the errors on these estimates each pass a normality test, my inclination is to treat the combined estimate as if it were the mean of n normal variables, and calculate the point estimate and a confidence interval for that mean accordingly.

Is this a reasonable methodology? Is there a better one?

This is my first venture into ensemble estimation, and I am pretty much ignorant of whatever results are in the literature. So do not hesitate to belabor the obvious

• note that, even with the normal approximation, you need to consider correlations between the predictors in calculating the variance of the average Commented May 25 at 9:28

That's not generally what people mean when they refer to an ensemble forecast.

What you have, abstractly, is a collection of predictive distributions for future paths, each conditional on the history, and on the correct model $$M$$ being that specific model:

$$p(Y_{t+1},...,Y_{t+m}|Y_1,...,Y_t; M = i), \quad i=1,...,n$$

What you're after is the predictive distribution but not conditioned on any particular model; you want to hedge your bets and not believe in only one model. To do this you need to integrate out the model index:

$$p(Y_{t+1},...,Y_{t+m}|Y_1,...,Y_t)=\sum_{i=1}^np(Y_{t+1},...,Y_{t+m}|Y_1,...,Y_t; M = i)\times p(M=i|Y_1,...,Y_t)$$

The ensemble predictive distribution is a weighted average of the predictive distributions for each model, where the weights are the model posteriors.

Now, this does imply that the ensemble predictive distribution's mean is the same weighted average of the means, so in a point forecasting exercise, you can indeed average the point forecasts. But it does not imply that this is true for the higher moments.

Computing results from the ensemble distribution is most easily done via simulation: draw an integer $$i$$ between 1 and $$n$$, then draw a sample from the $$i$$-th model. Repeat these two steps $$M$$ times. This yields a sample of $$M$$ paths from the ensemble.

Here's an example. Say you're forecasting $$m=1$$ periods ahead and you have two models, both with normal predictive distributions:

$$p(Y_{t+1}|Y_t, M=1) \sim \mathcal{N}(1, 1)$$ $$p(Y_{t+1}|Y_t, M=2) \sim \mathcal{N}(-1, 10)$$

Your approach will always yield a distribution that is normal in that case. But the actual ensemble distribution does not need to be, and here it is highly skewed (it may also be multimodal, fat tailed, etc...):

Ensemble methods can sometimes be used to enable the use of more complex models relatively cheaply by combining simple models that each represent different approximations of the more complex model. If you were limited to always producing normal predictive distributions it wouldn't be as useful.

Some other notes:

• You mention confidence intervals, but you probably mean prediction intervals.
• It can be difficult to get at the predictive distribution for each model. In a time series context it is often substituted by the plug-in estimate (i.e. collapse the parameter posterior at the point estimate $$\hat{\theta}$$):

\begin{align} p(Y_{t+1},...,Y_{t+m}|Y_1,...,Y_t; M = i)&=\int p(Y_{t+1},...,Y_{t+m}|Y_1,...,Y_t; M = i, \theta)p(\theta|Y_1,...,Y_t; M = i)d\theta \\&\approx p(Y_{t+1},...,Y_{t+m}|Y_1,...,Y_t; M = i, \hat{\theta})\end{align}

• It can also be difficult to estimate the model posterior weights reliably. In practice, in many contexts taking the weights to be uniform can work well.