File: dshw.r

package info (click to toggle)
r-cran-forecast 8.13-1
  • links: PTS, VCS
  • area: main
  • in suites: bullseye
  • size: 2,248 kB
  • sloc: cpp: 975; ansic: 648; sh: 13; makefile: 2
file content (279 lines) | stat: -rw-r--r-- 9,284 bytes parent folder | download
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)
}