1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279
|
####################################################################
## 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 \code{period1=48} for the daily period and \code{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 \code{\link{ets}} function.
#'
#' @param y Either an \code{\link{msts}} object with two seasonal periods or a
#' numeric vector.
#' @param period1 Period of the shorter seasonal period. Only used if \code{y}
#' is not an \code{\link{msts}} object.
#' @param period2 Period of the longer seasonal period. Only used if \code{y}
#' is not an \code{\link{msts}} object.
#' @param h Number of periods for forecasting.
#' @param alpha Smoothing parameter for the level. If \code{NULL}, the
#' parameter is estimated using least squares.
#' @param beta Smoothing parameter for the slope. If \code{NULL}, the parameter
#' is estimated using least squares.
#' @param gamma Smoothing parameter for the first seasonal period. If
#' \code{NULL}, the parameter is estimated using least squares.
#' @param omega Smoothing parameter for the second seasonal period. If
#' \code{NULL}, the parameter is estimated using least squares.
#' @param phi Autoregressive parameter. If \code{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
#'
#' @return An object of class "\code{forecast}" which is a list that includes 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{x}{The original time series.}
#' \item{residuals}{Residuals from the fitted model. That is x 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{dshw}.
#'
#' @author Rob J Hyndman
#' @seealso \code{\link[stats]{HoltWinters}}, \code{\link{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{http://www.exponentialsmoothing.net}.
#' @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,.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 <- deparse(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 (any(class(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 1: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)
}
return(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
return(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)]
}
return(si)
}
|