R/ets.R
In toppyy/forecastlight: Forecasting Functions for Time Series and Linear Models
Defines functions
is.ets
nobs.ets
logLik.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
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)
#' @param upper Upper bounds for the parameters (alpha, beta, gamma, phi)
#' @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
#'
#' @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 (class(model) == "ets") {
# 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)
fit <- 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, ...
)
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)
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 = "forecastlight"
)
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 = "forecastlight"
)
res
}
initparam <- function(alpha, beta, gamma, phi, trendtype, seasontype, damped, lower, upper, m) {
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(paste("Call:\n", deparse(x$call), "\n\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 = "forecastlight"
)
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
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")
}
toppyy/forecastlight documentation built on Feb. 4, 2022, 12:28 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions is.ets nobs.ets logLik.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 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)
#' @param upper Upper bounds for the parameters (alpha, beta, gamma, phi)
#' @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
#'
#' @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 (class(model) == "ets") {
# 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)
fit <- 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, ...
)
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)
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 = "forecastlight"
)
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 = "forecastlight"
)
res
}
initparam <- function(alpha, beta, gamma, phi, trendtype, seasontype, damped, lower, upper, m) {
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(paste("Call:\n", deparse(x$call), "\n\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 = "forecastlight"
)
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
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")
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.