# Forecast accuracy metric that involves prediction intervals

I'm in the process of generating a time series forecast for a company's product revenue and am looking for some way to show accuracy over time - e.g. after say 6 months they want to see how the actual revenue compared to the forecasted revenue generated 6 months ago.

I'm generating a forecast using the R ets() package and making predictions for each month over the next 6 months, including prediction intervals.

Are there any forecast accuracy metrics that take these prediction intervals into account?

I know of the standard MAPE, MASE etc but these all apply to point forecasts. What I'm looking for is a measure that also takes into account how accurate the prediction intervals are - e.g. if we're generating 95% prediction intervals, but the actual value only appears in them 10% of time, I want to be able to identify this.

• Another scenario might be that we're comparing two forecasting models. Model 1 has closer point forecasts than Model 2, but the prediction intervals from Model 1 are very narrow and hardly ever cover the actual value, while Model 2's do. All the accuracy measures that I know would point towards Model 1 as the best, but in practice Model 2 might be more appropriate. – Fergus Feb 9 '16 at 1:33
• It should be possible to compute the coverage probability. That value is a Bernoulli trial & you can compute CIs for a binomial to see if the coverage seems right (you might be interested in this: Confidence interval for the confidence interval?). – gung - Reinstate Monica Feb 9 '16 at 1:44
• What's wrong with MASE? – Aksakal Feb 9 '16 at 1:44
• Sounds a bit like you're wanting to check calibration. If you generate 25%, 50%, 75%, and 95% prediction intervals (say), you should find about 25%, 50%, 75%, and 95% of the actuals within those bounds. So I'd suggest searching on "model calibration" to get ideas. – Wayne Feb 9 '16 at 2:16
• @Aksakal: the Mean Absolute Scaled Error (MASE) is a point forecast accuracy measure. It does not evaluate intervals. – Stephan Kolassa Feb 9 '16 at 7:34

This is a good question. Unfortunately, while the academic forecasting literature is indeed (slowly) moving from an almost exclusive emphasis on point forecasts towards interval forecasts and predictive densities, there has been little work on evaluating interval forecasts. (EDIT: See the bottom of this answer for an update.)

As gung notes, whether or not a given 95% prediction interval contains the true actual is in principle a Bernoulli trial with $$p=0.95$$. Given enough PIs and realizations, you can in principle test the null hypothesis that the actual coverage probability is in fact at least 95%.

However, you will need to think about statistical power and sample sizes. It's probably best to decide beforehand what kind of deviation from the target coverage is still acceptable (is 92% OK? 90%?), then find the minimum sample size needed to detect a deviation that is this strong or stronger with a given probability, say $$\beta=80\%$$, which is standard. You can do this by a straightforward simulation: simulate $$n$$ Bernoulli trials with $$p=0.92$$, estimate $$\hat{p}$$ with a confidence interval for it, see whether it contains the value 95%, do this "often", and tweak $$n$$ until 95% is outside the CI in $$\beta=80\%$$ of your cases. Or use any Bernoulli power calculator.

OK, now that we have our sample size, you can batch your PIs and realizations in batches of this sample size, see how often your PIs contain the true realization, and start testing. Your batches can be the last $$n$$ PI/realizations of a single time series, or all the latest PIs/realizations of a large number of time series you are forecasting, or whatever.

This approach has the advantage of being rather easy to explain and to understand. Of course, if you have a large number of trials, even small deviations from the target coverage will be statistically significant, which is why you'll need to think about what deviation actually is significant from a business perspective, as per above.

Alternatively, quantile forecasts (say, a 2.5% and a 97.5% quantile forecast, to yield a 95% PI) arise naturally as optimal point forecasts under certain loss functions, which are parameterized based on the target quantile. This paper gives a nice overview. This may be an alternative to the Bernoulli tests above: find the correct loss function for your target upper and lower quantile, then evaluate the two endpoints of your PIs under these loss functions. However, the loss functions are rather abstract and not easily understood, especially for nontechnical audiences.

If you are comparing, say, multiple forecasting methods, you could first discard those whose PIs significantly underperform, based on Bernoulli hypothesis tests or loss functions, then assess the ones that passed this initial screening based on the width of their PIs. Among two PIs with the same correct coverage rate, the narrower one is usually better.

For a simple evaluation of PIs using null hypothesis significance tests, see this paper. There are also some far more elaborate schemes for evaluating PIs, which can also deal with serial dependence in deviations in coverage (maybe your financial PIs are good some part of the year, but bad at specific times), like this paper and that paper. Unfortunately, these require quite a large number of PI/realizations and so are likely only relevant for high-frequency financial data, like stock prices reported multiple times per day.

Finally, there has recently been some interest in going beyond PIs to the underlying predictive densities, which can be evaluated using (proper) scoring rules. Tilmann Gneiting has been very active in this area, and this paper of his gives a good introduction. However, even if you do decide to go deeper into predictive densities, scoring rules are again quite abstract and hard to communicate to a nontechnical audience.

EDIT - an update:

Your quality measure needs to balance coverage and length of the prediction intervals: yes, we want high coverage, but we also want short intervals.

There is a quality measure that does precisely this and has attractive properties: the interval score. Let $$\ell$$ and $$u$$ be the lower and the upper end of the prediction interval. The score is given by

$$S(\ell,u,h) = (u-\ell)+\frac{2}{\alpha}(\ell-h)1(h<\ell)+\frac{2}{\alpha}(h-u)1(h>u).$$

Here $$1$$ is the indicator function, and $$\alpha$$ is the coverage your algorithm is aiming for. (You will need to prespecify this, based on what you plan on doing with the prediction interval. It makes no sense to aim for $$\alpha=100\%$$ coverage, because the resulting intervals will be too wide to be useful for anything.)

You can then average the interval score over many predictions. The lower the average score, the better. See Gneiting & Raftery (2007, JASA)] for a discussion and pointers to further literature. A scaled version of this score was used, for instance, in assessing predictions intervals in the recent M4 forecasting competition.

(Full disclosure: this was shamelessly cribbed from this answer of mine.)

As nicely reflected in Stephan Kolassa's answer, there is a large academic literature on forecast evaluation. Let me add another idea from that literature, which I think matches your problem well.

Suppose your prediction interval has lower limit $p_l$ and upper limit $p_u$, and the goal is to design a central prediction interval at level $\alpha$ (i.e. you allow for $\alpha$/2 percent of the observations to fall below $p_l$, and for $\alpha$/2 percent to fall above $p_u$).

Given a realization $y$, a statistical score for the prediction interval is given by $$\text{Score} = (p_u - p_l) + \frac{2}{\alpha}(p_l - y)\mathbb{1}(y < p_l) + \frac{2}{\alpha}(y - p_u)\mathbb{1}(y > p_u),$$ where $\mathbb{1}$ is the indicator function; see Section 6.2 in Gneiting and Raftery ("Strictly Proper Scoring Rules, Prediction and Estimation", Journal of the American Statistical Association, 2007). Importantly, note that a smaller score corresponds to a better prediction interval.

The intuition behind the score formula is as follows: In the best of all worlds, we would want to have very short prediction intervals ($p_u$ only slightly larger than $p_l$), which nevertheless cover the realization $y$ ($p_l \le y \le p_u$). In this ideal world, the first summand of the score formula would be close to zero, whereas the second and third summands would be exactly zero.

Of course, in practice there is a tradeoff between producing short prediction intervals and actually covering the realization $y$. The score above strikes a balance between those two goals.

This seems like something you could accomplish with a dataframe--say you have your confidence intervals at each day--that is two columns, with day rows: an upper estimate, and a lower estimate. You could also just use the mean and std deviation of the gaussian with that day's confidence interval as your two columns for each day.

Then, you just store another column of actual measurements of product revenue.

Next, you load just the dates with estimation data and actual data, take the difference between the mean of the gaussian and the actual data, divide by the standard deviation, and use that number to 'rank' the effectiveness of your predictions.

If it is within ~2 standard deviations ($\frac{|predicted mean - actual|}{\sigma} \leq 2$), then it is within your confidence interval.

So the key problems are:

using a dataframe in R, storing the data (database, etc.), choosing the data properly...just the mechanics.