--- title: "Applying stlpro" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Applying stlpro} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>" ) ``` ```{r setup} library(stlpro) ``` # Introduction `stlpro` provides an extension of Seasonal-Trend decomposition using LOESS (STL, developed by Cleveland, Cleveland, McRae, and Terpenning, 1990). Like `stats::stl()`, it decomposes a time series into trend, seasonal, and remainder components. In addition, `stlpro` supports: - additive and multiplicative decompositions; - non-integer seasonal periods, such as `52.18` for weekly data; - handling of missing values, - cross-validation of the seasonal smoothing window. The additive decomposition is \[ Y_t = T_t + S_t + R_t, \] where \(T_t\) is the trend component, \(S_t\) is the seasonal component, and \(R_t\) is the remainder. For strictly positive data, `stlpro` can also fit a multiplicative decomposition, \[ Y_t = T_t \times S_t \times R_t. \] This is particularly useful when seasonal variation grows as the level of the series increases, which usually can be observed in economic time series. # Additive decomposition We begin with a simple monthly example. The data contain a smooth trend, a stable annual seasonal pattern, and random noise. ```{r additive-data, fig.width=8, fig.height=4} set.seed(2106) n <- 15 * 12 tt <- seq_len(n) trend <- 100 + 0.20 * tt season <- 12 * sin(2 * pi * tt / 12) + 4 * cos(2 * pi * tt / 12) noise <- rnorm(n, sd = 3) y <- trend + season + noise y_ts <- ts(y, start = c(2010, 1), frequency = 12) plot(y_ts, main = "Simulated monthly time series", ylab = "y") ``` Fit an additive STL decomposition. ```{r additive-fit} fit_additive <- stl( y_ts, period = 12, s_window = 13 ) lapply(fit_additive, head) ``` The fitted components can be inspected individually. ```{r additive-components, fig.width=8, fig.height=8} plot(fit_additive) ``` The estimated components, they can of course also be plotted manually. ```{r additive-manual-plot, eval = FALSE, fig.width=8, fig.height=8} par(mfrow = c(4, 1), mar = c(3, 4, 2, 1)) plot(y_ts, main = "Observed series", ylab = "") plot(fit_additive$trend, main = "Trend", ylab = "") plot(fit_additive$seasonal, main = "Seasonal component", ylab = "") plot(fit_additive$remainder, main = "Remainder", ylab = "") par(mfrow = c(1, 1)) ``` # Multiplicative decomposition An additive decomposition assumes that the magnitude of the seasonal variation is approximately constant over time. In many economic, business, energy, and demand series, however, larger values are accompanied by larger seasonal fluctuations. For such data a multiplicative decomposition is often more appropriate: \[ Y_t = T_t S_t R_t. \] The following simulated example has a growing level and a seasonal pattern whose absolute amplitude grows with the level. ```{r multiplicative-data, fig.width=8, fig.height=4} n <- 20 * 12 tt <- seq_len(n) trend_mult <- exp(4.5 + 0.005 * tt) season_mult <- exp( 0.20 * sin(2 * pi * tt / 12) + 0.08 * cos(2 * pi * tt / 12) ) remainder_mult <- exp(rnorm(n, sd = 0.06)) y_mult <- trend_mult * season_mult * remainder_mult y_mult_ts <- ts(y_mult, start = c(2000, 1), frequency = 12) plot( y_mult_ts, main = "Simulated multiplicative monthly series", ylab = "y" ) ``` Fit a multiplicative decomposition. ```{r multiplicative-fit} fit_mult <- stl( y_mult_ts, period = 12, s_window = 13, multiplicative = TRUE ) ``` ```{r multiplicative-plot, fig.width=8, fig.height=8} plot(fit_mult) ``` # Weekly data with a non-integer seasonal period A common problem with weekly data is that the annual seasonal cycle is not exactly 52 weeks long. The mean number of weeks in a year is \[ \frac{365.2425}{7} \approx 52.18. \] Using a period of exactly 52 weeks can gradually shift the estimated seasonal pattern relative to the calendar year. `stlpro` supports non-integer seasonal periods and can therefore use `period = 52.18`. Internally, the package constructs cycle subseries at positions such as \[ t,\quad t + 52.18,\quad t + 2(52.18),\ldots \] When one of these positions lies between observed weekly data points, its value is obtained by interpolation between the neighboring observations. ## Simulated weekly series The following example contains 25 years of weekly observations with an annual seasonal cycle of 52.18 weeks. ```{r weekly-data, fig.width=8, fig.height=4} n_years <- 25 period_weekly <- 52.18 n <- round(n_years * period_weekly) tt <- seq_len(n) trend_weekly <- 200 + 0.12 * tt season_weekly <- 25 * sin(2 * pi * tt / period_weekly) + 8 * cos(2 * pi * tt / period_weekly) noise_weekly <- rnorm(n, sd = 3) y_weekly <- trend_weekly + season_weekly + noise_weekly weekly_start <- as.Date("2015-01-05") weekly_dates <- weekly_start + 7 * (tt - 1) plot( weekly_dates, y_weekly, type = "l", xlab = "", ylab = "y", main = "Simulated weekly series with a 52.18-week seasonal cycle" ) ``` Fit the decomposition using a non-integer seasonal period. ```{r weekly-fit} fit_weekly <- stl( y_weekly, period = 52.18, s_window = 7 ) ``` The seasonal window `s_window = 7` means that the seasonal smoother uses seven neighbouring observations of the same period (here: pseudo-week) in the weighted local regression (loess). ```{r weekly-components, eval = TRUE, fig.width=8, fig.height=8} plot(fit_weekly) ``` ## Integer versus non-integer weekly periods For comparison, fit the same series with a fixed period of 52 weeks. ```{r weekly-compare} fit_weekly_52 <- stl( y_weekly, period = 52, s_window = 7 ) fit_weekly_5218 <- stl( y_weekly, period = 52.18, s_window = 7 ) ``` ```{r weekly-compare-plot, eval = TRUE, fig.width=8, fig.height=8} plot( weekly_dates, fit_weekly_52$seasonal, type = "l", col = "orange", lwd = 1.5, xlab = "", ylab = "Seasonal component", main = "Seasonal estimates using 52 and 52.18 week periods", ylim = c(-30, 35) ) lines( weekly_dates, fit_weekly_5218$seasonal, col = "darkblue", lwd = 1.5 ) legend( "topright", legend = c("period = 52", "period = 52.18"), col = c("orange", "darkblue"), lwd = 1.5, bty = "n" ) ``` The difference between `period = 52` and `period = 52.18` is usually modest over short series, but can become relevant for long weekly series because the fixed 52-week cycle slowly drifts relative to an annual calendar cycle. # Selecting the seasonal window by cross-validation The seasonal window `ns` controls the smoothness of the estimated seasonal component. - Smaller values of `ns` allow the seasonal pattern to change more quickly over time. - Larger values of `ns` impose a more stable seasonal pattern. - Very small windows can overfit short-term noise. - Very large windows can miss genuine changes in seasonality. `stlpro` can select or evaluate the seasonal window using leave-one-out cross-validation. In this approach, the observation being predicted is excluded from the local seasonal fit. The following example evaluates a set of candidate seasonal windows for the weekly series. ```{r cv-grid} candidate_ns <- seq(from = 7, to = 51, by = 2) ``` Then, ```{r cv-fit-direct, eval = TRUE} fit_weekly_cv <- stl( y_weekly, period = 52.18, s_window_candidates = candidate_ns ) fit_weekly_cv$cv_results$cv_choice knitr::kable(fit_weekly_cv$cv_results$error_metrics) ``` A cross-validation result typically contains the candidate window sizes and their prediction errors. ```{r cv-plot, eval = TRUE, fig.width=8, fig.height=4} plot( candidate_ns, fit_weekly_cv$cv_results$error_metrics, type = "b", pch = 19, xlab = "Seasonal window (ns)", ylab = "Cross-validation error", main = "Seasonal-window cross-validation" ) ``` # Robust decomposition and outliers Robust STL fitting reduces the influence of unusual observations. This can be particularly useful for retail demand, web traffic, energy data, and other series that may contain isolated spikes or drops. ```{r outlier-data, fig.width=8, fig.height=4} y_outlier <- y_weekly y_outlier[c(80, 180, 260)] <- y_outlier[c(80, 180, 260)] + c(80, -100, 120) plot( weekly_dates, y_outlier, type = "l", xlab = "", ylab = "y", main = "Weekly series with artificial outliers" ) ``` ```{r robust-fit} fit_nonrobust <- stl( y_outlier, period = 52.18, s_window = 7 ) fit_robust <- stl( y_outlier, period = 52.18, s_window = 7, robust = TRUE ) ``` ```{r robust-plot, eval = TRUE, fig.width=8, fig.height=4} plot( weekly_dates, fit_nonrobust$seasonal, type = "l", col = "grey40", xlab = "", ylab = "Seasonal component", main = "Effect of robust fitting" ) lines( weekly_dates, fit_robust$seasonal, col = "steelblue", lwd = 1.5 ) legend( "topright", legend = c("Non-robust", "Robust"), col = c("grey40", "steelblue"), lwd = c(1, 1.5), bty = "n" ) ``` The robust estimate should generally be less affected by isolated extreme observations. # Results comparison with `stats::stl()` When identical settings are used, the results from `stats::stl()` and `stlpro` generally differ by less than `1e-12`. Differences that may cause larger deviations are: - If cross-validation is used, the s_jump parameter is set to 1, while otherwise and the default in `stats::stl()` is that it depends on s_window: s_jump = max(1, as.integer(floor(s_window / 10 + 0.9))) - multiplicative decomposition is not possible in `stats::stl()`, approximating it by using the well-known combination of `log(series)` and `exp(result)` will not yield identical results ```{r comparison-data} n <- 12 * 13 tt <- seq_len(n) y_benchmark <- 100 + 0.15 * tt + 12 * sin(2 * pi * tt / 12) + rnorm(n, sd = 4) y_benchmark_ts <- ts(y_benchmark, frequency = 12) ``` ```{r comparison_results, eval = TRUE} stlpro_result <- stl( y_benchmark_ts, period = 12, s_window = 13, s_jump = 1, l_jump = 1, t_jump = 1, robust = FALSE ) stats_stl_result <- stats::stl( y_benchmark_ts, s.window = 13, s.jump = 1, t.jump = 1, l.jump = 1, robust = FALSE ) res <- data.frame(stlpro = stlpro_result$seasonal, stats_stl = stats_stl_result$time.series[,1], diff = stlpro_result$seasonal - stats_stl_result$time.series[,1]) res max(res[,3]) ``` # Speed comparison with `stats::stl()` `stats::stl()` is highly optimized and is an appropriate benchmark for ordinary integer-period STL decompositions. Since `stats::stl()` is designed for conventional `ts` data, the fairest comparison uses an integer period such as monthly data with period 12. The comparison below measures the elapsed time for repeated decompositions of a moderately long monthly series. ```{r benchmark-system-time, eval = TRUE} n_rep <- 1000 time_stlpro <- system.time({ for (i in seq_len(n_rep)) { stl( y_benchmark_ts, period = 12, s_window = 13, s_jump = 1, l_jump = 1, t_jump = 1, robust = FALSE ) } }) time_stats <- system.time({ for (i in seq_len(n_rep)) { stats::stl( y_benchmark_ts, s.window = 13, s.jump = 1, t.jump = 1, l.jump = 1, robust = FALSE ) } }) rbind( stlpro = time_stlpro, stats_stl = time_stats ) ``` Additionally, we can investigate the time it takes to handle a daily time series. Here, we focus on a single time series, because usually only few time series are adjusted in production, but we turn on cross-validation, i.e. leaving s_window = NULL. Cross-validation is particularly useful for higher-frequency time series, as they tend to be harder to assess using traditional tools to select the filter length, such as graphics contrasting the seasonal-irregular component with the estimated seasonal component(s). ```{r benchmark-daily, eval = TRUE} n <- 12 * 366 tt <- seq_len(n) y_benchmark <- 100 + 0.15 * tt + 20 * sin(2 * pi * tt / 365.2425) + rnorm(n, sd = 4) y_daily <- ts(y_benchmark, frequency = 365.2425) n_rep <- 1 time_stlpro_daily <- system.time({ for (i in seq_len(n_rep)) { stl( y_daily, period = 12 ) } }) time_stlpro_daily ``` ## Interpreting the timing comparison As we see, the stl implementation in the stats package is considerably faster than the implementation in stlpro. This is partly due to the fact, that while implementing stlpro, we focussed on interpretability of the code instead of speed optimization. This is because the additional capabilities of stlpro are so new that we want to make it possible for interested parties to evaluate them easily. With this said, even decomposing 1,000 monthly time series usually takes around a second, of course depending on the length of the series. Using cross-validation on a daily time series to select among 5 lengths of the s_window parameter takes only fractions of a second. # Summary `stlpro` extends STL decomposition to settings where the standard integer seasonal-period assumption is too restrictive. In particular, it can be used to: - decompose additive time series; - fit multiplicative decompositions for positive series; - model weekly annual seasonality using `period = 52.18`; - interpolate cycle-subseries values at fractional seasonal positions; - use robust fitting in the presence of outliers; - select the seasonal smoothing window using cross-validation.