#' Strength of trend and seasonality of a time series
#'
#' Computes various measures of trend and seasonality of a time series based on
#' an STL decomposition. The number of seasonal periods, and the length of the
#' seasonal periods are returned. Also, the strength of seasonality corresponding
#' to each period is estimated. The \code{\link[forecast]{mstl}} function is used
#' to do the decomposition.
#' @param x a univariate time series.
#' @param ... Other arguments are passed to \code{\link[forecast]{mstl}}.
#' @return A vector of numeric values.
#' @author Rob J Hyndman
#' @export
stl_features <- function(x, ...) {
if ("msts" %in% class(x)) {
msts <- attributes(x)$msts
nperiods <- length(msts)
}
else if ("ts" %in% class(x)) {
msts <- frequency(x)
nperiods <- msts > 1
if(length(x) <= 2*msts) {
warning("Insufficient data to compute STL decomposition")
x <- c(x)
}
season <- 0
}
else {
msts <- 1
nperiods <- 0L
season <- 0
}
if(NCOL(x) > 1){
stop("x must be a univariate time series.")
}
trend <- linearity <- curvature <- season <- spike <- peak <- trough <- acfremainder <- NA
# STL fits
stlfit <- forecast::mstl(x, ...)
trend0 <- stlfit[, "Trend"]
remainder <- stlfit[, "Remainder"]
seasonal <- stlfit[, grep("Season", colnames(stlfit)), drop = FALSE]
# When the maximum frequency is dropped
tsp(x) <- tsp(trend0)
# De-trended and de-seasonalized data
detrend <- x - trend0
deseason <- forecast::seasadj(stlfit)
fits <- x - remainder
# Summary stats
n <- length(x)
varx <- var(x, na.rm = TRUE)
vare <- var(remainder, na.rm = TRUE)
vardetrend <- var(detrend, na.rm = TRUE)
vardeseason <- var(deseason, na.rm = TRUE)
nseas <- NCOL(seasonal)
# Measure of trend strength
if(varx < .Machine$double.eps)
trend <- 0
else if (vardeseason / varx < 1e-10) {
trend <- 0
} else {
trend <- max(0, min(1, 1 - vare / vardeseason))
}
if (nseas > 0) {
# Measure of seasonal strength
season <- numeric(nseas)
for (i in seq(nseas))
season[i] <- max(0, min(1, 1 - vare / var(remainder + seasonal[, i], na.rm = TRUE)))
# Find time of peak and trough for each component
peak <- trough <- numeric(nseas)
for (i in seq(nseas))
{
startx <- start(x)[2L] - 1L
pk <- (startx + which.max(seasonal[, i])) %% msts[i]
th <- (startx + which.min(seasonal[, i])) %% msts[i]
peak[i] <- ifelse(pk == 0, msts[i], pk)
trough[i] <- ifelse(th == 0, msts[i], th)
}
}
# Compute measure of spikiness
d <- (remainder - mean(remainder, na.rm = TRUE))^2
varloo <- (vare * (n - 1) - d) / (n - 2)
spike <- var(varloo, na.rm = TRUE)
# Compute measures of linearity and curvature
tren.coef <- coef(lm(trend0 ~ poly(seq(n), degree = min(n-1, 2L))))[2L:3L]
linearity <- tren.coef[1L]
curvature <- tren.coef[2L]
# ACF of remainder
acfremainder <- unname(acf_features(remainder))
# Assemble results
output <- c(
nperiods = nperiods, seasonal_period = msts, trend = trend,
spike = spike, linearity = unname(linearity), curvature = unname(curvature),
e_acf1 = acfremainder[1L], e_acf10 = acfremainder[2L]
)
if (nseas > 0) {
output <- c(output, seasonal_strength = season, peak = peak, trough = trough)
}
return(output)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.