3
$\begingroup$

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?

$\endgroup$
4
$\begingroup$

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/

$\endgroup$
  • $\begingroup$ Thank you for sharing this! I have a small additional question about this method: What if I have data for the first 3 months of 2013? Should I go with y <- filter(x,rep(1,12), sides=1); fit <- auto.arima(y); forecast(fit,h=9)$mean[9] or y <- filter(x,rep(1,9), sides=1); fit <- auto.arima(y); forecast(fit,h=9)$mean[9]+sum(x[(length(x)-2):length(x)])? $\endgroup$ – Nikita Samoylov May 15 '13 at 14:02
  • 1
    $\begingroup$ The first option with the same filter command, but grabbing the ninth forecast. $\endgroup$ – Rob Hyndman May 15 '13 at 22:25
  • $\begingroup$ Does not the first option ignore the fact that some of the forecast (the first 3 months) has already eventuated? $\endgroup$ – Nikita Samoylov May 16 '13 at 0:33
  • $\begingroup$ No. It uses the actual values when they exist and adds in the forecast values when they don't. $\endgroup$ – Rob Hyndman May 16 '13 at 0:36
1
$\begingroup$

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)
)
$\endgroup$
  • $\begingroup$ Are you able to give a brief description of what the approach involves? Many people find it easier to comprehend what code is doing if they have an outline in mind first. $\endgroup$ – Glen_b -Reinstate Monica May 16 '13 at 8:09
  • 1
    $\begingroup$ The above method uses the patterns of "how the yearly total has been reached in previous years" (e.g from monthly data) as an estimation of current path. The method adjust the estimated total as data come in (at any sub time period). Essentially It has an intra year path with an estimated percentage over total at any sub time period. $\endgroup$ – D. Amberti Jun 3 '13 at 7:46

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.