In time series analysis and forecasting, it is often useful to identify the underlying patterns of time series data, to better understand the contributing phenomena and to enable better forecasting. Time series decomposition techniques have been introduced to decompose a time series into multiple components including trend, seasonality, and remainder. Such methods have important time series applications in seasonal adjustment procedures (Maravall, 2006; Thornton, 2013), forecasting (Koopman and Ooms, 2006; Hewamalage et al., 2021; Bandara et al., 2021a)2019, 2020).
Nowadays, with rapid growth of availability of sensors and data storage capabilities, there is a significant increase of time series with higher sampling rates (sub-hourly, hourly, daily). Compared with traditional time series data, the higher-frequency time series data may exhibit more complex properties, such as multiple seasonal cycles, non-integer seasonality, etc. These time series are commonly found in the utility demand industry (electricity and water usage), mainly due to the intricate usage patterns of humans. For example, Figure 1 illustrates the aggregated half-hourly energy demand time series in the state of Victoria, Australia that has two seasonal periods, a daily seasonality (period ) and a weekly seasonality (period ). It can be seen that the average energy consumption in the weekdays is relatively higher to those in the weekends. Furthermore, a longer version of this time series may even show a yearly seasonality (period
), with values systematically changing across the seasons such as summer and winter. In this scenario, daily and weekly consumption patterns are useful to estimate the short-term energy requirements; the yearly seasonal usage patterns are beneficial in long-term energy planning. Therefore, the accurate decomposition of time series with multiple seasonal cycles is useful to provide better grounds for decision-making in various circumstances.
There is a large body of literature available on time series decomposition techniques. Widely used traditional time series decomposition methods are Seasonal-Trend decomposition using Loess (STL, Cleveland et al., 1990), X-13-ARIMA-SEATS (Bell and Hillmer, 1984), and X-12-ARIMA (Findley et al., 1998). Although these methods have been heavily used in many real-world applications due to their robustness and efficiency, these techniques can handle only time series with a single seasonality. More recently, methods to decompose time series with multiple seasonal patterns have been introduced (Dokumentov and Hyndman, 2021; Wen et al., 2020). For example, Dokumentov and Hyndman (2021) introduced Seasonal-Trend decomposition by Regression (STR), a regression based, additive decomposition technique, which is also capable of modelling the influence of external factors towards the seasonal patterns in a time series. Wen et al. (2020) recently developed Fast-RobustSTL, a decomposition technique that accounts for multiple seasonal patterns and noise in time series. Several studies also support the use of forecasting models to extract seasonal patterns from time series (Bandara et al., 2021a, 2020, b). The idea is to fit a forecasting model to the time series that is capable of handling time series with multiple seasonal patterns, such as TBATS (Trigonometric Exponential Smoothing State Space model with Box-Cox transformation, ARMA errors, Trend and Seasonal Components) (De Livera et al., 2011), and Prophet (Taylor and Letham, 2021); thereafter, to extract the fitted time series components, i.e., trend, multiple seasonal components, from the fitted forecast model. Nonetheless, the objective function of these forecast models is to minimise the prediction error and therefore these methods have to limit themselves to only past data for decomposition, whereas dedicated decomposition methods can use both past and future information. Therefore, the use of prediction based approaches for time series decomposition can be inaccurate and may give not the best possible decomposition of a time series.
In this paper, we introduce Multiple STL Decomposition (MSTL), a fully automated, additive time series decomposition algorithm to handle time series with multiple seasonal cycles. The proposed MSTL algorithm is an extended version of the STL decomposition algorithm, where the STL procedure is applied iteratively to estimate the multiple seasonal components in a time series. This allows MSTL to control the smoothness of the change of seasonal components for each seasonal cycle extracted from the time series, and seamlessly separate their seasonal variations (e.g., deterministic and stochastic seasonalities). For non-seasonal time series, MSTL determines only the trend and remainder components of the time series.
Specifically, the MSTL algorithm initially determines the number of distinct seasonal patterns available in the time series. Often times, the multiple seasonal patterns are structurally unnested and interlacing together. As a result, during decomposition, the seasonal components relevant to a lower seasonal cycle can be excessively absorbed by a higher seasonal cycle. To minimise such seasonal confounding, as the second step, MSTL arranges the identified seasonal cycles in an ascending order. Then, if the time series is seasonal, MSTL applies the STL algorithm iteratively to each of the identified seasonal frequencies. Next, the trend component of the time series is computed using the last iteration of STL. On the other hand, if the time series is non-seasonal, MSTL uses the Friedman’s Super Smoother function, supsmu, available in R (R Core Team, 2020), to directly estimate the trend of the time series. Finally, to calculate the remainder part of seasonal time series, the trend component is subtracted from the seasonally adjusted time series. Whereas, for non-seasonal time series, the trend component is subtracted from the original time series to derive the remainder.
MSTL is a robust, accurate seasonal-trend decomposition algorithm that is designed to capture multiple seasonal patterns in a time series. Most importantly, compared with other decomposition alternatives, MSTL is an extremely fast, computationally efficient algorithm, which is scalable to increasing volumes of time series data. In R, the proposed MSTL algorithm is implemented in the
mstl function from the forecast package (Hyndman and Khandakar, 2008; Hyndman et al., 2021).
2 Model Overview
Similar to the STL algorithm, MSTL gives an additive decomposition of the time series. Given is the observation at time , the additive decomposition can be defined as follows:
Here, , , denote the seasonal, trend, and remainder components of the observation, respectively. MSTL extends Equation 1 to include multiple seasonal patterns in a time series as follows:
Here, represents the number of seasonal cycles present in .
To summarize, a scheme of the MSTL procedure is given in Algorithm 1
. At first, the frequencies of the seasonal patterns in the time series are identified and sorted in an ascending order. Here, the frequencies which are smaller than half of the length of the series are ignored, as those frequencies cannot exhibit any seasonal patterns. After identifying the seasonality, the missing values of the time series are imputed using the
na.interpfunction available from the forecast package (Hyndman et al., 2021). Next, if is given, a Box-Cox transformation is applied to the time series accordingly, using the
BoxCoxfunction available from the forecast package (Hyndman et al., 2021). Thereafter, to every selected seasonal cycle the STL decomposition is fitted. Here, the inner-loop repeats the STL procedure to extract the seasonal components from the time series. The STL model is implemented using the
stlfunction provided by the
statspackage in R (R Core Team, 2020). In STL, the rate of seasonal variation is controlled by the
s.windowparameter. Here, a smaller value of
is set if the seasonal pattern evolves quickly, whereas a higher value is used if the seasonal pattern is constant over time. For example, adjusting the s.window parameter to “periodic” limits the change in the seasonal components to zero, which can extract the deterministic seasonality from a time series. In MSTL, a vector of
s.windowcan be provided to control the variation of the seasonal components of each seasonal cycle. The outer-loop iterates the STL procedure multiple times to refine the extracted seasonal components. After executing the outer-loop, MSTL calculates the trend component of the time series using the final iteration of STL. In situations where a time series fails the seasonality test, MSTL applies the supsmu function to ascertain the trend of the time series. Finally, the remainder component is retrieved by deducting the trend component from the seasonally-adjusted time series. Also, since MSTL extends the STL algorithm, the other parameters of STL (e.g.,
l.window) are inherited by MSTL, and can be used for decomposition.
3 Experimental Setup
We evaluate the proposed MSTL decomposition algorithm on both simulated and a perturbed real-world time series dataset, where we know the true composition, i.e., trend, seasonality, and remainder, of time series.
We compare MSTL against a collection of current state-of-the-art techniques in decomposing time series with multiple seasonal cycles. This includes STR (Dokumentov and Hyndman, 2021) as a pure decomposition technique, TBATS (De Livera et al., 2011) and Prophet (Taylor and Letham, 2021) as forecasting techniques, where we use the decomposition of the time series.
STR: A regression based, additive decomposition technique. In R, the STR algorithm is available through the
STRfunction from the stR package (Dokumentov and Hyndman, 2018).
TBATS: A state-of-the-art technique to forecast time series with multiple seasonal cycles. TBATS uses trigonometric expression terms to model complex seasonal terms in a time series. In our experiments, we use the R implementation of the TBATS algorithm,
tbats, from the forecast package (Hyndman et al., 2021).
PROPHET: An automated forecasting framework, developed by Facebook, that can handle multiple seasonal patterns. Similar to STR and MSTL, Prophet is an additive decomposition technique. In our evaluation, we apply the Prophet algorithm available through the
Prophetpackage in R (Taylor and Letham, 2021).
3.2 Evaluation metric
The accuracy of the decomposition methods are evaluated using the root mean square error (RMSE). The RMSE is defined as follows:
Here, represents the actual decomposition value at time , and is the estimated decomposition value. Also, is the number of observations in the time series.
3.3 Simulated Data
We simulate daily and hourly data using four time series components. The additive decomposition of the daily time series () and the hourly time series () can be formulated as follows:
In Equation (4), is the trend, is the weekly seasonal component, is the yearly seasonal component, and is the remainder of . Whereas in Equation (5), , , , and corresponds to the trend, daily seasonal component, weekly seasonal component, and remainder component of respectively. Here, , , and are parameters which control the contribution of the components to and . Also, denotes the length of the daily time series and is the length of the hourly time series.
We simulate the time series data using two data generating processes (DGPs), the “Deterministic DGP” and the “Stochastic DGP””. Here, the Deterministic DGP generates time series that has deterministic components that are invariant to time, whereas the Stochastic DGP gives time series with time-varying components.
In our experiments, the deterministic trend components are generated using a quadratic trend function with random coefficients, where and
are independent N(0,1) random variables. The deterministic seasonal componentsare composed of five pairs of Fourier terms with random N(0,1) coefficients. Both and
are normalised to give mean zero and unit variance.
To generate stochastic trend components , we use an ARIMA(0,2,0) model with standard normal errors. With respect to the stochastic seasonal components , we introduce an additional error term N(0,) to to change the coefficients for the Fourier terms from one seasonal cycle to another. Here, the stochastic strength of the seasonality is controlled by the parameter value. In other words, the scenario of represents the .
Figure 2 illustrates the examples of simulated daily and hourly time series generated by the two DGPs. Here, we set the and for both the DGPs. The lengths of the daily and hourly time series are equivalent to 1096 days and 505 days respectively (one observation more than three seasonal cycles of the highest available seasonality, i.e., 365 * 3 + 1 for the daily data and 168 * 3 + 1 for the hourly data).
Table 1 summarises the parameter sets (, , , ) used in our experiments. We create 150 time series for each DGP, generating 300 time series of daily and hourly data. As we know the true values of the seasonal, trend, and remainder components of the simulated time series, we calculate the root mean square error (RMSE) for each component by averaging across all datasets and dates. The presence of statistical significance of differences within multiple decomposition methods is assessed using a linear model applied to the squared errors.
Moreover, for the deterministic DGPs we set the
s.window values of MSTL to , whereas for the stochastic DGPs, we use the default
s.window values of the MSTL, which were obtained through a simulation study (refer to A).
|Method||Trend RMSE||Weekly RMSE||Yearly RMSE||Remainder RMSE|
Table 2 shows the evaluation summary for the daily simulated time series. With respect to deterministic DGPs, it can be seen that on the Trend, Weekly and Yearly seasonal components, MSTL outperforms the TBATS method in the majority of cases. Whereas, MSTL achieves the best Weekly RMSE values for stochastic DGPs, outperforming all the benchmarks. Also, for stochastic DGPs, MSTL gives better results compared to TBATS on the Trend, Weekly and Yearly seasonal components.
|Method||Trend RMSE||Daily RMSE||Weekly RMSE||Remainder RMSE|
Table 3 shows the results of all the decomposition techniques for the hourly simulated time series. For the deterministic and stochastic DGPs, we see that the proposed MSTL algorithm gives better RMSE on all components, compared to the PROPHET method. Also, for the scenario on the Trend component, we that the MSTL achieves the best RMSE.
3.4 Perturbed real-world data
We also evaluate the performance of MSTL using real-world time series data. We use the moving block bootstrap (MBB) for real-world time series perturbing, following the procedure introduced in Bergmeir et al. (2016). Here, we first use MSTL to extract the seasonal, trend and remainder of components of a time series. Assuming these extracted components are the true decomposition of the time series, next we apply the MBB technique to the remainder component of the time series to generate multiple versions of the residual components. Finally, these bootstrapped residual components are added back together with the previously extracted seasonal and trend components to produce new perturbed versions of a time series, where we know the true composition of the series. In our experiments, we use the MBB implementation available through the
MBB function from the forecast package Hyndman and Khandakar (2008); Hyndman et al. (2021).
We select the half-hourly electricity consumption in the state of Victoria, Australia, extracted from the vic_elec dataset (O’Hara-Wild et al., 2021) to generate multiple versions of the time series. As illustrated in Figure 1, this time series has two seasonal patterns, the daily pattern and the weekly seasonal pattern. In our experiments, we first aggregate the half-hourly data to hourly data, and select 149 days starting from 01 January 2012. Figure 3 shows the application of MSTL to the hourly electricity demand in Victoria (3601 hourly observations).
As discussed earlier, we use the MBB technique to generate 100 bootstrapped versions of the hourly electricity demand time series, so we can assess the performance of MSTL on a real-world dataset. Figure 4 illustrates the snippet of original hourly electricity demand time series and the bootstrapped time series.
|Method||Trend RMSE||Daily RMSE||Weekly RMSE||Remainder RMSE|
Table 4 summarises the overall performance of MSTL and the benchmarks on the perturbed electricity demand time series. According to Table 4, MSTL significantly outperforms STR, TBATS, and PROPHET in estimating all components. Furthermore, Table 5 provides a summary of the computational cost of MSTL and the benchmarks over the 100 bootstrapped time series. The experiments are run on an Intel(R) i7 processor (1.8 GHz), with 2 threads per core, 8 cores in total, and 16GB of main memory. As shown in Table 5, MSTL has the lowest execution time compared to other benchmarks, highlighting the scalability of MSTL to the increasing volumes of time series data. When used for sub-daily time series data, i.e, hourly, half-hourly, the computational efficiency of a decomposition method can be important as those time series are generally longer and contain a higher number of observations.
The total computational cost of the decomposition methods for the electricity demand dataset, measured in seconds.
Time series datasets with multiple seasonalities are nowadays common in real-world applications. To better understand the variations of such time series, it is often important to decompose the time series into their subcomponents, such as trend, seasonality, and remainder. The existing techniques available for decomposing time series with multiple seasonal cycles are mostly based on complex procedures that can be computationally inefficient for long time series.
To this end, we have introduced MSTL, a fast time series decomposition algorithm that is capable of handling time series with multiple seasonal cycles. The proposed MSTL algorithm is an extension of the STL decomposition algorithm, which can only extract a single seasonality from a time series. Experimental results on both simulated data and perturbed real-world data have demonstrated that MSTL provides competitive results with lower computational cost in comparison with other state-of-the-art decomposition algorithms, such as STR, TBATS, and PROPHET.
The MSTL algorithm is implemented in the
mstl function in the forecast package (Hyndman et al., 2021), which is available on the Comprehensive R Archive Network (CRAN).
This research was supported by the Australian Research Council under grant DE190100045, and the Australian Centre of Excellence for Mathematical and Statistical Frontiers (Grant CE140100049).
Bandara et al. (2020)
Bandara, K., Bergmeir, C.,
Campbell, S., Scott, D.,
Lubman, D., 2020.
Towards accurate predictions and causal ‘what-if’ analyses for planning and policy-making: A case study in emergency medical services demand, in: 2020 International Joint Conference on Neural Networks (IJCNN), pp. 1–10.
- Bandara et al. (2021a) Bandara, K., Bergmeir, C., Hewamalage, H., 2021a. LSTM-MSNet: Leveraging forecasts on sets of related time series with multiple seasonal patterns. IEEE Transactions on Neural Networks and Learning Systems 32, 1586–1599.
- Bandara et al. (2021b) Bandara, K., Hewamalage, H., Liu, Y.H., Kang, Y., Bergmeir, C., 2021b. Improving the accuracy of global forecasting models using time series data augmentation. Pattern Recognittion , 108148.
- Bell and Hillmer (1984) Bell, W.R., Hillmer, S.C., 1984. Issues involved with the seasonal adjustment of economic time series. J. Bus. Econ. Stat. 2, 291–320.
- Bergmeir et al. (2016) Bergmeir, C., Hyndman, R.J., Benítez, J.M., 2016. Bagging exponential smoothing methods using STL decomposition and Box–Cox transformation. Int. J. Forecast. 32, 303–312.
- Cleveland et al. (1990) Cleveland, R.B., Cleveland, W.S., McRae, J.E., Terpenning, I.J., 1990. STL : A seasonal-trend decomposition procedure based on loess. Journal of Official Statistics 6, 3–73.
- De Livera et al. (2011) De Livera, A.M., Hyndman, R.J., Snyder, R., 2011. Forecasting time series with complex seasonal patterns using exponential smoothing. J. Am. Stat. Assoc. 106, 1513–1527.
- Dokumentov and Hyndman (2018) Dokumentov, A., Hyndman, R.J., 2018. stR: STR Decomposition. URL: https://CRAN.R-project.org/package=stR. R package version 0.4.
- Dokumentov and Hyndman (2021) Dokumentov, A., Hyndman, R.J., 2021. STR: A seasonal-trend decomposition procedure based on regression. INFORMS J on Data Science (in press).
- Findley et al. (1998) Findley, D.F., Monsell, B.C., Bell, W.R., Otto, M.C., Chen, B.C., 1998. New capabilities and methods of the X-12-ARIMA Seasonal-Adjustment program. J. Bus. Econ. Stat. 16, 127–152.
- Hewamalage et al. (2021) Hewamalage, H., Bergmeir, C., Bandara, K., 2021. Recurrent neural networks for time series forecasting: Current status and future directions. International Journal of Forecasting 37, 388–427.
- Hyndman et al. (2021) Hyndman, R.J., Athanasopoulos, G., Bergmeir, C., Caceres, G., Chhay, L., O’Hara-Wild, M., Petropoulos, F., Razbash, S., Wang, E., Yasmeen, F., R Core Team, Ihaka, R., Reid, D., Shaub, D., Tang, Y., Zhou, Z., 2021. forecast: Forecasting Functions for Time Series and Linear Models. URL: https://CRAN.R-project.org/package=forecast. R package version 8.15.
- Hyndman and Khandakar (2008) Hyndman, R.J., Khandakar, Y., 2008. Automatic time series forecasting: the forecast package for R. Journal of Statistical Software 26, 1–22.
- Koopman and Ooms (2006) Koopman, S.J., Ooms, M., 2006. Forecasting daily time series using periodic unobserved components time series models. Computational Statistics & Data Analysis 51, 885–903.
- Maravall (2006) Maravall, A., 2006. An application of the TRAMO-SEATS automatic procedure; direct versus indirect adjustment. Computational Statistics & Data Analysis 50, 2167–2190.
- O’Hara-Wild et al. (2021) O’Hara-Wild, M., Hyndman, R.J., Wang, E., 2021. tsibbledata: Diverse Datasets for ’tsibble’. URL: https://CRAN.R-project.org/package=tsibbledata. R package version 0.3.0.
- R Core Team (2020) R Core Team, 2020. R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing. Vienna, Austria. URL: https://www.R-project.org/.
- Taylor and Letham (2021) Taylor, S., Letham, B., 2021. prophet: Automatic Forecasting Procedure. URL: https://CRAN.R-project.org/package=prophet. R package version 1.0.
- Thornton (2013) Thornton, M.A., 2013. Removing seasonality under a changing regime: Filtering new car sales. Computational Statistics & Data Analysis 58, 4–14.
Wen et al. (2019)
Wen, Q., Gao, J., Song,
X., Sun, L., Xu, H.,
Zhu, S., 2019.
RobustSTL: A robust Seasonal-Trend decomposition
algorithm for long time series.
Proceedings of the AAAI Conference on Artificial Intelligence 33, 5409–5416.
- Wen et al. (2020) Wen, Q., Zhang, Z., Li, Y., Sun, L., 2020. Fast RobustSTL: Efficient and robust Seasonal-Trend decomposition for time series with complex patterns, in: Proceedings of the 26th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining, Association for Computing Machinery, New York, NY, USA. pp. 2203–2213.
Appendix A Results of seasonal window simulations
To determine the default parameter values for the
s.window parameter values of MSTL, we conduct a series of experiments using a simulation setup similar to Section 3.3. Table 1 summarises the parameter sets (, , , ) used in our experiments. We only use stochastic DGPs to generate time series for this simulation study as real-world time series often contain stochastic time series components.
As we use daily and hourly simulated data, we define the first
s.window value for the lowest seasonal cycle as S1.Window and the second
s.window value for the next highest seasonal cycle as S2.Window. For example, S1.Window and S2.window are the
s.window values for the weekly and yearly seasonal cycles of daily time series and the
s.window values for the daily and weekly seasonal cycles of hourly time series. Figure 1 shows the RMSE error distribution boxplots for the different types of S1.Window and S2.Window combinations on the weekly and hourly datasets. Here, the S1.Window and S2.Window pairs are chosen from the vector .
To identify default parameter values for the
s.window parameter, we extend this experiments to include more combinations of
s.window values that satisfy the S1.Window S2.Window premise. Figure 2 shows RMSE error distribution boxplots for the different combinations of S1.Window and S2.Window values on the weekly and hourly datasets. The S1.Window and S2.Window pairs are generated from , where and . Here represents the seasonal cycle number, i.e., for the first and second seasonality and
respectively. Also, we select the smallest odd value from, when choosing a
Overall, it can be seen that the S1.Window value of 11 and the S2.Window value of 15 gives the best median RMSE for weekly and hourly time series, which is the and instance of the above formula. Based on these results, we set the default
s.window parameter of MSTL using the formula, where and values are set to 7 and 4 respectively.
According to Figure 1, the
s.window combinations that meet the S1.Window S2.Window condition give the best median RMSE.