Forecast total for a year given monthly time series

Question

I have a monthly time series (for 2009-2012 non-stationary, with seasonality). I can use ARIMA (or ETS) to obtain point and interval forecasts for each month of 2013, but I am interested in forecasting the total for the whole year, including prediction intervals. Is there an easy way in R to obtain interval forecasts for the total for 2013?

Answer 1

Here is a trick I've used before, although I don't think I've ever published it. If x is your monthly time series, then you can construct annual totals as follows.

y <- filter(x,rep(1,12), sides=1) # Total of previous 12 months

To get the forecasts of the annual totals:

library(forecast)
fit <- auto.arima(y)
forecast(fit,h=12)

The last forecast is for the total of the next year.

An extended version of this answer is at http://robjhyndman.com/hyndsight/forecasting-annual-totals/

Answer 2

I suggest that you use what is described in

Mendoza, M. and E. de Alba (2006),
"Forecasting an Accumulated Series Based on Partial Accumulation II: A New Bayesian Method for Short Series With Stable Seasonal Patterns",
International Journal of Forecasting, 2006, vol. 22, issue 4, pages 781-798.

I wrote some R code for that a couple of years ago. I'd like to add this to Rob's forecast package, but I didn't have time, may be in the future (Rob, what do you think about this?). Please refer to me as the author of code below and note that the code is unfinished.

  pas <-  function(y){

  if (class(y) == "data.frame" | class(y) == "list" | class(y) == 
    "matrix" | is.element("mts", class(y))) 
    stop("y should be a univariate time series")

  y <- as.ts(y)
  freq <- frequency(y)
  if(freq <= 1)
    stop("function need y frequency > 1")

  l <- floor(length(y)/freq)
  u <- ceiling(length(y)/freq)

  if(l != u)
    stop("historical data should represent complete seasons")
  if(l < 2)
    stop("historical data should have at least two complete seasons")

  pattern <- function(x, method = "manhattan", ...){
    P <- colCumsumsMatrix(x)/matrix(
      colSums(x), 
      nrow(x), ncol(x), byrow= T)
    return(P)
  }

  colCumsumsMatrix <- function(x, na.rm = FALSE, ...){
    if(!is.matrix(x)) stop("x is not a matrix")
    if (na.rm)
      x <- na.omit(x)
    ans <- apply(x, 2, cumsum, ...)
    # special treatment when x has one row because apply returns a vector
    if (NROW(x) == 1)
      ans <- matrix(ans, nrow = 1, dimnames = dimnames(x))
    ans
  }

  x <- matrix(data = as.numeric(y),  ncol = l, byrow = FALSE)

  P <- pattern(x)
  P[P == 0] <- 2e-100 #because 1/0 == Inf and 1/2e-309 == Inf

  W <- 1/P-1 #1/0 == Inf

  Sd <- apply(log(W), 1, sd)
  Sd[length(Sd)] <- 0

  model <- list(w = W, sd = Sd, pattern = P, 
                frequency = freq, cumseries = ts(colSums(x), start= start(y), frequency=1))

  return(structure(model, class = "pas"))
}

plot.pas <- function(x, method = "manhattan", ...){

    P <- x$pattern
    D <- dist(t(P)/colSums(P), method = method)
    plot(as.ts(P), plot.type="single", 
         col = c("black", rainbow(ncol(P))),
         main ="Patterns", ylab = "%", ...)
    legend("topleft", legend=paste(1:ncol(P), 
                                   "dist", 
                                   c("ref", round(as.vector(D),3))),
           fill = c("black", rainbow(ncol(P))), ...)
    invisible(D)
}


forecast.pas <- function(object, newdata, level = c(5, 20, 80, 95)
                         #, onlylast = TRUE
                         ){

  if (class(newdata) == "data.frame" | class(newdata) == "list" | class(newdata) == 
    "matrix" | is.element("mts", class(newdata))) 
    stop("newdata should be a univariate time series")
  if(frequency(newdata) != object$frequency)
    stop("newdata should have the frequancy specified in object")

  if (min(level) > 0 & max(level) < 1) 
    level <- 100 * level
  else if (min(level) < 0 | max(level) > 99.99) 
    stop("Confidence limit out of range")

  level <- sort(unique(c(level, 50)))/100

  frcFun <- function(x, y, w, p, sd){
    sum(y[1:x]) + 
      sum(y[1:x])*
      ((prod(w[x,]))^(1/(ncol(w))))*
      (exp(qt(p, df = ncol(w) -1 )))^((1+(1/ncol(w)))^(1/2)*sd[x])
  }

#   if(onlylast == FALSE){
#     frc <- rep(as.double(NA), length(newdata) )  
#     frc <- sapply(1:length(newdata), frcFun, y = newdata, w = object$w, p = 0.5, sd = object$sd)    
#   } else {
     frc <- rep(as.double(NA), length(level))
     for(j in 1:length(level)){
       frc[j] <- sapply(length(newdata), frcFun, y = newdata, w = object$w, p = level[j], sd = object$sd)   
     }


#   }
  return(frc)

}


KwhIowa <- c(523, 502, 439, 420, 387, 453, 
             630, 637, 576, 411, 455, 512, 
             530, 507, 436, 407, 392, 531, 
             710, 658, 500, 414, 418, 520,
             535, 503, 464, 414, 383, 472, 
             676, 622, 652, 474, 422, 501)
KwhIowaTs <- ts(KwhIowa, frequency= 12)
Pas <- pas(y=window(x=KwhIowaTs, start =1, end =2.999))
plot(Pas)
library(forecast)
forecast(Pas, window(x=KwhIowaTs, start =3, end = 3.5)
         #, onlylast=T)
)