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