Formula that R predict function uses to calculate intervals for exponential smoothing

I have been trying to figure out the exact formula that the R predict function uses to calculate prediction intervals for simple exponential smoothing. The prediction interval formula seems to vary according to the software used (Gretl is different from Minitab is different from SAS). I have tried to look at the predict code by simply typing it in without parentheses but R won't display it.

I spent the better part of a day trying to figure out R's formula. Does anyone know the exact formula that R uses? Also, I'd appreciate a link showing the derivation of the formula, as there seems to be some disagreement in the literature about what the "right" formula to use is.

• Are you solely talking about predict using the default features? You can specify different prediction intervals if you specify the parameters yourself so which one is it? Also, are you using predict on a linear model or.....? There are several different predict commands in R.
– user25658
Sep 16 '13 at 15:50
– John
Sep 16 '13 at 16:23
• ... to expand on BabakP's comment, if, for example, you're using lm to estimate the model, the predict function is predict.lm (that's what you'd type in w/o parentheses to see the code.) If you're using "glm", then predict.glm ... and so forth. Sep 16 '13 at 16:36
• I was thinking about "predict" as it's used in this calculator from Wessa.net This calculator feeds a HoltWinters object without trend and season to the predict function: ' p <- predict(fit, par1, prediction.interval=TRUE) " Sep 16 '13 at 20:14

You can see the code for hidden functions using package:::function. In this case:

> stats:::predict.HoltWinters
function (object, n.ahead = 1L, prediction.interval = FALSE,
level = 0.95, ...)
{
f <- frequency(object$x) vars <- function(h) { psi <- function(j) object$alpha * (1 + j * object$beta) + (j%%f == 0) * object$gamma * (1 - object$alpha) var(residuals(object)) * if (object$seasonal == "additive")
sum(1, (h > 1) * sapply(1L:(h - 1), function(j) crossprod(psi(j))))
else {
rel <- 1 + (h - 1)%%f
sum(sapply(0:(h - 1), function(j) crossprod(psi(j) *
object$coefficients[2 + rel]/object$coefficients[2 +
(rel - j)%%f])))
}
}
fit <- rep(as.vector(object$coefficients[1L]), n.ahead) if (!is.logical(object$beta) || object$beta) fit <- fit + as.vector((1L:n.ahead) * object$coefficients[2L])
if (!is.logical(object$gamma) || object$gamma)
if (object$seasonal == "additive") fit <- fit + rep(object$coefficients[-(1L:(1 + (!is.logical(object$beta) || object$beta)))], length.out = length(fit))
else fit <- fit * rep(object$coefficients[-(1L:(1 + (!is.logical(object$beta) ||
object\$beta)))], length.out = length(fit))
if (prediction.interval)
int <- qnorm((1 + level)/2) * sqrt(sapply(1L:n.ahead,
vars))
ts(cbind(fit = fit, upr = if (prediction.interval)
fit + int, lwr = if (prediction.interval)
fit - int), start = end(lag(fitted(object)[, 1], k = -1)),
frequency = frequency(fitted(object)[, 1]))
}


The forecast variances for the non-seasonal cases are based on the equivalent ARIMA model. The forecast variances for the seasonal cases are computed using the Chatfield-Yar formulae (Chatfield & Yar, IJF, 1990; Yar & Chatfield, IJF 1991). There should be no dispute. These are the "right" formula assuming the errors are additive for the non-seasonal and additive HW cases, and multiplicative for multiplicative HW. More general results plus derivations are given in my Springer book (Hyndman, Koehler, Ord & Snyder, 2008, ch6).

• Thanks very much. You're right, there is no theoretical dispute. But I know some software (Gretl and Minitab) gave me different values for the same data. I'm still new to R, so I didn't know of the ::: trick, so thanks for that as well. Sep 17 '13 at 0:41
• Check if the differences are due to different parameter values or different formulas. Different parameter values can arise due to different ways of initializing, and convergence to local optima. Sep 17 '13 at 1:30