globalVariables("self")
specials_stl <- fabletools::new_specials(
trend = function(window, degree, jump){
args <- as.list(environment())
args <- args[!map_lgl(args, rlang::is_missing)]
if(isFALSE(is.finite(args$window)) && sign(args$window) == 1){
args$window <- NROW(self$data) + 1e7
}
if(length(args) > 0){
set_names(args, paste0("t.", names(args)))
}
},
season = function(period = NULL, window = NULL, degree, jump){
args <- as.list(environment())
args <- args[!map_lgl(args, rlang::is_missing)]
args <- args[names(args)!="period"]
if(isFALSE(is.finite(args$window)) && is.numeric(args$window)){
if(sign(args$window) == 1){
args$window <- "periodic"
}
}
args <- set_names(args, paste0("s.", names(args)))
nm <- period
period <- get_frequencies(period, self$data, .auto = "all")
if(is.null(names(period))) names(period) <- nm
period <- period[NROW(self$data)/period >= 2]
if(!is.null(period)){
map(period, function(.x) c(period = .x, args))
}
},
lowpass = function(window, degree, jump){
args <- as.list(environment())
args <- args[!map_lgl(args, rlang::is_missing)]
if(isFALSE(is.finite(args$window)) && sign(args$window) == 1){
args$window <- NROW(self$data) + 1e7
}
if(length(args) > 0){
set_names(args, paste0("l.", names(args)))
}
},
.required_specials = c("trend", "season", "lowpass")
)
estimate_stl <- function(y, trend.args, season.args, lowpass.args,
iterations = 2, ...){
if(any(is.na(y))){
abort("STL decomposition does not support series with missing values.")
}
deseas <- y
season.args <- Filter(function(x) x[["period"]] > 1, season.args)
period <- map_dbl(season.args, function(x) x[["period"]])
season.args <- map2(season.args, 7 + 4*order(period),
function(x, default_window){
x$s.window <- x$s.window %||% default_window
x
})
season.args <- season.args[order(period)]
seas <- set_names(as.list(rep(0, length(season.args))),
sprintf("season_%s", names(season.args)%||%map(season.args, function(x) x[["period"]])))
if(length(season.args) > 0){
for (j in seq_len(iterations))
{
for (i in seq_along(season.args))
{
deseas <- ts(deseas + seas[[i]], frequency = season.args[[i]][[1]])
fit <- eval_tidy(expr(stl(deseas, !!!c(trend.args, season.args[[i]][-1], lowpass.args), ...)))
seas[[i]] <- as.numeric(fit$time.series[, "seasonal"])
deseas <- deseas - seas[[i]]
}
}
trend <- fit$time.series[, "trend"]
}
else{
trend <- stats::supsmu(seq_along(y), y)$y
}
trend <- as.numeric(trend)
deseas <- as.numeric(deseas)
list2(trend = trend, !!!seas, remainder = deseas - trend, season_adjust = deseas)
}
train_stl <- function(.data, formula, specials, iterations = 2, ...){
stopifnot(is_tsibble(.data))
y <- .data[[measured_vars(.data)]]
trend.args <- specials$trend[[1]]
season.args <- unlist(specials$season, recursive = FALSE)
season.args <- Filter(function(x) x[["period"]] > 1, season.args)
lowpass.args <- specials$lowpass[[1]]
decomposition <- .data %>%
mutate(!!!estimate_stl(y, trend.args, season.args, lowpass.args, ...))
seas_cols <- sprintf("season_%s", names(season.args)%||%map(season.args, function(x) x[["period"]]))
seasonalities <- lapply(season.args, function(x){
x["base"] <- 0
x[c("period", "base")]
})
names(seasonalities) <- seas_cols
aliases <- list2(
!!measured_vars(.data) := reduce(syms(c("trend", seas_cols, "remainder")),
function(x,y) call2("+", x, y)),
season_adjust = call2("+", sym("trend"), sym("remainder"))
)
structure(
list(decomposition = decomposition,
response = measured_vars(.data), method = "STL",
seasons = seasonalities, aliases = aliases
),
class = "stl_decomposition"
)
}
#' @importFrom fabletools components as_dable
#' @export
components.stl_decomposition <- function(object, ...){
as_dable(object[["decomposition"]], response = !!sym(object[["response"]]),
method = object[["method"]], seasons = object[["seasons"]],
aliases = object[["aliases"]])
}
#' @export
fitted.stl_decomposition <- function(object, ...) {
object[["decomposition"]][[object[["response"]]]] - residuals(object)
}
#' @importFrom stats residuals
#' @export
residuals.stl_decomposition <- function(object, ...) {
object[["decomposition"]][["remainder"]]
}
MBB <- function (x, window_size) {
bx <- array(0, (floor(length(x)/window_size) + 2) * window_size)
for (i in 1:(floor(length(x)/window_size) + 2)) {
c <- sample(1:(length(x) - window_size + 1), 1)
bx[((i - 1) * window_size + 1):(i * window_size)] <- x[c:(c + window_size - 1)]
}
start_from <- sample(0:(window_size - 1), 1) + 1
bx[start_from:(start_from + length(x) - 1)]
}
#' Generate block bootstrapped series from an STL decomposition
#'
#' Produces new data with the same structure by resampling the residuals using
#' a block bootstrap procedure. This method can only generate within sample, and
#' any generated data out of the trained sample will produce NA simulations.
#'
#' @inheritParams fable::generate.ARIMA
#'
#' @examples
#' as_tsibble(USAccDeaths) %>%
#' model(STL(log(value))) %>%
#' generate(as_tsibble(USAccDeaths), times = 3)
#'
#' @references
#' Bergmeir, C., R. J. Hyndman, and J. M. Benitez (2016). Bagging Exponential Smoothing Methods using STL Decomposition and Box-Cox Transformation. International Journal of Forecasting 32, 303-312.
#'
#' @importFrom fabletools generate
#'
#' @export
generate.stl_decomposition <- function(x, new_data, specials = NULL, ...){
dcmp <- x$decomposition
# Match new_data index with dcmp index
pos <- vec_match(new_data[[index_var(new_data)]], dcmp[[index_var(dcmp)]])
if(!(".innov" %in% names(new_data))){
# Block bootstrap for each replicate
kr <- tsibble::key_rows(new_data)
# Get default bootstrap params
period <- max(vapply(x$seasons, `[[`, double(1L), "period"))
block_size <- ifelse(period > 1, 2 * period, min(8, floor(length(x)/2)))
# Block bootstrap
innov <- lapply(seq_along(kr), function(...) MBB(dcmp[["remainder"]], block_size))
innov <- mapply(function(i, e) e[pos[i]], kr, innov, SIMPLIFY = FALSE)
new_data$.innov <- vec_c(!!!innov)
}
dcmp <- as_tibble(dcmp)[pos,]
dcmp$remainder <- new_data$.innov
new_data[[".sim"]] <- as.numeric(eval_tidy(x$aliases[[x$response]], dcmp))
new_data[[".innov"]] <- NULL
new_data
}
#' @importFrom fabletools model_sum
#' @export
model_sum.stl_decomposition <- function(x){
"STL"
}
#' Multiple seasonal decomposition by Loess
#'
#' @inherit forecast::mstl
#'
#' @inheritParams classical_decomposition
#' @param iterations Number of iterations to use to refine the seasonal component.
#' @param ... Other arguments passed to [stats::stl()].
#'
#' @return A [`fabletools::dable()`] containing the decomposed trend, seasonality
#' and remainder from the STL decomposition.
#'
#' @section Specials:
#'
#' \subsection{trend}{
#' The `trend` special is used to specify the trend extraction parameters.
#' \preformatted{
#' trend(window, degree, jump)
#' }
#'
#' \tabular{ll}{
#' `window` \tab The span (in lags) of the loess window, which should be odd. If NULL, the default, nextodd(ceiling((1.5*period) / (1-(1.5/s.window)))), is taken.\cr
#' `degree` \tab The degree of locally-fitted polynomial. Should be zero or one. \cr
#' `jump` \tab Integers at least one to increase speed of the respective smoother. Linear interpolation happens between every `jump`th value.
#' }
#' }
#'
#' \subsection{season}{
#' The `season` special is used to specify the season extraction parameters.
#' \preformatted{
#' season(period = NULL, window = NULL, degree, jump)
#' }
#'
#' \tabular{ll}{
#' `period` \tab The periodic nature of the seasonality. This can be either a number indicating the number of observations in each seasonal period, or text to indicate the duration of the seasonal window (for example, annual seasonality would be "1 year").\cr
#' `window` \tab The span (in lags) of the loess window, which should be odd. If the `window` is set to `"periodic"` or `Inf`, the seasonal pattern will be fixed. The window size should be odd and at least 7, according to Cleveland et al. The default (NULL) will choose an appropriate default, for a dataset with one seasonal pattern this would be 11, the second larger seasonal window would be 15, then 19, 23, ... onwards.\cr
#' `degree` \tab The degree of locally-fitted polynomial. Should be zero or one. \cr
#' `jump` \tab Integers at least one to increase speed of the respective smoother. Linear interpolation happens between every `jump`th value.
#' }
#' }
#'
#' \subsection{lowpass}{
#' The `lowpass` special is used to specify the low-pass filter parameters.
#' \preformatted{
#' lowpass(window, degree, jump)
#' }
#'
#' \tabular{ll}{
#' `window` \tab The span (in lags) of the loess window of the low-pass filter used for each subseries. Defaults to the smallest odd integer greater than or equal to the seasonal `period` which is recommended since it prevents competition between the trend and seasonal components. If not an odd integer its given value is increased to the next odd one. \cr
#' `degree` \tab The degree of locally-fitted polynomial. Must be zero or one. \cr
#' `jump` \tab Integers at least one to increase speed of the respective smoother. Linear interpolation happens between every `jump`th value.
#' }
#' }
#'
#' @references
#' R. B. Cleveland, W. S. Cleveland, J.E. McRae, and I. Terpenning (1990) STL: A Seasonal-Trend Decomposition Procedure Based on Loess. Journal of Official Statistics, 6, 3–73.
#'
#' @examples
#' as_tsibble(USAccDeaths) %>%
#' model(STL(value ~ trend(window = 10))) %>%
#' components()
#'
#' @importFrom stats ts stl
#' @importFrom fabletools new_model_class new_model_definition
#' @export
STL <- function(formula, iterations = 2, ...){
dcmp <- new_model_class("STL",
train = train_stl, specials = specials_stl,
check = all_tsbl_checks)
new_model_definition(dcmp, !!enquo(formula), iterations = iterations, ...)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.