1# File src/library/stats/R/HoltWinters.R 2# Part of the R package, https://www.R-project.org 3# 4# Copyright (C) 2002-2020 The R Core Team 5# 6# This program is free software; you can redistribute it and/or modify 7# it under the terms of the GNU General Public License as published by 8# the Free Software Foundation; either version 2 of the License, or 9# (at your option) any later version. 10# 11# This program is distributed in the hope that it will be useful, 12# but WITHOUT ANY WARRANTY; without even the implied warranty of 13# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 14# GNU General Public License for more details. 15# 16# A copy of the GNU General Public License is available at 17# https://www.R-project.org/Licenses/ 18 19# Originally contributed by David Meyer 20 21HoltWinters <- 22 function (x, 23 # smoothing parameters 24 alpha = NULL, # level 25 beta = NULL, # trend 26 gamma = NULL, # seasonal component 27 seasonal = c("additive", "multiplicative"), 28 start.periods = 2, 29 30 # starting values 31 l.start = NULL, # level 32 b.start = NULL, # trend 33 s.start = NULL, # seasonal components vector of length `period' 34 35 # starting values for optim 36 optim.start = c(alpha = 0.3, beta = 0.1, gamma = 0.1), 37 optim.control = list() 38 ) 39{ 40 x <- as.ts(x) 41 seasonal <- match.arg(seasonal) 42 f <- frequency(x) 43 44 if(!is.null(alpha) && (alpha == 0)) 45 stop ("cannot fit models without level ('alpha' must not be 0 or FALSE)") 46 if(!is.null(abg <- c(alpha, beta, gamma)) && any(abg < 0 | abg > 1)) 47 stop ("'alpha', 'beta' and 'gamma' must be within the unit interval") 48 if(is.null(gamma) || gamma > 0) { 49 if (seasonal == "multiplicative" && any(x == 0)) 50 stop ("data must be non-zero for multiplicative Holt-Winters") 51 if (start.periods < 2) 52 stop ("need at least 2 periods to compute seasonal start values") 53 } 54 55 ## initialization 56 if(!is.null(gamma) && is.logical(gamma) && !gamma) { 57 ## non-seasonal Holt-Winters 58 expsmooth <- !is.null(beta) && is.logical(beta) && !beta 59 if(is.null(l.start)) 60 l.start <- if(expsmooth) x[1L] else x[2L] 61 if(is.null(b.start)) 62 if(is.null(beta) || !is.logical(beta) || beta) 63 b.start <- x[2L] - x[1L] 64 start.time <- 3 - expsmooth 65 s.start <- 0 66 } else { 67 ## seasonal Holt-Winters 68 start.time <- f + 1 69 wind <- start.periods * f 70 71 ## decompose series 72 st <- decompose(ts(x[1L:wind], start = start(x), frequency = f), 73 seasonal) 74 75 if (is.null(l.start) || is.null(b.start)) { 76 ## level & intercept 77 dat <- na.omit(st$trend) 78 cf <- coef(.lm.fit(x=cbind(1,seq_along(dat)), y=dat)) 79 if (is.null(l.start)) l.start <- cf[1L] 80 if (is.null(b.start)) b.start <- cf[2L] 81 } 82 if (is.null(s.start)) s.start <- st$figure 83 } 84 85 ## Call to filtering loop 86 lenx <- as.integer(length(x)) 87 if (is.na(lenx)) stop("invalid length(x)") 88 89 len <- lenx - start.time + 1 90 hw <- function(alpha, beta, gamma) 91 .C(C_HoltWinters, 92 as.double(x), 93 lenx, 94 as.double(max(min(alpha, 1), 0)), 95 as.double(max(min(beta, 1), 0)), 96 as.double(max(min(gamma, 1), 0)), 97 as.integer(start.time), 98 ## no idea why this is so: same as seasonal != "multiplicative" 99 as.integer(! + (seasonal == "multiplicative")), 100 as.integer(f), 101 as.integer(!is.logical(beta) || beta), 102 as.integer(!is.logical(gamma) || gamma), 103 104 a = as.double(l.start), 105 b = as.double(b.start), 106 s = as.double(s.start), 107 108 ## return values 109 SSE = as.double(0), 110 level = double(len + 1L), 111 trend = double(len + 1L), 112 seasonal = double(len + f) 113 ) 114 115 ## if alpha and/or beta and/or gamma are omitted, use optim to find the 116 ## values minimizing the squared prediction error 117 if (is.null(gamma)) { 118 ## optimize gamma 119 if (is.null(alpha)) { 120 ## optimize alpha 121 if (is.null(beta)) { 122 ## optimize beta 123 ## --> optimize alpha, beta, and gamma 124 error <- function (p) hw(p[1L], p[2L], p[3L])$SSE 125 sol <- optim(optim.start, error, method = "L-BFGS-B", 126 lower = c(0, 0, 0), upper = c(1, 1, 1), 127 control = optim.control) 128 if(sol$convergence || any(sol$par < 0 | sol$par > 1)) { 129 if (sol$convergence > 50) { 130 warning(gettextf("optimization difficulties: %s", 131 sol$message), domain = NA) 132 } else stop("optimization failure") 133 } 134 alpha <- sol$par[1L] 135 beta <- sol$par[2L] 136 gamma <- sol$par[3L] 137 } else { 138 ## !optimize beta 139 ## --> optimize alpha and gamma 140 error <- function (p) hw(p[1L], beta, p[2L])$SSE 141 sol <- optim(c(optim.start["alpha"], optim.start["gamma"]), 142 error, method = "L-BFGS-B", 143 lower = c(0, 0), upper = c(1, 1), 144 control = optim.control) 145 if(sol$convergence || any(sol$par < 0 | sol$par > 1)) { 146 if (sol$convergence > 50) { 147 warning(gettextf("optimization difficulties: %s", 148 sol$message), domain = NA) 149 } else stop("optimization failure") 150 } 151 alpha <- sol$par[1L] 152 gamma <- sol$par[2L] 153 } 154 } else { 155 ## !optimize alpha 156 if (is.null(beta)) { 157 ## optimize beta 158 ## --> optimize beta and gamma 159 error <- function (p) hw(alpha, p[1L], p[2L])$SSE 160 sol <- optim(c(optim.start["beta"], optim.start["gamma"]), 161 error, method = "L-BFGS-B", 162 lower = c(0, 0), upper = c(1, 1), 163 control = optim.control) 164 if(sol$convergence || any(sol$par < 0 | sol$par > 1)) { 165 if (sol$convergence > 50) { 166 warning(gettextf("optimization difficulties: %s", 167 sol$message), domain = NA) 168 } else stop("optimization failure") 169 } 170 beta <- sol$par[1L] 171 gamma <- sol$par[2L] 172 } else { 173 ## !optimize beta 174 ## --> optimize gamma 175 error <- function (p) hw(alpha, beta, p)$SSE 176 gamma <- optimize(error, lower = 0, upper = 1)$minimum 177 } 178 } 179 } else { 180 ## !optimize gamma 181 if (is.null(alpha)) { 182 ## optimize alpha 183 if (is.null(beta)) { 184 ## optimize beta 185 ## --> optimize alpha and beta 186 error <- function (p) hw(p[1L], p[2L], gamma)$SSE 187 sol <- optim(c(optim.start["alpha"], optim.start["beta"]), 188 error, method = "L-BFGS-B", 189 lower = c(0, 0), upper = c(1, 1), 190 control = optim.control) 191 if(sol$convergence || any(sol$par < 0 | sol$par > 1)) { 192 if (sol$convergence > 50) { 193 warning(gettextf("optimization difficulties: %s", 194 sol$message), domain = NA) 195 } else stop("optimization failure") 196 } 197 alpha <- sol$par[1L] 198 beta <- sol$par[2L] 199 } else { 200 ## !optimize beta 201 ## --> optimize alpha 202 error <- function (p) hw(p, beta, gamma)$SSE 203 alpha <- optimize(error, lower = 0, upper = 1)$minimum 204 } 205 } else { 206 ## !optimize alpha 207 if(is.null(beta)) { 208 ## optimize beta 209 ## --> optimize beta 210 error <- function (p) hw(alpha, p, gamma)$SSE 211 beta <- optimize(error, lower = 0, upper = 1)$minimum 212 } ## else optimize nothing! 213 } 214 } 215 216 ## get (final) results 217 final.fit <- hw(alpha, beta, gamma) 218 219 ## return fitted values and estimated coefficients along with parameters used 220 fitted <- ts(cbind(xhat = final.fit$level[-len-1], 221 level = final.fit$level[-len-1], 222 trend = if (!is.logical(beta) || beta) 223 final.fit$trend[-len-1], 224 season = if (!is.logical(gamma) || gamma) 225 final.fit$seasonal[1L:len]), 226 start = start(lag(x, k = 1 - start.time)), 227 frequency = frequency(x) 228 ) 229 if (!is.logical(beta) || beta) fitted[,1] <- fitted[,1] + fitted[,"trend"] 230 if (!is.logical(gamma) || gamma) 231 fitted[,1] <- if (seasonal == "multiplicative") 232 fitted[,1] * fitted[,"season"] 233 else 234 fitted[,1] + fitted[,"season"] 235 structure(list(fitted = fitted, 236 x = x, 237 alpha = alpha, 238 beta = beta, 239 gamma = gamma, 240 coefficients = c(a = final.fit$level[len + 1], 241 b = if (!is.logical(beta) || beta) final.fit$trend[len + 1], 242 s = if (!is.logical(gamma) || gamma) final.fit$seasonal[len + 1L:f]), 243 seasonal = seasonal, 244 SSE = final.fit$SSE, 245 call = match.call() 246 ), 247 class = "HoltWinters" 248 ) 249} 250 251## Predictions, optionally with prediction intervals 252predict.HoltWinters <- 253 function (object, n.ahead = 1L, prediction.interval = FALSE, 254 level = 0.95, ...) 255{ 256 f <- frequency(object$x) 257 258 vars <- function(h) { 259 psi <- function(j) 260 object$alpha * (1 + j * object$beta) + 261 (j %% f == 0) * object$gamma * (1 - object$alpha) 262 var(residuals(object)) * if (object$seasonal == "additive") 263 sum(1, (h > 1) * sapply(1L:(h-1), function(j) crossprod(psi(j)))) 264 else { 265 rel <- 1 + (h - 1) %% f 266 sum(sapply(0:(h-1), function(j) crossprod (psi(j) * object$coefficients[2 + rel] / object$coefficients[2 + (rel - j) %% f]))) 267 } 268 } 269 270 ## compute predictions 271 # level 272 fit <- rep(as.vector(object$coefficients[1L]) ,n.ahead) 273 # trend 274 if (!is.logical(object$beta) || object$beta) 275 fit <- fit + as.vector((1L:n.ahead)*object$coefficients[2L]) 276 # seasonal component 277 if (!is.logical(object$gamma) || object$gamma) 278 if (object$seasonal == "additive") 279 fit <- fit + rep(object$coefficients[-(1L:(1+(!is.logical(object$beta) || object$beta)))], 280 length.out=length(fit)) 281 else 282 fit <- fit * rep(object$coefficients[-(1L:(1+(!is.logical(object$beta) || object$beta)))], 283 length.out=length(fit)) 284 285 ## compute prediction intervals 286 if (prediction.interval) 287 int <- qnorm((1 + level) / 2) * sqrt(sapply(1L:n.ahead,vars)) 288 ts( 289 cbind(fit = fit, 290 upr = if(prediction.interval) fit + int, 291 lwr = if(prediction.interval) fit - int 292 ), 293 start = end(lag(fitted(object)[,1], k = -1)), 294 frequency = frequency(fitted(object)[,1]) 295 ) 296} 297 298residuals.HoltWinters <- function (object, ...) object$x - object$fitted[,1] 299 300plot.HoltWinters <- 301 function (x, predicted.values = NA, intervals = TRUE, separator = TRUE, 302 col = 1, col.predicted = 2, col.intervals = 4, col.separator = 1, 303 lty = 1, lty.predicted = 1, lty.intervals = 1, lty.separator = 3, 304 ylab = "Observed / Fitted", main = "Holt-Winters filtering", 305 ylim = NULL, ...) 306{ 307 if (is.null(ylim)) 308 ylim <- range(na.omit(c(fitted(x)[,1], x$x, predicted.values))) 309 310 preds <- length(predicted.values) > 1 || !is.na(predicted.values) 311 312 dev.hold(); on.exit(dev.flush()) 313 ## plot fitted/predicted values 314 plot(ts(c(fitted(x)[,1], if(preds) predicted.values[,1]), 315 start = start(fitted(x)[,1]), 316 frequency = frequency(fitted(x)[,1])), 317 col = col.predicted, 318 ylim = ylim, 319 ylab = ylab, main = main, 320 lty = lty.predicted, 321 ... 322 ) 323 324 ## plot prediction interval 325 if(preds && intervals && ncol(predicted.values) > 1) { 326 lines(predicted.values[,2], col = col.intervals, lty = lty.intervals) 327 lines(predicted.values[,3], col = col.intervals, lty = lty.intervals) 328 } 329 330 ## plot observed values 331 lines(x$x, col = col, lty = lty) 332 333 ## plot separator 334 if (separator && preds) 335 abline (v = time(x$x)[length(x$x)], lty = lty.separator, col = col.separator) 336} 337 338## print function 339print.HoltWinters <- function (x, ...) 340{ 341 cat("Holt-Winters exponential smoothing", 342 if (is.logical(x$beta) && !x$beta) "without" else "with", "trend and", 343 if (is.logical(x$gamma) && !x$gamma) "without" else 344 paste0(if (is.logical(x$beta) && !x$beta) "with ", x$seasonal), 345 "seasonal component.") 346 cat("\n\nCall:\n", deparse (x$call), "\n\n", sep = "") 347 cat("Smoothing parameters:\n") 348 cat(" alpha: ", x$alpha, "\n", sep = "") 349 cat(" beta : ", x$beta, "\n", sep = "") 350 cat(" gamma: ", x$gamma, "\n\n", sep = "") 351 352 cat("Coefficients:\n") 353 print(t(t(x$coefficients))) 354 invisible(x) 355} 356 357# decompose additive/multiplicative series into trend/seasonal figures/noise 358decompose <- 359function (x, type = c("additive", "multiplicative"), filter = NULL) 360{ 361 type <- match.arg(type) 362 l <- length(x) 363 f <- frequency(x) 364 if (f <= 1 || length(na.omit(x)) < 2 * f) 365 stop("time series has no or less than 2 periods") 366 367 ## filter out seasonal components 368 if (is.null(filter)) 369 filter <- if (!f %% 2) 370 c(0.5, rep_len(1, f - 1), 0.5) / f 371 else 372 rep_len(1, f) / f 373 trend <- filter(x, filter) 374 375 ## compute seasonal components 376 season <- if (type == "additive") 377 x - trend 378 else 379 x / trend 380 381 ## average seasonal figures 382 periods <- l %/% f 383 index <- seq.int(1L, l, by = f) - 1L 384 figure <- numeric(f) 385 for (i in 1L:f) 386 figure[i] <- mean(season[index + i], na.rm = TRUE) 387 388 ## normalize figure 389 figure <- if (type == "additive") 390 figure - mean(figure) 391 else figure / mean(figure) 392 393 seasonal <- ts(rep(figure, periods+1)[seq_len(l)], 394 start = start(x), frequency = f) 395 396 ## return values 397 structure(list(x = x, 398 seasonal = seasonal, 399 trend = trend, 400 random = if (type == "additive") 401 x - seasonal - trend 402 else 403 x / seasonal / trend, 404 figure = figure, 405 type = type), 406 class = "decomposed.ts") 407} 408 409plot.decomposed.ts <- function(x, ...) 410{ 411 xx <- x$x %||% # added in 2.14.0 412 with(x, if (type == "additive") random + trend + seasonal 413 else random * trend * seasonal) 414 plot(cbind(observed = xx, 415 trend = x$trend, 416 seasonal = x$seasonal, 417 random = x$random 418 ), 419 main = paste("Decomposition of", x$type, "time series"), 420 ...) 421} 422 423