title: Forecasting with Multiple Seasonality
authors: Tianyang Xie Jie Ding
text: 
  
  
  
  
  
  Forecasting with Multiple Seasonality
      Tianyang Xie
  School of Statistics 
  
  University of Minnesota Twin Cities
  
  Minneapolis, US 
  
  tianyang.xie.xty@gmail.com
  Jie Ding
  School of Statistics 
  
  University of Minnesota Twin Cities
  
  Minneapolis, US 
  
  dingj@umn.edu
      December 30, 2023
  ==============================================================================================================================================================================================================================
  
  
  
  
  An emerging number of modern applications involve forecasting time series data that exhibit both short-time dynamics and long-time seasonality. Specifically, time series with multiple seasonality is a difficult task with comparatively fewer discussions. In this paper, we propose a two-stage method for time series with multiple seasonality, which does not require pre-determined seasonality periods. In the first stage, we generalize the classical seasonal autoregressive moving average (ARMA) model in multiple seasonality regime. In the second stage, we utilize an appropriate criterion for lag order selection. Simulation and empirical studies show the excellent predictive performance of our method, especially compared to a recently popular `Facebook Prophet' model for time series.
  
  
  
  Time series, Model selection, Multiple seasonality, Nowcasting. 
  
  
  
  
  § INTRODUCTION
  
  In time series, seasonality is defined as the presence of variations that occur at specific regular intervals. Forecasting on time series with multiple seasonality that has different lengths of seasonality cycles is usually considered a difficult task. Therefore, detection and accommodation of the seasonality effect play an important role in time series forecasting.
  
  Among all the forecasting techniques, the seasonal ARIMA model <cit.> and exponential smoothing technique  <cit.> are the classical approaches. The basic forms of these methods are only suited in modeling single seasonality, and unable to account for multiple seasonal patterns. Many studies have been conducted to extend the classical statistical forecasting models to accommodate multiple seasonal patterns. Representative methods include the double seasonal ARIMA model and exponential smoothing technique adapted from the simple Holt-Winters method <cit.>, and the hidden Markov model with multiple seasonality <cit.>. However, existing techniques tend to be sensitive to over-parameterization and optimization issues, and are also unable to model complex seasonal patterns in a time series, as pointed out in <cit.>. A more flexible method from the state space perspective was developed in <cit.>, which is the state-of-the-art method for multiple seasonality. 
  
  Recently, a new model named Facebook Prophet <cit.> is introduced that utilizes the generalized additive model <cit.>, and succeeds in forecasting Facebook's business data. In parallel to this development, neural networks (NNs) have been advocated as a strong alternative for forecasting time series with multiple seasonality. The recurrent neural network (RNN) <cit.> based deep learning techniques such as long short-term memory (LSTM) <cit.>, gated recurrent unit (GRU) <cit.> achieve promising results in many real data applications <cit.>. 
  In addition to the development in neural networks, hybridization of classical methods and deep learning has also been advocated. For example, in the winning submission of the M4 forecasting competition <cit.>, a hybrid approach called exponential smoothing-recurrent neural network (ES-RNN) succeeds in forecasting hourly time series data. 
  
  Most of the classical parametric methods require the pre-specification of seasonality periods. For example, we usually fix the seasonality period to be 24 for hourly data, then we train model or employ seasonal adjustment on the belief of that. However, in multiple seasonality situation, the pre-specification can be risky when the data doesn't exhibit clear seasonal patterns. An example concerning the total sunspot number is included in Section <ref>. Unfortunately, for methods trained in state-space innovative algorithms, model selection for seasonality periods is computationally prohibited. 
  On the other hand, though deep learning techniques are free from pre-specification, they require a sophisticated tuning process and large datasets. 
  Moreover, it is not clear how to extract the information of seasonality from existing deep learning architectures. 
  
  This paper introduces a simple yet powerful method named Multiple Seasonality (MS) modeling procedure. The proposed method is to forecast centered or detrended time series with multiple seasonality by two stages. In the first stage, it detects seasonal components and fits generalized seasonal ARMA models. In the second stage, it selects an appropriate model and provides forecasting results. Instead of requiring the pre-specification of seasonality periods, it only needs to know how many seasonality components to include in the model. The procedure has a reasonable number of hyper-parameters, enabling easy interpretation and tuning compared with existing approaches. Simulation and empirical study show the superior performance of our method over the Facebook Prophet model and TBATS (which requires pre-specified seasonality periods). We have developed a comprehensive R package with documentation and plan to publish on CRAN.
  
  The paper is outlined below. In Section II, we generalize the classical seasonal ARMA model and propose the MS modeling procedure followed by discussions. In Section III, we compare MS modeling procedure with other methods by various simulations and real-data case studies and show its suitability in modeling multiple seasonality.
  
  
  
  § METHOD
  
  
  
  
   §.§ Metric
  
  
  Observing a time series X_t,t=1,…,N where N is the sample size, our goal is to perform n-step ahead prediction X̂_N+n with minimized cumulative mean squared error (CMSE).
  
  
      CMSE_n = 1/n∑_j=1^n (X_N+j - X̂_N+j)^2
  
  
  
  
   §.§ Model Formulation
  
  
  We first generalize the classical multiplicative seasonal ARMA model <cit.> with normal distributed noise. The model has pre-determined constants: r, the number of seasonality components; τ, the maximum neighborhood length of each seasonality component; p̅, the maximum AR components allowed; q̅, the maximum MA components allowed. It also has parameters p,q,C which need to be selected by an information criterion. The model is formed as following:
  
      X_t = Φ' X_t-1^t-p +    Ψ' ϵ_t-1^t-q + ∑_i=1^r (Γ_i' X_S_i + Λ_i' ϵ_S_i) + ϵ_t
  
  
  The observed data is X_t ∈ℝ,  where t=1,…,N, and auto-regressive component is defined as X_t-1^t-p = (X_t-1,…,X_t-p)', with coefficients Φ = (ϕ_1,…,ϕ_p)' ∈ℝ^p; The noise random variable is assumed as normal distributed ϵ_t ∼ ^i.i.d N(0,σ^2) ∈ℝ, and moving average component is defined as ϵ_t-1^t-q=(ϵ_t-1,…,ϵ_t-q)', with coefficients Ψ = (ψ_1,…,ψ_q)' ∈ℝ^q; Seasonality periods' collection C is a collection of sets C={S_i,  for 1≤ i ≤ r}, where r=|C| (the number of seasonality components) is pre-determined. S_i is defined as set of neighboring integers where S_i = {t-S_i,1,…,t-S_i,τ}, τ=|S_i| (the neighborhood length of each seasonality component) is pre-determined for all 1≤ i ≤ r. Accordingly, we let 
  
      X_S_i = (X_t-S_i,1,…,X_t-S_i,τ)' ∈ℝ^τ
         Γ_i = (γ_i,1,…,γ_i,τ)' ∈ℝ^τ
         ϵ_S_i=(ϵ_t-S_i,1,…,ϵ_t-S_i,τ)' ∈ℝ^τ
         Λ_i = (λ_i,1,…,λ_i,τ)' ∈ℝ^τ
  
  To exemplify the formulation, we take the selected model in Facebook Events Data case study in Section <ref> as an example. Three seasonality components are assumed in the model (r=3), since more than two seasonality are suspected in the original paper that include the data <cit.>. The neighborhood length is allowed to be 13 (τ=13). Parameters are selected as p=3,q=4,C={(116,…,128),(278,…,290),(360,…,372)}. The model is written as:
  
      X_t    = Φ' X_t-1^t-3 + Ψ' ϵ_t-1^t-4 + (Γ_1' X_t-116^t-128 + Λ_1' ϵ_t-116^t-128) 
            + (Γ_2' X_t-278^t-290 + Λ_2' ϵ_t-278^t-290) 
            + (Γ_3' X_t-360^t-372 + Λ_3' ϵ_t-360^t-372) + ϵ_t
  
  More details of this model will be discussed in Section <ref>. 
  
  
  
  For estimating coefficients Ω = {Φ,Ψ,σ,Γ_i,Λ_i for  1 ≤ i ≤ r}, maximum likelihood estimation (MLE) is utilized. In our package, the default optimization algorithm is BFGS <cit.>. The formula of log-likelihood function and its gradient is given below.
  
  Objective:        Minimize  l    = ∑_t=1^N ϵ_t^2    
  
         ϵ_t = X_t - Φ' X_t-1^t-p    - Ψ' ϵ_t-1^t-q - ∑_i=1^r (Γ_i' X_S_i + Λ_i' ϵ_S_i)    
  
  Gradient:      ∂l/∂Φ = 2∑_t=1^N ϵ_t ∂ϵ_t/∂Φ     , 
  ∂l/∂Ψ = 2∑_t=1^N ϵ_t ∂ϵ_t/∂Ψ    
  
        ∂l/∂Γ = 2∑_t=1^N ϵ_t ∂ϵ_t/∂Γ     , 
  ∂l/∂Λ = 2∑_t=1^N ϵ_t ∂ϵ_t/∂Λ    
  
  Where      ∂ϵ_t/∂Φ = - X_t-1^t-p -    Ψ' ∂ϵ_t-1^t-q/∂Φ - ∑_i=1^r Λ_i' ∂ϵ_S_i/∂Φ    
  
        ∂ϵ_t/∂Ψ = - ϵ_t-1^t-q -    Ψ' ∂ϵ_t-1^t-q/∂Ψ - ∑_i=1^r Λ_i' ∂ϵ_S_i/∂Ψ    
  
        ∂ϵ_t/∂Γ_i = - X_S_i -    Ψ' ∂ϵ_t-1^t-q/∂Γ_i - ∑_i=1^r Λ_i' ∂ϵ_S_i/∂Γ_i    
  
        ∂ϵ_t/∂Λ_i = - ϵ_S_i -    Ψ' ∂ϵ_t-1^t-q/∂Λ_i - ∑_i=1^r Λ_i' ∂ϵ_S_i/∂Λ_i    
  
  
  
  
   §.§ Modeling Procedure
  
  
  Based on the generalized seasonal ARIMA model, for a given r,τ, p̅, q̅, we introduce the following Multiple Seasonality (MS) modeling procedure.
  
  
  	
    * Detect seasonality candidates S_i's by spectrum analysis of X_t. The time series will be decomposed into (N/2)'s trigonometric components by Discrete Fourier Transformation (DFT) <cit.>. 
  	
      X_t = a_0 + ∑_j=1^N/2 [a_j cos(2π tj/N) + b_j sin(2π tj/N)]
  
  	Each trigonometric component is defined on an unique frequency 2π j/N assigned with a signal strength P(j/N)=a_j^2 + b_j^2. r+2's frequencies with the largest signal strength will be selected as candidates in C.
  	
    * At each parameter combination {p,q,C}, fit multiple seasonality model by MLE. Calculate the information criterion for each combination. Bridge criterion (BC) <cit.> is employed as the default information criterion in the package.
  	Then the combination of parameters is selected by optimizing the information criterion. Based on optimal model, perform forecast by X̂_N+j = E[X_N+j|X_1^N].
  
  
  
  
  The MS procedure's performance depends on the pre-specification of r,τ, p̅, q̅. In practice, τ should be considerably large to allow flexibility in fitting. The default value for τ in our developed package (and experiments) is 6. p̅, q̅ shouldn't be too large. 
  In practice, allowing overly large p̅, q̅ will overweight the importance of short-term dynamics and neglect potential seasonal effects. From various experiments, we found that the number of seasonality components r is the most influential parameter to the model, and it deserves sophisticated tuning. 
  
  
  
  
   §.§ Discussion
  
  In the sequel, we elaborate on some other aspects of the above model and procedure.
  
  	
    * Overfitting. Since our method has a linear ARMA structure and allows limited ranges of ARMA components, it is unlikely to have an overfitting issue.
  	
    * Detection & forecasting? The idea of our work is to propose a new modeling procedure to perform prediction, which takes into account multiple seasonality. The detection of seasonality is part of the method. We select seasonality by comparing the information criterion. While many state-of-the-art statistical methods such as TBATS, the Prophet model relies on pre-specification of seasonality, our method detects and suggests appropriate seasonality over a range of potential seasonality orders. As a result, our method performs well for multiple seasonality data scenarios where the seasonality pattern is not visually straightforward.
  	
    * Relation to ARIMA. The ARIMA model is an extension of the ARMA model when the data series are taken differencing operations. A usual ARIMA model can be expanded to ARMA expression. The point of introducing the ARIMA structure in the first place is to maintain the stationarity assumption of the ARMA model. Since we do not assume stationarity in our method, we only discussed the relationship between our method to the ARMA model.
  	
    * Stationarity. Since the essence of seasonality is the recurrence of long-term signals, 
  	a seasonal time series model is generally non-stationary, e.g. the model X_t = X_t-12 + ϵ_t. 
  	
    * Theory derivation. Because we do not assume stationarity in the time series, it is typically difficult to analyze the model's theoretical property unless the time series is assumed to be decomposed into a periodic part and stationary part or by assuming local stationarity <cit.>.
  	
    * Identifiability. As τ,p̅, q̅ increase, the identifiability of the model will be questionable. For example, the true model X_t = X_t-12 + ϵ_t 
  	may be also written as X_t = X_t-24 + ϵ_t-12 +ϵ_t. In general, we require that τ > max{p̅, q̅}.
  	
    * Information criterion. As stated in the last subsection, BC is taken as the default information criterion for the model selection, because of its theoretical advantages comparing to standard AIC or BIC  <cit.>. Nevertheless, we repeat our numerical and empirical studies with AIC and BIC as information criterion as well. They show similar results.
  	
    * “Detrended” Data. Currently, MS modeling is only suitable for detrended data. As future work, a component accommodating trend effects will be included in the formulation.
  	
    * Ultra-long term prediction. As a consequence of recurrent linear structure, the n-step ahead prediction will either converge or vibrate drastically when n is very large. A suitable correction technique may be used for ultra long term prediction <cit.>. From our experiments, this is not an issue if n<100 for N=1000.
  	
    * Computational complexity for long series. Since the MLE involves the computation of gradient recursively, the computational costs cumulate as the length of series increases. The computational complexity follows the order of O(N), where N is the sample size. The track of the computational time of MS procedure in single seasonality and double seasonality with trigonometric components in Section <ref> is shown below.
  	
  	
  	
  
  
  
  
  
  
  
  § EXPERIMENTAL STUDIES
  
  In this section, we evaluate the forecasting capabilities of the proposed method through various experimental studies.
  
  
  
   §.§ Numerical Examples
   
  
  In the following experiment studies, we will compare our MS modeling procedure to benchmark methods, including the vanilla ARMA model selected by BC, the Facebook Prophet model, the TBATS model (both with correct pre-determined seasonality orders) and a two-layer LSTM with 50 neurons on each layer. The LSTM model is a simple LSTM model recommended in the Pytorch tutorial for sin-wave time sequence prediction. We use this LSTM model because we try to compare our method with benchmark methods with similar complexity and computational cost. We train it individually in each repetition on each series by Adam optimizer, with learning rate 0.01 and 600 iterations. 
  
  The simulation includes four different data setups. In each setup, 30 parallel univariate sequences are generated, and the models are evaluated 50 times on each sequence. Specifically, in the i-th evaluation where 1≤ i ≤ 50, models will be trained on the first 650 + 5*(i-1) observations from the sequence, and predict the following 100. The mean value and standard error of CMSE for n-step ahead prediction is reported.
  
  In the first data setup, we consider single seasonality in trigonometric form. The reason to start with single seasonality is that most literature of forecasting time series with seasonality focus on single seasonality. We need to ensure that our MS procedure has satisfactory performance in a simple situation. The data is generated as an aggregation of a trigonometric component and a short-term random effect, which is the standard way to simulate time series with seasonality.
  
  The first column refers to different prediction length n. For each value of n, the mean value of cumulative mean squared error is recorded in the first row, and the standard error is recorded in parenthesis. Mean value with distance to smallest one in a row less than the summation of their standard errors, will be highlighted.
  
  
      X_t    = 10 sin(2π/50t) + Z_t 
      
      Z_t    = 0.8 Z_t-1 - 0.3 Z_t-2 + 0.5 ϵ_t-1 + ϵ_t 
      ϵ_t    ∼^i.i.d N(0,σ^2)  where σ=2
  
  
  
  The results show that MS has the best performance for 1-step ahead prediction. It also provides comparable results to TBATS in long-term prediction, where n=5,15. MS also outperforms Prophet when n=1,5,15.
  
  In the second data setup, we consider double seasonality in trigonometric form. Since double seasonality is the simplest situation in multiple seasonality regimes, we will emphasize simulation in a double seasonality situation. The data is generated as an aggregation of two trigonometric components with different frequencies and a short-term random effect.
  
  
      X_t    = 10 sin(2π/50t) + 5 sin(2π/15t)+ Z_t 
      
      Z_t    = 0.8 Z_t-1 - 0.3 Z_t-2 + 0.5 ϵ_t-1 + ϵ_t 
      ϵ_t    ∼^i.i.d N(0,σ^2)  where σ=2
  
  
  
  The results show that TBATS is better than MS for most of the prediction length. The results are conceivable because TBATS utilize trigonometric components to accommodate seasonality. Naturally, it is suitable for modeling seasonality generated in trigonometric form. Furthermore, seasonality orders of TBATS are correctly pre-determined. However, the Prophet model performs poorly due to the disturbance from short term random effect, although it also utilizes trigonometric components for seasonality. MS outperforms the Prophet model at all prediction lengths.
  
  Therefore, to further explore the methods' capability of modeling seasonality in a complex form, we consider non-trigonometric double seasonality in the third data setup. The seasonality component comes from the repetition of exogenous white-noise series instead of trigonometric functions. Then the observed data is generated as an aggregation of two non-trigonometric components with different frequencies and a short-term random effect.
  
  
      X_t    = A_t + B_t + Z_t 
      
      Z_t    = 0.8 Z_t-1 - 0.3 Z_t-2 + 0.5 ϵ_t-1 + ϵ_t 
      ϵ_t    ∼^i.i.d N(0,σ^2)  where σ=2
  
  
      A_t = 
      a_t ∼^i.i.d N(0,1)     for  0 < t < 100
      
      A_t-50    for  100≤ t
      
      
      B_t = 
      b_t ∼^i.i.d N(0,1)     for  0 < t < 100
      
      B_t-15    for  100≤ t
  
  
  
  The results show that MS have superiority when n=1,5,15, and TBATS model has best performance when n=50,100. MS is better than the Prophet at all prediction length. Overall, comparing to modeling trigonometric seasonality, MS has more potential for modeling non-trigonometric seasonality.
  
  In the fourth data setup, we add up both trigonometric and non-trigonometric components for double seasonality. The complexity of seasonality structure increased one more time.
  
  
      X_t    = 10 sin(2π/50t) + 5 sin(2π/15t)+ A_t + B_t + Z_t 
      
      Z_t    = 0.8 Z_t-1 - 0.3 Z_t-2 + 0.5 ϵ_t-1 + ϵ_t 
      ϵ_t    ∼^i.i.d N(0,σ^2)  where σ=2
  
  
      A_t = 
      a_t ∼^i.i.d N(0,1)     for  0 < t < 100
      
      A_t-50    for  100≤ t
      
      
      B_t = 
      b_t ∼^i.i.d N(0,1)     for  0 < t < 100
      
      B_t-15    for  100≤ t
  
  
  The results show that MS outperforms TBATS and Prophet at all prediction lengths, especially for short-term prediction. MS is suitable for modeling multiple seasonality in complex form.
  
  Generally speaking, in the simulation study above, MS modeling procedure shows comparable performance to the benchmark method such as TBATS and outperforms method like Prophet in most cases. MS also shows its suitability for modeling multiple seasonality in non-trigonometric form, which might be potentially useful in real data application, since the seasonality structure is usually more complicated than trigonometric functions in real data. Therefore, in the next subsection, we will evaluate MS procedure and other methods on four real-data case study, to further show its capability of modeling multiple seasonality.
  
  
  
  
   §.§ Empirical Studies
  
  
  
  
    §.§.§ PJM Load Forecasting
  
  In many countries worldwide, electricity is now traded under market rules using spot and derivative contracts <cit.>. At the corporate level, electricity load and price forecasts have become a fundamental input to energy companies’ decision making mechanisms. The costs of over- or undercontracting and then selling or buying power in the balancing market are typically so high that they can lead to huge financial losses and bankruptcy in the extreme case. <cit.> The risk goes up when pass on to consumers. A failure in forecasting electricity utility can eventually lead to disastrous power outage. (e.g. Manhattan blackout in July 2019)
  
  This case study focuses on forecasting electricity load for PJM, a regional transmission organization (RTO) that coordinates the movement of wholesale electricity in all or parts of 13 states and District of Columbia in US. The data is the hourly load of electricity across three regional company under PJM from 9/30/2018 to 10/1/2019. (PEPCO : Potomac Electric Power Company; PE : Pennsylvania Electric Company; COMED : ComEd) Forecasting will be performed weekly for the last five weeks on rolling basis. In this case, daily (period as 24) and weekly (period as 24*7) seasonality terms are specified in Prophet and TBATS model. Two seasonality components are allowed in MS procedure (r=2).
  
  
  
  
  
  The MS procedure indicates the seasonality terms at periods 21-27, 163-171, matching daily and weekly seasonality. It agrees with common knowledge of electricity load pattern. As it shows from the figure and table, MS procedure outperforms Prophet and TBATS model for most prediction lengths.
  
  
  
    §.§.§ Total Sun Spot Number
  
  
  This is a case study of monthly total sun spot number from 1749 to 2019. Existence of seasonality pattern in this univariate time series is well known. The seasonality period range from 9 to 14 years in history, and has averaged length as 11 years for last decades. MS, Prophet and TBATS model are evaluated 5 times, forecasting the following 150 steps. The seasonality period for Prophet and TBATS model is specified as 132 (11 years). Two seasonality components are allowed in MS (r=2).
  
  
  
  
  
  MS in this example suggests a range of seasonality period at 8-12 years instead of picking a single value. As it also shows from the figure and table, MS outperforms other two for most prediction lengths. This case reveals suitability of MS for modeling time series with unclear seasonality pattern. By allowing a range instead of fixing seasonality period, the model is able to adjust periodic pattern along the time series, creating more flexibility.
  
  
  
    §.§.§ Air Quality (CO Level)
  
  The data <cit.> contains the responses of a gas multi-sensor device deployed on the field in an Italian city. Hourly responses averages are recorded along with gas concentrations references from a certified analyzer. In this case study, we focus on forecasting CO level. MS, Prophet and TBATS are evaluated 10 times for forecasting the following 50 steps. Seasonality period as 24 is specified for Prophet and TBATS. Two seasonality components are allowed in MS (r=2).
  
  
  
  
  
  
  The MS procedure indicates seasonality periods at 18-30, 37-45. The former one corresponds to daily seasonality period, since the median of 18-30 is 24. The second one can be referring to the effect from two days ago. For this single seasonality data situation, MS provide satisfactory results, and outperforms other two for n=1,48.
  
  
  
    §.§.§ Facebook Events Data
  
  The data comes from the original paper of Facebook Prophet model <cit.>. It shows daily data for the number of events created on Facebook during the dates from 12/10/2007 to 01/20/2016. Weekly and yearly seasonality are strongly suspected. In this case study, MS,Prophet and TBATS are evaluated 7 times for forecasting the following 200 steps. Weekly and yearly seasonality periods are specified for Prophet and TBATS model, and three seasonality components are allowed in MS (r=3).
  
  
  
  
  
  The MS procedure indicates seasonality period range at 116-128, 278-290, 360-372, where the last one refer to yearly seasonality. MS provides comparable results to Prophet in terms of long term prediction, and has better performance for short term prediction. TBATS performs poorly for long-term prediction.
  
  
  
  § CONCLUSION
  
  In this paper, we have introduced MS modeling procedure, a simple yet powerful forecasting procedure that combines seasonality detection, estimation, model selection and forecasting for time series with multiple seasonality. The procedure doesn't require any pre-determined seasonality periods, and has a limited number of parameters to specify. We reveal its competence of forecasting by experimental and empirical studies. In most cases, MS outperforms the state-of-the-art such as the Facebook Prophet model. 
  
  
  
  
  IEEEtran
  
  
  