R/forecast2.R
In robjhyndman/forecast: Forecasting Functions for Time Series and Linear Models
Defines functions
croston2
croston
forecast.HoltWinters
forecast.StructTS
InvBoxCoxf
InvBoxCox
BoxCox
meanf
Documented in
BoxCox
croston
forecast.HoltWinters
forecast.StructTS
InvBoxCox
meanf
# Mean forecast
#' Mean Forecast
#'
#' Returns forecasts and prediction intervals for an iid model applied to y.
#'
#' The iid model is \deqn{Y_t=\mu + Z_t}{Y[t]=mu + Z[t]} where \eqn{Z_t}{Z[t]}
#' is a normal iid error. Forecasts are given by \deqn{Y_n(h)=\mu}{Y[n+h]=mu}
#' where \eqn{\mu}{mu} is estimated by the sample mean.
#'
#' @param y a numeric vector or time series of class \code{ts}
#' @param h Number of periods for forecasting
#' @param level Confidence levels for prediction intervals.
#' @param fan If TRUE, level is set to seq(51,99,by=3). This is suitable for
#' fan plots.
#' @param bootstrap If TRUE, use a bootstrap method to compute prediction intervals.
#' Otherwise, assume a normal distribution.
#' @param npaths Number of bootstrapped sample paths to use if \code{bootstrap==TRUE}.
#' @param x Deprecated. Included for backwards compatibility.
#' @inheritParams forecast.ts
#'
#' @return An object of class "\code{forecast}".
#'
#' The function \code{summary} is used to obtain and print a summary of the
#' results, while the function \code{plot} produces a plot of the forecasts and
#' prediction intervals.
#'
#' The generic accessor functions \code{fitted.values} and \code{residuals}
#' extract useful features of the value returned by \code{meanf}.
#'
#' An object of class \code{"forecast"} is a list containing at least the
#' following elements: \item{model}{A list containing information about the
#' fitted model} \item{method}{The name of the forecasting method as a
#' character string} \item{mean}{Point forecasts as a time series}
#' \item{lower}{Lower limits for prediction intervals} \item{upper}{Upper
#' limits for prediction intervals} \item{level}{The confidence values
#' associated with the prediction intervals} \item{x}{The original time series
#' (either \code{object} itself or the time series used to create the model
#' stored as \code{object}).} \item{residuals}{Residuals from the fitted model.
#' That is x minus fitted values.} \item{fitted}{Fitted values (one-step
#' forecasts)}
#' @author Rob J Hyndman
#' @seealso \code{\link{rwf}}
#' @keywords ts
#' @examples
#' nile.fcast <- meanf(Nile, h=10)
#' plot(nile.fcast)
#'
#' @export
meanf <- function(y, h=10, level=c(80, 95), fan=FALSE, lambda=NULL, biasadj=FALSE,
bootstrap=FALSE, npaths=5000, x=y) {
n <- length(x)
if (!is.null(lambda)) {
origx <- x
x <- BoxCox(x, lambda)
lambda <- attr(x, "lambda")
}
meanx <- mean(x, na.rm = TRUE)
fits <- rep(meanx, length(x))
res <- x - fits
f <- rep(meanx, h)
if (fan) {
level <- seq(51, 99, by = 3)
} else {
if (min(level) > 0 && max(level) < 1) {
level <- 100 * level
} else if (min(level) < 0 || max(level) > 99.99) {
stop("Confidence limit out of range")
}
}
nconf <- length(level)
s <- sd(x, na.rm = TRUE)
if (bootstrap) {
e <- na.omit(res) - mean(res, na.rm = TRUE)
sim <- matrix(sample(e, size = npaths * h, replace = TRUE), ncol = npaths, nrow = h)
sim <- sweep(sim, 1, f, "+")
lower <- t(apply(sim, 1, quantile, prob = .5 - level / 200))
upper <- t(apply(sim, 1, quantile, prob = .5 + level / 200))
}
else {
lower <- upper <- matrix(NA, nrow = h, ncol = nconf)
for (i in 1:nconf)
{
if (n > 1) {
tfrac <- qt(0.5 - level[i] / 200, n - 1)
} else {
tfrac <- -Inf
}
w <- -tfrac * s * sqrt(1 + 1 / n)
lower[, i] <- f - w
upper[, i] <- f + w
}
}
colnames(lower) <- colnames(upper) <- paste(level, "%", sep = "")
if (is.ts(x)) {
fits <- copy_msts(x, fits)
res <- copy_msts(x, res)
f <- future_msts(x, f)
lower <- future_msts(x, lower)
upper <- future_msts(x, upper)
}
if (!is.null(lambda)) {
fits <- InvBoxCox(fits, lambda)
x <- origx
f <- InvBoxCox(f, lambda, biasadj, list(level = level, upper = upper, lower = lower))
lower <- InvBoxCox(lower, lambda)
upper <- InvBoxCox(upper, lambda)
}
out <- list(
method = "Mean", level = level, x = x, series = deparse(substitute(y)), mean = f, lower = lower, upper = upper,
model = structure(list(mu = f[1], mu.se = s / sqrt(length(x)), sd = s, bootstrap = bootstrap), class = "meanf"), lambda = lambda, fitted = fits, residuals = res
)
out$model$call <- match.call()
return(structure(out, class = "forecast"))
}
#' Box Cox Transformation
#'
#' BoxCox() returns a transformation of the input variable using a Box-Cox
#' transformation. InvBoxCox() reverses the transformation.
#'
#' The Box-Cox transformation (as given by Bickel & Doksum 1981) is given by \deqn{f_\lambda(x) =(sign(x)|x|^\lambda -
#' 1)/\lambda}{f(x;lambda)=(sign(x)|x|^lambda - 1)/lambda} if \eqn{\lambda\ne0}{lambda
#' is not equal to 0}. For \eqn{\lambda=0}{lambda=0},
#' \deqn{f_0(x)=\log(x)}{f(x;0)=log(x)}.
#'
#' @param x a numeric vector or time series of class \code{ts}.
#' @param lambda transformation parameter. If \code{lambda = "auto"}, then
#' the transformation parameter lambda is chosen using BoxCox.lambda (with a lower bound of -0.9)
#' @param biasadj Use adjusted back-transformed mean for Box-Cox
#' transformations. If transformed data is used to produce forecasts and fitted values,
#' a regular back transformation will result in median forecasts. If biasadj is TRUE,
#' an adjustment will be made to produce mean forecasts and fitted values.
#' @param fvar Optional parameter required if biasadj=TRUE. Can either be the
#' forecast variance, or a list containing the interval \code{level}, and the
#' corresponding \code{upper} and \code{lower} intervals.
#' @return a numeric vector of the same length as x.
#' @author Rob J Hyndman & Mitchell O'Hara-Wild
#' @seealso \code{\link{BoxCox.lambda}}
#' @references Box, G. E. P. and Cox, D. R. (1964) An analysis of
#' transformations. \emph{JRSS B} \bold{26} 211--246.
#' Bickel, P. J. and Doksum K. A. (1981) An Analysis of Transformations Revisited. \emph{JASA} \bold{76} 296-311.
#' @keywords ts
#' @examples
#'
#' lambda <- BoxCox.lambda(lynx)
#' lynx.fit <- ar(BoxCox(lynx,lambda))
#' plot(forecast(lynx.fit,h=20,lambda=lambda))
#'
#' @export
BoxCox <- function(x, lambda) {
if (lambda == "auto") {
lambda <- BoxCox.lambda(x, lower = -0.9)
}
if (lambda < 0) {
x[x < 0] <- NA
}
if (lambda == 0) {
out <- log(x)
} else {
out <- (sign(x) * abs(x) ^ lambda - 1) / lambda
}
if (!is.null(colnames(x))) {
colnames(out) <- colnames(x)
}
attr(out, "lambda") <- lambda
return(out)
}
#' @rdname BoxCox
#' @export
InvBoxCox <- function(x, lambda, biasadj=FALSE, fvar=NULL) {
if (lambda < 0) {
x[x > -1 / lambda] <- NA
}
if (lambda == 0) {
out <- exp(x)
} else {
xx <- x * lambda + 1
out <- sign(xx) * abs(xx) ^ (1 / lambda)
}
if (!is.null(colnames(x))) {
colnames(out) <- colnames(x)
}
if (is.null(biasadj)) {
biasadj <- attr(lambda, "biasadj")
}
if (!is.logical(biasadj)) {
warning("biasadj information not found, defaulting to FALSE.")
biasadj <- FALSE
}
if (biasadj) {
if (is.null(fvar)) {
stop("fvar must be provided when biasadj=TRUE")
}
if (is.list(fvar)) { # Create fvar from forecast interval
level <- max(fvar$level)
if (NCOL(fvar$upper) > 1 && NCOL(fvar$lower)) {
i <- match(level, fvar$level)
fvar$upper <- fvar$upper[, i]
fvar$lower <- fvar$lower[, i]
}
if (level > 1) {
level <- level / 100
}
level <- mean(c(level, 1))
# Note: Use BoxCox transformed upper and lower values
fvar <- as.numeric((fvar$upper - fvar$lower) / stats::qnorm(level) / 2) ^ 2
}
if (NCOL(fvar) > 1) {
fvar <- diag(fvar)
}
out <- out * (1 + 0.5 * as.numeric(fvar) * (1 - lambda) / (out) ^ (2 * lambda))
}
return(out)
}
# Deprecated
InvBoxCoxf <- function(x=NULL, fvar=NULL, lambda=NULL) {
message("Deprecated, use InvBoxCox instead")
if (is.null(lambda)) {
stop("Must specify lambda using lambda=numeric(1)")
}
if (is.null(fvar)) {
level <- max(x$level)
if (NCOL(x$upper) > 1 && NCOL(x$lower)) {
i <- match(level, x$level)
x$upper <- x$upper[, i]
x$lower <- x$lower[, i]
}
if (level > 1) {
level <- level / 100
}
level <- mean(c(level, 1))
# Note: Use BoxCox transformed upper and lower values
fvar <- ((x$upper - x$lower) / stats::qnorm(level) / 2) ^ 2
}
else {
x <- list(mean = x)
}
if ("matrix" %in% class(fvar)) {
fvar <- diag(fvar)
}
return(x$mean * (1 + 0.5 * fvar * (1 - lambda) / (x$mean) ^ (2 * lambda)))
}
#' Forecasting using Structural Time Series models
#'
#' Returns forecasts and other information for univariate structural time
#' series models.
#'
#' This function calls \code{predict.StructTS} and constructs an object of
#' class "\code{forecast}" from the results.
#'
#' @param object An object of class "\code{StructTS}". Usually the result of a
#' call to \code{\link[stats]{StructTS}}.
#' @param h Number of periods for forecasting
#' @param level Confidence level for prediction intervals.
#' @param fan If TRUE, level is set to seq(51,99,by=3). This is suitable for
#' fan plots.
#' @param ... Other arguments.
#' @inheritParams forecast.ts
#'
#' @return An object of class "\code{forecast}".
#'
#' The function \code{summary} is used to obtain and print a summary of the
#' results, while the function \code{plot} produces a plot of the forecasts and
#' prediction intervals.
#'
#' The generic accessor functions \code{fitted.values} and \code{residuals}
#' extract useful features of the value returned by \code{forecast.StructTS}.
#'
#' An object of class \code{"forecast"} is a list containing at least the
#' following elements: \item{model}{A list containing information about the
#' fitted model} \item{method}{The name of the forecasting method as a
#' character string} \item{mean}{Point forecasts as a time series}
#' \item{lower}{Lower limits for prediction intervals} \item{upper}{Upper
#' limits for prediction intervals} \item{level}{The confidence values
#' associated with the prediction intervals} \item{x}{The original time series
#' (either \code{object} itself or the time series used to create the model
#' stored as \code{object}).} \item{residuals}{Residuals from the fitted model.
#' That is x minus fitted values.} \item{fitted}{Fitted values (one-step
#' forecasts)}
#' @author Rob J Hyndman
#' @seealso \code{\link[stats]{StructTS}}.
#' @keywords ts
#' @examples
#' fit <- StructTS(WWWusage,"level")
#' plot(forecast(fit))
#'
#' @export
forecast.StructTS <- function(object, h=ifelse(object$coef["epsilon"] > 1e-10, 2 * object$xtsp[3], 10), level=c(80, 95), fan=FALSE, lambda=NULL, biasadj=NULL, ...) {
x <- object$data
pred <- predict(object, n.ahead = h)
if (fan) {
level <- seq(51, 99, by = 3)
} else {
if (min(level) > 0 && max(level) < 1) {
level <- 100 * level
} else if (min(level) < 0 || max(level) > 99.99) {
stop("Confidence limit out of range")
}
}
nint <- length(level)
upper <- lower <- matrix(NA, ncol = nint, nrow = length(pred$pred))
for (i in 1:nint)
{
qq <- qnorm(0.5 * (1 + level[i] / 100))
lower[, i] <- pred$pred - qq * pred$se
upper[, i] <- pred$pred + qq * pred$se
}
colnames(lower) <- colnames(upper) <- paste(level, "%", sep = "")
if (is.element("seas", names(object$coef))) {
method <- "Basic structural model"
} else if (is.element("slope", names(object$coef))) {
method <- "Local linear structural model"
} else {
method <- "Local level structural model"
}
# Compute fitted values and residuals
sigma2 <- c(predict(object, n.ahead=1)$se)
res <- residuals(object) * sigma2
fits <- x - res
if (!is.null(lambda)) {
fits <- InvBoxCox(fits, lambda)
x <- InvBoxCox(x, lambda)
pred$pred <- InvBoxCox(pred$pred, lambda, biasadj, list(level = level, upper = upper, lower = lower))
lower <- InvBoxCox(lower, lambda)
upper <- InvBoxCox(upper, lambda)
}
mean <- future_msts(x, pred$pred)
lower <- future_msts(x, lower)
upper <- future_msts(x, upper)
fits <- copy_msts(x, fits)
res <- copy_msts(x, res)
return(structure(
list(
method = method, model = object, level = level, mean = pred$pred,
lower = lower, upper = upper, x = x, series = object$series, fitted = fits, residuals = res
),
class = "forecast"
))
}
#' Forecasting using Holt-Winters objects
#'
#' Returns forecasts and other information for univariate Holt-Winters time
#' series models.
#'
#' This function calls \code{\link[stats]{predict.HoltWinters}} and constructs
#' an object of class "\code{forecast}" from the results.
#'
#' It is included for completeness, but the \code{\link{ets}} is recommended
#' for use instead of \code{\link[stats]{HoltWinters}}.
#'
#' @param object An object of class "\code{HoltWinters}". Usually the result of
#' a call to \code{\link[stats]{HoltWinters}}.
#' @param h Number of periods for forecasting
#' @param level Confidence level for prediction intervals.
#' @param fan If TRUE, level is set to seq(51,99,by=3). This is suitable for
#' fan plots.
#' @param ... Other arguments.
#' @inheritParams forecast.ts
#'
#' @return An object of class "\code{forecast}".
#'
#' The function \code{summary} is used to obtain and print a summary of the
#' results, while the function \code{plot} produces a plot of the forecasts and
#' prediction intervals.
#'
#' The generic accessor functions \code{fitted.values} and \code{residuals}
#' extract useful features of the value returned by
#' \code{forecast.HoltWinters}.
#'
#' An object of class \code{"forecast"} is a list containing at least the
#' following elements: \item{model}{A list containing information about the
#' fitted model} \item{method}{The name of the forecasting method as a
#' character string} \item{mean}{Point forecasts as a time series}
#' \item{lower}{Lower limits for prediction intervals} \item{upper}{Upper
#' limits for prediction intervals} \item{level}{The confidence values
#' associated with the prediction intervals} \item{x}{The original time series
#' (either \code{object} itself or the time series used to create the model
#' stored as \code{object}).} \item{residuals}{Residuals from the fitted
#' model.} \item{fitted}{Fitted values (one-step forecasts)}
#' @author Rob J Hyndman
#' @seealso \code{\link[stats]{predict.HoltWinters}},
#' \code{\link[stats]{HoltWinters}}.
#' @keywords ts
#' @examples
#' fit <- HoltWinters(WWWusage,gamma=FALSE)
#' plot(forecast(fit))
#'
#' @export
forecast.HoltWinters <- function(object, h=ifelse(frequency(object$x) > 1, 2 * frequency(object$x), 10),
level=c(80, 95), fan=FALSE, lambda=NULL, biasadj=NULL, ...) {
x <- object$x
if (!is.null(object$exponential)) {
if (object$exponential) {
stop("Forecasting for exponential trend not yet implemented.")
}
}
if (fan) {
level <- seq(51, 99, by = 3)
} else {
if (min(level) > 0 && max(level) < 1) {
level <- 100 * level
} else if (min(level) < 0 || max(level) > 99.99) {
stop("Confidence limit out of range")
}
}
nint <- length(level)
pred <- predict(object, n.ahead = h, prediction.interval = TRUE, level = level[1] / 100)
pmean <- pred[, 1]
upper <- lower <- matrix(NA, ncol = nint, nrow = length(pred[, 1]))
se <- (pred[, 2] - pred[, 3]) / (2 * qnorm(0.5 * (1 + level[1] / 100)))
for (i in 1:nint)
{
qq <- qnorm(0.5 * (1 + level[i] / 100))
lower[, i] <- pmean - qq * se
upper[, i] <- pmean + qq * se
}
colnames(lower) <- colnames(upper) <- paste(level, "%", sep = "")
if (!is.null(lambda)) {
fitted <- InvBoxCox(object$fitted[, 1], lambda)
x <- InvBoxCox(x, lambda)
pmean <- InvBoxCox(pmean, lambda, biasadj, list(level = level, upper = upper, lower = lower))
lower <- InvBoxCox(lower, lambda)
upper <- InvBoxCox(upper, lambda)
}
else {
fitted <- object$fitted[, 1]
}
# Pad fitted values with NAs
nf <- length(fitted)
n <- length(x)
fitted <- ts(c(rep(NA, n - nf), fitted))
fitted <- copy_msts(object$x, fitted)
pmean <- future_msts(object$x, pmean)
lower <- future_msts(object$x, lower)
upper <- future_msts(object$x, upper)
return(structure(
list(
method = "HoltWinters", model = object, level = level,
mean = pmean, lower = lower, upper = upper, x = x, series = deparse(object$call$x),
fitted = fitted, residuals = x - fitted
),
class = "forecast"
))
}
## CROSTON
#' Forecasts for intermittent demand using Croston's method
#'
#' Returns forecasts and other information for Croston's forecasts applied to
#' y.
#'
#' Based on Croston's (1972) method for intermittent demand forecasting, also
#' described in Shenstone and Hyndman (2005). Croston's method involves using
#' simple exponential smoothing (SES) on the non-zero elements of the time
#' series and a separate application of SES to the times between non-zero
#' elements of the time series. The smoothing parameters of the two
#' applications of SES are assumed to be equal and are denoted by \code{alpha}.
#'
#' Note that prediction intervals are not computed as Croston's method has no
#' underlying stochastic model.
#'
#' @param y a numeric vector or time series of class \code{ts}
#' @param h Number of periods for forecasting.
#' @param alpha Value of alpha. Default value is 0.1.
#' @param x Deprecated. Included for backwards compatibility.
#' @return An object of class \code{"forecast"} is a list containing at least
#' the following elements: \item{model}{A list containing information about the
#' fitted model. The first element gives the model used for non-zero demands.
#' The second element gives the model used for times between non-zero demands.
#' Both elements are of class \code{forecast}.} \item{method}{The name of the
#' forecasting method as a character string} \item{mean}{Point forecasts as a
#' time series} \item{x}{The original time series (either \code{object} itself
#' or the time series used to create the model stored as \code{object}).}
#' \item{residuals}{Residuals from the fitted model. That is y minus fitted
#' values.} \item{fitted}{Fitted values (one-step forecasts)}
#'
#' The function \code{summary} is used to obtain and print a summary of the
#' results, while the function \code{plot} produces a plot of the forecasts.
#'
#' The generic accessor functions \code{fitted.values} and \code{residuals}
#' extract useful features of the value returned by \code{croston} and
#' associated functions.
#' @author Rob J Hyndman
#' @seealso \code{\link{ses}}.
#' @references Croston, J. (1972) "Forecasting and stock control for
#' intermittent demands", \emph{Operational Research Quarterly}, \bold{23}(3),
#' 289-303.
#'
#' Shenstone, L., and Hyndman, R.J. (2005) "Stochastic models underlying
#' Croston's method for intermittent demand forecasting". \emph{Journal of
#' Forecasting}, \bold{24}, 389-402.
#' @keywords ts
#' @examples
#' y <- rpois(20,lambda=.3)
#' fcast <- croston(y)
#' plot(fcast)
#'
#' @export
croston <- function(y, h=10, alpha=0.1, x=y) {
if (sum(x < 0) > 0) {
stop("Series should not contain negative values")
}
out <- croston2(x, h, alpha)
out$x <- x
if (!is.null(out$fitted)) {
out$residuals <- x - out$fitted
}
out$method <- "Croston's method"
out$series <- deparse(substitute(y))
return(structure(out, class = "forecast"))
}
croston2 <- function(x, h=10, alpha=0.1, nofits=FALSE) {
x <- as.ts(x)
y <- x[x > 0]
tsp.x <- tsp(x)
freq.x <- tsp.x[3]
start.f <- tsp.x[2] + 1 / freq.x
if (length(y) == 0) # All historical values are equal to zero
{
fc <- ts(rep(0, h), start = start.f, frequency = freq.x)
if (nofits) {
return(fc)
} else {
return(list(mean = fc, fitted = ts(x * 0, start = tsp.x[1], frequency = freq.x)))
}
}
tt <- diff(c(0, (1:length(x))[x > 0])) # Times between non-zero observations
if (length(y) == 1 && length(tt) == 1) # Only one non-zero observation
{
y.f <- list(mean = ts(rep(y, h), start = start.f, frequency = freq.x))
p.f <- list(mean = ts(rep(tt, h), start = start.f, frequency = freq.x))
}
else if (length(y) <= 1 || length(tt) <= 1) { # length(tt)==0 but length(y)>0. How does that happen?
return(list(mean = ts(rep(NA, h), start = start.f, frequency = freq.x)))
} else {
y.f <- ses(y, alpha = alpha, initial = "simple", h = h, PI = FALSE)
p.f <- ses(tt, alpha = alpha, initial = "simple", h = h, PI = FALSE)
}
ratio <- ts(y.f$mean / p.f$mean, start = start.f, frequency = freq.x)
if (nofits) {
return(ratio)
} else {
n <- length(x)
fits <- x * NA
if (n > 1) {
for (i in 1:(n - 1))
fits[i + 1] <- croston2(x[1:i], h = 1, alpha = alpha, nofits = TRUE)
}
ratio <- future_msts(x, ratio)
fits <- copy_msts(x, fits)
return(list(mean = ratio, fitted = fits, model = list(demand = y.f, period = p.f)))
}
}
robjhyndman/forecast documentation built on March 14, 2024, 11:18 p.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 croston2 croston forecast.HoltWinters forecast.StructTS InvBoxCoxf InvBoxCox BoxCox meanf
Documented in BoxCox croston forecast.HoltWinters forecast.StructTS InvBoxCox meanf
# Mean forecast
#' Mean Forecast
#'
#' Returns forecasts and prediction intervals for an iid model applied to y.
#'
#' The iid model is \deqn{Y_t=\mu + Z_t}{Y[t]=mu + Z[t]} where \eqn{Z_t}{Z[t]}
#' is a normal iid error. Forecasts are given by \deqn{Y_n(h)=\mu}{Y[n+h]=mu}
#' where \eqn{\mu}{mu} is estimated by the sample mean.
#'
#' @param y a numeric vector or time series of class \code{ts}
#' @param h Number of periods for forecasting
#' @param level Confidence levels for prediction intervals.
#' @param fan If TRUE, level is set to seq(51,99,by=3). This is suitable for
#' fan plots.
#' @param bootstrap If TRUE, use a bootstrap method to compute prediction intervals.
#' Otherwise, assume a normal distribution.
#' @param npaths Number of bootstrapped sample paths to use if \code{bootstrap==TRUE}.
#' @param x Deprecated. Included for backwards compatibility.
#' @inheritParams forecast.ts
#'
#' @return An object of class "\code{forecast}".
#'
#' The function \code{summary} is used to obtain and print a summary of the
#' results, while the function \code{plot} produces a plot of the forecasts and
#' prediction intervals.
#'
#' The generic accessor functions \code{fitted.values} and \code{residuals}
#' extract useful features of the value returned by \code{meanf}.
#'
#' An object of class \code{"forecast"} is a list containing at least the
#' following elements: \item{model}{A list containing information about the
#' fitted model} \item{method}{The name of the forecasting method as a
#' character string} \item{mean}{Point forecasts as a time series}
#' \item{lower}{Lower limits for prediction intervals} \item{upper}{Upper
#' limits for prediction intervals} \item{level}{The confidence values
#' associated with the prediction intervals} \item{x}{The original time series
#' (either \code{object} itself or the time series used to create the model
#' stored as \code{object}).} \item{residuals}{Residuals from the fitted model.
#' That is x minus fitted values.} \item{fitted}{Fitted values (one-step
#' forecasts)}
#' @author Rob J Hyndman
#' @seealso \code{\link{rwf}}
#' @keywords ts
#' @examples
#' nile.fcast <- meanf(Nile, h=10)
#' plot(nile.fcast)
#'
#' @export
meanf <- function(y, h=10, level=c(80, 95), fan=FALSE, lambda=NULL, biasadj=FALSE,
bootstrap=FALSE, npaths=5000, x=y) {
n <- length(x)
if (!is.null(lambda)) {
origx <- x
x <- BoxCox(x, lambda)
lambda <- attr(x, "lambda")
}
meanx <- mean(x, na.rm = TRUE)
fits <- rep(meanx, length(x))
res <- x - fits
f <- rep(meanx, h)
if (fan) {
level <- seq(51, 99, by = 3)
} else {
if (min(level) > 0 && max(level) < 1) {
level <- 100 * level
} else if (min(level) < 0 || max(level) > 99.99) {
stop("Confidence limit out of range")
}
}
nconf <- length(level)
s <- sd(x, na.rm = TRUE)
if (bootstrap) {
e <- na.omit(res) - mean(res, na.rm = TRUE)
sim <- matrix(sample(e, size = npaths * h, replace = TRUE), ncol = npaths, nrow = h)
sim <- sweep(sim, 1, f, "+")
lower <- t(apply(sim, 1, quantile, prob = .5 - level / 200))
upper <- t(apply(sim, 1, quantile, prob = .5 + level / 200))
}
else {
lower <- upper <- matrix(NA, nrow = h, ncol = nconf)
for (i in 1:nconf)
{
if (n > 1) {
tfrac <- qt(0.5 - level[i] / 200, n - 1)
} else {
tfrac <- -Inf
}
w <- -tfrac * s * sqrt(1 + 1 / n)
lower[, i] <- f - w
upper[, i] <- f + w
}
}
colnames(lower) <- colnames(upper) <- paste(level, "%", sep = "")
if (is.ts(x)) {
fits <- copy_msts(x, fits)
res <- copy_msts(x, res)
f <- future_msts(x, f)
lower <- future_msts(x, lower)
upper <- future_msts(x, upper)
}
if (!is.null(lambda)) {
fits <- InvBoxCox(fits, lambda)
x <- origx
f <- InvBoxCox(f, lambda, biasadj, list(level = level, upper = upper, lower = lower))
lower <- InvBoxCox(lower, lambda)
upper <- InvBoxCox(upper, lambda)
}
out <- list(
method = "Mean", level = level, x = x, series = deparse(substitute(y)), mean = f, lower = lower, upper = upper,
model = structure(list(mu = f[1], mu.se = s / sqrt(length(x)), sd = s, bootstrap = bootstrap), class = "meanf"), lambda = lambda, fitted = fits, residuals = res
)
out$model$call <- match.call()
return(structure(out, class = "forecast"))
}
#' Box Cox Transformation
#'
#' BoxCox() returns a transformation of the input variable using a Box-Cox
#' transformation. InvBoxCox() reverses the transformation.
#'
#' The Box-Cox transformation (as given by Bickel & Doksum 1981) is given by \deqn{f_\lambda(x) =(sign(x)|x|^\lambda -
#' 1)/\lambda}{f(x;lambda)=(sign(x)|x|^lambda - 1)/lambda} if \eqn{\lambda\ne0}{lambda
#' is not equal to 0}. For \eqn{\lambda=0}{lambda=0},
#' \deqn{f_0(x)=\log(x)}{f(x;0)=log(x)}.
#'
#' @param x a numeric vector or time series of class \code{ts}.
#' @param lambda transformation parameter. If \code{lambda = "auto"}, then
#' the transformation parameter lambda is chosen using BoxCox.lambda (with a lower bound of -0.9)
#' @param biasadj Use adjusted back-transformed mean for Box-Cox
#' transformations. If transformed data is used to produce forecasts and fitted values,
#' a regular back transformation will result in median forecasts. If biasadj is TRUE,
#' an adjustment will be made to produce mean forecasts and fitted values.
#' @param fvar Optional parameter required if biasadj=TRUE. Can either be the
#' forecast variance, or a list containing the interval \code{level}, and the
#' corresponding \code{upper} and \code{lower} intervals.
#' @return a numeric vector of the same length as x.
#' @author Rob J Hyndman & Mitchell O'Hara-Wild
#' @seealso \code{\link{BoxCox.lambda}}
#' @references Box, G. E. P. and Cox, D. R. (1964) An analysis of
#' transformations. \emph{JRSS B} \bold{26} 211--246.
#' Bickel, P. J. and Doksum K. A. (1981) An Analysis of Transformations Revisited. \emph{JASA} \bold{76} 296-311.
#' @keywords ts
#' @examples
#'
#' lambda <- BoxCox.lambda(lynx)
#' lynx.fit <- ar(BoxCox(lynx,lambda))
#' plot(forecast(lynx.fit,h=20,lambda=lambda))
#'
#' @export
BoxCox <- function(x, lambda) {
if (lambda == "auto") {
lambda <- BoxCox.lambda(x, lower = -0.9)
}
if (lambda < 0) {
x[x < 0] <- NA
}
if (lambda == 0) {
out <- log(x)
} else {
out <- (sign(x) * abs(x) ^ lambda - 1) / lambda
}
if (!is.null(colnames(x))) {
colnames(out) <- colnames(x)
}
attr(out, "lambda") <- lambda
return(out)
}
#' @rdname BoxCox
#' @export
InvBoxCox <- function(x, lambda, biasadj=FALSE, fvar=NULL) {
if (lambda < 0) {
x[x > -1 / lambda] <- NA
}
if (lambda == 0) {
out <- exp(x)
} else {
xx <- x * lambda + 1
out <- sign(xx) * abs(xx) ^ (1 / lambda)
}
if (!is.null(colnames(x))) {
colnames(out) <- colnames(x)
}
if (is.null(biasadj)) {
biasadj <- attr(lambda, "biasadj")
}
if (!is.logical(biasadj)) {
warning("biasadj information not found, defaulting to FALSE.")
biasadj <- FALSE
}
if (biasadj) {
if (is.null(fvar)) {
stop("fvar must be provided when biasadj=TRUE")
}
if (is.list(fvar)) { # Create fvar from forecast interval
level <- max(fvar$level)
if (NCOL(fvar$upper) > 1 && NCOL(fvar$lower)) {
i <- match(level, fvar$level)
fvar$upper <- fvar$upper[, i]
fvar$lower <- fvar$lower[, i]
}
if (level > 1) {
level <- level / 100
}
level <- mean(c(level, 1))
# Note: Use BoxCox transformed upper and lower values
fvar <- as.numeric((fvar$upper - fvar$lower) / stats::qnorm(level) / 2) ^ 2
}
if (NCOL(fvar) > 1) {
fvar <- diag(fvar)
}
out <- out * (1 + 0.5 * as.numeric(fvar) * (1 - lambda) / (out) ^ (2 * lambda))
}
return(out)
}
# Deprecated
InvBoxCoxf <- function(x=NULL, fvar=NULL, lambda=NULL) {
message("Deprecated, use InvBoxCox instead")
if (is.null(lambda)) {
stop("Must specify lambda using lambda=numeric(1)")
}
if (is.null(fvar)) {
level <- max(x$level)
if (NCOL(x$upper) > 1 && NCOL(x$lower)) {
i <- match(level, x$level)
x$upper <- x$upper[, i]
x$lower <- x$lower[, i]
}
if (level > 1) {
level <- level / 100
}
level <- mean(c(level, 1))
# Note: Use BoxCox transformed upper and lower values
fvar <- ((x$upper - x$lower) / stats::qnorm(level) / 2) ^ 2
}
else {
x <- list(mean = x)
}
if ("matrix" %in% class(fvar)) {
fvar <- diag(fvar)
}
return(x$mean * (1 + 0.5 * fvar * (1 - lambda) / (x$mean) ^ (2 * lambda)))
}
#' Forecasting using Structural Time Series models
#'
#' Returns forecasts and other information for univariate structural time
#' series models.
#'
#' This function calls \code{predict.StructTS} and constructs an object of
#' class "\code{forecast}" from the results.
#'
#' @param object An object of class "\code{StructTS}". Usually the result of a
#' call to \code{\link[stats]{StructTS}}.
#' @param h Number of periods for forecasting
#' @param level Confidence level for prediction intervals.
#' @param fan If TRUE, level is set to seq(51,99,by=3). This is suitable for
#' fan plots.
#' @param ... Other arguments.
#' @inheritParams forecast.ts
#'
#' @return An object of class "\code{forecast}".
#'
#' The function \code{summary} is used to obtain and print a summary of the
#' results, while the function \code{plot} produces a plot of the forecasts and
#' prediction intervals.
#'
#' The generic accessor functions \code{fitted.values} and \code{residuals}
#' extract useful features of the value returned by \code{forecast.StructTS}.
#'
#' An object of class \code{"forecast"} is a list containing at least the
#' following elements: \item{model}{A list containing information about the
#' fitted model} \item{method}{The name of the forecasting method as a
#' character string} \item{mean}{Point forecasts as a time series}
#' \item{lower}{Lower limits for prediction intervals} \item{upper}{Upper
#' limits for prediction intervals} \item{level}{The confidence values
#' associated with the prediction intervals} \item{x}{The original time series
#' (either \code{object} itself or the time series used to create the model
#' stored as \code{object}).} \item{residuals}{Residuals from the fitted model.
#' That is x minus fitted values.} \item{fitted}{Fitted values (one-step
#' forecasts)}
#' @author Rob J Hyndman
#' @seealso \code{\link[stats]{StructTS}}.
#' @keywords ts
#' @examples
#' fit <- StructTS(WWWusage,"level")
#' plot(forecast(fit))
#'
#' @export
forecast.StructTS <- function(object, h=ifelse(object$coef["epsilon"] > 1e-10, 2 * object$xtsp[3], 10), level=c(80, 95), fan=FALSE, lambda=NULL, biasadj=NULL, ...) {
x <- object$data
pred <- predict(object, n.ahead = h)
if (fan) {
level <- seq(51, 99, by = 3)
} else {
if (min(level) > 0 && max(level) < 1) {
level <- 100 * level
} else if (min(level) < 0 || max(level) > 99.99) {
stop("Confidence limit out of range")
}
}
nint <- length(level)
upper <- lower <- matrix(NA, ncol = nint, nrow = length(pred$pred))
for (i in 1:nint)
{
qq <- qnorm(0.5 * (1 + level[i] / 100))
lower[, i] <- pred$pred - qq * pred$se
upper[, i] <- pred$pred + qq * pred$se
}
colnames(lower) <- colnames(upper) <- paste(level, "%", sep = "")
if (is.element("seas", names(object$coef))) {
method <- "Basic structural model"
} else if (is.element("slope", names(object$coef))) {
method <- "Local linear structural model"
} else {
method <- "Local level structural model"
}
# Compute fitted values and residuals
sigma2 <- c(predict(object, n.ahead=1)$se)
res <- residuals(object) * sigma2
fits <- x - res
if (!is.null(lambda)) {
fits <- InvBoxCox(fits, lambda)
x <- InvBoxCox(x, lambda)
pred$pred <- InvBoxCox(pred$pred, lambda, biasadj, list(level = level, upper = upper, lower = lower))
lower <- InvBoxCox(lower, lambda)
upper <- InvBoxCox(upper, lambda)
}
mean <- future_msts(x, pred$pred)
lower <- future_msts(x, lower)
upper <- future_msts(x, upper)
fits <- copy_msts(x, fits)
res <- copy_msts(x, res)
return(structure(
list(
method = method, model = object, level = level, mean = pred$pred,
lower = lower, upper = upper, x = x, series = object$series, fitted = fits, residuals = res
),
class = "forecast"
))
}
#' Forecasting using Holt-Winters objects
#'
#' Returns forecasts and other information for univariate Holt-Winters time
#' series models.
#'
#' This function calls \code{\link[stats]{predict.HoltWinters}} and constructs
#' an object of class "\code{forecast}" from the results.
#'
#' It is included for completeness, but the \code{\link{ets}} is recommended
#' for use instead of \code{\link[stats]{HoltWinters}}.
#'
#' @param object An object of class "\code{HoltWinters}". Usually the result of
#' a call to \code{\link[stats]{HoltWinters}}.
#' @param h Number of periods for forecasting
#' @param level Confidence level for prediction intervals.
#' @param fan If TRUE, level is set to seq(51,99,by=3). This is suitable for
#' fan plots.
#' @param ... Other arguments.
#' @inheritParams forecast.ts
#'
#' @return An object of class "\code{forecast}".
#'
#' The function \code{summary} is used to obtain and print a summary of the
#' results, while the function \code{plot} produces a plot of the forecasts and
#' prediction intervals.
#'
#' The generic accessor functions \code{fitted.values} and \code{residuals}
#' extract useful features of the value returned by
#' \code{forecast.HoltWinters}.
#'
#' An object of class \code{"forecast"} is a list containing at least the
#' following elements: \item{model}{A list containing information about the
#' fitted model} \item{method}{The name of the forecasting method as a
#' character string} \item{mean}{Point forecasts as a time series}
#' \item{lower}{Lower limits for prediction intervals} \item{upper}{Upper
#' limits for prediction intervals} \item{level}{The confidence values
#' associated with the prediction intervals} \item{x}{The original time series
#' (either \code{object} itself or the time series used to create the model
#' stored as \code{object}).} \item{residuals}{Residuals from the fitted
#' model.} \item{fitted}{Fitted values (one-step forecasts)}
#' @author Rob J Hyndman
#' @seealso \code{\link[stats]{predict.HoltWinters}},
#' \code{\link[stats]{HoltWinters}}.
#' @keywords ts
#' @examples
#' fit <- HoltWinters(WWWusage,gamma=FALSE)
#' plot(forecast(fit))
#'
#' @export
forecast.HoltWinters <- function(object, h=ifelse(frequency(object$x) > 1, 2 * frequency(object$x), 10),
level=c(80, 95), fan=FALSE, lambda=NULL, biasadj=NULL, ...) {
x <- object$x
if (!is.null(object$exponential)) {
if (object$exponential) {
stop("Forecasting for exponential trend not yet implemented.")
}
}
if (fan) {
level <- seq(51, 99, by = 3)
} else {
if (min(level) > 0 && max(level) < 1) {
level <- 100 * level
} else if (min(level) < 0 || max(level) > 99.99) {
stop("Confidence limit out of range")
}
}
nint <- length(level)
pred <- predict(object, n.ahead = h, prediction.interval = TRUE, level = level[1] / 100)
pmean <- pred[, 1]
upper <- lower <- matrix(NA, ncol = nint, nrow = length(pred[, 1]))
se <- (pred[, 2] - pred[, 3]) / (2 * qnorm(0.5 * (1 + level[1] / 100)))
for (i in 1:nint)
{
qq <- qnorm(0.5 * (1 + level[i] / 100))
lower[, i] <- pmean - qq * se
upper[, i] <- pmean + qq * se
}
colnames(lower) <- colnames(upper) <- paste(level, "%", sep = "")
if (!is.null(lambda)) {
fitted <- InvBoxCox(object$fitted[, 1], lambda)
x <- InvBoxCox(x, lambda)
pmean <- InvBoxCox(pmean, lambda, biasadj, list(level = level, upper = upper, lower = lower))
lower <- InvBoxCox(lower, lambda)
upper <- InvBoxCox(upper, lambda)
}
else {
fitted <- object$fitted[, 1]
}
# Pad fitted values with NAs
nf <- length(fitted)
n <- length(x)
fitted <- ts(c(rep(NA, n - nf), fitted))
fitted <- copy_msts(object$x, fitted)
pmean <- future_msts(object$x, pmean)
lower <- future_msts(object$x, lower)
upper <- future_msts(object$x, upper)
return(structure(
list(
method = "HoltWinters", model = object, level = level,
mean = pmean, lower = lower, upper = upper, x = x, series = deparse(object$call$x),
fitted = fitted, residuals = x - fitted
),
class = "forecast"
))
}
## CROSTON
#' Forecasts for intermittent demand using Croston's method
#'
#' Returns forecasts and other information for Croston's forecasts applied to
#' y.
#'
#' Based on Croston's (1972) method for intermittent demand forecasting, also
#' described in Shenstone and Hyndman (2005). Croston's method involves using
#' simple exponential smoothing (SES) on the non-zero elements of the time
#' series and a separate application of SES to the times between non-zero
#' elements of the time series. The smoothing parameters of the two
#' applications of SES are assumed to be equal and are denoted by \code{alpha}.
#'
#' Note that prediction intervals are not computed as Croston's method has no
#' underlying stochastic model.
#'
#' @param y a numeric vector or time series of class \code{ts}
#' @param h Number of periods for forecasting.
#' @param alpha Value of alpha. Default value is 0.1.
#' @param x Deprecated. Included for backwards compatibility.
#' @return An object of class \code{"forecast"} is a list containing at least
#' the following elements: \item{model}{A list containing information about the
#' fitted model. The first element gives the model used for non-zero demands.
#' The second element gives the model used for times between non-zero demands.
#' Both elements are of class \code{forecast}.} \item{method}{The name of the
#' forecasting method as a character string} \item{mean}{Point forecasts as a
#' time series} \item{x}{The original time series (either \code{object} itself
#' or the time series used to create the model stored as \code{object}).}
#' \item{residuals}{Residuals from the fitted model. That is y minus fitted
#' values.} \item{fitted}{Fitted values (one-step forecasts)}
#'
#' The function \code{summary} is used to obtain and print a summary of the
#' results, while the function \code{plot} produces a plot of the forecasts.
#'
#' The generic accessor functions \code{fitted.values} and \code{residuals}
#' extract useful features of the value returned by \code{croston} and
#' associated functions.
#' @author Rob J Hyndman
#' @seealso \code{\link{ses}}.
#' @references Croston, J. (1972) "Forecasting and stock control for
#' intermittent demands", \emph{Operational Research Quarterly}, \bold{23}(3),
#' 289-303.
#'
#' Shenstone, L., and Hyndman, R.J. (2005) "Stochastic models underlying
#' Croston's method for intermittent demand forecasting". \emph{Journal of
#' Forecasting}, \bold{24}, 389-402.
#' @keywords ts
#' @examples
#' y <- rpois(20,lambda=.3)
#' fcast <- croston(y)
#' plot(fcast)
#'
#' @export
croston <- function(y, h=10, alpha=0.1, x=y) {
if (sum(x < 0) > 0) {
stop("Series should not contain negative values")
}
out <- croston2(x, h, alpha)
out$x <- x
if (!is.null(out$fitted)) {
out$residuals <- x - out$fitted
}
out$method <- "Croston's method"
out$series <- deparse(substitute(y))
return(structure(out, class = "forecast"))
}
croston2 <- function(x, h=10, alpha=0.1, nofits=FALSE) {
x <- as.ts(x)
y <- x[x > 0]
tsp.x <- tsp(x)
freq.x <- tsp.x[3]
start.f <- tsp.x[2] + 1 / freq.x
if (length(y) == 0) # All historical values are equal to zero
{
fc <- ts(rep(0, h), start = start.f, frequency = freq.x)
if (nofits) {
return(fc)
} else {
return(list(mean = fc, fitted = ts(x * 0, start = tsp.x[1], frequency = freq.x)))
}
}
tt <- diff(c(0, (1:length(x))[x > 0])) # Times between non-zero observations
if (length(y) == 1 && length(tt) == 1) # Only one non-zero observation
{
y.f <- list(mean = ts(rep(y, h), start = start.f, frequency = freq.x))
p.f <- list(mean = ts(rep(tt, h), start = start.f, frequency = freq.x))
}
else if (length(y) <= 1 || length(tt) <= 1) { # length(tt)==0 but length(y)>0. How does that happen?
return(list(mean = ts(rep(NA, h), start = start.f, frequency = freq.x)))
} else {
y.f <- ses(y, alpha = alpha, initial = "simple", h = h, PI = FALSE)
p.f <- ses(tt, alpha = alpha, initial = "simple", h = h, PI = FALSE)
}
ratio <- ts(y.f$mean / p.f$mean, start = start.f, frequency = freq.x)
if (nofits) {
return(ratio)
} else {
n <- length(x)
fits <- x * NA
if (n > 1) {
for (i in 1:(n - 1))
fits[i + 1] <- croston2(x[1:i], h = 1, alpha = alpha, nofits = TRUE)
}
ratio <- future_msts(x, ratio)
fits <- copy_msts(x, fits)
return(list(mean = ratio, fitted = fits, model = list(demand = y.f, period = p.f)))
}
}
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.