R/etsforecast.R

In ttnsdcn/forecast-package: Forecasting functions for time series

forecast.ets <- function(object, h=ifelse(object$m>1, 2*object$m, 10),
    level=c(80,95), fan=FALSE, simulate=FALSE, bootstrap=FALSE, npaths=5000, PI=TRUE, lambda=object$lambda, ...)
{
    # Check inputs
    if(h>2000 | h<=0)
        stop("Forecast horizon out of bounds")
    if(!PI)
    {
      simulate <- bootstrap <- fan <- FALSE
      npaths <- 2 # Just to avoid errors
      level <- 90
    }
    if(fan)
        level <- seq(51,99,by=3)
    else
    {
        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")
    }

    n <- length(object$x)
    damped <- as.logical(object$components[4])

    if(simulate)
        f <- pegelsfcast.C(h,object,level=level,bootstrap=bootstrap,npaths=npaths)
    else if(object$components[1]=="A" & is.element(object$components[2],c("A","N")) & is.element(object$components[3],c("N","A")))
        f <- class1(h,object$states[n+1,],object$components[2],object$components[3],damped,object$m,object$sigma2,object$par)
    else if(object$components[1]=="M" & is.element(object$components[2],c("A","N")) & is.element(object$components[3],c("N","A")))
        f <- class2(h,object$states[n+1,],object$components[2],object$components[3],damped,object$m,object$sigma2,object$par)
    else if(object$components[1]=="M" & object$components[3]=="M" & object$components[2]!="M")
        f <- class3(h,object$states[n+1,],object$components[2],object$components[3],damped,object$m,object$sigma2,object$par)
    else
        f <- pegelsfcast.C(h,object,level=level,bootstrap=bootstrap,npaths=npaths)
    if(!PI)
      f$var <- f$lower <- f$upper <- NULL

    tsp.x <- tsp(object$x)
    if(!is.null(tsp.x))
        start.f <- tsp(object$x)[2] + 1/object$m
    else
        start.f <- length(object$x)+1
    out <- list(model=object,mean=ts(f$mu,f=object$m,s=start.f),level=level,x=object$x)
    if(!is.null(f$var))
    {
        out$lower <- out$upper <- matrix(NA,ncol=length(level),nrow=h)
        for(i in 1:length(level))
        {
            marg.error <- sqrt(f$var) * abs(qnorm((100-level[i])/200))
            out$lower[,i] <- out$mean - marg.error
            out$upper[,i] <- out$mean + marg.error
        }
    }
    else if(!is.null(f$lower))
    {
        out$lower <- ts(f$lower,f=object$m,s=start.f)
        out$upper <- ts(f$upper,f=object$m,s=start.f)
    }
    else if(PI)
        warning("No prediction intervals for this model")
		
  out$fitted <- fitted(object)
  out$method <- object$method
  out$residuals <- residuals(object)

  if(!is.null(lambda))
  {
	#out$x <- InvBoxCox(object$x,lambda)
	#out$fitted <- InvBoxCox(out$fitted,lambda)
    out$mean <- InvBoxCox(out$mean,lambda)
    out$lower <- InvBoxCox(out$lower,lambda)
    out$upper <- InvBoxCox(out$upper,lambda)
  }
	
    return(structure(out,class="forecast"))
}

pegelsfcast.C <- function(h,obj,npaths,level,bootstrap)
{
    y.paths <- matrix(NA,nrow=npaths,ncol=h)
    for(i in 1:npaths)
        y.paths[i,] <- simulate.ets(obj, h, future=TRUE, bootstrap=bootstrap)
    y.f <- .C("etsforecast",
            as.double(obj$state[length(obj$x)+1,]),
            as.integer(obj$m),
            as.integer(switch(obj$components[2],"N"=0,"A"=1,"M"=2)),
            as.integer(switch(obj$components[3],"N"=0,"A"=1,"M"=2)),
            as.double(ifelse(obj$components[4]=="FALSE",1,obj$par["phi"])),
            as.integer(h),
            as.double(numeric(h)),
        PACKAGE="forecast")[[7]]
    if(abs(y.f[1]+99999) < 1e-7)
        stop("Problem with multiplicative damped trend")

    lower <- apply(y.paths,2,quantile,0.5 - level/200,type=8)
    upper <- apply(y.paths,2,quantile,0.5 + level/200,type=8)
    if(length(level)>1)
    {
        lower <- t(lower)
        upper <- t(upper)
    }
    return(list(mu=y.f,lower=lower,upper=upper))
}

class1 <- function(h,last.state,trendtype,seasontype,damped,m,sigma2,par)
{
    p <- length(last.state)
    H <- matrix(c(1,rep(0,p-1)),nrow=1)
    if(seasontype=="A")
        H[1,p] <- 1
    if(trendtype=="A")
        H[1,2] <- 1
    F <- matrix(0,p,p)
    F[1,1] <- 1
    if(trendtype=="A")
    {
        F[1,2] <- 1
        if(damped)
            F[2,2] <- par["phi"]
        else
            F[2,2] <- 1
    }
    if(seasontype=="A")
    {
        F[p-m+1,p] <- 1
        F[(p-m+2):p,(p-m+1):(p-1)] <- diag(m-1)
    }
    G <- matrix(0,nrow=p,ncol=1)
    G[1,1] <- par["alpha"]
    if(trendtype=="A")
        G[2,1] <- par["beta"]
    if(seasontype=="A")
        G[3,1] <- par["gamma"]
    mu <- numeric(h)
    Fj <- diag(p)
    cj <- numeric(h-1)
    if(h>1)
    {
        for(i in 1:(h-1))
        {
            mu[i] <- H %*% Fj %*% last.state
            Fj <- Fj %*% F
            cj[i] <- H %*% Fj %*% G
        }
        cj2 <- cumsum(cj^2)
        var <- sigma2 * c(1,1+cj2)
    }
    else
        var <- sigma2
    mu[h] <- H %*% Fj %*% last.state

    return(list(mu=mu,var=var,cj=cj))
}

class2 <- function(h,last.state,trendtype,seasontype,damped,m,sigma2,par)
{
    tmp <- class1(h,last.state,trendtype,seasontype,damped,m,sigma2,par)
    theta <- numeric(h)
    theta[1] <- tmp$mu[1]^2
    if(h>1)
    {
        for(j in 2:h)
            theta[j] <- tmp$mu[j]^2 + sigma2 * sum(tmp$cj[1:(j-1)]^2*theta[(j-1):1])
    }
    var <- (1+sigma2)*theta - tmp$mu^2
    return(list(mu=tmp$mu,var=var))
}

class3 <- function(h,last.state,trendtype,seasontype,damped,m,sigma2,par)
{
    p <- length(last.state)
    H1 <- matrix(rep(1,1+(trendtype!="N")),nrow=1)
    H2 <- matrix(c(rep(0,m-1),1),nrow=1)
    if(trendtype=="N")
    {
        F1 <- 1
        G1 <- par["alpha"]
    }
    else
    {
        F1 <- rbind(c(1,1),c(0,ifelse(damped,par["phi"],1)))
        G1 <- rbind(c(par["alpha"],par["alpha"]),c(par["beta"],par["beta"]))
    }
    F2 <- rbind(c(rep(0,m-1),1),cbind(diag(m-1),rep(0,m-1)))

    G2 <- matrix(0,m,m)
    G2[1,m] <- par["gamma"]
    Mh <- matrix(last.state[1:(p-m)]) %*% matrix(last.state[(p-m+1):p],nrow=1)
    Vh <- matrix(0,length(Mh),length(Mh))
    H21 <- H2 %x% H1
    F21 <- F2 %x% F1
    G21 <- G2 %x% G1
    K <- (G2 %x% F1) + (F2 %x% G1)
    mu <- var <- numeric(h)
    for(i in 1:h)
    {
        mu[i] <- H1 %*% Mh %*% t(H2)
        var[i] <- (1+sigma2) * H21 %*% Vh %*% t(H21) + sigma2*mu[i]^2
        vecMh <- c(Mh)
        Vh <- F21 %*% Vh %*% t(F21) + sigma2 * (F21 %*% Vh %*% t(G21) + G21 %*% Vh %*% t(F21) +
            K %*% (Vh + vecMh %*% t(vecMh)) %*% t(K) + sigma2 * G21 %*% (3*Vh + 2*vecMh%*%t(vecMh))%*%t(G21))
        Mh <- F1 %*% Mh %*% t(F2) + G1 %*% Mh %*% t(G2) * sigma2
    }
    return(list(mu=mu,var=var))
}

ses <- function(x,h=10,level=c(80,95),fan=FALSE,...)
{
    fcast <- forecast(ets(x,"ANN"),h,level=level,fan=fan,...)
    fcast$method <- "Simple exponential smoothing"
    fcast$model$call <- match.call()
    return(fcast)
}

holt <- function(x,h=10, damped=FALSE, level=c(80,95), fan=FALSE, ...)
{
    junk <- forecast(ets(x,"AAN",damped=damped),h,level=level,fan=fan,...)
    if(damped)
        junk$method <- "Damped Holt's method"
    else
        junk$method <- "Holt's method"
    junk$model$call <- match.call()
    return(junk)
}

hw <- function(x,h=2*frequency(x),seasonal="additive",damped=FALSE,level=c(80,95), fan=FALSE, ...)
{
    if(seasonal=="additive")
    {
        junk <- forecast(ets(x,"AAA",damped=damped),h,level=level,fan=fan,...)
        junk$method <- "Holt-Winters' additive method"
    }
    else
    {
        junk <- forecast(ets(x,"MAM",damped=damped),h,level=level,fan=fan,...)
        junk$method <- "Holt-Winters' multiplicative method"
    }
    junk$model$call <- match.call()
    return(junk)
}