Applying stlpro
library(stlpro)
#>
#> Attache Paket: 'stlpro'
#> Das folgende Objekt ist maskiert 'package:stats':
#>
#> stl
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.
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.
fit_additive <- stl(
y_ts,
period = 12,
s_window = 13
)
lapply(fit_additive, head)
#> $seasonal
#> Jan Feb Mar Apr May Jun
#> 2010 8.905377 12.864770 11.068121 8.059292 5.147047 -4.441502
#>
#> $trend
#> Jan Feb Mar Apr May Jun
#> 2010 99.87488 100.15369 100.43251 100.71132 101.01336 101.31539
#>
#> $remainder
#> Jan Feb Mar Apr May Jun
#> 2010 -1.2754633 -1.5520801 4.2781431 -0.7908847 3.8848175 0.4738356
#>
#> $sa
#> Jan Feb Mar Apr May Jun
#> 2010 98.59942 98.60161 104.71065 99.92044 104.89818 101.78923
#>
#> $weights
#> Jan Feb Mar Apr May Jun
#> 2010 1 1 1 1 1 1
#>
#> $call
#> $call$period
#> [1] 12
#>
#> $call$multiplicative
#> [1] FALSE
#>
#> $call$s_window
#> [1] 13
#>
#> $call$s_window_used
#> [1] 13
#>
#> $call$t_window
#> [1] 21
#>
#> $call$l_window
#> [1] 13
#>
#>
#> $cv_results
#> NULL
The fitted components can be inspected individually.
The estimated components, they can of course also be plotted
manually.
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.
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.
fit_mult <- stl(
y_mult_ts,
period = 12,
s_window = 13,
multiplicative = TRUE
)
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.
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.
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).
Integer versus non-integer weekly periods
For comparison, fit the same series with a fixed period of 52
weeks.
fit_weekly_52 <- stl(
y_weekly,
period = 52,
s_window = 7
)
fit_weekly_5218 <- stl(
y_weekly,
period = 52.18,
s_window = 7
)
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.
candidate_ns <- seq(from = 7, to = 51, by = 2)
Then,
fit_weekly_cv <- stl(
y_weekly,
period = 52.18,
s_window_candidates = candidate_ns
)
fit_weekly_cv$cv_results$cv_choice
#> [1] 51
knitr::kable(fit_weekly_cv$cv_results$error_metrics)
| 7 |
3.373198 |
| 9 |
3.264886 |
| 11 |
3.206524 |
| 13 |
3.172830 |
| 15 |
3.151831 |
| 17 |
3.137699 |
| 19 |
3.129897 |
| 21 |
3.125588 |
| 23 |
3.124037 |
| 25 |
3.123155 |
| 27 |
3.122733 |
| 29 |
3.122247 |
| 31 |
3.121734 |
| 33 |
3.121246 |
| 35 |
3.120815 |
| 37 |
3.120452 |
| 39 |
3.120157 |
| 41 |
3.119923 |
| 43 |
3.119740 |
| 45 |
3.119600 |
| 47 |
3.119495 |
| 49 |
3.119416 |
| 51 |
3.119359 |
A cross-validation result typically contains the candidate window
sizes and their prediction errors.
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.
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"
)
fit_nonrobust <- stl(
y_outlier,
period = 52.18,
s_window = 7
)
fit_robust <- stl(
y_outlier,
period = 52.18,
s_window = 7,
robust = TRUE
)
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
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)
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
#> stlpro stats_stl diff
#> 1 5.2500188 5.2500188 -1.598721e-14
#> 2 11.1966214 11.1966214 -1.065814e-14
#> 3 12.9520989 12.9520989 -3.552714e-15
#> 4 10.9615021 10.9615021 1.776357e-15
#> 5 6.0933260 6.0933260 6.217249e-15
#> 6 0.7089599 0.7089599 4.773959e-15
#> 7 -8.3840355 -8.3840355 -2.309264e-14
#> 8 -11.4492858 -11.4492858 8.881784e-15
#> 9 -12.0215842 -12.0215842 -1.776357e-15
#> 10 -9.8670891 -9.8670891 1.421085e-14
#> 11 -5.9876855 -5.9876855 0.000000e+00
#> 12 0.5111988 0.5111988 1.132427e-14
#> 13 5.2570666 5.2570666 -1.865175e-14
#> 14 11.2193484 11.2193484 -1.065814e-14
#> 15 12.9331006 12.9331006 0.000000e+00
#> 16 11.0418702 11.0418702 5.329071e-15
#> 17 6.0826157 6.0826157 6.217249e-15
#> 18 0.6480595 0.6480595 5.440093e-15
#> 19 -8.3608638 -8.3608638 -2.131628e-14
#> 20 -11.4448196 -11.4448196 1.065814e-14
#> 21 -12.0018359 -12.0018359 -3.552714e-15
#> 22 -9.9174954 -9.9174954 1.776357e-14
#> 23 -5.9899376 -5.9899376 8.881784e-16
#> 24 0.4912013 0.4912013 1.088019e-14
#> 25 5.2638132 5.2638132 -1.687539e-14
#> 26 11.2569190 11.2569190 -1.421085e-14
#> 27 12.9083339 12.9083339 3.552714e-15
#> 28 11.1173724 11.1173724 0.000000e+00
#> 29 6.0796601 6.0796601 3.552714e-15
#> 30 0.5694223 0.5694223 5.551115e-15
#> 31 -8.3320446 -8.3320446 -2.131628e-14
#> 32 -11.4452714 -11.4452714 3.552714e-15
#> 33 -11.9768974 -11.9768974 -3.552714e-15
#> 34 -9.9786683 -9.9786683 1.776357e-14
#> 35 -5.9876301 -5.9876301 1.776357e-15
#> 36 0.4845221 0.4845221 1.099121e-14
#> 37 5.2830955 5.2830955 -1.776357e-14
#> 38 11.3113183 11.3113183 -1.065814e-14
#> 39 12.8752440 12.8752440 1.776357e-15
#> 40 11.1881002 11.1881002 1.776357e-15
#> 41 6.0611481 6.0611481 2.664535e-15
#> 42 0.4461417 0.4461417 6.439294e-15
#> 43 -8.2928327 -8.2928327 -2.131628e-14
#> 44 -11.4246726 -11.4246726 1.776357e-15
#> 45 -11.9612949 -11.9612949 -5.329071e-15
#> 46 -10.0649865 -10.0649865 1.776357e-14
#> 47 -5.9679792 -5.9679792 8.881784e-16
#> 48 0.5187698 0.5187698 1.054712e-14
#> 49 5.3386781 5.3386781 -1.598721e-14
#> 50 11.3687248 11.3687248 -1.065814e-14
#> 51 12.8220848 12.8220848 0.000000e+00
#> 52 11.2875189 11.2875189 3.552714e-15
#> 53 5.9793484 5.9793484 3.552714e-15
#> 54 0.2274499 0.2274499 7.382983e-15
#> 55 -8.2308624 -8.2308624 -2.131628e-14
#> 56 -11.3333725 -11.3333725 7.105427e-15
#> 57 -11.9651857 -11.9651857 -7.105427e-15
#> 58 -10.1916349 -10.1916349 1.598721e-14
#> 59 -5.9346826 -5.9346826 0.000000e+00
#> 60 0.6139717 0.6139717 8.548717e-15
#> 61 5.4701478 5.4701478 -1.509903e-14
#> 62 11.4069858 11.4069858 -8.881784e-15
#> 63 12.6874321 12.6874321 8.881784e-15
#> 64 11.5366401 11.5366401 1.776357e-15
#> 65 5.7758663 5.7758663 3.552714e-15
#> 66 -0.1349596 -0.1349596 1.013079e-14
#> 67 -8.1194415 -8.1194415 -1.953993e-14
#> 68 -11.1143747 -11.1143747 5.329071e-15
#> 69 -11.9648441 -11.9648441 -1.065814e-14
#> 70 -10.3958280 -10.3958280 7.105427e-15
#> 71 -5.9331848 -5.9331848 8.881784e-16
#> 72 0.7306297 0.7306297 6.883383e-15
#> 73 5.6677206 5.6677206 -1.421085e-14
#> 74 11.3411282 11.3411282 -3.552714e-15
#> 75 12.3677968 12.3677968 7.105427e-15
#> 76 12.1202582 12.1202582 0.000000e+00
#> 77 5.5787713 5.5787713 2.664535e-15
#> 78 -0.4348531 -0.4348531 1.187939e-14
#> 79 -7.9052680 -7.9052680 -2.131628e-14
#> 80 -10.9117899 -10.9117899 0.000000e+00
#> 81 -11.6048659 -11.6048659 -1.598721e-14
#> 82 -10.6314679 -10.6314679 7.105427e-15
#> 83 -6.0788900 -6.0788900 8.881784e-16
#> 84 0.5766264 0.5766264 2.442491e-15
#> 85 6.2328112 6.2328112 3.552714e-15
#> 86 11.2328290 11.2328290 7.105427e-15
#> 87 11.5691576 11.5691576 2.131628e-14
#> 88 11.3844438 11.3844438 5.329071e-15
#> 89 5.5095033 5.5095033 1.776357e-15
#> 90 -0.6076414 -0.6076414 2.153833e-14
#> 91 -7.3100461 -7.3100461 -2.575717e-14
#> 92 -10.2417339 -10.2417339 3.552714e-15
#> 93 -11.1607780 -11.1607780 -3.019807e-14
#> 94 -10.7686821 -10.7686821 -1.776357e-14
#> 95 -6.2950157 -6.2950157 2.664535e-14
#> 96 0.7823270 0.7823270 -1.576517e-14
#> 97 6.4487779 6.4487779 1.154632e-14
#> 98 11.2514304 11.2514304 1.953993e-14
#> 99 11.2289064 11.2289064 2.486900e-14
#> 100 10.7491447 10.7491447 1.421085e-14
#> 101 5.4283970 5.4283970 -8.881784e-16
#> 102 -0.6574570 -0.6574570 3.252953e-14
#> 103 -7.0125264 -7.0125264 -3.641532e-14
#> 104 -9.9371669 -9.9371669 5.329071e-15
#> 105 -11.1103513 -11.1103513 -6.039613e-14
#> 106 -10.7644630 -10.7644630 -6.039613e-14
#> 107 -6.3235989 -6.3235989 5.329071e-14
#> 108 0.8518071 0.8518071 -5.406786e-14
#> 109 6.5447989 6.5447989 2.575717e-14
#> 110 11.3056840 11.3056840 3.552714e-14
#> 111 11.0807809 11.0807809 3.730349e-14
#> 112 10.4248766 10.4248766 2.842171e-14
#> 113 5.3931336 5.3931336 3.552714e-15
#> 114 -0.6962780 -0.6962780 4.318768e-14
#> 115 -6.8721742 -6.8721742 -3.996803e-14
#> 116 -9.8304623 -9.8304623 1.776357e-15
#> 117 -11.0828648 -11.0828648 -7.815970e-14
#> 118 -10.7735339 -10.7735339 -8.526513e-14
#> 119 -6.3117723 -6.3117723 6.927792e-14
#> 120 0.9022232 0.9022232 -7.482903e-14
#> 121 6.6196360 6.6196360 3.730349e-14
#> 122 11.3415934 11.3415934 3.907985e-14
#> 123 10.9716161 10.9716161 4.618528e-14
#> 124 10.2386807 10.2386807 3.197442e-14
#> 125 5.3703568 5.3703568 4.440892e-15
#> 126 -0.7440089 -0.7440089 4.951595e-14
#> 127 -6.7846490 -6.7846490 -3.996803e-14
#> 128 -9.7673297 -9.7673297 3.552714e-15
#> 129 -11.0285237 -11.0285237 -8.704149e-14
#> 130 -10.7919983 -10.7919983 -9.592327e-14
#> 131 -6.3086635 -6.3086635 7.815970e-14
#> 132 0.9549577 0.9549577 -8.626433e-14
#> 133 6.6969943 6.6969943 4.263256e-14
#> 134 11.3524405 11.3524405 4.618528e-14
#> 135 10.8590943 10.8590943 4.796163e-14
#> 136 10.0967317 10.0967317 3.375078e-14
#> 137 5.3491710 5.3491710 4.440892e-15
#> 138 -0.7913790 -0.7913790 5.195844e-14
#> 139 -6.7084975 -6.7084975 -4.263256e-14
#> 140 -9.6994688 -9.6994688 1.776357e-15
#> 141 -10.9640934 -10.9640934 -8.881784e-14
#> 142 -10.8112957 -10.8112957 -1.048051e-13
#> 143 -6.3176147 -6.3176147 8.437695e-14
#> 144 1.0120331 1.0120331 -9.303669e-14
#> 145 6.7729826 6.7729826 4.440892e-14
#> 146 11.3512856 11.3512856 4.973799e-14
#> 147 10.7474872 10.7474872 5.151435e-14
#> 148 9.9590510 9.9590510 4.085621e-14
#> 149 5.3245377 5.3245377 6.217249e-15
#> 150 -0.8297047 -0.8297047 5.495604e-14
#> 151 -6.6342305 -6.6342305 -4.174439e-14
#> 152 -9.6195256 -9.6195256 5.329071e-15
#> 153 -10.9077016 -10.9077016 -9.592327e-14
#> 154 -10.8267821 -10.8267821 -1.101341e-13
#> 155 -6.3309961 -6.3309961 8.704149e-14
#> 156 1.0656927 1.0656927 -9.836576e-14
max(res[,3])
#> [1] 8.704149e-14
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.
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
)
#> user.self sys.self elapsed user.child sys.child
#> stlpro 0.83 0.01 0.86 NA NA
#> stats_stl 0.27 0.02 0.28 NA NA
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).
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
#> User System verstrichen
#> 0.04 0.02 0.05
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.