Hybrid Model of Singular Spectrum Analysis and ARIMA for Seasonal Time Series Data


CAUCHY –Jurnal Matematika Murni dan Aplikasi 
Volume 7(2) (2022), Pages 302-315 
p-ISSN: 2086-0382; e-ISSN: 2477-3344 

Submitted: December 01, 2021 Reviewed: December 10, 2021 Accepted: December 23, 2021 
DOI: http://dx.doi.org/10.18860/ca.v7i1.14136 

Hybrid Model of Singular Spectrum Analysis and ARIMA for 
Seasonal Time Series Data 

Gumgum Darmawan1,2,*, Dedi Rosadi1, Budi N Ruchjana2  

1Gadjah Mada University, Yogyakarta, Indonesia 
2Padjadjaran University, Bandung, Indonesia 

*Corresponding Author  
Email: gumgum.darmawan@gmail.ugm.ac.id*,  

dedirosadi@gadjahmada.edu, budi.nurani@unpad.ac.id 

ABSTRACT  

Hybrid models between Singular Spectrum Analysis (SSA) and Autoregressive Integrated Moving 
Average (ARIMA) have been developed by several researchers. In the SSA-ARIMA hybrid model, 
SSA is used in the decomposition and reconstruction process, while forecasting is done through 
the ARIMA model. In this paper, hybrid SSA-ARIMA uses two auto grouping models. The purpose 
of this paper is to analyze seasonal data using the SSA-ARIMA hybrid by auto grouping. The first 
model namely the Alexandrov method and the second method is alternative auto grouping with 
long memory approach. The two hybrid models were tested for two types of seasonal pattern, 
multiplicative and additive seasonal time series data. The analysis results using both methods give 
accurate result; as seen from the MAPE generated the 12 observations for future, the value is 
below 5%. For additive seasonal pattern, The hybrid SSA-ARIMA method with Alexandrov auto 
grouping is more accurate (MAPE= 0.13%) than the hybrid SSA-ARIMA method with Alternative 
method but for multiplicative seasonal pattern  the hybrid SSA-ARIMA with alternative auto 
grouping is more accurate (MAPE = 3.63%) than the hybrid SSA-ARIMA method with Alexandrov 
method.  
 
Keywords:  ARIMA; Automatic Grouping; Long Memory Effect; Seasonal Pattern, Singular 
Spectrum Analysis 

INTRODUCTION 

Singular Spectrum Analysis (SSA) is a relatively new non-parametric method that 
has proved its capability in various time series types. Solving all these problems 
correspond to the so-called basic capabilities of SSA. Besides, the method has several 
extensions. First, the multivariate version of the method permits the simultaneous 
expansion of several time series data; see, for example, [1]. Second, the SSA ideas lead to 
several forecasting procedures for time series; see [2]. Third, SSA has been utilized for 
change-point detection in time series. The SSA technique has been used as a filtering 
method in [3]. Fifth, a family of the causality test based on the multivariate SSA technique 
has been introduced in [4]. Sixth, SSA can be applied for missing value imputation [5].  

SSA can be applied in various disciplines, from mathematics and physics to 
economics and financial mathematics, meteorology and oceanography, to social sciences. 

http://dx.doi.org/10.18860/ca.v7i1.14136
mailto:gumgum.darmawan@gmail.ugm.ac.id
mailto:dedirosadi@gadjahmada.edu
mailto:budi.nurani@unpad.com


Hybrid Model of Singular Spectrum Analysis and ARIMA for Seasonal Time Series Data 

G.Darmawan 303 

For instance, in climatology ([6], [7], [8]) and biomedical data time series analysis [9]. 
Hybrid modeling of SSA in time series data has been carried out by many researchers. The 
hybrid model is carried out so that the advantages of two or more models make a positive 
contribution to the forecasting results. 

SSA hybrid model with other time series models includes ARIMA, Neural network, 
ARIMAX, PAR, VARIMAX, and others. [10], performed hybrid SSA with Neural 
Network.[11] perform the hybrid SSA-Algorithm Firefly-BP Neural Network process. [12] 
carried out a hybrid SSA model with ARMAX. [13] Combining the SSA model with PAR(p), 
this model was applied to wind speed data.[14], built the SSA-VARIMAX hybrid model and 
used it for climate data. The ARIMA model is often used as a comparison for the SSA model, 
such as [15], comparing SSA, ARIMA, and other time series models for tourism cases in 
various countries in Europe. The result has indicated that there is no good time series 
model for all tourism data. [16] compared SSA and ARIMA for predicting ambulance 
demand. 
 The SSA-ARIMA hybrid model studied by [17] was applied to the annual Runoff 
data. [18], the SSA-ARIMA hybrid model was compared with the basic SSA and ARIMA 
models. The result showed that the SSA-ARIMA hybrid model was the most accurate. 
However, many of these papers do not discuss specific data forms (e.g., seasonal patterns), 
so we consider it necessary to examine this hybrid model for seasonal data. In this study, 
the SSA and the ARIMA were employed collectively to forecast two types of time series 
data. Both models run to get fast and accurate computation. In SSA, there are two methods 
of automatic grouping (Alexandrov and Alternative). The forecasting performance of the 
hybrid SSA-ARIMA model was compared between the two methods (alternative vs. 
Alexandrov). This paper contributes to the analysis of the seasonal patterns (additive and 
multiplicative) by the SSA-ARIMA hybrid. The purpose of this paper was to analyze 
seasonal data using the SSA-ARIMA hybrid by auto grouping for two types of seasonal 
patterns.  

This paper was organized as follows: The current section was an introduction 
where we briefly outlined the use of SSA and introduced our study. In the next section, the 
methods section, we described the detailed methodology of SSA and ARIMA, briefly  
outlined forecasting using a linear recurrent formula, identification of fractional 
differencing parameter, identification of hidden periodicities based on Periodogram and 
Automatic grouping on Alexandrov Method ([19], [20]) also alternative Automatic 
grouping [21]. This section also included a proposed algorithm for automatic hybrid SSA-
ARIMA. In the results and discussion section, we demonstrated the abilities of hybrid SSA-
ARIMA in real-time series data. In this part, we also investigated three types of time series 
data: Seasonal with no trend, multiplicative seasonal with the trend, and Additive 
seasonal with the trend. This section also discussed the comparison result between 
hybrid SSA-ARIMA with the Alexandrov method and hybrid SSA-ARIMA with an 
alternative method for real data analysis.  

 

METHODS  

Singular Spectrum Analysis  

The (non-parametric) SSA method has received a fair amount of attention in the 
literature. The first phase of SSA is the decomposition, where the time series are broken 
down into four components: trend, seasonal, cyclical, and noise. This phase consists of the 
Embedding and Singular Value Decomposition steps. The second phase, namely the 
Reconstruction phase, consists of Grouping and Diagonal Average process. The 



Hybrid Model of Singular Spectrum Analysis and ARIMA for Seasonal Time Series Data 

G.Darmawan 304 

forecasting process can be done once the four stages have been completed. For the 
completeness of presentation of our method, we presented the complete phase of the SSA 
algorithm in the following section. 
 

Embedding 

 The embedding step will transform one-dimensional time series  1 2 TX = x , x , ....., x  into 

multi-dimensional series 1 2 KX , X , ..., X with vectors 𝑋 = (𝑋𝑖,𝑋𝑖+1,𝑋𝑖+2, . . ,𝑋𝑖+𝐿−1)
𝑇 ∈ 𝑅𝐿  

, where 𝑖 = 1,2,…,𝐾, 𝐾 = 𝑇 − 𝐿 + 1. The parameter window length L defines the 
embedding process, where 2 ≤ 𝐿 ≤ 𝑇 − 1 [22]. If we need to emphasize the size 
(dimension) of the vectors Xi, then we shall call them L-lagged vectors. The L-trajectory 
matrix (or simply the trajectory matrix) of the series X is defined as 

𝑋 = [

𝑥1 𝑥2 ⋯ 𝑥𝐾
𝑥2
⋮

𝑥3
⋮

…
⋯

𝑥𝐾+1
⋮

𝑥𝐿 𝑥𝐿+1 … 𝑥𝑇

] (1) 

The lagged vectors Xi are the columns of the trajectory matrix X. Both the rows and column 

of X are sub-series of the original series. The (i,j) element of matrix X is 𝑥𝑖𝑗 = 𝑥𝑖+𝑗−1 which 
yields that X has equal elements on the ‘antidiagonals’ i+j=const. Hence the trajectory matrix 

is a Hankel matrix.  

Singular Value Decomposition 

The second step, the SVD step, makes the singular value decomposition of the 
trajectory matrix X and represents it as a sum of rank-one bi-orthogonal elementary 
matrices. Set 𝑆 = 𝑋𝑋𝑇 and denoted by 𝜆1,𝜆2,…,𝜆𝐿the eigenvalues of S taken in the 
decreasing order of magnitude  (𝜆1 ≥ 𝜆2 ≥ ⋯ ≥ 𝜆𝐿 ≥ 0)and by U1, U2,…., UL  the 
orthonormal system of the eigenvectors of the matrix S corresponding to these 
eigenvalues. 

𝑑 = 𝑚𝑎𝑥{𝑖,𝑠𝑢𝑐ℎ 𝑡ℎ𝑎𝑡 𝜆𝑖 > 0} = 𝑟𝑎𝑛𝑘 𝑋  

If we denote  𝑉𝑖 =
𝑋𝑇𝑈𝑖

√𝜆𝑖
, then the SVD of the trajectory matrix can be written as  

𝑋 = 𝑋1 + 𝑋2 + ⋯+ 𝑋𝑑,  

 

where eigenvector Ui, Eigenvalues iλ  form matrix  𝑉𝑖
𝑇𝑋. The three elements of SVD 

forming are called eigen triple.  

Grouping 

The purpose of this step is to appropriately identify the trend, the oscillatory components 

with different periods and noise. This step can be skipped if one does not want to extract hidden 

information by regrouping and filtering components precisely. The grouping procedure 

partitions the set of indices 1,2,…., L into m disjoint subsets 𝐼 = 𝐼1, 𝐼2,…,𝐼𝑚, so the elementary 

matrix in equation (2) is regrouped into m groups. Let 𝐼 = {𝑖1, 𝑖2,…, 𝑖𝑝}. Then the resultant 

matrix Xi corresponding to the group i is defined as 𝑋𝑖 = 𝑋𝑖1 + 𝑋𝑖2+.. .+𝑋𝑖𝑝. The matrices are 

computed for I1, I2,…Im, and substituted into equation (2) to obtain the new expansion. 

The grouping process is the phase when the LxK matrix is grouped into several sub-groups, 

namely trend patterns, seasonal or periodic, and noise patterns. Here, in this paper, the patterns 

are identified by Fourier series analysis and long-memory analysis. Fourier series analysis is 



Hybrid Model of Singular Spectrum Analysis and ARIMA for Seasonal Time Series Data 

G.Darmawan 305 

used to identify a seasonal pattern, and long memory series analysis is used to identify the 

differencing parameter of data. We use the GPH method [23] to identify the differencing 

parameter of time series.  

 

Diagonal Averaging 

The next step in Basic SSA transforms each resultant matrix of the grouped 
decomposition (3) into a new one-dimensional series of length N and is called diagonal 
averaging. Let Y denote a matrix with orde (LxK), with the elements    ijy ,1 i L,1 j K , 

and define L* = min(L, K), K*=max(L, K), and T=L+K-1. Let 𝑦𝑖𝑗
∗ = 𝑦𝑖𝑗 if L<K and

*

ij ji
y y  

otherwise. Diagonal averaging transfers matrix Y into a series 1 2 Ny , y , ...., y by the formula. 

𝑦𝑘 =

{
 
 
 
 

 
 
 
 1

𝑘
∑ 𝑦𝑚,𝑘−𝑚+1

∗

𝑘

𝑚=1

𝑓𝑜𝑟 1 ≤ 𝑘 < 𝐿

1

𝐿
∑ 𝑦𝑚,𝑘−𝑚+1

∗

𝐿∗

𝑚=1

𝑓𝑜𝑟 𝐿∗ ≤ 𝑘 < 𝐾∗

1

𝑇 − 𝑘 + 1
∑ 𝑦𝑚,𝑘−𝑚+1

∗

𝑇−𝐾∗+1

𝑚=𝑘−𝐾∗+1

𝑓𝑜𝑟 𝐾∗ < 𝑘 ≤ 𝑇

(2) 

 
          Using equation (3), the resultant matrix 

iI
X , i = 1,2, ..., m  will form the series

        k k k ktY y y y1 2, ,..., ; therefore, the original sequence will be decomposed into the 
sum of the m series. 

 

Forecasting Using Linear Recurrent Formula 

The forecasting method used in this paper is the Linear Recurrent Formula (LRF). 
In the above notation, the recurrent forecasting algorithm (briefly, R-Forecasting) can be 
formulated as follows. 

1. The time series 𝑌𝑇+𝑀 = (𝑦1,𝑦2,…𝑦𝑇+𝑀)is defined by  

𝑦𝑖 = {
�̃�𝑖 𝑓𝑜𝑟 𝑖 = 1,2,…. ,𝑇

∑ 𝑎𝑗𝑦𝑖−𝑗
𝐿−1
𝑗=1 𝑓𝑜𝑟 𝑖 = 𝑇 + 1,…,𝑇 + 𝑀

(3)

                                

2. The numbers 𝑦𝑇+1,𝑦𝑇+2,…𝑦𝑇+𝑀 form the M terms of the recurrent forecast 

Thus, the R-forecasting is formed by the direct use of the LRF with 
coefficients {𝑎𝑗, 𝑗 = 1,…,𝐿 − 1}. 

Autoregressive Integrating Moving Average 

A time series  tx is said to follow the Autoregressive Integrated Moving Average 

model if the differenced that is 𝑊𝑡 = ∇
𝑑𝑥𝑡then this process is ARMA processes. If 𝑊𝑡 tw is 

ARMA(p,q), then {𝑋𝑡}  is ARIMA(p,d,q) process. The mathematic model of ARIMA(p,d,q) is 
as follow; 



Hybrid Model of Singular Spectrum Analysis and ARIMA for Seasonal Time Series Data 

G.Darmawan 306 

∅(𝐵)(1 − 𝐵)𝑑(𝑋𝑡 − 𝜇) = 𝜃(𝐵)𝑎𝑡 (4)          
   

Where 
t = index of observation  

d = differencing parameter (1,2,3,..) 

 = mean of observation  

𝑎𝑡~𝐼𝐼𝐷(0,𝜎𝑎
2) 

 

∅(𝐵) =  1 − ∅1𝐵 − ∅2𝐵
2 − ⋯− ∅𝑝𝐵

𝑝 is a polynomial of AR(p)    

𝜃(𝐵) =  1 − 𝜃1𝐵 − 𝜃2𝐵
2 − ⋯− 𝜃𝑝𝐵

𝑝 is polynomial of MA(q) 

Identification of Fraction differencing parameter  

The fractional differencing operator is described as a limitless binomial series 
growth in powers of the backward-shift operator [24]. One of the methods to estimate the 
value of the fractional differencing parameter is the Geweke and Porter-Hudak method 
(GPH). According to the GPH method, the differencing parameter of time series data can 
be estimated by the least square method of the ARFIMA (autoregressive fractionally 
integrated moving average) model. The GPH procedure can be described as follows. For 
the ARFIMA model given in (4), let 𝑊𝑡 = (1 − 𝐵)

𝛿𝑍𝑡 and let 𝑓𝑊(𝜔)and 𝑓𝑍(𝜔)be the 
spectral density function of {𝑊𝑡} and {𝑍𝑡}, respectively. Then,   

𝑓𝑧(𝜔𝑗) = 𝑓𝑊(𝜔𝑗){2𝑠𝑖𝑛(
𝜔𝑗

2
⁄ )}

−2𝛿

,𝜔𝑗 ∈ (−𝜋,𝜋) (5) 

 
where; 

𝑓𝑧(𝜔𝑗) =
𝜎𝑎
2

2𝜋
|
𝜃𝑞(𝑒

−𝑖𝜔𝑗)

𝜙𝑝(𝑒
−𝑖𝜔𝑗)

|

2

,   𝜔𝑗 ∈ (−𝜋,𝜋) is the spectral density of a regular ARMA(p,q) 

model. 

 ln[𝑓𝑧(𝜔𝑗)] = 𝛿ln{2𝑠𝑖𝑛(
𝜔𝑗

2
⁄ )}

−2

+ 𝑙𝑛[𝑓𝑊(𝜔𝑗)], 

ln[𝑓𝑧(𝜔𝑗)] = 𝑙𝑛[𝑓𝑊(0)] + 𝛿 ln{2𝑠𝑖𝑛(
𝜔𝑗

2
⁄ )}

−2

+ 𝑙𝑛[
𝑓𝑊(𝜔𝑗)

𝑓𝑊(0)
]. 

Replacing 𝜔𝑗 by the Fourier frequencies 𝜔𝑗 =
2𝜋𝑗

𝑇
⁄  and adding 𝑙𝑛[𝐼𝑍(𝜔𝑗)] to both 

sides, the periodogram{𝑍𝑡} gives 

ln[𝐼𝑍(𝜔𝑗)] = 𝑙𝑛[𝑓𝑊(0)] + 𝛿ln{2𝑠𝑖𝑛(
𝜔𝑗

2
⁄ )}

−2

+ 𝑙𝑛[
𝐼𝑍(𝜔𝑗)

𝑓𝑍(𝜔𝑗)
] + 𝑙𝑛[

𝑓𝑊(𝜔𝑗)

𝑓𝑊(0)
]. 

For  𝜔𝑗 closes to zero, i.e. for 𝑗 = 1,2,…,
𝑇
2⁄  such that 

𝑚
𝑇⁄ → ∞ we have 𝑙𝑛[

𝑓𝑊(𝜔𝑗)

𝑓𝑊(0)
] ≈ 0, 

thus, 

𝑌𝑗 = 𝛽0 + 𝛽1𝑋𝑗 + 𝑒𝑗,    𝑗 = 1,2,…,𝑚 (6) 

𝑌𝑗 = 𝑙𝑛[𝐼𝑍(𝜔𝑗)],   𝑋𝑗 = 𝑙𝑛(
1

4𝑠𝑖𝑛2 (
𝜔𝑗

2
⁄ )

) 



Hybrid Model of Singular Spectrum Analysis and ARIMA for Seasonal Time Series Data 

G.Darmawan 307 

�̂�0 = 𝑙𝑛[𝑓𝑊(0)],𝐼𝑍 =
1

2𝜋
[𝛾0 + ∑𝛾𝑡𝑐𝑜𝑠(𝜔𝑗. 𝑡)

𝑇−1

𝑡=0

]    

 
The least-square estimator of long memory coefficient is given by 

 

�̂�1 = 𝛿 =
∑ (𝑋𝑗 − �̅�)
𝑚
𝑗=1 (𝑌𝑗 − �̅�)

∑ (𝑋𝑗 − �̅�)
2𝑚

𝑗=1

;𝑚 = 𝑔(𝑇) = 𝑇𝛼,0 < 𝛼 < 1 (7) 

 

Identification of hidden Periodicities based on Periodogram  

Identification of hidden periodicities based on Periodogram analysis can be described 
as follows; 

1. Fit the data to equation (7): 

𝑋𝑡 = ∑(𝑎𝑘𝑐𝑜𝑠𝜔𝑘𝑡 + 𝑏𝑘𝑠𝑖𝑛𝜔𝑘𝑡),

[𝑇 2⁄ ]

𝑘=0

(8) 

where k = 0,1,…,[T/2] and k  
is Fourier frequency determined by 𝜔𝑘 =

2𝜋𝑘
𝑇⁄ . 

2. Compute ka and kb  by (8) and (9) : 

𝑎𝑘 =

{
 
 

 
 1

𝑇
∑𝑥𝑡𝑐𝑜𝑠𝜔𝑘𝑡

𝑇

𝑡=1

𝑓𝑜𝑟 𝑘 = 0 𝑎𝑛𝑑 𝑘 =
𝑇

2
 𝑖𝑓 𝑇 𝑒𝑣𝑒𝑛

2

𝑇
∑𝑥𝑡𝑐𝑜𝑠𝜔𝑘𝑡

𝑇

𝑡=1

𝑓𝑜𝑟  𝑘 = 1,2,…,
(𝑇 − 1)

2

(9) 

𝑏𝑡 =
2

𝑇
∑𝑥𝑡𝑠𝑖𝑛𝜔𝑘𝑡,

𝑇

𝑡=1

                     𝑘 = 1,2,…
(𝑇 − 1)

2
(10) 

3. Compute the values of ordinate ( )kI   by formula (10) ; 

𝐼(𝜔𝑘) =

{
 
 

 
 

𝑇𝑎0
2 𝑘 = 0

𝑇

2
(𝑎𝑘

2 + 𝑏𝑘
2) 𝑘 = 1,2,…,[

𝑇 − 1

2
]

𝑇𝑎𝑇
2⁄

2 𝑘 = [
𝑇

2
]

(11) 

      
4. Test for significance of every Fourier frequency as follow; 

0
H : 0    (data do not have a seasonal pattern)   

1
H : 0   or 0   (data have a seasonal pattern)  

Statistic test: 

𝑇 =
𝐼(1)(𝜔(1))

∑ 𝐼(𝜔𝑘)
[𝑇 2⁄ ]

𝑘=1

(12) 



Hybrid Model of Singular Spectrum Analysis and ARIMA for Seasonal Time Series Data 

G.Darmawan 308 

 

   

     
where, 

(1)
(1)I (ω ) : maximum ordinate of periodogram of Fourier  frequency 

kI(ω ) 
 : Value of periodogram ordinate at k-th’ Fourier frequency.  

Test criteria: 

Reject 0H  if T >gα with α = significant level. The value of gα can be seen in table 

Fisher [25]. 

Automatic Grouping in SSA 

In this section, we explain two methods SSA automatic grouping based on trend 
extraction based namely Alexandrov, and SSA automatic grouping based on long memory 
approach, namely Alternative method. First, we introduce the periodogram of a time 
series. Let us consider the Fourier representation of the elements of a time series  𝑋 =
(𝑥1,𝑥2,…,𝑥𝑇)of length T. 

𝑋𝑡 = 𝐶0 + ∑ (𝑐𝑘𝑐𝑜𝑠(
2𝜋𝑡𝑘

𝑇⁄ )) + 𝑠𝑘𝑐𝑜𝑠(
2𝜋𝑡𝑘

𝑇⁄ ) +
(−1)𝑡𝐶𝑇

2⁄
,

1≤𝑘≤
𝑁−1

2

(13)
 

Where ,  k 0 t T - 1  and 
T 2c = 0if  T is an odd number. The periodogram of X 

at the frequencies  
N 2

k=0
ω k N  is defined as 

𝐼𝑋
𝑁(𝑘 𝑇⁄ ) =

𝑇

2
{
 

 
2𝑐0

2 𝑘 = 0

(𝑐𝑘
2 + 𝑠𝑘

2) 0 < 𝑘 < 𝑇 2⁄

2𝑐𝑇
2⁄

2 𝑘 = 𝑇 2⁄ ,𝑇 = 𝑒𝑣𝑒𝑛

(14) 

Define the kth element of the discrete Fourier transform of X as 

𝐹𝑘(𝑋) = ∑ 𝑒
−𝑖2𝜋𝑡𝑘 𝑇⁄ 𝑥𝑡

𝑇−1
𝑡=0 .  

The periodogram  TXI  at the frequencies  
2

0

T

k
k T

  


 can be simplified as 

𝐼𝑋
𝑇(𝑘 𝑇⁄ ) =

1

𝑇
{
2|𝐹𝑘(𝑋)|

2 𝑓𝑜𝑟 0 < 𝑘 < 𝑇/2

|𝐹𝑘(𝑋)|
2 𝑓𝑜𝑟  𝑘 = 0 𝑜𝑟 𝑇 𝑖𝑠 𝑒𝑣𝑒𝑛 𝑎𝑛𝑑 𝑘 = 𝑇/2

 

All frequency in the interval (0,0.5) is multiplied by two. This is done to ensure the 
following property;  

‖𝑋‖2
2 = ∑ 𝑥𝑡

2𝑇−1
𝑡=0 = ∑ 𝐼𝑋

𝑁(𝑘 𝑇⁄ )
⌊𝑇 2⁄ ⌋

𝑘=0
        

Where 𝜋𝑋
𝑇(𝜔) = ∑ 𝐼𝑇

𝑋(𝑘 𝑇⁄ ),𝜔 ∈
[0,0.5]𝑘:0≤𝑘≤𝜔  is the cumulative contribution of 

the frequencies 0, ? Then, for given  0 0, 0.5  , [19] define the contribution of low 

frequencies from the interval  00,  to 
T

X R as 

𝐶(𝑋,𝜔0) = 𝜋𝑋
𝑇(𝜔0) 𝜋𝑋

𝑇(0.5)⁄ (15) 

Then, given parameters  0ω 0,0.5 and  0ω 0,1 ,[19] proposed to select those 
SVD components whose eigenvectors satisfy the following criteria: 

𝐶(𝑈𝑗,𝜔0) ≤ 𝐶0 (16) 



Hybrid Model of Singular Spectrum Analysis and ARIMA for Seasonal Time Series Data 

G.Darmawan 309 

where Uj is a corresponding jth eigenvector. 

The second procedure (alternative method) of the automatic Grouping Algorithm 

can be described using the following steps: 

1. Denote the original series as 𝑋 = (𝑥1,𝑥2,…,𝑥𝑇), the univariate time series data 

with length T 

2. Define Window Length (L) 2 ≤ 𝐿 ≤ 𝑇 2⁄ .. 

3. Determine the shape of the trajectory matrix based on L as in equation (3). 

4. Find the Eigen Triple  i i iλ ,V ,U , i = 1,2,….,d (d is the number of positive 
eigenvalues) of the XXT matrix, where X is the Trajectory matrix. 

5. Identify the ̂ (Order difference) of each series Ui by equation (6) in section 2.3. 

6. Make a differencing for all Ui based on the value of ̂  (step 5) to obtain the UDi 

series. Using this step, then we obtain data with the stationary pattern.  

7. Specify the maximum value of seasonal order of Pi of  UDi; Here, the  Pi values 

are between 0 and 24.  

8. Group UDi series based on Pi values obtained from the previous step 7, a 

grouping of UDi is based on periods, 3,6,12, and 24, UDi with  Pi=0 is grouped 

as non-seasonal data, where when Pi>24 is considered as cycle pattern. 

9. Create diagonal averaging for grouping components from step 8; see equation 

(3).  

10. Forecasting (Linear Recurrent Method) with equation (4). 

11. Determine MAPE values based on the number of our sample data. 

The MAPE of the h- step ahead forecast is defined as 

𝑀𝐴𝑃𝐸 = (
1

ℎ
∑|

𝑦𝑛+ℎ − �̂�𝑛+ℎ
�̂�𝑛+ℎ

|

ℎ

𝑖=1

) × 100% (17) 

 

where n denotes the last data used in the in-sample data and ˆ
n h

y


denotes the 

forecast values. 

12. Go to step 2 until all of the L choices have been calculated. The best model will 

give the smallest MAPE value. 

 Algorithm for hybrid SSA-ARIMA 

In this section, we proposed our new automatic algorithm on Singular Spectrum 

Analysis. The procedure of our automatic Grouping Algorithm can be described 

using the following steps: 

1. Denote the original series as 𝑋 = (𝑥1,𝑥2,…,𝑥𝑇), the univariate time series data 

with length T. 

2. Define Window Length (L) 2 ≤ 𝐿 ≤ 𝑇 2⁄ . 

3. Determine the shape of the trajectory matrix based on L as in equation (3). 

4. Find the Eigen Triple  i i iλ ,V ,U , i = 1,2,….,d (d is the number of positive 
eigenvalues) of the XXT matrix, where X is the Trajectory matrix. 



Hybrid Model of Singular Spectrum Analysis and ARIMA for Seasonal Time Series Data 

G.Darmawan 310 

5. Use one method of automatic grouping procedures (section 2.4). This procedure 

takes the best parameter (window length (L)).   

6. Build Model ARIMA of all grouping (RC’s) by Auto-ARIMA procedures. 

7. Forecast a new series of all groups by model ARIMA. 

8. Sum all results of forecast to get predicted of new time series.   

9. Determine MAPE values based on the number of our sample data (17). 

This procedure is described in figure 1. 

 

Figure 1. Flowchart of Hybrid SSA-ARIMA 

RESULTS AND DISCUSSION  

For our empirical study, we analyzed AirPassengers, and CO2 data, available in R 
package datasets. Here, we consider two seasonal data, namely AirPassengers, which has 
a trend and multiplicative seasonal pattern; and CO2, which contains an additive seasonal 
pattern with the trend.  



Hybrid Model of Singular Spectrum Analysis and ARIMA for Seasonal Time Series Data 

G.Darmawan 311 

  
Figure 2. AirPassengers Figure 3. CO2 

 
Figure 2 reveals that there has been a gradual increase in several passengers with a 
multiplicative seasonal pattern from 1949 to 1960. Figure 3 shows that there has been a 
sharp increase in Mauna Loa Atmospheric CO2 Concentration with an additive seasonal 
pattern. 

 
Table 1. Grouping and Contribution of Data CO2 

RC 

Hybrid-SSA-ARIMA 
(Alexandrov Grouping) 

Hybrid-SSA-ARIMA 
(Alternative Grouping) 

Principle Component (PC) Contribution Principle Component (PC) Contribution 

RC1 
PC1,PC2,PC3,PC6 

PC7,PC8 
99.49% 

PC1,PC8,PC9,PC16,PC17 
PC21,PC23,PC25 

98.40% 

RC2 
PC4,PC5,PC12,PC15,PC16 

PC17 
0.28% 

PC2,PC3,PC10,PC12,PC13 
PC14,PC15,PC18,PC19,PC20 

0.93% 

RC3 
PC10,PC11,PC18,PC19 

PC22,PC28,PC29 
0.07% PC24 0.01% 

RC4 PC13,PC14,PC26,PC30,PC31 0.06% 
PC4,PC5,PC6,PC22,PC28 

PC31,PC32,PC33,PC34,PC35 
0.40% 

RC5 
PC20,PC21,PC23,PC24,PC25 

PC27 
0.06% 

PC11,PC26,PC27,PC37,PC39 
PC40,PC41 

0.06% 

CO2 data had been successfully decomposed into five groups through SSA, each group 
consisting of 5 or 6 Principle Components. PC1 to PC31 were PCs that had successfully 
identified the pattern, while PC32 and so on were PCs with random patterns. RC1 was the 
largest contribution of 99.49%; this showed that the trend pattern was very dominant in 
this data. 
On the right side of this table, CO2 data had been successfully decomposed into five groups 
through SSA; the number of PCs for each RC is uneven. RC3 only consists of 1 PC. PC1 to 
PC41 were PCs that had successfully identified the pattern, while PC42 and so on were 
PCs with random patterns. RC1 had the largest contribution of 98.49%; this showed that 
the trend pattern was very dominant in this data. 

 
Table 2. Separability of  CO2 

RC 
Hybrid-SSA-ARIMA 

 (Alexandrov Grouping) 
Hybrid-SSA-ARIMA  

(Alternative Grouping) 
RC1 RC2 RC3 RC4 RC5 RC1 RC2 RC3 RC4 RC5 

RC1 1 0.021 0.023 0.000 0.006 1 0.010 0.006 0.351* 0.000 
RC2 0.021 1 0,000 0,001 0,001 0.010 1 0.000 0.000 0.001 
RC3 0.023 0,000 1 0,001 0.001 0.006 0.000 1 0.000 0.002 
RC4 0.000 0,001 0,001 1 0.001 0.351* 0.000 0.000 1 -0.002 
RC5 0.00 0,001 0.001 0.001 1 0.000 0.001 0.002 -0.002 1 

Time

co
2

1960 1970 1980 1990

3
2

0
3

3
0

3
4

0
3

5
0

3
6

0

Time

A
ir

P
a

s
s
e

n
g

e
rs

1950 1952 1954 1956 1958 1960

1
0

0
2

0
0

3
0

0
4

0
0

5
0

0
6

0
0



Hybrid Model of Singular Spectrum Analysis and ARIMA for Seasonal Time Series Data 

G.Darmawan 312 

This table shows the correlation values between groups from RC1 to RC5. This best model 
turned out to not have a very good separability which can be seen from the correlation 
between RC1 and RC4, which was quite large. However, the correlation values for the 
others were very small. 
The right side of this table shows the correlation values between groups from RC1 to RC5. 
This best model turned out to not have a very good separability. It can be seen from the 
correlation between RC1 and RC4, which is quite large. However, the correlation values 
for the others were very small.  
 

Table 3. Grouping and  Contribution of AirPassengers 

RC 
Hybrid-SSA-ARIMA 

(Alexandrov Grouping) 
Hybrid-SSA-ARIMA 

(Alternative Grouping) 
Principle Component (PC) Contribution Principle Component (PC) Contribution 

RC1 PC1,PC2,PC3,PC6,PC11 
PC14,PC15 

83.58% PC15,PC19 0.80% 

RC2 PC4,PC5,PC16,PC17,PC23 
PC24,PC29 

7.20% PC1,PC2,PC3,PC16,PC17,PC18
PC21,PC22,PC23,PC28,PC29 

81.51% 

RC3 PC9,PC10,PC19,PC20 
PC33,PC35,PC36 

3.26% PC4,PC5,PC6,PC20,PC24,PC26,
PC30 

8.54% 

RC4 PC7,PC8,PC25,PC26,PC30 
PC31 

3.13% PC7,PC8,PC9,PC10,PC11 
PC33,PC34,PC36,PC37 

5.98% 

RC5 PC12,PC13,PC18,PC21,PC22 
PC27,PC28,PC34 

2.91% PC31,PC32 0.31% 

 

AirPassengers data had been successfully decomposed into five groups through SSA with 
each group consisted of 6 to 8 Principle Components. PC1 to PC36 were PCs that had 
successfully identified the pattern, while PC37 and so on were PCs with random patterns. 
RC1 had the largest contribution of 83.589% which indicated that the trend pattern 
dominated in this data.  
AirPassengers data had been successfully decomposed into five groups via SSA, RC1 and 
RC5 consisted of 2 PCs, RC2 consists of 11 PCs, RC3 consisted of 7 PCs, and RC4 9PCs. PC1 
to PC36 were PCs whose patterns had been identified, while PC38 and so on were PCs 
with random patterns. RC2 had the most considerable contribution of 81, 51%, which 
indicated that seasonal patterns dominated in this data. 
 

Table 4. Separability of AirPassengers 

RC 
Hybrid-SSA-ARIMA 

(Alexandrov Grouping) 
Hybrid-SSA-ARIMA 

(Alternative Grouping) 
RC1 RC2 RC3 RC4 RC5 RC1 RC2 RC3 RC4 RC5 

RC1 1 -0.053 -0.003 0.008 -0.015 1 -0.017 0.000 0.000 0.000 
RC2 -0.053 1 0.005 0.001 0.000 -0.017 1 -0.133 0.033 -0.001 
RC3 -0.003 0.005 1 -0.003 -0.002 0.000 -0.133 1 0.000 -0.002 
RC4 0.008 0.001 -0.003 1 -0.001 0.000 0.033 0.000 1 0.000 
RC5 -0.015 0.000 -0.002 -0.001 1 0.000 -0.001 -0.002 0.000 1 

 

 
This table showed the correlation values between groups from RC1 to RC5. This best 
model displayed good separability, which could be seen from RC's very small correlation 
value. The Hybrid-SSA-ARIMA (Alternative Grouping) model showed the correlation 
value between groups from RC1 to RC5. This best model showed good separability which 
could be seen from the very small correlation value between RC. 

 

 



Hybrid Model of Singular Spectrum Analysis and ARIMA for Seasonal Time Series Data 

G.Darmawan 313 

Table 5. MAPE of Two Methods from AirPassengers dan CO2  
Data Hybrid-SSA-ARIMA 

(Alternative Grouping) 
Hybrid-SSA-ARIMA 

(Alexandrov-Grouping) 
AirPassengers 3.63% 4.64% 
CO2 0.48% 0.13% 

 
This table showed that the AirPassengers data using the Hybrid-SSA-ARIMA (Alternative 
Grouping) method was more accurate than the hybrid-SSA-ARIMA (Alexandrov grouping) 
method. The MAPE value of the Hybrid-SSA-ARIMA (Alternative Grouping) method was 
3.63%, while the the value of Hybrid-SSA-ARIMA (Alexandrov grouping) method was 
4.64%. On the CO2 data, on the other hand, the Hybrid SSA-ARIMA (Alexandrov grouping) 
method was more accurate than the Hybrid SSA-ARIMA (Alternative Grouping) method 
because the MAPE value for CO2 data using the hybrid-SSA-ARIMA (Alexandrov grouping) 
method was 0.13%, while for the hybrid-SSA-ARIMA (Alternative Grouping) method, the 
value was 0.48%. 
 

CONCLUSIONS 

Based on the analysis results, it appeared that the SSA-ARIMA method with 
Alexandrov auto grouping was better for data with additive seasonal patterns such as CO2 
data. The MAPE value for this method was 0.13% for CO2 data while the SSA-ARIMA 
method with alternative auto grouping was 0.48% for CO2 data. On the other hand, for 
data with multiplicative seasonal patterns such as AirPassengers data, the SSA-ARIMA 
hybrid method with alternative auto grouping produced more accurate forecasts with the 
MAPE value of 3.63%. In comparison, the SSA-ARIMA hybrid method with Alexandrov 
auto grouping containing the same data gave a MAPE value of 4.64%. 

In general, the separability value had good separability as the correlation value 
was close to zero - only one significant correlation value is 0.351. The correlation value 
between RC1 and RC4 for CO2 data generated from the SSA-ARIMA hybrid method with 
alternative auto grouping. As for the contribution, in general, it is dominated by RC1, 
except for AirPassengers data which was analyzed through the SSA-ARIMA auto hybrid 
method with alternative auto grouping; the contribution is dominated by RC2, amounting 
to 81.51%.   In general, both the SSA-ARIMA  with Alexadarov auto grouping and the SSA-
ARIMA with alternative auto grouping provide satisfactory forecast results. The MAPE 
value of the two methods was below 5% for the two different data patterns. 

ACKNOWLEDGMENTS  

The author thanks the Rector of Universitas Padjadjaran who has provided financial 
assistance for dissemination lecturer and student research through the Academic 
Leadership Grant through contract number:1959/UN6.3.1/PT.00/2021. 
 

REFERENCES 

[1]  K. Patterson, H. Hassani, S. Heravi and A. Zhigljavsky, "Multivariate Singular 
Spectrum Analysis for forecasting revisions to real-time series data," 
Journal.Application.Stat,vol.38,no.10,pp.2183-2211 2011.  

[2]   H. Hassani, E. S. Silva, N. Antonakakis, G. Filis, and R. Gupta, Forecasting Accuracy 
evaluation of tourist arrivals, ann. Tour. Res, no.63.pp.112-127, 2017.  



Hybrid Model of Singular Spectrum Analysis and ARIMA for Seasonal Time Series Data 

G.Darmawan 314 

[3]  T. J. Harris and H. Yuan, “Filtering and frequency interpretations of Singular 
Spectrum Analysis,” Phys. D Nonlinear Phenom, vol. 239, no. 20–22, pp. 1958–1967, 
2010.  

[4]  R. Menezes, A. Dionísio, and H. Hassani, "On the globalization of stock markets: An 
application of Vector Error Correction Model, Mutual Information and Singular 
Spectrum Analysis to the G7 countries," Q. Rev. Econ. Finance., vol. 52, no. 4, pp. 369–
384, 2012.  

[5]  P. C. Rodrigues and M. de Carvalho, "Spectral modeling of time series with missing 
data," Appl. Math. Model., vol. 37, no. 7, pp. 4676–4684, 2013.  

[6]  C. L. Wu and K. W. Chau, "Rainfall-runoff modeling using artificial neural network 
coupled with singular spectrum analysis," J. Hydrology., vol. 399, no. 3–4, pp. 394–
409, 2011. 

[7] C. Sivapragasam, S. Liong, and M. Pasha, “Rainfall and runoff forecasting with SSA-
SVM approach,” J. Hydroinformatics, no. April 2016, pp. 141–152, 2001. 
 

[8] R. Benzi, R. Deidda, and M. Marrocu, “Characterization of temperature and 
precipitation fields over Sardinia with principal component analysis and singular 
spectrum analysis,” Int. J. Climatol., vol. 17, no. 11, pp. 1231–1262, 1997. 

[9] S. Hassani Hossein; Sanei, Singular spectrum analysis of biomedical signals. 2016. 
 

[10] F.Lisi, O.Nicolas, M.Sandri, "Combining Singular Spectrum Analysis and Neural 
Networks for time series letters, vol.2, No. forecasting", Neural processing letters, 
vol.2, No.4, 6-10, 1995. 

[11] Y, Gao, C, Qu and K, Zhang, "A Hybrid Method Based on Singular Spectrum Analysis, 
Firefly Algorithm and BP Neural Network for Short-Term Wind Speed Forecasting," 
Energies 9, 757. 

[12] K. Laksmi, A.Rao, and N.Gopalakhrisnan, "Singular Spectrum Analysis combined 
with ARMAX model for structural damage detection," structure Control Health 
Monit,1-21, 2016. 

[13] M.L, Menezes, R.C, Sousa and Pessanha, "Combining Singular spectrum Analysis and 
PAR(p) Structures to Model Wind Speed Time Series," J Syst Sci Complex,27,29-46, 
2014. 

[14] Y.Feliks, Y.Groth, A.W, Robertson and M, Gill, "Oscillatory climate modes in the 
Indian monsoon, North Atlantic, tropical Pacific" Journal of Climate,vol26, 9528- 
9544,2013. 

[15] 
 

H.Hassani, E.Sirimal,N. Antonakakis and R.Gupta," Forecasting Accuracy Evaluation 
of Tourism Arrivals," Annals of Tourism Research,63,112-127, 2017. 

[16] 
 

J.Vile, J.W, Gillard, P.R, Harper and V.A.Knight," Predicting ambulance demand using 
Singular Spectrum Analysis, "Journal of Operation Research Society,63,1556-
1565,2012. 

[17] 
 

O.Zhang,B.D,Wang,B.He, Y.Peng and M.LRen” Singular Spectrum Analysis and 
ARIMA hybrid Model for Annual Runoff Forecasting, "Water Resource 
Manage,25:2683-2703,2011. 

[18] 
 
 

G, Liu, G.Zhang and T.Zhang ”Software Reliability Forecasting: Singular Spectrum 
Analysis and ARIMA Hybrid Model” 2015 International Symposium on Theoretical 
Aspects of Software Engineering, IEEE Computer Society,111-118, 2015. 



Hybrid Model of Singular Spectrum Analysis and ARIMA for Seasonal Time Series Data 

G.Darmawan 315 

[19] T. Alexandrov, “A Method of Trend Extraction Using Singular Spectrum Analysis,” 
pp. 1–23, 2008 

[20] T. Alexandrov and N. Golyandina, “Automatic extraction and forecast of time series 
cyclic components within the framework of SSA,” Proc. St. Petersburg. Work. Simul., 
no. September, pp. 45–50, 2005 

[21] G.Darmawan, D.Rosadi, and B.N, Ruchjana, "Daily Forecast for COVID-19 During    
Ramadhan by Singular Spectrum Analysis in Indonesia”, Jour of Adv Research in 
Dynamical & Control Systems, Vol. 12, Issue-02, 2020. 

[22] N. Golyandina, A. Zhigljavsky, G. N., and Z. A., Singular Spectrum Analysis for Time 
Series. 2013. 

[23] S. Geweke J and Porter-Hudak, “The Estimation and Application of Long memory 
Time Series Models,” vol. 4, pp. 221–238, 1983. 

[24] J. R.. Hosking, "Fractional differencing," vol. 68, no. 1, pp. 165–176, 1981. 

[25] W. Wei, Time Series Analysis Univariate and Multivariate Methods. California: 
Addison-Wesley Publishing Company, 2006