R/ets.R

Defines functions is.ets nobs.ets logLik.ets hfitted.ets fitted.ets coef.ets print.summary.ets summary.ets plot.ets admissible pegelsresid.C print.ets lik initstate check.param initparam etsTargetFunctionInit etsmodel as.character.ets ets

Documented in as.character.ets coef.ets ets fitted.ets is.ets plot.ets print.ets summary.ets

#' Exponential smoothing state space model
#'
#' Returns ets model applied to \code{y}.
#'
#' Based on the classification of methods as described in Hyndman et al (2008).
#'
#' The methodology is fully automatic. The only required argument for ets is
#' the time series. The model is chosen automatically if not specified. This
#' methodology performed extremely well on the M3-competition data. (See
#' Hyndman, et al, 2002, below.)
#'
#' @aliases print.ets summary.ets as.character.ets coef.ets tsdiag.ets
#'
#' @param y a numeric vector or time series of class \code{ts}
#' @param model Usually a three-character string identifying method using the
#' framework terminology of Hyndman et al. (2002) and Hyndman et al. (2008).
#' The first letter denotes the error type ("A", "M" or "Z"); the second letter
#' denotes the trend type ("N","A","M" or "Z"); and the third letter denotes
#' the season type ("N","A","M" or "Z"). In all cases, "N"=none, "A"=additive,
#' "M"=multiplicative and "Z"=automatically selected. So, for example, "ANN" is
#' simple exponential smoothing with additive errors, "MAM" is multiplicative
#' Holt-Winters' method with multiplicative errors, and so on.
#'
#' It is also possible for the model to be of class \code{"ets"}, and equal to
#' the output from a previous call to \code{ets}. In this case, the same model
#' is fitted to \code{y} without re-estimating any smoothing parameters. See
#' also the \code{use.initial.values} argument.
#' @param damped If TRUE, use a damped trend (either additive or
#' multiplicative). If NULL, both damped and non-damped trends will be tried
#' and the best model (according to the information criterion \code{ic})
#' returned.
#' @param alpha Value of alpha. If NULL, it is estimated.
#' @param beta Value of beta. If NULL, it is estimated.
#' @param gamma Value of gamma. If NULL, it is estimated.
#' @param phi Value of phi. If NULL, it is estimated.
#' @param additive.only If TRUE, will only consider additive models. Default is
#' FALSE.
#' @param lambda Box-Cox transformation parameter. If \code{lambda="auto"},
#' then a transformation is automatically selected using \code{BoxCox.lambda}.
#' The transformation is ignored if NULL. Otherwise,
#' data transformed before model is estimated. When \code{lambda} is specified,
#' \code{additive.only} is set to \code{TRUE}.
#' @param lower Lower bounds for the parameters (alpha, beta, gamma, phi). Ignored if \code{bounds=="admissible"}.
#' @param upper Upper bounds for the parameters (alpha, beta, gamma, phi). Ignored if \code{bounds=="admissible"}.
#' @param opt.crit Optimization criterion. One of "mse" (Mean Square Error),
#' "amse" (Average MSE over first \code{nmse} forecast horizons), "sigma"
#' (Standard deviation of residuals), "mae" (Mean of absolute residuals), or
#' "lik" (Log-likelihood, the default).
#' @param nmse Number of steps for average multistep MSE (1<=\code{nmse}<=30).
#' @param bounds Type of parameter space to impose: \code{"usual" } indicates
#' all parameters must lie between specified lower and upper bounds;
#' \code{"admissible"} indicates parameters must lie in the admissible space;
#' \code{"both"} (default) takes the intersection of these regions.
#' @param ic Information criterion to be used in model selection.
#' @param restrict If \code{TRUE} (default), the models with infinite variance
#' will not be allowed.
#' @param allow.multiplicative.trend If \code{TRUE}, models with multiplicative
#' trend are allowed when searching for a model. Otherwise, the model space
#' excludes them. This argument is ignored if a multiplicative trend model is
#' explicitly requested (e.g., using \code{model="MMN"}).
#' @param use.initial.values If \code{TRUE} and \code{model} is of class
#' \code{"ets"}, then the initial values in the model are also not
#' re-estimated.
#' @param na.action A function which indicates what should happen when the data
#' contains NA values. By default, the largest contiguous portion of the
#' time-series will be used.
#' @param ... Other undocumented arguments.
#' @inheritParams forecast.ts
#'
#' @return An object of class "\code{ets}".
#'
#' The generic accessor functions \code{fitted.values} and \code{residuals}
#' extract useful features of the value returned by \code{ets} and associated
#' functions.
#' @author Rob J Hyndman
#' @seealso \code{\link[stats]{HoltWinters}}, \code{\link{rwf}},
#' \code{\link{Arima}}.
#' @references Hyndman, R.J., Koehler, A.B., Snyder, R.D., and Grose, S. (2002)
#' "A state space framework for automatic forecasting using exponential
#' smoothing methods", \emph{International J. Forecasting}, \bold{18}(3),
#' 439--454.
#'
#' Hyndman, R.J., Akram, Md., and Archibald, B. (2008) "The admissible
#' parameter space for exponential smoothing models". \emph{Annals of
#' Statistical Mathematics}, \bold{60}(2), 407--426.
#'
#' Hyndman, R.J., Koehler, A.B., Ord, J.K., and Snyder, R.D. (2008)
#' \emph{Forecasting with exponential smoothing: the state space approach},
#' Springer-Verlag. \url{http://www.exponentialsmoothing.net}.
#' @keywords ts
#' @examples
#' fit <- ets(USAccDeaths)
#' plot(forecast(fit))
#'
#' @export
ets <- function(y, model="ZZZ", damped=NULL,
                alpha=NULL, beta=NULL, gamma=NULL, phi=NULL, additive.only=FALSE, lambda=NULL, biasadj=FALSE,
                lower=c(rep(0.0001, 3), 0.8), upper=c(rep(0.9999, 3), 0.98),
                opt.crit=c("lik", "amse", "mse", "sigma", "mae"), nmse=3, bounds=c("both", "usual", "admissible"),
                ic=c("aicc", "aic", "bic"), restrict=TRUE, allow.multiplicative.trend=FALSE,
                use.initial.values=FALSE, na.action = c("na.contiguous", "na.interp", "na.fail"), ...) {
  # dataname <- substitute(y)
  opt.crit <- match.arg(opt.crit)
  bounds <- match.arg(bounds)
  ic <- match.arg(ic)
  if(!is.function(na.action)){
    na.fn_name <- match.arg(na.action)
    na.action <- get(na.fn_name)
  }

  seriesname <- deparse(substitute(y))

  if (any(class(y) %in% c("data.frame", "list", "matrix", "mts"))) {
    stop("y should be a univariate time series")
  }
  y <- as.ts(y)

  # Check if data is constant
  if (missing(model) && is.constant(y)) {
    return(ses(y, alpha = 0.99999, initial = "simple")$model)
  }

  # Remove missing values near ends
  ny <- length(y)
  y <- na.action(y)
  if (ny != length(y) && na.fn_name == "na.contiguous") {
    warning("Missing values encountered. Using longest contiguous portion of time series")
    ny <- length(y)
  }

  orig.y <- y
  if (identical(class(model), "ets") && is.null(lambda)) {
    lambda <- model$lambda
  }
  if (!is.null(lambda)) {
    y <- BoxCox(y, lambda)
    lambda <- attr(y, "lambda")
    additive.only <- TRUE
  }

  if (nmse < 1 || nmse > 30) {
    stop("nmse out of range")
  }
  m <- frequency(y)

  if (any(upper < lower)) {
    stop("Lower limits must be less than upper limits")
  }

  # If model is an ets object, re-fit model to new data
  if ("ets" %in% class(model)) {
    # Prevent alpha being zero (to avoid divide by zero in the C code)
    alpha <- max(model$par["alpha"], 1e-10)
    beta <- model$par["beta"]
    if (is.na(beta)) {
      beta <- NULL
    }
    gamma <- model$par["gamma"]
    if (is.na(gamma)) {
      gamma <- NULL
    }
    phi <- model$par["phi"]
    if (is.na(phi)) {
      phi <- NULL
    }
    modelcomponents <- paste(model$components[1], model$components[2], model$components[3], sep = "")
    damped <- (model$components[4] == "TRUE")
    if (use.initial.values) {
      errortype <- substr(modelcomponents, 1, 1)
      trendtype <- substr(modelcomponents, 2, 2)
      seasontype <- substr(modelcomponents, 3, 3)

      # Recompute errors from pegelsresid.C
      e <- pegelsresid.C(y, m, model$initstate, errortype, trendtype, seasontype, damped, alpha, beta, gamma, phi, nmse)

      # Compute error measures
      np <- length(model$par) + 1
      model$loglik <- -0.5 * e$lik
      model$aic <- e$lik + 2 * np
      model$bic <- e$lik + log(ny) * np
      model$aicc <- model$aic + 2 * np * (np + 1) / (ny - np - 1)
      model$mse <- e$amse[1]
      model$amse <- mean(e$amse)

      # Compute states, fitted values and residuals
      tsp.y <- tsp(y)
      model$states <- ts(e$states, frequency = tsp.y[3], start = tsp.y[1] - 1 / tsp.y[3])
      colnames(model$states)[1] <- "l"
      if (trendtype != "N") {
        colnames(model$states)[2] <- "b"
      }
      if (seasontype != "N") {
        colnames(model$states)[(2 + (trendtype != "N")):ncol(model$states)] <- paste("s", 1:m, sep = "")
      }
      if (errortype == "A") {
        model$fitted <- ts(y - e$e, frequency = tsp.y[3], start = tsp.y[1])
      } else {
        model$fitted <- ts(y / (1 + e$e), frequency = tsp.y[3], start = tsp.y[1])
      }
      model$residuals <- ts(e$e, frequency = tsp.y[3], start = tsp.y[1])
      model$sigma2 <- sum(model$residuals ^ 2, na.rm = TRUE) / (ny - np)
      model$x <- orig.y
      model$series <- seriesname
      if (!is.null(lambda)) {
        model$fitted <- InvBoxCox(model$fitted, lambda, biasadj, var(model$residuals))
        attr(lambda, "biasadj") <- biasadj
      }
      model$lambda <- lambda

      # Return model object
      return(model)
    }
    else {
      model <- modelcomponents
      if (missing(use.initial.values)) {
        message("Model is being refit with current smoothing parameters but initial states are being re-estimated.\nSet 'use.initial.values=TRUE' if you want to re-use existing initial values.")
      }
    }
  }

  errortype <- substr(model, 1, 1)
  trendtype <- substr(model, 2, 2)
  seasontype <- substr(model, 3, 3)

  if (!is.element(errortype, c("M", "A", "Z"))) {
    stop("Invalid error type")
  }
  if (!is.element(trendtype, c("N", "A", "M", "Z"))) {
    stop("Invalid trend type")
  }
  if (!is.element(seasontype, c("N", "A", "M", "Z"))) {
    stop("Invalid season type")
  }

  if (m < 1 || length(y) <= m) {
    # warning("I can't handle data with frequency less than 1. Seasonality will be ignored.")
    seasontype <- "N"
  }
  if (m == 1) {
    if (seasontype == "A" || seasontype == "M") {
      stop("Nonseasonal data")
    } else {
      substr(model, 3, 3) <- seasontype <- "N"
    }
  }
  if (m > 24) {
    if (is.element(seasontype, c("A", "M"))) {
      stop("Frequency too high")
    } else if (seasontype == "Z") {
      warning("I can't handle data with frequency greater than 24. Seasonality will be ignored. Try stlf() if you need seasonal forecasts.")
      substr(model, 3, 3) <- seasontype <- "N"
      # m <- 1
    }
  }

  # Check inputs
  if (restrict) {
    if ((errortype == "A" && (trendtype == "M" || seasontype == "M")) |
      (errortype == "M" && trendtype == "M" && seasontype == "A") ||
      (additive.only && (errortype == "M" || trendtype == "M" || seasontype == "M"))) {
      stop("Forbidden model combination")
    }
  }

  data.positive <- (min(y) > 0)

  if (!data.positive && errortype == "M") {
    stop("Inappropriate model for data with negative or zero values")
  }

  if (!is.null(damped)) {
    if (damped && trendtype == "N") {
      stop("Forbidden model combination")
    }
  }

  n <- length(y)
  # Check we have enough data to fit a model
  npars <- 2L # alpha + l0
  if (trendtype == "A" || trendtype == "M") {
    npars <- npars + 2L
  } # beta + b0
  if (seasontype == "A" || seasontype == "M") {
    npars <- npars + m
  } # gamma + s
  if (!is.null(damped)) {
    npars <- npars + as.numeric(damped)
  }

  # Produce something non-optimized for tiny data sets
  if (n <= npars + 4L) {
    if (!is.null(damped)) {
      if (damped) {
        warning("Not enough data to use damping")
      }
    }
    if (seasontype == "A" || seasontype == "M") {
      fit <- try(HoltWintersZZ(
        orig.y,
        alpha = alpha, beta = beta, gamma = gamma, phi = phi,
        exponential = (trendtype == "M"),
        seasonal = ifelse(seasontype != "A", "multiplicative", "additive"),
        lambda = lambda, biasadj = biasadj, warnings = FALSE
      ), silent = TRUE)
      if (!("try-error" %in% class(fit))) {
        fit$call <- match.call()
        fit$method <- as.character(fit)
        fit$series <- deparse(substitute(y))
        return(fit)
      }
      else {
        warning("Seasonal component could not be estimated")
      }
    }
    if (trendtype == "A" || trendtype == "M") {
      fit <- try(HoltWintersZZ(
        orig.y,
        alpha = alpha, beta = beta, gamma = FALSE, phi = phi,
        exponential = (trendtype == "M"),
        lambda = lambda, biasadj = biasadj, warnings = FALSE
      ), silent = TRUE)
      if (!("try-error" %in% class(fit))) {
        fit$call <- match.call()
        fit$method <- as.character(fit)
        fit$series <- deparse(substitute(y))
        return(fit)
      }
      else {
        warning("Trend component could not be estimated")
      }
    }
    if (trendtype == "N" && seasontype == "N") {
      fit <- try(HoltWintersZZ(
        orig.y,
        alpha = alpha, beta = FALSE, gamma = FALSE,
        lambda = lambda, biasadj = biasadj, warnings = FALSE
      ), silent = TRUE)
      if (!("try-error" %in% class(fit))) {
        fit$call <- match.call()
        fit$method <- as.character(fit)
        fit$series <- deparse(substitute(y))
        return(fit)
      }
    }
    # Try holt and ses and return best
    fit1 <- try(HoltWintersZZ(
      orig.y,
      alpha = alpha, beta = beta, gamma = FALSE, phi = phi,
      exponential = (trendtype == "M"),
      lambda = lambda, biasadj = biasadj, warnings = FALSE
    ), silent = TRUE)
    fit2 <- try(HoltWintersZZ(
      orig.y,
      alpha = alpha, beta = FALSE, gamma = FALSE, phi = phi,
      exponential = (trendtype == "M"),
      lambda = lambda, biasadj = biasadj, warnings = FALSE
    ), silent = TRUE)
    if ("try-error" %in% class(fit1)) {
      fit <- fit2
    } else if (fit1$sigma2 < fit2$sigma2) {
      fit <- fit1
    } else {
      fit <- fit2
    }
    if("try-error" %in% class(fit))
      stop("Unable to estimate a model.")
    fit$call <- match.call()
    fit$method <- as.character(fit)
    fit$series <- deparse(substitute(y))
    return(fit)
  }

  # Fit model (assuming only one nonseasonal model)
  if (errortype == "Z") {
    errortype <- c("A", "M")
  }
  if (trendtype == "Z") {
    if (allow.multiplicative.trend) {
      trendtype <- c("N", "A", "M")
    } else {
      trendtype <- c("N", "A")
    }
  }
  if (seasontype == "Z") {
    seasontype <- c("N", "A", "M")
  }
  if (is.null(damped)) {
    damped <- c(TRUE, FALSE)
  }
  best.ic <- Inf
  for (i in 1:length(errortype))
  {
    for (j in 1:length(trendtype))
    {
      for (k in 1:length(seasontype))
      {
        for (l in 1:length(damped))
        {
          if (trendtype[j] == "N" && damped[l]) {
            next
          }
          if (restrict) {
            if (errortype[i] == "A" && (trendtype[j] == "M" || seasontype[k] == "M")) {
              next
            }
            if (errortype[i] == "M" && trendtype[j] == "M" && seasontype[k] == "A") {
              next
            }
            if (additive.only && (errortype[i] == "M" || trendtype[j] == "M" || seasontype[k] == "M")) {
              next
            }
          }
          if (!data.positive && errortype[i] == "M") {
            next
          }
          fit <- try(etsmodel(
            y, errortype[i], trendtype[j], seasontype[k], damped[l], alpha, beta, gamma, phi,
            lower = lower, upper = upper, opt.crit = opt.crit, nmse = nmse, bounds = bounds, ...
          ), silent=TRUE)
          if(is.element("try-error", class(fit)))
            fit.ic <- Inf
          else
            fit.ic <- switch(ic, aic = fit$aic, bic = fit$bic, aicc = fit$aicc)
          if (!is.na(fit.ic)) {
            if (fit.ic < best.ic) {
              model <- fit
              best.ic <- fit.ic
              best.e <- errortype[i]
              best.t <- trendtype[j]
              best.s <- seasontype[k]
              best.d <- damped[l]
            }
          }
        }
      }
    }
  }
  if (best.ic == Inf) {
    stop("No model able to be fitted")
  }

  model$m <- m
  model$method <- paste("ETS(", best.e, ",", best.t, ifelse(best.d, "d", ""), ",", best.s, ")", sep = "")
  model$series <- seriesname
  model$components <- c(best.e, best.t, best.s, best.d)
  model$call <- match.call()
  model$initstate <- model$states[1, ]
  np <- length(model$par)
  model$sigma2 <- sum(model$residuals^2, na.rm = TRUE) / (ny - np)
  model$x <- orig.y
  if (!is.null(lambda)) {
    model$fitted <- InvBoxCox(model$fitted, lambda, biasadj, model$sigma2)
    attr(lambda, "biasadj") <- biasadj
  }

  model$lambda <- lambda
  # model$call$data <- dataname

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

#' @export
as.character.ets <- function(x, ...) {
  paste(
    "ETS(",
    x$components[1], ",",
    x$components[2],
    ifelse(x$components[4], "d", ""), ",",
    x$components[3], ")",
    sep = ""
  )
}

# myRequire <- function(libName) {

#   req.suc <- require(libName, quietly=TRUE, character.only=TRUE)
#   if(!req.suc) stop("The ",libName," package is not available.")

#   req.suc
# }

# getNewBounds <- function(par, lower, upper, nstate) {

#   myLower <- NULL
#   myUpper <- NULL

#   if("alpha" %in% names(par)) {
#     myLower <- c(myLower, lower[1])
#     myUpper <- c(myUpper, upper[1])
#   }
#   if("beta" %in% names(par)) {
#     myLower <- c(myLower, lower[2])
#     myUpper <- c(myUpper, upper[2])
#   }
#   if("gamma" %in% names(par)) {
#     myLower <- c(myLower, lower[3])
#     myUpper <- c(myUpper, upper[3])
#   }
#   if("phi" %in% names(par)) {
#     myLower <- c(myLower, lower[4])
#     myUpper <- c(myUpper, upper[4])
#   }

#   myLower <- c(myLower,rep(-1e8,nstate))
#   myUpper <- c(myUpper,rep(1e8,nstate))

#   list(lower=myLower, upper=myUpper)
# }

etsmodel <- function(y, errortype, trendtype, seasontype, damped,
                     alpha=NULL, beta=NULL, gamma=NULL, phi=NULL,
                     lower, upper, opt.crit, nmse, bounds, maxit=2000, control=NULL, seed=NULL, trace=FALSE) {
  tsp.y <- tsp(y)
  if (is.null(tsp.y)) {
    tsp.y <- c(1, length(y), 1)
  }
  if (seasontype != "N") {
    m <- tsp.y[3]
  } else {
    m <- 1
  }

  # Modify limits if alpha, beta or gamma have been specified.
  if (!is.null(alpha)) {
    upper[2] <- min(alpha, upper[2])
    upper[3] <- min(1 - alpha, upper[3])
  }
  if (!is.null(beta)) {
    lower[1] <- max(beta, lower[1])
  }
  if (!is.null(gamma)) {
    upper[1] <- min(1 - gamma, upper[1])
  }

  # Initialize smoothing parameters
  par <- initparam(alpha, beta, gamma, phi, trendtype, seasontype, damped, lower, upper, m, bounds)
  names(alpha) <- names(beta) <- names(gamma) <- names(phi) <- NULL
  par.noopt <- c(alpha = alpha, beta = beta, gamma = gamma, phi = phi)
  if (!is.null(par.noopt)) {
    par.noopt <- c(na.omit(par.noopt))
  }
  if (!is.na(par["alpha"])) {
    alpha <- par["alpha"]
  }
  if (!is.na(par["beta"])) {
    beta <- par["beta"]
  }
  if (!is.na(par["gamma"])) {
    gamma <- par["gamma"]
  }
  if (!is.na(par["phi"])) {
    phi <- par["phi"]
  }

  #    if(errortype=="M" | trendtype=="M" | seasontype=="M")
  #        bounds="usual"
  if (!check.param(alpha, beta, gamma, phi, lower, upper, bounds, m)) {
    print(paste("Model: ETS(", errortype, ",", trendtype, ifelse(damped, "d", ""), ",", seasontype, ")", sep = ""))
    stop("Parameters out of range")
  }

  # Initialize state
  init.state <- initstate(y, trendtype, seasontype)
  nstate <- length(init.state)
  par <- c(par, init.state)
  lower <- c(lower, rep(-Inf, nstate))
  upper <- c(upper, rep(Inf, nstate))

  np <- length(par)
  if (np >= length(y) - 1) { # Not enough data to continue
    return(list(aic = Inf, bic = Inf, aicc = Inf, mse = Inf, amse = Inf, fit = NULL, par = par, states = init.state))
  }

  # -------------------------------------------------

  #  if(is.null(seed)) seed <- 1000*runif(1)

  # if(solver=="malschains" || solver=="malschains_c") {

  #   malschains <- NULL
  #   if(!myRequire("Rmalschains"))
  #     stop("malschains optimizer unavailable")

  #   func <- NULL
  #   #env <- NULL

  #   if(solver=="malschains") {

  #     func <- function(myPar) {
  #       names(myPar) <- names(par)
  #       res <- lik(myPar,y=y,nstate=nstate, errortype=errortype, trendtype=trendtype,
  #           seasontype=seasontype, damped=damped, par.noopt=par.noopt, lowerb=lower, upperb=upper,
  #           opt.crit=opt.crit, nmse=nmse, bounds=bounds, m=m,pnames=names(par),pnames2=names(par.noopt))
  #       res
  #     }

  #     env <- new.env()

  #   } else {

  #     env <- etsTargetFunctionInit(par=par, y=y, nstate=nstate, errortype=errortype, trendtype=trendtype,
  #         seasontype=seasontype, damped=damped, par.noopt=par.noopt, lowerb=lower, upperb=upper,
  #         opt.crit=opt.crit, nmse=nmse, bounds=bounds, m=m,pnames=names(par),pnames2=names(par.noopt))

  #     func <- .Call("etsGetTargetFunctionRmalschainsPtr", PACKAGE="forecast")

  #   }

  #   myBounds <- getNewBounds(par, lower, upper, nstate)

  #   if(is.null(control)) {
  #     control <- Rmalschains::malschains.control(ls="simplex", lsOnly=TRUE)
  #   }

  #   control$optimum <- if(opt.crit=="lik") -1e12 else 0

  #   fredTmp <- Rmalschains::malschains(func, env=env, lower=myBounds$lower, upper=myBounds$upper,
  #       maxEvals=maxit, seed=seed, initialpop=par, control=control)

  #   fred <- NULL
  #   fred$par <- fredTmp$sol

  #   fit.par <- fred$par

  #   names(fit.par) <- names(par)

  #  } else if (solver=="Rdonlp2") {
  #
  #    donlp2 <- NULL
  #    myRequire("Rdonlp2")
  #
  #    env <- etsTargetFunctionInit(par=par, y=y, nstate=nstate, errortype=errortype, trendtype=trendtype,
  #        seasontype=seasontype, damped=damped, par.noopt=par.noopt, lowerb=lower, upperb=upper,
  #        opt.crit=opt.crit, nmse=nmse, bounds=bounds, m=m,pnames=names(par),pnames2=names(par.noopt))
  #
  #    func <- .Call("etsGetTargetFunctionRdonlp2Ptr", PACKAGE="forecast")
  #
  #    myBounds <- getNewBounds(par, lower, upper, nstate)
  #
  #    fred <- donlp2(par, func, env=env, par.lower=myBounds$lower, par.upper=myBounds$upper)#, nlin.lower=c(-1), nlin.upper=c(1)) #nlin.lower=c(0,-Inf, -Inf, -Inf), nlin.upper=c(0,0,0,0))
  #
  #    fit.par <- fred$par
  #
  #    names(fit.par) <- names(par)

  #  } else if(solver=="optim_c"){

  env <- etsTargetFunctionInit(
    par = par, y = y, nstate = nstate, errortype = errortype, trendtype = trendtype,
    seasontype = seasontype, damped = damped, par.noopt = par.noopt, lowerb = lower, upperb = upper,
    opt.crit = opt.crit, nmse = as.integer(nmse), bounds = bounds, m = m, pnames = names(par), pnames2 = names(par.noopt)
  )

  fred <- .Call(
    "etsNelderMead", par, env, -Inf,
    sqrt(.Machine$double.eps), 1.0, 0.5, 2.0, trace, maxit, PACKAGE = "forecast"
  )

  fit.par <- fred$par

  names(fit.par) <- names(par)

  # } else { #if(solver=="optim")

  #   # Optimize parameters and state
  #   if(length(par)==1)
  #     method <- "Brent"
  #   else
  #   	method <- "Nelder-Mead"

  #   fred <- optim(par,lik,method=method,y=y,nstate=nstate, errortype=errortype, trendtype=trendtype,
  #       seasontype=seasontype, damped=damped, par.noopt=par.noopt, lowerb=lower, upperb=upper,
  #       opt.crit=opt.crit, nmse=nmse, bounds=bounds, m=m,pnames=names(par),pnames2=names(par.noopt),
  #       control=list(maxit=maxit))

  #   fit.par <- fred$par
  #   names(fit.par) <- names(par)
  # }

  # -------------------------------------------------

  init.state <- fit.par[(np - nstate + 1):np]
  # Add extra state
  if (seasontype != "N") {
    init.state <- c(init.state, m * (seasontype == "M") - sum(init.state[(2 + (trendtype != "N")):nstate]))
  }

  if (!is.na(fit.par["alpha"])) {
    alpha <- fit.par["alpha"]
  }
  if (!is.na(fit.par["beta"])) {
    beta <- fit.par["beta"]
  }
  if (!is.na(fit.par["gamma"])) {
    gamma <- fit.par["gamma"]
  }
  if (!is.na(fit.par["phi"])) {
    phi <- fit.par["phi"]
  }

  e <- pegelsresid.C(y, m, init.state, errortype, trendtype, seasontype, damped, alpha, beta, gamma, phi, nmse)

  np <- np + 1
  ny <- length(y)
  aic <- e$lik + 2 * np
  bic <- e$lik + log(ny) * np
  aicc <- aic + 2 * np * (np + 1) / (ny - np - 1)

  mse <- e$amse[1]
  amse <- mean(e$amse)

  states <- ts(e$states, frequency = tsp.y[3], start = tsp.y[1] - 1 / tsp.y[3])
  colnames(states)[1] <- "l"
  if (trendtype != "N") {
    colnames(states)[2] <- "b"
  }
  if (seasontype != "N") {
    colnames(states)[(2 + (trendtype != "N")):ncol(states)] <- paste("s", 1:m, sep = "")
  }

  tmp <- c("alpha", rep("beta", trendtype != "N"), rep("gamma", seasontype != "N"), rep("phi", damped))
  fit.par <- c(fit.par, par.noopt)
  #    fit.par <- fit.par[order(names(fit.par))]
  if (errortype == "A") {
    fits <- y - e$e
  } else {
    fits <- y / (1 + e$e)
  }

  return(list(
    loglik = -0.5 * e$lik, aic = aic, bic = bic, aicc = aicc, mse = mse, amse = amse, fit = fred, residuals = ts(e$e, frequency = tsp.y[3], start = tsp.y[1]), fitted = ts(fits, frequency = tsp.y[3], start = tsp.y[1]),
    states = states, par = fit.par
  ))
}

etsTargetFunctionInit <- function(par, y, nstate, errortype, trendtype, seasontype, damped, par.noopt, lowerb, upperb,
                                  opt.crit, nmse, bounds, m, pnames, pnames2) {
  names(par) <- pnames
  names(par.noopt) <- pnames2
  alpha <- c(par["alpha"], par.noopt["alpha"])["alpha"]
  if (is.na(alpha)) {
    stop("alpha problem!")
  }
  if (trendtype != "N") {
    beta <- c(par["beta"], par.noopt["beta"])["beta"]
    if (is.na(beta)) {
      stop("beta Problem!")
    }
  }
  else {
    beta <- NULL
  }
  if (seasontype != "N") {
    gamma <- c(par["gamma"], par.noopt["gamma"])["gamma"]
    if (is.na(gamma)) {
      stop("gamma Problem!")
    }
  }
  else {
    m <- 1
    gamma <- NULL
  }
  if (damped) {
    phi <- c(par["phi"], par.noopt["phi"])["phi"]
    if (is.na(phi)) {
      stop("phi Problem!")
    }
  }
  else {
    phi <- NULL
  }

  # determine which values to optimize and which ones are given by the user/not needed
  optAlpha <- !is.null(alpha)
  optBeta <- !is.null(beta)
  optGamma <- !is.null(gamma)
  optPhi <- !is.null(phi)

  givenAlpha <- FALSE
  givenBeta <- FALSE
  givenGamma <- FALSE
  givenPhi <- FALSE

  if (!is.null(par.noopt["alpha"])) {
    if (!is.na(par.noopt["alpha"])) {
      optAlpha <- FALSE
      givenAlpha <- TRUE
    }
  }
  if (!is.null(par.noopt["beta"])) {
    if (!is.na(par.noopt["beta"])) {
      optBeta <- FALSE
      givenBeta <- TRUE
    }
  }
  if (!is.null(par.noopt["gamma"])) {
    if (!is.na(par.noopt["gamma"])) {
      optGamma <- FALSE
      givenGamma <- TRUE
    }
  }
  if (!is.null(par.noopt["phi"])) {
    if (!is.na(par.noopt["phi"])) {
      optPhi <- FALSE
      givenPhi <- TRUE
    }
  }

  if (!damped) {
    phi <- 1
  }
  if (trendtype == "N") {
    beta <- 0
  }
  if (seasontype == "N") {
    gamma <- 0
  }

  #  cat("alpha: ", alpha)
  #  cat(" beta: ", beta)
  #  cat(" gamma: ", gamma)
  #  cat(" phi: ", phi, "\n")
  #
  #  cat("useAlpha: ", useAlpha)
  #  cat(" useBeta: ", useBeta)
  #  cat(" useGamma: ", useGamma)
  #  cat(" usePhi: ", usePhi, "\n")

  env <- new.env()

  res <- .Call(
    "etsTargetFunctionInit", y = y, nstate = nstate, errortype = switch(errortype, "A" = 1, "M" = 2),
    trendtype = switch(trendtype, "N" = 0, "A" = 1, "M" = 2), seasontype = switch(seasontype, "N" = 0, "A" = 1, "M" = 2),
    damped = damped, lowerb = lowerb, upperb = upperb,
    opt.crit = opt.crit, nmse = as.integer(nmse), bounds = bounds, m = m,
    optAlpha, optBeta, optGamma, optPhi,
    givenAlpha, givenBeta, givenGamma, givenPhi,
    alpha, beta, gamma, phi, env, PACKAGE = "forecast"
  )
  res
}

initparam <- function(alpha, beta, gamma, phi, trendtype, seasontype, damped, lower, upper, m, bounds) {
  if(bounds == "admissible") {
    lower[1L:3L] <- lower[1L:3L]*0
    upper[1L:3L] <- upper[1L:3L]*0 + 1e-3
  } else if (any(lower > upper)) {
    stop("Inconsistent parameter boundaries")
  }

  # Select alpha
  if (is.null(alpha)) {
    alpha <- lower[1] + 0.2 * (upper[1] - lower[1]) / m
    if (alpha > 1 || alpha < 0) {
      alpha <- lower[1] + 2e-3
    }
    par <- c(alpha = alpha)
  }
  else {
    par <- numeric(0)
  }

  # Select beta
  if (trendtype != "N" && is.null(beta)) {
    # Ensure beta < alpha
    upper[2] <- min(upper[2], alpha)
    beta <- lower[2] + 0.1 * (upper[2] - lower[2])
    if (beta < 0 || beta > alpha) {
      beta <- alpha - 1e-3
    }
    par <- c(par, beta = beta)
  }

  # Select gamma
  if (seasontype != "N" && is.null(gamma)) {
    # Ensure gamma < 1-alpha
    upper[3] <- min(upper[3], 1 - alpha)
    gamma <- lower[3] + 0.05 * (upper[3] - lower[3])
    if (gamma < 0 || gamma > 1 - alpha) {
      gamma <- 1 - alpha - 1e-3
    }
    par <- c(par, gamma = gamma)
  }

  # Select phi
  if (damped && is.null(phi)) {
    phi <- lower[4] + .99 * (upper[4] - lower[4])
    if (phi < 0 || phi > 1) {
      phi <- upper[4] - 1e-3
    }
    par <- c(par, phi = phi)
  }

  return(par)
}

check.param <- function(alpha, beta, gamma, phi, lower, upper, bounds, m) {
  if (bounds != "admissible") {
    if (!is.null(alpha)) {
      if (alpha < lower[1] || alpha > upper[1]) {
        return(0)
      }
    }
    if (!is.null(beta)) {
      if (beta < lower[2] || beta > alpha || beta > upper[2]) {
        return(0)
      }
    }
    if (!is.null(phi)) {
      if (phi < lower[4] || phi > upper[4]) {
        return(0)
      }
    }
    if (!is.null(gamma)) {
      if (gamma < lower[3] || gamma > 1 - alpha || gamma > upper[3]) {
        return(0)
      }
    }
  }
  if (bounds != "usual") {
    if (!admissible(alpha, beta, gamma, phi, m)) {
      return(0)
    }
  }
  return(1)
}

initstate <- function(y, trendtype, seasontype) {
  if (seasontype != "N") {
    # Do decomposition
    m <- frequency(y)
    n <- length(y)
    if (n < 4) {
      stop("You've got to be joking (not enough data).")
    } else if (n < 3 * m) # Fit simple Fourier model.
    {
      fouriery <- fourier(y, 1)
      fit <- tslm(y ~ trend + fouriery)
      if (seasontype == "A") {
        y.d <- list(seasonal = y - fit$coefficients[1] - fit$coefficients[2] * (1:n))
      } else { # seasontype=="M". Biased method, but we only need a starting point
        y.d <- list(seasonal = y / (fit$coefficients[1] + fit$coefficients[2] * (1:n)))
      }
    }
    else { # n is large enough to do a decomposition
      y.d <- decompose(y, type = switch(seasontype, A = "additive", M = "multiplicative"))
    }

    init.seas <- rev(y.d$seasonal[2:m]) # initial seasonal component
    names(init.seas) <- paste("s", 0:(m - 2), sep = "")
    # Seasonally adjusted data
    if (seasontype == "A") {
      y.sa <- y - y.d$seasonal
    } else {
      init.seas <- pmax(init.seas, 1e-2) # We do not want negative seasonal indexes
      if (sum(init.seas) > m) {
        init.seas <- init.seas / sum(init.seas + 1e-2)
      }
      y.sa <- y / pmax(y.d$seasonal, 1e-2)
    }
  }
  else # non-seasonal model
  {
    m <- 1
    init.seas <- NULL
    y.sa <- y
  }

  maxn <- min(max(10, 2 * m), length(y.sa))

  if (trendtype == "N") {
    l0 <- mean(y.sa[1:maxn])
    b0 <- NULL
  }
  else # Simple linear regression on seasonally adjusted data
  {
    fit <- lsfit(1:maxn, y.sa[1:maxn])
    if (trendtype == "A") {
      l0 <- fit$coefficients[1]
      b0 <- fit$coefficients[2]
      # If error type is "M", then we don't want l0+b0=0.
      # So perturb just in case.
      if (abs(l0 + b0) < 1e-8) {
        l0 <- l0 * (1 + 1e-3)
        b0 <- b0 * (1 - 1e-3)
      }
    }
    else # if(trendtype=="M")
    {
      l0 <- fit$coefficients[1] + fit$coefficients[2] # First fitted value
      if (abs(l0) < 1e-8) {
        l0 <- 1e-7
      }
      b0 <- (fit$coefficients[1] + 2 * fit$coefficients[2]) / l0 # Ratio of first two fitted values
      l0 <- l0 / b0 # First fitted value divided by b0
      if (abs(b0) > 1e10) { # Avoid infinite slopes
        b0 <- sign(b0) * 1e10
      }
      if (l0 < 1e-8 || b0 < 1e-8) # Simple linear approximation didn't work.
      {
        l0 <- max(y.sa[1], 1e-3)
        b0 <- max(y.sa[2] / y.sa[1], 1e-3)
      }
    }
  }

  names(l0) <- "l"
  if (!is.null(b0)) {
    names(b0) <- "b"
  }
  return(c(l0, b0, init.seas))
}

lik <- function(par, y, nstate, errortype, trendtype, seasontype, damped, par.noopt, lowerb, upperb,
                opt.crit, nmse, bounds, m, pnames, pnames2) {

  # browser()

  # cat("par: ", par, "\n")

  names(par) <- pnames
  names(par.noopt) <- pnames2
  alpha <- c(par["alpha"], par.noopt["alpha"])["alpha"]
  if (is.na(alpha)) {
    stop("alpha problem!")
  }
  if (trendtype != "N") {
    beta <- c(par["beta"], par.noopt["beta"])["beta"]
    if (is.na(beta)) {
      stop("beta Problem!")
    }
  }
  else {
    beta <- NULL
  }
  if (seasontype != "N") {
    gamma <- c(par["gamma"], par.noopt["gamma"])["gamma"]
    if (is.na(gamma)) {
      stop("gamma Problem!")
    }
  }
  else {
    m <- 1
    gamma <- NULL
  }
  if (damped) {
    phi <- c(par["phi"], par.noopt["phi"])["phi"]
    if (is.na(phi)) {
      stop("phi Problem!")
    }
  }
  else {
    phi <- NULL
  }

  if (!check.param(alpha, beta, gamma, phi, lowerb, upperb, bounds, m)) {
    return(Inf)
  }

  np <- length(par)

  init.state <- par[(np - nstate + 1):np]
  # Add extra state
  if (seasontype != "N") {
    init.state <- c(init.state, m * (seasontype == "M") - sum(init.state[(2 + (trendtype != "N")):nstate]))
  }
  # Check states
  if (seasontype == "M") {
    seas.states <- init.state[-(1:(1 + (trendtype != "N")))]
    if (min(seas.states) < 0) {
      return(Inf)
    }
  }

  e <- pegelsresid.C(y, m, init.state, errortype, trendtype, seasontype, damped, alpha, beta, gamma, phi, nmse)

  if (is.na(e$lik)) {
    return(Inf)
  }
  if (e$lik < -1e10) { # Avoid perfect fits
    return(-1e10)
  }

  #      cat("lik: ", e$lik, "\n")
  #    points(alpha,e$lik,col=2)

  if (opt.crit == "lik") {
    return(e$lik)
  } else if (opt.crit == "mse") {
    return(e$amse[1])
  } else if (opt.crit == "amse") {
    return(mean(e$amse))
  } else if (opt.crit == "sigma") {
    return(mean(e$e ^ 2))
  } else if (opt.crit == "mae") {
    return(mean(abs(e$e)))
  }
}

#' @export
print.ets <- function(x, ...) {
  cat(paste(x$method, "\n\n"))
  if(!is.null(x$call)) {
    cat("Call:", deparse(x$call), "", sep = "\n")
  }
  ncoef <- length(x$initstate)
  if (!is.null(x$lambda)) {
    cat("  Box-Cox transformation: lambda=", round(x$lambda, 4), "\n\n")
  }

  cat("  Smoothing parameters:\n")
  cat(paste("    alpha =", round(x$par["alpha"], 4), "\n"))
  if (x$components[2] != "N") {
    cat(paste("    beta  =", round(x$par["beta"], 4), "\n"))
  }
  if (x$components[3] != "N") {
    cat(paste("    gamma =", round(x$par["gamma"], 4), "\n"))
  }
  if (x$components[4] != "FALSE") {
    cat(paste("    phi   =", round(x$par["phi"], 4), "\n"))
  }

  cat("\n  Initial states:\n")
  cat(paste("    l =", round(x$initstate[1], 4), "\n"))
  if (x$components[2] != "N") {
    cat(paste("    b =", round(x$initstate[2], 4), "\n"))
  } else {
    x$initstate <- c(x$initstate[1], NA, x$initstate[2:ncoef])
    ncoef <- ncoef + 1
  }
  if (x$components[3] != "N") {
    cat("    s = ")
    if (ncoef <= 8) {
      cat(round(x$initstate[3:ncoef], 4))
    } else {
      cat(round(x$initstate[3:8], 4))
      cat("\n           ")
      cat(round(x$initstate[9:ncoef], 4))
    }
    cat("\n")
  }

  cat("\n  sigma:  ")
  cat(round(sqrt(x$sigma2), 4))
  if (!is.null(x$aic)) {
    stats <- c(x$aic, x$aicc, x$bic)
    names(stats) <- c("AIC", "AICc", "BIC")
    cat("\n\n")
    print(stats)
  }
  #    cat("\n  AIC:    ")
  #    cat(round(x$aic,4))
  #    cat("\n  AICc:   ")
  #    cat(round(x$aicc,4))
  #    cat("\n  BIC:    ")
  #    cat(round(x$bic,4))
}

pegelsresid.C <- function(y, m, init.state, errortype, trendtype, seasontype, damped, alpha, beta, gamma, phi, nmse) {
  n <- length(y)
  p <- length(init.state)
  x <- numeric(p * (n + 1))
  x[1:p] <- init.state
  e <- numeric(n)
  lik <- 0
  if (!damped) {
    phi <- 1
  }
  if (trendtype == "N") {
    beta <- 0
  }
  if (seasontype == "N") {
    gamma <- 0
  }

  amse <- numeric(nmse)

  Cout <- .C(
    "etscalc",
    as.double(y),
    as.integer(n),
    as.double(x),
    as.integer(m),
    as.integer(switch(errortype, "A" = 1, "M" = 2)),
    as.integer(switch(trendtype, "N" = 0, "A" = 1, "M" = 2)),
    as.integer(switch(seasontype, "N" = 0, "A" = 1, "M" = 2)),
    as.double(alpha),
    as.double(beta),
    as.double(gamma),
    as.double(phi),
    as.double(e),
    as.double(lik),
    as.double(amse),
    as.integer(nmse),
    PACKAGE = "forecast"
  )
  if (!is.na(Cout[[13]])) {
    if (abs(Cout[[13]] + 99999) < 1e-7) {
      Cout[[13]] <- NA
    }
  }
  tsp.y <- tsp(y)
  e <- ts(Cout[[12]])
  tsp(e) <- tsp.y

  return(list(lik = Cout[[13]], amse = Cout[[14]], e = e, states = matrix(Cout[[3]], nrow = n + 1, ncol = p, byrow = TRUE)))
}

admissible <- function(alpha, beta, gamma, phi, m) {
  if (is.null(phi)) {
    phi <- 1
  }
  if (phi < 0 || phi > 1 + 1e-8) {
    return(0)
  }
  if (is.null(gamma)) {
    if (alpha < 1 - 1 / phi || alpha > 1 + 1 / phi) {
      return(0)
    }
    if (!is.null(beta)) {
      if (beta < alpha * (phi - 1) || beta > (1 + phi) * (2 - alpha)) {
        return(0)
      }
    }
  }
  else if (m > 1) # Seasonal model
  {
    if (is.null(beta)) {
      beta <- 0
    }
    if (gamma < max(1 - 1 / phi - alpha, 0) || gamma > 1 + 1 / phi - alpha) {
      return(0)
    }
    if (alpha < 1 - 1 / phi - gamma * (1 - m + phi + phi * m) / (2 * phi * m)) {
      return(0)
    }
    if (beta < -(1 - phi) * (gamma / m + alpha)) {
      return(0)
    }

    # End of easy tests. Now use characteristic equation
    P <- c(phi * (1 - alpha - gamma), alpha + beta - alpha * phi + gamma - 1, rep(alpha + beta - alpha * phi, m - 2), (alpha + beta - phi), 1)
    roots <- polyroot(P)

    # cat("maxpolyroots: ", max(abs(roots)), "\n")

    if (max(abs(roots)) > 1 + 1e-10) {
      return(0)
    }
  }
  # Passed all tests
  return(1)
}

### PLOT COMPONENTS

#' Plot components from ETS model
#'
#' Produces a plot of the level, slope and seasonal components from an ETS
#' model.
#'
#' \code{autoplot} will produce an equivalent plot as a ggplot object.
#'
#' @param x Object of class \dQuote{ets}.
#' @param object Object of class \dQuote{ets}. Used for ggplot graphics (S3
#' method consistency).
#' @param range.bars Logical indicating if each plot should have a bar at its
#' right side representing relative size. If NULL, automatic selection takes
#' place.
#' @param ... Other plotting parameters to affect the plot.
#' @return None. Function produces a plot
#' @author Rob J Hyndman & Mitchell O'Hara-Wild
#' @seealso \code{\link{ets}}
#' @keywords hplot
#' @examples
#'
#' fit <- ets(USAccDeaths)
#' plot(fit)
#' plot(fit,plot.type="single",ylab="",col=1:3)
#'
#' library(ggplot2)
#' autoplot(fit)
#'
#' @export
plot.ets <- function(x, ...) {
  if (!is.null(x$lambda)) {
    y <- BoxCox(x$x, x$lambda)
  } else {
    y <- x$x
  }
  if (x$components[3] == "N" && x$components[2] == "N") {
    plot(
      cbind(observed = y, level = x$states[, 1]),
      main = paste("Decomposition by", x$method, "method"), ...
    )
  }
  else if (x$components[3] == "N") {
    plot(
      cbind(observed = y, level = x$states[, 1], slope = x$states[, "b"]),
      main = paste("Decomposition by", x$method, "method"), ...
    )
  }
  else if (x$components[2] == "N") {
    plot(
      cbind(observed = y, level = x$states[, 1], season = x$states[, "s1"]),
      main = paste("Decomposition by", x$method, "method"), ...
    )
  }
  else {
    plot(
      cbind(
        observed = y, level = x$states[, 1], slope = x$states[, "b"],
        season = x$states[, "s1"]
      ),
      main = paste("Decomposition by", x$method, "method"), ...
    )
  }
}

#' @export
summary.ets <- function(object, ...) {
  class(object) <- c("summary.ets", class(object))
  object
}

#' @export
print.summary.ets <- function(x, ...) {
  NextMethod()
  cat("\nTraining set error measures:\n")
  print(accuracy(x))
}

#' @export
coef.ets <- function(object, ...) {
  object$par
}

#' @rdname fitted.Arima
#' @export
fitted.ets <- function(object, h=1, ...) {
  if (h == 1) {
    return(object$fitted)
  }
  else {
    return(hfitted(object = object, h = h, FUN = "ets", ...))
  }
}

#' @export
hfitted.ets <- function(object, h=1, ...) {
  n <- length(object$x)
  out <- rep(NA_real_, n)
  for(i in seq_len(n-h+1)) {
    out[i+h-1] <- .C(
      "etsforecast",
      as.double(object$states[i, ]),
      as.integer(object$m),
      as.integer(switch(object$components[2], "N" = 0, "A" = 1, "M" = 2)),
      as.integer(switch(object$components[3], "N" = 0, "A" = 1, "M" = 2)),
      as.double(ifelse(object$components[4] == "FALSE", 1, object$par["phi"])),
      as.integer(h),
      as.double(numeric(h)),
      PACKAGE = "forecast"
    )[[7]][h]
  }
  out
}

#' @export
logLik.ets <- function(object, ...) {
  structure(object$loglik, df = length(object$par) + 1, class = "logLik")
}

#' @export
nobs.ets <- function(object, ...) {
  length(object$x)
}


#' Is an object a particular model type?
#'
#' Returns true if the model object is of a particular type
#'
#' @param x object to be tested
#' @export
is.ets <- function(x) {
  inherits(x, "ets")
}

Try the forecast package in your browser

Any scripts or data that you put into this service are public.

forecast documentation built on June 22, 2024, 9:20 a.m.