Prediction intervals for mixture models for time series forecasting - is it really an average of the prediction intervals of the averaged models?

Question

I'm trying to find out how to do forecasting with a mixture model (averaging the forecasts of an ets, an arima and an stlf model). I do not have a huge amount of statistics experience and so I'm struggling with finding out how to do it.

The point forecasts will just be the average of the point forecasts of the three methods, no problem.

The problem is how to calculate the prediction intervals.

I have found an R script with an attempt to do it, but the mixture prediction intervals are just calculated as an average of the prediction intervals of the models, and I am pretty sceptical about this approach - is it really that easy?

If not, how do I go about calculating them?

Answer 1

You are correct that you cannot simply average the quantiles of the constituent distributions.

As a little example, let's assume that we have a mixture of two normal distributions with means 0 and 1, both with standard deviations of 1, and with equal probability. Below is a histogram, the true 90% quantile (the vertical green line) and the equally weighted average of the separate 90% quantiles for the two normals. As you see, the two vertical lines are very different indeed.

means <- c(0,3)
sds <- c(1,1)
probs <- c(0.5,0.5)
nn <- 1e4

set.seed(1)
sims <- unlist(mapply(function(pp,mm,ss)rnorm(floor(pp*nn),mm,ss),probs,means,sds))

qq <- 0.9
quantile(sims,qq)
# 3.854995
sum(probs*qnorm(qq,means,sds))
# 2.781552

hist(sims,las=1,xlab="")
abline(v=c(quantile(sims,qq),sum(probs*qnorm(qq,means,sds))),col=c("green","red"),lwd=2)

So, what can we do? In the case of a gaussian-mixture as here, we should be able to use the qmixnorm() function in the KScorrect package for R (as per the earlier thread Compute quantile function from a mixture of Normal distribution, which I very much recommend):

library(KScorrect)
qmixnorm(qq,means,sds,probs)
# 3.841839

This is close enough to what we got empirically that I trust it.

However, this currently gives a warning which is passed through from stats::spline(). I already submitted a bug report here and asked the maintainer to leave a comment to this thread if/when this is addressed.

Alternatively, you could get quantiles of Gaussian mixtures by a bisection search on the convex combination of CDFs per this earlier thread. Or you could draw a large number of samples from the mixture and take empirical quantiles.

The other problem is that the forecasting tools you mention apparently use $t$ distributions for their predictive densities, not normal ones. So what you are really looking for is quantiles of mixtures of $t$ distributions. There doesn't seem to be anything analytical on that. So it's probably again easiest to do the bisection search, using the $t$ CDF, or the sampling approach.

Answer 2

The warning isn't an error. The warning is caused by the spline approximation function used within 'qmixnorm', which appears to have had a change in the default behavior an an additional internal function, 'regularize.values', that triggers this warning. The behavior shouldn't alter the way 'qmixnorm' works.

Interestingly, using a smaller, non-default 'expand' argument value (which changes how the spline approximation handles probabilities close to 0 or 1) doesn't trigger the warning.

qmixnorm(0.9, c(0, 3), c(1, 1), c(0.3, 0.7), expand = 0.5)

You can get qmixnorm not trigger the warning in this case by modifying the 'expand' argument to a smaller value.

Here's an example:

library(KScorrect)
qmixnorm(qq, means, sds, probs, expand = 0.5)