1# Mean forecast 2 3 4#' Mean Forecast 5#' 6#' Returns forecasts and prediction intervals for an iid model applied to y. 7#' 8#' The iid model is \deqn{Y_t=\mu + Z_t}{Y[t]=mu + Z[t]} where \eqn{Z_t}{Z[t]} 9#' is a normal iid error. Forecasts are given by \deqn{Y_n(h)=\mu}{Y[n+h]=mu} 10#' where \eqn{\mu}{mu} is estimated by the sample mean. 11#' 12#' @param y a numeric vector or time series of class \code{ts} 13#' @param h Number of periods for forecasting 14#' @param level Confidence levels for prediction intervals. 15#' @param fan If TRUE, level is set to seq(51,99,by=3). This is suitable for 16#' fan plots. 17#' @param bootstrap If TRUE, use a bootstrap method to compute prediction intervals. 18#' Otherwise, assume a normal distribution. 19#' @param npaths Number of bootstrapped sample paths to use if \code{bootstrap==TRUE}. 20#' @param x Deprecated. Included for backwards compatibility. 21#' @inheritParams forecast 22#' 23#' @return An object of class "\code{forecast}". 24#' 25#' The function \code{summary} is used to obtain and print a summary of the 26#' results, while the function \code{plot} produces a plot of the forecasts and 27#' prediction intervals. 28#' 29#' The generic accessor functions \code{fitted.values} and \code{residuals} 30#' extract useful features of the value returned by \code{meanf}. 31#' 32#' An object of class \code{"forecast"} is a list containing at least the 33#' following elements: \item{model}{A list containing information about the 34#' fitted model} \item{method}{The name of the forecasting method as a 35#' character string} \item{mean}{Point forecasts as a time series} 36#' \item{lower}{Lower limits for prediction intervals} \item{upper}{Upper 37#' limits for prediction intervals} \item{level}{The confidence values 38#' associated with the prediction intervals} \item{x}{The original time series 39#' (either \code{object} itself or the time series used to create the model 40#' stored as \code{object}).} \item{residuals}{Residuals from the fitted model. 41#' That is x minus fitted values.} \item{fitted}{Fitted values (one-step 42#' forecasts)} 43#' @author Rob J Hyndman 44#' @seealso \code{\link{rwf}} 45#' @keywords ts 46#' @examples 47#' nile.fcast <- meanf(Nile, h=10) 48#' plot(nile.fcast) 49#' 50#' @export 51meanf <- function(y, h=10, level=c(80, 95), fan=FALSE, lambda=NULL, biasadj=FALSE, 52 bootstrap=FALSE, npaths=5000, x=y) { 53 n <- length(x) 54 if (!is.null(lambda)) { 55 origx <- x 56 x <- BoxCox(x, lambda) 57 lambda <- attr(x, "lambda") 58 } 59 meanx <- mean(x, na.rm = TRUE) 60 fits <- rep(meanx, length(x)) 61 res <- x - fits 62 f <- rep(meanx, h) 63 if (fan) { 64 level <- seq(51, 99, by = 3) 65 } else { 66 if (min(level) > 0 && max(level) < 1) { 67 level <- 100 * level 68 } else if (min(level) < 0 || max(level) > 99.99) { 69 stop("Confidence limit out of range") 70 } 71 } 72 nconf <- length(level) 73 s <- sd(x, na.rm = TRUE) 74 if (bootstrap) { 75 e <- na.omit(res) - mean(res, na.rm = TRUE) 76 sim <- matrix(sample(e, size = npaths * h, replace = TRUE), ncol = npaths, nrow = h) 77 sim <- sweep(sim, 1, f, "+") 78 lower <- t(apply(sim, 1, quantile, prob = .5 - level / 200)) 79 upper <- t(apply(sim, 1, quantile, prob = .5 + level / 200)) 80 } 81 else { 82 lower <- upper <- matrix(NA, nrow = h, ncol = nconf) 83 for (i in 1:nconf) 84 { 85 if (n > 1) { 86 tfrac <- qt(0.5 - level[i] / 200, n - 1) 87 } else { 88 tfrac <- -Inf 89 } 90 w <- -tfrac * s * sqrt(1 + 1 / n) 91 lower[, i] <- f - w 92 upper[, i] <- f + w 93 } 94 } 95 colnames(lower) <- colnames(upper) <- paste(level, "%", sep = "") 96 if (is.ts(x)) { 97 fits <- ts(fits) 98 res <- ts(res) 99 tsp(fits) <- tsp(res) <- tsp(x) 100 freq <- frequency(x) 101 f <- ts(f, start = tsp(x)[2] + 1 / freq, frequency = freq) 102 lower <- ts(lower, start = tsp(x)[2] + 1 / freq, frequency = freq) 103 upper <- ts(upper, start = tsp(x)[2] + 1 / freq, frequency = freq) 104 } 105 106 if (!is.null(lambda)) { 107 fits <- InvBoxCox(fits, lambda) 108 x <- origx 109 f <- InvBoxCox(f, lambda, biasadj, list(level = level, upper = upper, lower = lower)) 110 lower <- InvBoxCox(lower, lambda) 111 upper <- InvBoxCox(upper, lambda) 112 } 113 114 out <- list( 115 method = "Mean", level = level, x = x, series = deparse(substitute(y)), mean = f, lower = lower, upper = upper, 116 model = structure(list(mu = f[1], mu.se = s / sqrt(length(x)), sd = s, bootstrap = bootstrap), class = "meanf"), lambda = lambda, fitted = fits, residuals = res 117 ) 118 out$model$call <- match.call() 119 120 return(structure(out, class = "forecast")) 121} 122 123 124#' Box Cox Transformation 125#' 126#' BoxCox() returns a transformation of the input variable using a Box-Cox 127#' transformation. InvBoxCox() reverses the transformation. 128#' 129#' The Box-Cox transformation (as given by Bickel & Doksum 1981) is given by \deqn{f_\lambda(x) =sign(x)(|x|^\lambda - 130#' 1)/\lambda}{f(x;lambda)=sign(x)(|x|^lambda - 1)/lambda} if \eqn{\lambda\ne0}{lambda 131#' is not equal to 0}. For \eqn{\lambda=0}{lambda=0}, 132#' \deqn{f_0(x)=\log(x)}{f(x;0)=log(x)}. 133#' 134#' @param x a numeric vector or time series of class \code{ts}. 135#' @param lambda transformation parameter. If \code{lambda = "auto"}, then 136#' the transformation parameter lambda is chosen using BoxCox.lambda (with a lower bound of -0.9) 137#' @param biasadj Use adjusted back-transformed mean for Box-Cox 138#' transformations. If transformed data is used to produce forecasts and fitted values, 139#' a regular back transformation will result in median forecasts. If biasadj is TRUE, 140#' an adjustment will be made to produce mean forecasts and fitted values. 141#' @param fvar Optional parameter required if biasadj=TRUE. Can either be the 142#' forecast variance, or a list containing the interval \code{level}, and the 143#' corresponding \code{upper} and \code{lower} intervals. 144#' @return a numeric vector of the same length as x. 145#' @author Rob J Hyndman & Mitchell O'Hara-Wild 146#' @seealso \code{\link{BoxCox.lambda}} 147#' @references Box, G. E. P. and Cox, D. R. (1964) An analysis of 148#' transformations. \emph{JRSS B} \bold{26} 211--246. 149#' Bickel, P. J. and Doksum K. A. (1981) An Analysis of Transformations Revisited. \emph{JASA} \bold{76} 296-311. 150#' @keywords ts 151#' @examples 152#' 153#' lambda <- BoxCox.lambda(lynx) 154#' lynx.fit <- ar(BoxCox(lynx,lambda)) 155#' plot(forecast(lynx.fit,h=20,lambda=lambda)) 156#' 157#' @export 158BoxCox <- function(x, lambda) { 159 if (lambda == "auto") { 160 lambda <- BoxCox.lambda(x, lower = -0.9) 161 } 162 if (lambda < 0) { 163 x[x < 0] <- NA 164 } 165 if (lambda == 0) { 166 out <- log(x) 167 } else { 168 out <- (sign(x) * abs(x) ^ lambda - 1) / lambda 169 } 170 if (!is.null(colnames(x))) { 171 colnames(out) <- colnames(x) 172 } 173 attr(out, "lambda") <- lambda 174 return(out) 175} 176 177#' @rdname BoxCox 178#' @export 179InvBoxCox <- function(x, lambda, biasadj=FALSE, fvar=NULL) { 180 if (lambda < 0) { 181 x[x > -1 / lambda] <- NA 182 } 183 if (lambda == 0) { 184 out <- exp(x) 185 } else { 186 xx <- x * lambda + 1 187 out <- sign(xx) * abs(xx) ^ (1 / lambda) 188 } 189 if (!is.null(colnames(x))) { 190 colnames(out) <- colnames(x) 191 } 192 193 if (is.null(biasadj)) { 194 biasadj <- attr(lambda, "biasadj") 195 } 196 if (!is.logical(biasadj)) { 197 warning("biasadj information not found, defaulting to FALSE.") 198 biasadj <- FALSE 199 } 200 if (biasadj) { 201 if (is.null(fvar)) { 202 stop("fvar must be provided when biasadj=TRUE") 203 } 204 if (is.list(fvar)) { # Create fvar from forecast interval 205 level <- max(fvar$level) 206 if (NCOL(fvar$upper) > 1 && NCOL(fvar$lower)) { 207 i <- match(level, fvar$level) 208 fvar$upper <- fvar$upper[, i] 209 fvar$lower <- fvar$lower[, i] 210 } 211 if (level > 1) { 212 level <- level / 100 213 } 214 level <- mean(c(level, 1)) 215 # Note: Use BoxCox transformed upper and lower values 216 fvar <- as.numeric((fvar$upper - fvar$lower) / stats::qnorm(level) / 2) ^ 2 217 } 218 if (NCOL(fvar) > 1) { 219 fvar <- diag(fvar) 220 } 221 out <- out * (1 + 0.5 * as.numeric(fvar) * (1 - lambda) / (out) ^ (2 * lambda)) 222 } 223 return(out) 224} 225 226# Deprecated 227InvBoxCoxf <- function(x=NULL, fvar=NULL, lambda=NULL) { 228 message("Deprecated, use InvBoxCox instead") 229 if (is.null(lambda)) { 230 stop("Must specify lambda using lambda=numeric(1)") 231 } 232 if (is.null(fvar)) { 233 level <- max(x$level) 234 if (NCOL(x$upper) > 1 && NCOL(x$lower)) { 235 i <- match(level, x$level) 236 x$upper <- x$upper[, i] 237 x$lower <- x$lower[, i] 238 } 239 if (level > 1) { 240 level <- level / 100 241 } 242 level <- mean(c(level, 1)) 243 # Note: Use BoxCox transformed upper and lower values 244 fvar <- ((x$upper - x$lower) / stats::qnorm(level) / 2) ^ 2 245 } 246 else { 247 x <- list(mean = x) 248 } 249 if ("matrix" %in% class(fvar)) { 250 fvar <- diag(fvar) 251 } 252 253 return(x$mean * (1 + 0.5 * fvar * (1 - lambda) / (x$mean) ^ (2 * lambda))) 254} 255 256 257 258#' Forecasting using Structural Time Series models 259#' 260#' Returns forecasts and other information for univariate structural time 261#' series models. 262#' 263#' This function calls \code{predict.StructTS} and constructs an object of 264#' class "\code{forecast}" from the results. 265#' 266#' @param object An object of class "\code{StructTS}". Usually the result of a 267#' call to \code{\link[stats]{StructTS}}. 268#' @param h Number of periods for forecasting 269#' @param level Confidence level for prediction intervals. 270#' @param fan If TRUE, level is set to seq(51,99,by=3). This is suitable for 271#' fan plots. 272#' @param ... Other arguments. 273#' @inheritParams forecast 274#' 275#' @return An object of class "\code{forecast}". 276#' 277#' The function \code{summary} is used to obtain and print a summary of the 278#' results, while the function \code{plot} produces a plot of the forecasts and 279#' prediction intervals. 280#' 281#' The generic accessor functions \code{fitted.values} and \code{residuals} 282#' extract useful features of the value returned by \code{forecast.StructTS}. 283#' 284#' An object of class \code{"forecast"} is a list containing at least the 285#' following elements: \item{model}{A list containing information about the 286#' fitted model} \item{method}{The name of the forecasting method as a 287#' character string} \item{mean}{Point forecasts as a time series} 288#' \item{lower}{Lower limits for prediction intervals} \item{upper}{Upper 289#' limits for prediction intervals} \item{level}{The confidence values 290#' associated with the prediction intervals} \item{x}{The original time series 291#' (either \code{object} itself or the time series used to create the model 292#' stored as \code{object}).} \item{residuals}{Residuals from the fitted model. 293#' That is x minus fitted values.} \item{fitted}{Fitted values (one-step 294#' forecasts)} 295#' @author Rob J Hyndman 296#' @seealso \code{\link[stats]{StructTS}}. 297#' @keywords ts 298#' @examples 299#' fit <- StructTS(WWWusage,"level") 300#' plot(forecast(fit)) 301#' 302#' @export 303forecast.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, ...) { 304 x <- object$data 305 pred <- predict(object, n.ahead = h) 306 if (fan) { 307 level <- seq(51, 99, by = 3) 308 } else { 309 if (min(level) > 0 && max(level) < 1) { 310 level <- 100 * level 311 } else if (min(level) < 0 || max(level) > 99.99) { 312 stop("Confidence limit out of range") 313 } 314 } 315 nint <- length(level) 316 upper <- lower <- matrix(NA, ncol = nint, nrow = length(pred$pred)) 317 for (i in 1:nint) 318 { 319 qq <- qnorm(0.5 * (1 + level[i] / 100)) 320 lower[, i] <- pred$pred - qq * pred$se 321 upper[, i] <- pred$pred + qq * pred$se 322 } 323 colnames(lower) <- colnames(upper) <- paste(level, "%", sep = "") 324 if (is.element("seas", names(object$coef))) { 325 method <- "Basic structural model" 326 } else if (is.element("slope", names(object$coef))) { 327 method <- "Local linear structural model" 328 } else { 329 method <- "Local level structural model" 330 } 331 332 # Compute fitted values and residuals 333 sigma2 <- c(predict(object, n.ahead=1)$se) 334 res <- residuals(object) * sigma2 335 fits <- x - res 336 337 if (!is.null(lambda)) { 338 fits <- InvBoxCox(fits, lambda) 339 x <- InvBoxCox(x, lambda) 340 pred$pred <- InvBoxCox(pred$pred, lambda, biasadj, list(level = level, upper = upper, lower = lower)) 341 lower <- InvBoxCox(lower, lambda) 342 upper <- InvBoxCox(upper, lambda) 343 } 344 345 return(structure( 346 list( 347 method = method, model = object, level = level, mean = pred$pred, 348 lower = lower, upper = upper, x = x, series = object$series, fitted = fits, residuals = res 349 ), 350 class = "forecast" 351 )) 352} 353 354#' Forecasting using Holt-Winters objects 355#' 356#' Returns forecasts and other information for univariate Holt-Winters time 357#' series models. 358#' 359#' This function calls \code{\link[stats]{predict.HoltWinters}} and constructs 360#' an object of class "\code{forecast}" from the results. 361#' 362#' It is included for completeness, but the \code{\link{ets}} is recommended 363#' for use instead of \code{\link[stats]{HoltWinters}}. 364#' 365#' @param object An object of class "\code{HoltWinters}". Usually the result of 366#' a call to \code{\link[stats]{HoltWinters}}. 367#' @param h Number of periods for forecasting 368#' @param level Confidence level for prediction intervals. 369#' @param fan If TRUE, level is set to seq(51,99,by=3). This is suitable for 370#' fan plots. 371#' @param ... Other arguments. 372#' @inheritParams forecast 373#' 374#' @return An object of class "\code{forecast}". 375#' 376#' The function \code{summary} is used to obtain and print a summary of the 377#' results, while the function \code{plot} produces a plot of the forecasts and 378#' prediction intervals. 379#' 380#' The generic accessor functions \code{fitted.values} and \code{residuals} 381#' extract useful features of the value returned by 382#' \code{forecast.HoltWinters}. 383#' 384#' An object of class \code{"forecast"} is a list containing at least the 385#' following elements: \item{model}{A list containing information about the 386#' fitted model} \item{method}{The name of the forecasting method as a 387#' character string} \item{mean}{Point forecasts as a time series} 388#' \item{lower}{Lower limits for prediction intervals} \item{upper}{Upper 389#' limits for prediction intervals} \item{level}{The confidence values 390#' associated with the prediction intervals} \item{x}{The original time series 391#' (either \code{object} itself or the time series used to create the model 392#' stored as \code{object}).} \item{residuals}{Residuals from the fitted 393#' model.} \item{fitted}{Fitted values (one-step forecasts)} 394#' @author Rob J Hyndman 395#' @seealso \code{\link[stats]{predict.HoltWinters}}, 396#' \code{\link[stats]{HoltWinters}}. 397#' @keywords ts 398#' @examples 399#' fit <- HoltWinters(WWWusage,gamma=FALSE) 400#' plot(forecast(fit)) 401#' 402#' @export 403forecast.HoltWinters <- function(object, h=ifelse(frequency(object$x) > 1, 2 * frequency(object$x), 10), 404 level=c(80, 95), fan=FALSE, lambda=NULL, biasadj=NULL, ...) { 405 x <- object$x 406 if (!is.null(object$exponential)) { 407 if (object$exponential) { 408 stop("Forecasting for exponential trend not yet implemented.") 409 } 410 } 411 412 if (fan) { 413 level <- seq(51, 99, by = 3) 414 } else { 415 if (min(level) > 0 && max(level) < 1) { 416 level <- 100 * level 417 } else if (min(level) < 0 || max(level) > 99.99) { 418 stop("Confidence limit out of range") 419 } 420 } 421 nint <- length(level) 422 423 pred <- predict(object, n.ahead = h, prediction.interval = TRUE, level = level[1] / 100) 424 pmean <- pred[, 1] 425 upper <- lower <- matrix(NA, ncol = nint, nrow = length(pred[, 1])) 426 se <- (pred[, 2] - pred[, 3]) / (2 * qnorm(0.5 * (1 + level[1] / 100))) 427 for (i in 1:nint) 428 { 429 qq <- qnorm(0.5 * (1 + level[i] / 100)) 430 lower[, i] <- pmean - qq * se 431 upper[, i] <- pmean + qq * se 432 } 433 colnames(lower) <- colnames(upper) <- paste(level, "%", sep = "") 434 435 436 if (!is.null(lambda)) { 437 fitted <- InvBoxCox(object$fitted[, 1], lambda) 438 x <- InvBoxCox(x, lambda) 439 pmean <- InvBoxCox(pmean, lambda, biasadj, list(level = level, upper = upper, lower = lower)) 440 lower <- InvBoxCox(lower, lambda) 441 upper <- InvBoxCox(upper, lambda) 442 } 443 else { 444 fitted <- object$fitted[, 1] 445 } 446 447 # Pad fitted values with NAs 448 nf <- length(fitted) 449 n <- length(x) 450 fitted <- ts(c(rep(NA, n - nf), fitted)) 451 tsp(fitted) <- tsp(object$x) 452 453 return(structure( 454 list( 455 method = "HoltWinters", model = object, level = level, 456 mean = pmean, lower = lower, upper = upper, x = x, series = deparse(object$call$x), 457 fitted = fitted, residuals = x - fitted 458 ), 459 class = "forecast" 460 )) 461} 462 463 464## CROSTON 465 466#' Forecasts for intermittent demand using Croston's method 467#' 468#' Returns forecasts and other information for Croston's forecasts applied to 469#' y. 470#' 471#' Based on Croston's (1972) method for intermittent demand forecasting, also 472#' described in Shenstone and Hyndman (2005). Croston's method involves using 473#' simple exponential smoothing (SES) on the non-zero elements of the time 474#' series and a separate application of SES to the times between non-zero 475#' elements of the time series. The smoothing parameters of the two 476#' applications of SES are assumed to be equal and are denoted by \code{alpha}. 477#' 478#' Note that prediction intervals are not computed as Croston's method has no 479#' underlying stochastic model. 480#' 481#' @param y a numeric vector or time series of class \code{ts} 482#' @param h Number of periods for forecasting. 483#' @param alpha Value of alpha. Default value is 0.1. 484#' @param x Deprecated. Included for backwards compatibility. 485#' @return An object of class \code{"forecast"} is a list containing at least 486#' the following elements: \item{model}{A list containing information about the 487#' fitted model. The first element gives the model used for non-zero demands. 488#' The second element gives the model used for times between non-zero demands. 489#' Both elements are of class \code{forecast}.} \item{method}{The name of the 490#' forecasting method as a character string} \item{mean}{Point forecasts as a 491#' time series} \item{x}{The original time series (either \code{object} itself 492#' or the time series used to create the model stored as \code{object}).} 493#' \item{residuals}{Residuals from the fitted model. That is y minus fitted 494#' values.} \item{fitted}{Fitted values (one-step forecasts)} 495#' 496#' The function \code{summary} is used to obtain and print a summary of the 497#' results, while the function \code{plot} produces a plot of the forecasts. 498#' 499#' The generic accessor functions \code{fitted.values} and \code{residuals} 500#' extract useful features of the value returned by \code{croston} and 501#' associated functions. 502#' @author Rob J Hyndman 503#' @seealso \code{\link{ses}}. 504#' @references Croston, J. (1972) "Forecasting and stock control for 505#' intermittent demands", \emph{Operational Research Quarterly}, \bold{23}(3), 506#' 289-303. 507#' 508#' Shenstone, L., and Hyndman, R.J. (2005) "Stochastic models underlying 509#' Croston's method for intermittent demand forecasting". \emph{Journal of 510#' Forecasting}, \bold{24}, 389-402. 511#' @keywords ts 512#' @examples 513#' y <- rpois(20,lambda=.3) 514#' fcast <- croston(y) 515#' plot(fcast) 516#' 517#' @export 518croston <- function(y, h=10, alpha=0.1, x=y) { 519 if (sum(x < 0) > 0) { 520 stop("Series should not contain negative values") 521 } 522 out <- croston2(x, h, alpha) 523 out$x <- x 524 if (!is.null(out$fitted)) { 525 out$residuals <- x - out$fitted 526 } 527 out$method <- "Croston's method" 528 out$series <- deparse(substitute(y)) 529 return(structure(out, class = "forecast")) 530} 531 532croston2 <- function(x, h=10, alpha=0.1, nofits=FALSE) { 533 x <- as.ts(x) 534 y <- x[x > 0] 535 tsp.x <- tsp(x) 536 freq.x <- tsp.x[3] 537 start.f <- tsp.x[2] + 1 / freq.x 538 if (length(y) == 0) # All historical values are equal to zero 539 { 540 fc <- ts(rep(0, h), start = start.f, frequency = freq.x) 541 if (nofits) { 542 return(fc) 543 } else { 544 return(list(mean = fc, fitted = ts(x * 0, start = tsp.x[1], frequency = freq.x))) 545 } 546 } 547 tt <- diff(c(0, (1:length(x))[x > 0])) # Times between non-zero observations 548 if (length(y) == 1 && length(tt) == 1) # Only one non-zero observation 549 { 550 y.f <- list(mean = ts(rep(y, h), start = start.f, frequency = freq.x)) 551 p.f <- list(mean = ts(rep(tt, h), start = start.f, frequency = freq.x)) 552 } 553 else if (length(y) <= 1 || length(tt) <= 1) { # length(tt)==0 but length(y)>0. How does that happen? 554 return(list(mean = ts(rep(NA, h), start = start.f, frequency = freq.x))) 555 } else { 556 y.f <- ses(y, alpha = alpha, initial = "simple", h = h, PI = FALSE) 557 p.f <- ses(tt, alpha = alpha, initial = "simple", h = h, PI = FALSE) 558 } 559 ratio <- ts(y.f$mean / p.f$mean, start = start.f, frequency = freq.x) 560 if (nofits) { 561 return(ratio) 562 } else { 563 n <- length(x) 564 fits <- x * NA 565 if (n > 1) { 566 for (i in 1:(n - 1)) 567 fits[i + 1] <- croston2(x[1:i], h = 1, alpha = alpha, nofits = TRUE) 568 } 569 fits <- ts(fits) 570 tsp(fits) <- tsp.x 571 return(list(mean = ratio, fitted = fits, model = list(demand = y.f, period = p.f))) 572 } 573} 574