aTSA: Alternative Time Series Analysis

#' Holt's Two-parameter Exponential Smoothing
#' @description Performs Holt's two-parameter exponential smoothing for linear
#' trend or damped trend.
#' @param x a numeric vector or univariate time series.
#' @param type the type of interaction between the level and the linear trend. See 
#' details.
#' @param alpha the parameter for the level smoothing. The default is \code{0.2}.
#' @param beta the parameter for the trend smoothing. The default is \code{0.1057}.
#' @param lead the number of steps ahead for which prediction is required. 
#' The default is \code{0}.
#' @param damped a logical value indicating a damped trend. See details. The default is
#' \code{FALSE}.
#' @param phi a smoothing parameter for damped trend. The default is \code{0.98}, only valid
#' for \code{damped = TRUE}.
#' @param plot a logical value indicating to print the plot of original data v.s smoothed 
#' data. The default is \code{TRUE}.
#' @details Holt's two parameter is used to forecast a time series with trend, but 
#' wihtout seasonal pattern. For the additive model (\code{type = "additive"}), the 
#' \eqn{h}-step-ahead forecast is given by \eqn{hat{x}[t+h|t] = level[t] + h*b[t]},
#'  where
#' \deqn{level[t] = \alpha *x[t] + (1-\alpha)*(b[t-1] + level[t-1]),}
#' \deqn{b[t] = \beta*(level[t] - level[t-1]) + (1-\beta)*b[t-1],}
#' in which \eqn{b[t]} is the trend component.
#' For the multiplicative (\code{type = "multiplicative"}) model, the 
#' \eqn{h}-step-ahead forecast is given by \eqn{hat{x}[t+h|t] = level[t] + h*b[t]},
#'  where
#' \deqn{level[t] = \alpha *x[t] + (1-\alpha)*(b[t-1] * level[t-1]),}
#' \deqn{b[t] = \beta*(level[t] / level[t-1]) + (1-\beta)*b[t-1].} 
#' 
#' Compared with the Holt's linear trend that displays a constant increasing or 
#' decreasing, the damped trend generated by exponential smoothing method shows a 
#' exponential growth or decline, which is a situation between simple exponential
#' smoothing (with 0 increasing or decreasing rate) and Holt's two-parameter smoothing.
#' If \code{damped = TRUE}, the additive model becomes
#' \deqn{hat{x}[t+h|t] = level[t] + (\phi + \phi^{2} + ... + \phi^{h})*b[t],}
#' \deqn{level[t] = \alpha *x[t] + (1-\alpha)*(\phi*b[t-1] + level[t-1]),}
#' \deqn{b[t] = \beta*(level[t] - level[t-1]) + (1-\beta)*\phi*b[t-1].}
#' The multiplicative model becomes
#' \deqn{hat{x}[t+h|t] = level[t] *b[t]^(\phi + \phi^{2} + ... + \phi^{h}),}
#' \deqn{level[t] = \alpha *x[t] + (1-\alpha)*(b[t-1]^{\phi} * level[t-1]),}
#' \deqn{b[t] = \beta*(level[t] / level[t-1]) + (1-\beta)*b[t-1]^{\phi}.}
#' See Chapter 7.4 for more details in R. J. Hyndman and G. Athanasopoulos (2013).
#' @note Missing values are removed before analysis. 
#' @return A list with class "\code{Holt}" containing the following components:
#' \item{estimate}{the estimate values.}
#' \item{alpha}{the smoothing parameter used for level.}
#' \item{beta}{the smoothing parameter used for trend.}
#' \item{phi}{the smoothing parameter used for damped trend.}
#' \item{pred}{the predicted values, only available for \code{lead} > 0.}
#' \item{accurate}{the accurate measurements.}
#' @author Debin Qiu
#' @references R. J. Hyndman and G. Athanasopoulos, "Forecasting: principles and
#' practice," 2013. [Online]. Available: \url{http://otexts.org/fpp/}.
#' @seealso \code{\link{HoltWinters}}, \code{\link{expsmooth}}, \code{\link{Winters}}
#' @examples x <- (1:100)/100
#' y <- 2 + 1.2*x + rnorm(100)
#' 
#' ho0 <- Holt(y) # with additive interaction
#' ho1 <- Holt(y,damped = TRUE) # with damped trend
#' 
#' # multiplicative model for AirPassengers data, 
#' # although seasonal pattern exists.
#' ho2 <- Holt(AirPassengers,type = "multiplicative")
#' @importFrom stats is.ts
#' @importFrom stats ts
#' @importFrom graphics plot
#' @importFrom graphics lines
#' @export 
Holt <- function(x,type = c("additive","multiplicative"),
                 alpha = 0.2,beta = 0.1057,lead = 0,
                 damped = FALSE,phi = 0.98,plot = TRUE)
{
  if (NCOL(x) > 1)
    stop("'x' must be a numeric vector or univariate time series ")
  if (any(c(alpha,beta) > 1) || any(c(alpha,beta) < 0))
    stop("'alpha' or 'beta' must be between 0 and 1")
  if (phi > 1 || phi < 0)
    stop("'phi' must be between 0 and 1")
  type <- match.arg(type)
  if (any(!is.finite(x)))
    warning(paste("missing values exist at time",which(!is.finite(x)),
                  "and will be removed."))
  if (is.ts(x))
    x <- ts(x[is.finite(x)],start = time(x)[1],frequency = frequency(x))
  else
    x <- x[is.finite(x)]
  n <- length(x)
  if (n < 1L)
    stop("invalid length of 'x'")
  level <- c(x[1],numeric(n-1))
  trend <- numeric(n)
  trend[1] <- switch(type,additive = (x[n] - x[1])/n,
                     multiplicative = x[2]/x[1])
  x.hat <- c(x[1],numeric(n-1))
  phi <- ifelse(damped,phi,1)
  for (i in 2:n) {
    if (type == "additive") {
      level[i] <- alpha*x[i] + (1 - alpha)*(level[i-1] + phi*trend[i-1])
      trend[i] <- beta*(level[i] - level[i-1]) + (1 - beta)*phi*trend[i-1]
      x.hat[i] <- level[i-1] + phi*trend[i-1]
    }
    else {
      level[i] <- alpha*x[i] + (1 - alpha)*(level[i-1] * trend[i-1]^phi)
      trend[i] <- beta*(level[i] / level[i-1]) + (1 - beta)*trend[i-1]^phi
      x.hat[i] <- level[i-1] * trend[i-1]^phi
    }
  }
  if (is.ts(x))
    x.hat <- ts(x.hat,start = time(x)[1],frequency = frequency(x))
  result <- list(estimate = x.hat,alpha = alpha, beta = beta)
  if (lead > 0) {
    if (lead < 0 || lead%%1 != 0) 
      stop("'lead' must be a positive integer")
    l.s <- 1:lead
    x.pred <- switch(type, additive = level[n] + cumsum(phi^l.s)*trend[n],
                     multiplicative = level[n] + trend[n]^(cumsum(phi^l.s)))
    names(x.pred) <- n + 1:lead
    result <- c(result,list(phi = phi, pred = x.pred))
  }
  if (plot) {
    plot(x,main = "original v.s smoothed data",type = "l")
    lines(x.hat,col = 2)
  }
  k <- ifelse(damped,3,2)
  result <- c(result,list(accurate = accurate(x,x.hat,k,output = FALSE)))
  class(result) <- "Holt"
  return(result)
}