R/dshw.R

In forecast: Forecasting Functions for Time Series and Linear Models

####################################################################
## Double Seasonal Holt Winters method as per Taylor (2003)
## Periods must be nested.
## y can be an msts object, or periods can be passed explicitly.
####################################################################

#' Double-Seasonal Holt-Winters Forecasting
#'
#' Returns forecasts using Taylor's (2003) Double-Seasonal Holt-Winters method.
#'
#' Taylor's (2003) double-seasonal Holt-Winters method uses additive trend and
#' multiplicative seasonality, where there are two seasonal components which
#' are multiplied together. For example, with a series of half-hourly data, one
#' would set `period1 = 48` for the daily period and `period2 = 336` for
#' the weekly period. The smoothing parameter notation used here is different
#' from that in Taylor (2003); instead it matches that used in Hyndman et al
#' (2008) and that used for the [ets()] function.
#'
#' @param y Either an [msts()] object with two seasonal periods or a
#' numeric vector.
#' @param period1 Period of the shorter seasonal period. Only used if `y`
#' is not an [msts()] object.
#' @param period2 Period of the longer seasonal period.  Only used if `y`
#' is not an [msts()] object.
#' @param h Number of periods for forecasting.
#' @param alpha Smoothing parameter for the level. If `NULL`, the
#' parameter is estimated using least squares.
#' @param beta Smoothing parameter for the slope. If `NULL`, the parameter
#' is estimated using least squares.
#' @param gamma Smoothing parameter for the first seasonal period. If
#' `NULL`, the parameter is estimated using least squares.
#' @param omega Smoothing parameter for the second seasonal period. If
#' `NULL`, the parameter is estimated using least squares.
#' @param phi Autoregressive parameter. If `NULL`, the parameter is
#' estimated using least squares.
#' @param armethod If `TRUE`, the forecasts are adjusted using an AR(1)
#' model for the errors.
#' @param model If it's specified, an existing model is applied to a new data set.
#' @inheritParams forecast.ts
#' @inheritParams BoxCox
#' @return An object of class `forecast`.
#' @inherit forecast.ts format
#' @author Rob J Hyndman
#' @seealso [stats::HoltWinters()], [ets()].
#' @references Taylor, J.W. (2003) Short-term electricity demand forecasting
#' using double seasonal exponential smoothing. \emph{Journal of the
#' Operational Research Society}, \bold{54}, 799-805.
#'
#' 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{https://robjhyndman.com/expsmooth/}.
#' @keywords ts
#' @examples
#'
#' \dontrun{
#' fcast <- dshw(taylor)
#' plot(fcast)
#'
#' t <- seq(0, 5, by = 1 / 20)
#' x <- exp(sin(2 * pi * t) + cos(2 * pi * t * 4) + rnorm(length(t), 0, 0.1))
#' fit <- dshw(x, 20, 5)
#' plot(fit)
#' }
#'
#' @export
dshw <- function(
  y,
  period1 = NULL,
  period2 = NULL,
  h = 2 * max(period1, period2),
  alpha = NULL,
  beta = NULL,
  gamma = NULL,
  omega = NULL,
  phi = NULL,
  lambda = NULL,
  biasadj = FALSE,
  armethod = TRUE,
  model = NULL
) {
  if (min(y, na.rm = TRUE) <= 0) {
    stop("dshw not suitable when data contain zeros or negative numbers")
  }
  seriesname <- deparse1(substitute(y))
  if (!is.null(model) && model$method == "DSHW") {
    period1 <- model$period1
    period2 <- model$period2
  } else if (inherits(y, "msts") && (length(attr(y, "msts")) == 2)) {
    period1 <- as.integer(sort(attr(y, "msts"))[1])
    period2 <- as.integer(sort(attr(y, "msts"))[2])
  } else if (is.null(period1) || is.null(period2)) {
    stop(
      "Error in dshw(): y must either be an msts object with two seasonal periods OR the seasonal periods should be specified with period1= and period2="
    )
  } else {
    if (period1 > period2) {
      tmp <- period2
      period2 <- period1
      period1 <- tmp
    }
  }
  if (!inherits(y, "msts")) {
    y <- msts(y, c(period1, period2))
  }

  if (length(y) < 2 * max(period2)) {
    stop("Insufficient data to estimate model")
  }

  if (!armethod) {
    phi <- 0
  }

  if (period1 < 1 || period1 == period2) {
    stop("Inappropriate periods")
  }
  ratio <- period2 / period1
  if (ratio - trunc(ratio) > 1e-10) {
    stop("Seasonal periods are not nested")
  }

  if (!is.null(model)) {
    lambda <- model$model$lambda
  }

  if (!is.null(lambda)) {
    origy <- y
    y <- BoxCox(y, lambda)
    lambda <- attr(y, "lambda")
  }

  if (!is.null(model)) {
    pars <- model$model
    alpha <- pars$alpha
    beta <- pars$beta
    gamma <- pars$gamma
    omega <- pars$omega
    phi <- pars$phi
  } else {
    pars <- rep(NA, 5)
    if (!is.null(alpha)) {
      pars[1] <- alpha
    }
    if (!is.null(beta)) {
      pars[2] <- beta
    }
    if (!is.null(gamma)) {
      pars[3] <- gamma
    }
    if (!is.null(omega)) {
      pars[4] <- omega
    }
    if (!is.null(phi)) {
      pars[5] <- phi
    }
  }

  # Estimate parameters
  if (sum(is.na(pars)) > 0) {
    pars <- par_dshw(y, period1, period2, pars)
    alpha <- pars[1]
    beta <- pars[2]
    gamma <- pars[3]
    omega <- pars[4]
    phi <- pars[5]
  }

  ## Allocate space
  n <- length(y)
  yhat <- numeric(n)

  ## Starting values
  I <- seasindex(y, period1)
  wstart <- seasindex(y, period2)
  wstart <- wstart / rep(I, ratio)
  w <- wstart
  x <- c(0, diff(y[1:period2]))
  t <- t.start <- mean(
    ((y[1:period2] - y[(period2 + 1):(2 * period2)]) / period2) + x
  ) /
    2
  s <- s.start <- (mean(y[1:(2 * period2)]) - (period2 + 0.5) * t)

  ## In-sample fit
  for (i in seq_len(n)) {
    yhat[i] <- (s + t) * I[i] * w[i]
    snew <- alpha * (y[i] / (I[i] * w[i])) + (1 - alpha) * (s + t)
    tnew <- beta * (snew - s) + (1 - beta) * t
    I[i + period1] <- gamma * (y[i] / (snew * w[i])) + (1 - gamma) * I[i]
    w[i + period2] <- omega * (y[i] / (snew * I[i])) + (1 - omega) * w[i]
    s <- snew
    t <- tnew
  }

  # Forecasts
  fcast <- (s + (1:h) * t) *
    rep(I[n + (1:period1)], h / period1 + 1)[1:h] *
    rep(w[n + (1:period2)], h / period2 + 1)[1:h]
  fcast <- msts(fcast, c(period1, period2), start = tsp(y)[2] + 1 / tsp(y)[3])

  # Calculate MSE and MAPE
  yhat <- ts(yhat)
  tsp(yhat) <- tsp(y)
  yhat <- msts(yhat, c(period1, period2))
  e <- y - yhat
  e <- msts(e, c(period1, period2))
  if (armethod) {
    yhat <- yhat + phi * c(0, e[-n])
    fcast <- fcast + phi^(1:h) * e[n]
    e <- y - yhat
  }
  mse <- mean(e^2)
  mape <- mean(abs(e) / y) * 100

  end.y <- end(y)
  if (end.y[2] == frequency(y)) {
    end.y[1] <- end.y[1] + 1
    end.y[2] <- 1
  } else {
    end.y[2] <- end.y[2] + 1
  }

  fcast <- msts(fcast, c(period1, period2))

  if (!is.null(lambda)) {
    y <- origy
    fcast <- InvBoxCox(fcast, lambda, biasadj, var(e))
    attr(lambda, "biasadj") <- biasadj
    # Does this also need a biasadj backtransform?
    yhat <- InvBoxCox(yhat, lambda)
  }

  structure(
    list(
      mean = fcast,
      method = "DSHW",
      x = y,
      residuals = e,
      fitted = yhat,
      series = seriesname,
      model = list(
        mape = mape,
        mse = mse,
        alpha = alpha,
        beta = beta,
        gamma = gamma,
        omega = omega,
        phi = phi,
        lambda = lambda,
        l0 = s.start,
        b0 = t.start,
        s10 = wstart,
        s20 = I
      ),
      period1 = period1,
      period2 = period2
    ),
    class = "forecast"
  )
}

### Double Seasonal Holt-Winters smoothing parameter optimization

par_dshw <- function(y, period1, period2, pars) {
  start <- c(0.1, 0.01, 0.001, 0.001, 0.0)[is.na(pars)]
  out <- optim(
    start,
    dshw.mse,
    y = y,
    period1 = period1,
    period2 = period2,
    pars = pars
  )
  pars[is.na(pars)] <- out$par
  pars
}

dshw.mse <- function(par, y, period1, period2, pars) {
  pars[is.na(pars)] <- par
  if (max(pars) > 0.99 || min(pars) < 0 || pars[5] > .9) {
    return(Inf)
  } else {
    return(
      dshw(
        y,
        period1,
        period2,
        h = 1,
        pars[1],
        pars[2],
        pars[3],
        pars[4],
        pars[5],
        armethod = (abs(pars[5]) > 1e-7)
      )$model$mse
    )
  }
}

### Calculating seasonal indexes
seasindex <- function(y, p) {
  n <- length(y)
  n2 <- 2 * p
  shorty <- y[1:n2]
  average <- numeric(n)
  simplema <- zoo::rollmean.default(shorty, p)
  if (identical(p %% 2, 0)) {
    # Even order
    centeredma <- zoo::rollmean.default(simplema[1:(n2 - p + 1)], 2)
    average[p / 2 + 1:p] <- shorty[p / 2 + 1:p] / centeredma[1:p]
    si <- average[c(p + (1:(p / 2)), (1 + p / 2):p)]
  } else {
    # Odd order
    average[(p - 1) / 2 + 1:p] <- shorty[(p - 1) / 2 + 1:p] / simplema[1:p]
    si <- average[c(p + (1:((p - 1) / 2)), (1 + (p - 1) / 2):p)]
  }
  si
}