Applying stlpro

library(stlpro)  
#> 
#> Attaching package: 'stlpro'
#> The following object is masked from '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.

plot(fit_additive)

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
)  
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.

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).

plot(fit_weekly)

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)
x
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  8.881784e-16
#> 2    11.1966214  11.1966214 -1.065814e-14
#> 3    12.9520989  12.9520989  5.329071e-15
#> 4    10.9615021  10.9615021 -3.552714e-15
#> 5     6.0933260   6.0933260  2.664535e-15
#> 6     0.7089599   0.7089599 -4.440892e-16
#> 7    -8.3840355  -8.3840355 -5.329071e-15
#> 8   -11.4492858 -11.4492858  1.598721e-14
#> 9   -12.0215842 -12.0215842 -5.329071e-15
#> 10   -9.8670891  -9.8670891  1.243450e-14
#> 11   -5.9876855  -5.9876855 -7.105427e-15
#> 12    0.5111988   0.5111988 -6.661338e-16
#> 13    5.2570666   5.2570666  0.000000e+00
#> 14   11.2193484  11.2193484 -1.421085e-14
#> 15   12.9331006  12.9331006  7.105427e-15
#> 16   11.0418702  11.0418702 -1.776357e-15
#> 17    6.0826157   6.0826157  2.664535e-15
#> 18    0.6480595   0.6480595  3.330669e-16
#> 19   -8.3608638  -8.3608638 -3.552714e-15
#> 20  -11.4448196 -11.4448196  1.776357e-14
#> 21  -12.0018359 -12.0018359 -7.105427e-15
#> 22   -9.9174954  -9.9174954  1.776357e-14
#> 23   -5.9899376  -5.9899376 -7.105427e-15
#> 24    0.4912013   0.4912013 -5.551115e-16
#> 25    5.2638132   5.2638132  0.000000e+00
#> 26   11.2569190  11.2569190 -1.243450e-14
#> 27   12.9083339  12.9083339  5.329071e-15
#> 28   11.1173724  11.1173724  0.000000e+00
#> 29    6.0796601   6.0796601  8.881784e-16
#> 30    0.5694223   0.5694223 -6.661338e-16
#> 31   -8.3320446  -8.3320446 -5.329071e-15
#> 32  -11.4452714 -11.4452714  1.243450e-14
#> 33  -11.9768974 -11.9768974 -8.881784e-15
#> 34   -9.9786683  -9.9786683  1.776357e-14
#> 35   -5.9876301  -5.9876301 -5.329071e-15
#> 36    0.4845221   0.4845221  9.436896e-16
#> 37    5.2830955   5.2830955 -3.552714e-15
#> 38   11.3113183  11.3113183 -1.421085e-14
#> 39   12.8752440  12.8752440  5.329071e-15
#> 40   11.1881002  11.1881002  1.776357e-15
#> 41    6.0611481   6.0611481  0.000000e+00
#> 42    0.4461417   0.4461417  3.330669e-16
#> 43   -8.2928327  -8.2928327 -5.329071e-15
#> 44  -11.4246726 -11.4246726  1.421085e-14
#> 45  -11.9612949 -11.9612949 -7.105427e-15
#> 46  -10.0649865 -10.0649865  1.598721e-14
#> 47   -5.9679792  -5.9679792 -4.440892e-15
#> 48    0.5187698   0.5187698  3.441691e-15
#> 49    5.3386781   5.3386781 -8.881784e-16
#> 50   11.3687248  11.3687248 -1.776357e-14
#> 51   12.8220848  12.8220848  3.552714e-15
#> 52   11.2875189  11.2875189 -1.776357e-15
#> 53    5.9793484   5.9793484  1.776357e-15
#> 54    0.2274499   0.2274499  1.471046e-15
#> 55   -8.2308624  -8.2308624 -1.776357e-15
#> 56  -11.3333725 -11.3333725  1.598721e-14
#> 57  -11.9651857 -11.9651857 -8.881784e-15
#> 58  -10.1916349 -10.1916349  1.776357e-14
#> 59   -5.9346826  -5.9346826 -1.776357e-15
#> 60    0.6139717   0.6139717  5.884182e-15
#> 61    5.4701478   5.4701478 -6.217249e-15
#> 62   11.4069858  11.4069858 -1.776357e-14
#> 63   12.6874321  12.6874321  8.881784e-15
#> 64   11.5366401  11.5366401 -3.552714e-15
#> 65    5.7758663   5.7758663  8.881784e-16
#> 66   -0.1349596  -0.1349596  2.997602e-15
#> 67   -8.1194415  -8.1194415 -3.552714e-15
#> 68  -11.1143747 -11.1143747  1.598721e-14
#> 69  -11.9648441 -11.9648441 -8.881784e-15
#> 70  -10.3958280 -10.3958280  1.421085e-14
#> 71   -5.9331848  -5.9331848  8.881784e-16
#> 72    0.7306297   0.7306297  8.437695e-15
#> 73    5.6677206   5.6677206 -7.993606e-15
#> 74   11.3411282  11.3411282 -2.131628e-14
#> 75   12.3677968  12.3677968  3.552714e-15
#> 76   12.1202582  12.1202582 -3.552714e-15
#> 77    5.5787713   5.5787713 -8.881784e-16
#> 78   -0.4348531  -0.4348531  1.498801e-15
#> 79   -7.9052680  -7.9052680 -3.552714e-15
#> 80  -10.9117899 -10.9117899  8.881784e-15
#> 81  -11.6048659 -11.6048659 -1.953993e-14
#> 82  -10.6314679 -10.6314679  7.105427e-15
#> 83   -6.0788900  -6.0788900 -3.552714e-15
#> 84    0.5766264   0.5766264  5.107026e-15
#> 85    6.2328112   6.2328112  1.243450e-14
#> 86   11.2328290  11.2328290 -3.552714e-15
#> 87   11.5691576  11.5691576  2.486900e-14
#> 88   11.3844438  11.3844438  1.776357e-15
#> 89    5.5095033   5.5095033  0.000000e+00
#> 90   -0.6076414  -0.6076414  1.021405e-14
#> 91   -7.3100461  -7.3100461 -7.993606e-15
#> 92  -10.2417339 -10.2417339  1.243450e-14
#> 93  -11.1607780 -11.1607780 -3.019807e-14
#> 94  -10.7686821 -10.7686821 -1.598721e-14
#> 95   -6.2950157  -6.2950157  1.953993e-14
#> 96    0.7823270   0.7823270 -1.032507e-14
#> 97    6.4487779   6.4487779  2.131628e-14
#> 98   11.2514304  11.2514304  1.065814e-14
#> 99   11.2289064  11.2289064  2.664535e-14
#> 100  10.7491447  10.7491447  8.881784e-15
#> 101   5.4283970   5.4283970 -1.776357e-15
#> 102  -0.6574570  -0.6574570  2.220446e-14
#> 103  -7.0125264  -7.0125264 -2.220446e-14
#> 104  -9.9371669  -9.9371669  1.243450e-14
#> 105 -11.1103513 -11.1103513 -6.750156e-14
#> 106 -10.7644630 -10.7644630 -6.039613e-14
#> 107  -6.3235989  -6.3235989  4.529710e-14
#> 108   0.8518071   0.8518071 -4.585221e-14
#> 109   6.5447989   6.5447989  3.463896e-14
#> 110  11.3056840  11.3056840  2.664535e-14
#> 111  11.0807809  11.0807809  4.440892e-14
#> 112  10.4248766  10.4248766  1.953993e-14
#> 113   5.3931336   5.3931336  1.776357e-15
#> 114  -0.6962780  -0.6962780  3.275158e-14
#> 115  -6.8721742  -6.8721742 -2.398082e-14
#> 116  -9.8304623  -9.8304623  8.881784e-15
#> 117 -11.0828648 -11.0828648 -8.171241e-14
#> 118 -10.7735339 -10.7735339 -8.171241e-14
#> 119  -6.3117723  -6.3117723  6.039613e-14
#> 120   0.9022232   0.9022232 -6.639134e-14
#> 121   6.6196360   6.6196360  4.529710e-14
#> 122  11.3415934  11.3415934  3.375078e-14
#> 123  10.9716161  10.9716161  5.329071e-14
#> 124  10.2386807  10.2386807  2.486900e-14
#> 125   5.3703568   5.3703568  3.552714e-15
#> 126  -0.7440089  -0.7440089  3.808065e-14
#> 127  -6.7846490  -6.7846490 -2.664535e-14
#> 128  -9.7673297  -9.7673297  1.065814e-14
#> 129 -11.0285237 -11.0285237 -8.704149e-14
#> 130 -10.7919983 -10.7919983 -9.592327e-14
#> 131  -6.3086635  -6.3086635  6.750156e-14
#> 132   0.9549577   0.9549577 -7.727152e-14
#> 133   6.6969943   6.6969943  5.062617e-14
#> 134  11.3524405  11.3524405  3.907985e-14
#> 135  10.8590943  10.8590943  5.684342e-14
#> 136  10.0967317  10.0967317  2.664535e-14
#> 137   5.3491710   5.3491710  3.552714e-15
#> 138  -0.7913790  -0.7913790  4.096723e-14
#> 139  -6.7084975  -6.7084975 -2.842171e-14
#> 140  -9.6994688  -9.6994688  7.105427e-15
#> 141 -10.9640934 -10.9640934 -9.059420e-14
#> 142 -10.8112957 -10.8112957 -1.012523e-13
#> 143  -6.3176147  -6.3176147  7.283063e-14
#> 144   1.0120331   1.0120331 -8.282264e-14
#> 145   6.7729826   6.7729826  5.240253e-14
#> 146  11.3512856  11.3512856  4.440892e-14
#> 147  10.7474872  10.7474872  5.861978e-14
#> 148   9.9590510   9.9590510  3.197442e-14
#> 149   5.3245377   5.3245377  5.329071e-15
#> 150  -0.8297047  -0.8297047  4.363176e-14
#> 151  -6.6342305  -6.6342305 -2.842171e-14
#> 152  -9.6195256  -9.6195256  1.065814e-14
#> 153 -10.9077016 -10.9077016 -9.769963e-14
#> 154 -10.8267821 -10.8267821 -1.065814e-13
#> 155  -6.3309961  -6.3309961  7.638334e-14
#> 156   1.0656927   1.0656927 -8.837375e-14

max(res[,3])
#> [1] 7.638334e-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.593    0.017   0.610          0         0
#> stats_stl     0.273    0.000   0.273          0         0

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 elapsed 
#>   0.029   0.000   0.029

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.