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

,


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. 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 T X = x , x ,....., x into multi-dimensional series 1 2 K X , 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) 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 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 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     ,..., ; 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.

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

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, where; is the spectral density of a regular ARMA(p,q) model.
Replacing by the Fourier frequencies = 2 ⁄ and adding [ ( )] to both sides, the periodogram{ } gives For closes to zero, i.e. for = 1,2, … , 2 The least-square estimator of long memory coefficient is given by

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): Compute k a and k b by (8) and (9) : 3.
Compute the values of ordinate () k I  by formula (10) ;  H 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.
Where ,    k 0 t T -1 and T2 c = 0if T is an odd number. The periodogram of X at the frequencies Define the kth element of the discrete Fourier transform of X as can be simplified as All frequency in the interval (0,0.5) is multiplied by two. This is done to ensure the following property;  [19] proposed to select those SVD components whose eigenvectors satisfy the following criteria: 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). The MAPE of the h-step ahead forecast is defined as 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.   (17).
This procedure is described in figure 1.

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.   Figure 3 shows that there has been a sharp increase in Mauna Loa Atmospheric CO2 Concentration with an additive seasonal pattern. 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. 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. 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.  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.