Improving forecasting accuracy of time series data using a new ARIMA-ANN hybrid method and empirical mode decomposition

TL;DR

This paper introduces a new ARIMA-ANN hybrid method combined with empirical mode decomposition to improve time series forecasting accuracy across various applications.

Contribution

It proposes a more general hybrid framework that enhances forecasting performance by optimizing data decomposition and model combination strategies.

Findings

Hybrid method outperforms traditional approaches

Data decomposition strategies are crucial for accuracy

Effective in diverse time series applications

Abstract

Many applications in different domains produce large amount of time series data. Making accurate forecasting is critical for many decision makers. Various time series forecasting methods exist which use linear and nonlinear models separately or combination of both. Studies show that combining of linear and nonlinear models can be effective to improve forecasting performance. However, some assumptions that those existing methods make, might restrict their performance in certain situations. We provide a new Autoregressive Integrated Moving Average (ARIMA)-Artificial Neural Network(ANN) hybrid method that work in a more general framework. Experimental results show that strategies for decomposing the original data and for combining linear and nonlinear models throughout the hybridization process are key factors in the forecasting performance of the methods. By using appropriate strategies, our hybrid method can be an effective way to improve forecasting accuracy obtained by traditional hybrid methods and also either of the individual methods used separately.

Tables

Table 1. Table 1: Performance Comparison for All Datasets

Datasets	Methods Metrics	ANN	ARIMA	Zhang’s Method	Khashei-Bijari’s Method	Babu-Reddy’s Method	Proposed Method
Sunspot	MAE	14.23	13.37	13.14	10.62	11.39	10.48
	MSE	353.12	306.97	289.31	205.08	239.90	194.29
	MASE	0.629	0.591	0.581	0.470	0.504	0.463
Lynx	MAE	0.1249	0.1198	0.1003	0.1025	0.1102	0.1013
	MSE	0.0241	0.0231	0.0173	0.0175	0.0189	0.0162
	MASE	0.6185	0.5932	0.4966	0.50757	0.5457	0.5016
Gbp/Usd	MAE	428.55	435.72	429.52	406.22	436.34	404.90
	MSE	3.4681	3.5272	3.4496	3.1053	3.5053	2.9538
	MASE	1.085	1.103	1.087	1.028	1.104	1.025
Intraday	MAE	20.10	20.22	19.16	19.79	19.50	18.81
	MSE	617.46	652.72	594.09	600.93	619.67	581.38
	MASE	1.048	1.054	0.999	1.031	1.016	0.980

Table 2. Table 2: Performance Comparison of methods using EMD for All Datasets

Datasets	Methods Metrics	ANN	ARIMA	Zhang’s Method	Khashei-Bijari’s Method	Babu-Reddy’s Method	Proposed Method
Sunspot	MAE	8.33	7.72	7.46	7.76	7.92	7.28
	MSE	120.14	99.07	90.04	99.17	100.64	87.86
	MASE	0.368	0.341	0.330	0.343	0.350	0.322
Lynx	MAE	0.0912	0.0772	0.0764	0.0782	0.0788	0.0760
	MSE	0.01318	0.00996	0.00995	0.00998	0.01014	0.00923
	MASE	0.4516	0.3822	0.3783	0.3872	0.3902	0.3763
Gbp/Usd	MAE	190.08	146.03	142.92	146.11	147.18	141.38
	MSE	0.5610	0.3578	0.3479	0.3581	0.3633	0.3285
	MASE	0.4812	0.3697	0.3618	0.3699	0.3726	0.3579
Intraday	MAE	10.96	10.79	10.37	9.88	10.05	8.93
	MSE	182.86	158.64	158.88	162.55	165.41	130.79
	MASE	0.571	0.562	0.540	0.515	0.524	0.456

Table 3. Table 3: Percentage improvement in the methods when EMD is used

Datasets	Methods Metrics	ANN (%)	ARIMA (%)	Zhang’s Method (%)	Khashei-Bijari’s Method (%)	Babu-Reddy’s Method (%)	Proposed Method (%)
Sunspot	MAE	41.1	42.2	43.2	26.9	30.4	30.5
	MSE	65.9	67.7	68.8	51.6	58.0	54.7
	MASE	41.1	42.2	43.2	26.9	30.4	30.5
Lynx	MAE	26.9	35.5	25.1	23.7	28.4	28.4
	MSE	45.6	57.1	42.7	43.4	46.5	43.2
	MASE	26.9	35.5	25.1	23.7	28.4	28.4
Gbp/Usd	MAE	55.6	66.4	62.3	64.0	66.2	65.0
	MSE	83.8	90.0	90.1	88.7	89.7	89.1
	MASE	55.6	66.4	62.3	64.0	66.2	65.0
Intraday	MAE	45.4	46.6	44.4	50.0	48.4	52.5
	MSE	70.3	75.6	73.2	72.9	73.3	77.5
	MASE	45.4	46.6	44.4	50.0	48.4	52.5

Equations

$$\begin{array}{r}{y}_t={\alpha }_1{y}_{t-1}+{\alpha }_2{y}_{t-2}+...+{\alpha }_p{y}_{t-p}\\ +{ϵ}_t-{\theta }_1{ϵ}_{t-1}-{\theta }_2{ϵ}_{t-2}-...-{\theta }_q{ϵ}_{t-q}\end{array}$$

$$y_t=w_0+\sum _{j=1}^H{w}_{j}f(w_{0j}+\sum _{i=1}^N{w}_{ij}{y}_{t-i})+e_t$$

$$y_t=L_t+N_t$$

$$y_t=f(L_t,N_t)$$

$$l_t=\frac{1}{m}\sum _{i=t-m+1}^t{y}_i$$

$$r_t=y_t-l_t$$

$$\widehat{L_t}=g(l_{t-1},l_{t-2},...,l_{t-p},{ϵ}_{t-1},{ϵ}_{t-2},...,{ϵ}_{t-q})$$

$$\begin{array}{rl}{S}_t^1& =(y_{t-1},y_{t-2},...,y_{t-a})\\ S_t^2& =(r_{t-1},r_{t-2},...,r_{t-b})\\ \widehat{y_t}& =f(S_t^1,\widehat{L_t},S_t^2)\\ \widehat{y_t}& =f(y_{t-1},y_{t-2},\dots ,y_{t-a},\widehat{L_t},r_{t-1},r_{t-2},...,r_{t-b})\end{array}$$

$$\text{centering}\begin{array}{rl}\text{MAE}& =\frac{1}{n}\sum _{t=1}^n|e_t|\\ \text{MSE}& =\frac{1}{n}\sum _{t=1}^n{e}_t^2\\ \text{MASE}& =\frac{n-1}{n}\frac{{\sum }_{t=1}^n|e_t|}{{\sum }_{t=2}^n|y_t-y_{t-1}|}\end{array}\text{@add@centering}$$

$$y(t)=\sum _{i=1}^{n}IM{F}_i(t)+R_n(t)$$

Full text

1 Introduction and related work

Time series include data points listed in time order. It is generally a sequence of discrete-time data which consists of points equally spaced in time. In time series forecasting, we try to predict the future points by analyzing observed points in the series. It has been widely used in various applications of science, engineering, and business fields. However, time series data might show different characteristics and show increasing or decreasing trends. Some time series has seasonal trends in which variations are specific to a particular time range, e.g coat and boot sales increase in winter season while decrease in summer season. On the other hand, some time series are not seasonal, such as stock market data. Moreover, time series data might show different level of volatility. While USD/EUR (United States dollar/Euro) exchange rate shows high volatility, growth of an animal, plant, or human being show a linear change. Over past several decades, a considerable effort has been devoted to develop and improve time series forecasting models [1]. In the literature, various forecasting methods have been proposed which use linear and nonlinear models separately or combination of both. In this paper, we propose a hybrid algorithm of linear and nonlinear methods where we choose Autoregressive Integrated Moving Average (ARIMA) as a linear method and Artificial Neural Networks (ANNs) as a nonlinear method.

ARIMA is widely used linear time series forecasting method that is used in numerous applications including finance [2], engineering [3], social sciences [4], and agriculture [5]. ARIMA models are integration of Autoregressive models (AR) and Moving Average models (MA). The model building process of ARIMA depends on Box-Jenkins methodology [6] where it provides a step-by-step procedure for AR and MA analysis. ARIMA models give good accuracy in forecasting relatively stationary time series data. However it makes a strong assumption that the future data values are linearly dependent on the current and past data values. In addition, it also assumes that there is no missing data in the given time series [7]. However, many real world time series data presents complex nonlinear patterns which might not be modeled by ARIMA effectively.

For the nonlinear time series modeling, Artificial Neural Networks (ANNs) are one of the most widely used algorithms [8] in many fields, such as finance  [9], energy [10], hydrology [11], and network communications [12]. ANNs have several advantages over ARIMA and other forecasting models. Firstly, ANNs are capable of fitting a complex nonlinear function, therefore it does not need to make data stationary as ARIMA does. This ability helps ANNs to approximate any continuous measurable function with arbitrarily desired accuracy [13], [14]. Moreover, ANNs are adaptively data-driven in nature which means ANN models can be adaptively formed based on the features of time series data [15].

In the literature, there are studies which show the success of linear and nonlinear methods over each other. For example, [16, 17, 18] report that statistical and linear models give better results than ANNs. On the other hand, [19, 20, 21] report that ANN performs better than linear models when data exhibit high volatility and multicollinearity.

In short, each type of model outperform in different domains. It is apparent that there is no universal model which is suitable for all circumstances. In order to overcome this limitation, various hybrid techniques have been proposed which aim to take advantage of the unique strength of each different type of models.

The common practice in hybrid techniques is to decompose time series data into its linear and nonlinear forms, then use appropriate type of models on them separately. A hybrid ANN-ARIMA model proposed by Zhang [15] achieves more accurate forecast results in time series data as compared to using individual models, in applications such as electricity price forecasting [22] and water quality prediction [23]. Another successful hybrid ARIMA-ANN technique is presented by Khashei and Bijari [24] which defines functional relationship between components. Moreover, Babu and Reddy [25] offers a solution to volatility problem in time series data by smoothing out dataset with moving average filter.

Each hybrid method in the literature bring different perspectives to time series forecasting problem. However strong assumptions that these methods make might degenerate their performances if the contrary situations occur. [15] and [24] do not decompose time series data into linear and nonlinear components. Rather, they assume that linear component of the data is the output of the ARIMA model. Moreover, [15] and [25] assume that the output of their hybrid methods is linear combination of components. However, different datasets can suggest different type of relationships between output and the components.

In this study, we propose a novel hybrid method for time series forecasting which aims to overcome the limitations of the traditional hybrid methods by eliminating the need to make strong assumptions. In this method, nature of nonlinearity is first characterized by the help of moving-average (MA) filter, then ARIMA is applied to the linear component. In the final step, ANN is used to combine the output of ARIMA, the nonlinear component, and the original data. Three benchmark datasets, the Wolf’s sunspot data, the Canadian lynx data, and the British pound/US dollar exchange rate data and an additional public dataset, Turkey Intraday Electricity Market Price data are used in order to show the effectiveness of the proposed method in time series forecasting.

In addition, we propose an improvement to all hybrid methods mentioned in this paper, including ours, by adding Empirical mode decomposition (EMD) technique [26] to the models. When accuracy results of the hybrid methods with different types of datasets are compared, we observe that accuracy performance gets better with the increasing level of linearity in time series. Then, time series data can be considered as a merge of sub-series which each of them demonstrates more linearity. For this purpose, a well-known multiscale decomposition technique, EMD is used in the proposed hybrid method. The components achieved by EMD are relatively stationary and have simpler frequency range which make them strongly correlated in themselves. Thus, more accurate predictions can be obtained through the models [27, 28].

The rest of the paper is organized as follows. In the next section, we present existing different time series forecasting methods. In section 3 we present our proposed model. In Section 4, we show the evaluation results of our model and present comparison results with the other methods. In section 5, we present the positive effect of using EMD technique with the hybrid methods and our new method. We show that our method consistently outperforms other linear, nonlinear and even hybrid methods.

2 Time series forecasting methods

In order to give overall review before we present our method, we want to give a brief information about other time series forecasting methods, such as ARIMA, ANN, and the well-known hybrid methods.

2.1 Autoregressive Integrated Moving Average Method

ARIMA is a linear method which means future value of a variable to be forecasted is assumed to be linear function of past observations. As a consequence, time series data that is fed to ARIMA is expected to be linear and stationary.

The method is composed of three main steps in order: Identification, estimation parameters and forecasting steps. In the identification step, stationarity check is performed on given time series data. A stationary time series’ statistical characteristics such as mean and autocorrelation patterns constant over time. When a trend or heteroscedasticity is observed on a time serie data, differencing or power transformation methods are generally applied to remove the trend and make variance stabilize. If stationary property is not satisfied after first attempt, differencing (or power transformation) method is performed continuously until non-stationary property is disposed. If this operation is performed $d$ times, the integration order of ARIMA model is set to be $d$. Thereafter, an autoregressive moving average(ARMA) is applied on the resultant data as follows:

Let the actual data value is $y_{t}$ and random error $\epsilon_{t}$ at any given time $t$. This actual value $y_{t}$ is considered as a linear function of the past $p$ observation values, say $y_{t-1},y_{t-2},...,y_{t-p}$ and $q$ random erros, say $\epsilon_{t},\epsilon_{t-1},...,\epsilon_{t-q}$. The corresponding ARMA equation is given in the following equation:

[TABLE]

In Equation 1, the coefficients from $\alpha_{1}$ to $\alpha_{p}$ are Autoregression coefficients, $\theta_{1}$ to $\theta_{q}$ are Moving Average coefficients. Note that random errors $\epsilon_{t}$ are identically distributed with a mean of zero and a constant variance.

Similar to the $d$ parameter, $p$ and $q$ coefficients are referred to as the orders of the model. When $q$ equals to zero, the model is reduced to AR model of order $p$. If $p$ is equal to zero, the model becomes MA model of order $q$. The main issue in ARIMA modeling is to determine the appropriate model orders $(p,d,q)$. In order to estimate order of the ARIMA model, Box and Jenkins [6] proposed to use correlation analyses tools, such as the autocorrelation function (ACF) and the partial autocorrelation function (PACF). When model coefficient estimation is finalized, the future values of the time series data are forecasted using available past data values and estimated model coefficients.

2.2 Artificial Neural Networks Method

ANN provides flexible computation framework for nonlinear modeling in wide range of applications. Due to its flexible architecture, number of layers and the neurons at each layer can be easily varied. In addition, ANN does not require any prior assumption, such as data stationarity, in model building process. Therefore, the network model is largely determined by the characteristics of the data.

The architecture of the most widely used ANN model in time series forecasting, which is also called as multilayer perceptrons, contains three-layers. The neurons of the processing units are acyclically linked. In order to model time series data using such a network, nonlinear function $f$ of $y_{t}$ sequence from $y_{t-1}$ to $y_{t-N}$ is constructed as shown in the following equation:

[TABLE]

where, at any given time $t$, $w_{ij}$ and $w_{j}$ are model weights and $H$ and $N$ are the number of hidden and input nodes, respectively. In this equation, $e_{t}$ corresponds to a noise or error term. The transfer function of the hidden layers $f$ in ANN architecture is generally a sigmoid function.

The power of ANN comes from its flexibility to approximate any continuous function by changing the number of layers $N$ and hidden nodes $H$. The choice of number of layers and the nodes at each of them play important role in ANNs’ forecasting performance. Large numbers of $N$ and $H$ can give very high training accuracies, but since it tends to memorize the training data, it suffers from overfitting. On the other hand, a too simple network of ANN leads to poor generalization. Unfortunately, there is no systematic set of rules to decide the value of these parameters. Thus, extensive number of experiments are required to tune functions and the parameters.

2.3 Zhang’s Hybrid Method

For time series forecasting, Zhang proposed a hybrid ARIMA-ANN model[15]. According to this model, it is assumed that time series data is a sum of linear and nonlinear components, given in the form of:

[TABLE]

where $L_{t}$ denotes the linear and $N_{t}$ denotes the nonlinear component.Firstly, ARIMA is used with the given time series data and linear forecasts are obtained. Residuals from linear component is assumed to contain only nonlinear relationship. This method uses ARIMA to make forecast from the linear component and ANN from the nonlinear component. Then, these models are combined to improve overall forecasting performance. This method gives better foreceasting accuracy than using ARIMA and ANN methods individually, as seen in the experimental results in three well-known real data sets - the Wolf’s sunspot data, the Canadian lynx data, and the British pound/US dollar exchange rate data.

2.4 Khashei and Bijari’s Hybrid Method

For time series forecasting, Khashei and Bijari proposed another hybrid ARIMA-ANN [24]. Similar to Zhang’s model, this model also assumes that any time series data is composed of linear and nonlinear components. Likewise, ARIMA is used to extract linear component and make forecast on it and residuals, which are nonlinear components, are fed into ANN along with the original data, and linear forecast of ARIMA output. The difference from the Zhang’s model is to avoid the assumption that the relationship between linear and nonlinear components is additive. Rather, this method builds functional relationship between the components as shown in the following equation:

[TABLE]

where $L_{t}$ is the linear and $N_{t}$ is the nonlinear component.

In addition, one may not guarantee that the residuals of the linear component may comprise valid nonlinear patterns. Khashei and Bijari suggest that residuals should not put into ANN as an input alone.

2.5 Babu and Reddy’s Hybrid Method

The hybrid ARIMA-ANN method proposed by Babu and Reddy [25] integrates moving average filter into hybrid ANN-ARIMA model. Like other methods, this model also assumes that any time series data is composed of linear and nonlinear components. However, this study emphasizes that neither Zhang nor Khashei and Bijari’s methods decompose original time series data into its linear and nonlinear components; instead, they use a linear ARIMA model to extract the linear component and the error sequences is assumed to be nonlinear component. On the other hand, this study separates the linear and nonlinear components, then feed them into appropriate methods.

This method tries to fix Moving Average (MA) filter length until kurtosis value of the data becomes approximately 3. Kurtosis is a shape of a probability distribution which measure thickness or heaviness of the tails of a distribution. The kurtosis value is 3 if the data has normal distribution. Shortly, the method aims to find out normal distributed component, which shows low volatility, in time series data by using kurtosis value. When the low-volatile component is separated from the original data, high-volatile component, which is assumed to be nonlinear, is achieved. In the final step, like in Zhang’s method, the decomposed components are fed into ARIMA and ANN accordingly and forecast results are summed up to achieve final forecast.

3 Proposed Method

Many decision processes need high forecasting accuracies in time series applications. Although there are numerous available time series models, none of them consistently gives the best results in various situation. There are two main challenges for making an accurate forecast. The first challenge is that underlying data generating process of time series cannot easily identified [29]. The second one is that non-hybrid individual models are generally insufficient to determine all the characteristics of the time series [15]. Many researches in time series forecasting literature show that hybrid models improve the forecasting performances [30]. By taking the advantage of each individual method in a combined model, error risk of using an inappropriate method is reduced and more accurate results are obtained.

Each hybrid method mentioned in this paper bring different perspectives into time series forecasting problem. However, strong assumptions these methods make might degenerate their performances in certain circumstances. For example, Zhang’s and Babu-Reddy’s methods assume that the relationship between linear and nonlinear components is additive. If linear and nonlinear components are not additively associated and the relation is different (i.e., it can be multiplicative), the possible complex relationship between components might be overlooked and the forecasting performance might be degenerated. Another assumption that might not always hold is that residuals might not show nonlinear pattern property. Additionally, Zhang’s and Khashei-Bijari’s methods do not actually decompose data into linear and nonlinear components, but they assume that linear component can be extracted by ARIMA and error sequences shows nonlinear pattern. As a result; such assumptions may lead to low forecasting performances when unexpected scenarios happen.

In this study, we propose to novel hybrid method for time series forecasting which aims to overcome the limitations of traditional hybrid methods by eliminating strong assumptions. The architecture of the proposed hybrid method is shown in Figure 1.

The algorithm starts with data decomposition. In this method time series data $y_{t}$ is considered as a function of linear $L_{t}$ and nonlinear $N_{t}$ components in the same way as given in Equation 4.

These two components are separated from the original data by using moving average (MA) filter with the length of $m$, as given in Equation 5. While the linear component $l_{t}$ has low volatility, the residual $r_{t}$, which is the difference between the original data and the decomposed linear data in Equation 6, shows high fluctuation.

[TABLE]

[TABLE]

In order for a proper decomposition, the length of the MA filter $m$ has to be adjusted. Augmented Dickey Fuller (ADF) test which is unit root test can be performed to determine whether a given data series is stationary or not. The existence of a unit root on a given dataset indicates that there is an unpredictable systematic pattern. The more negative ADF test result means the stronger rejection of the existence of unit root for a given time series. Therefore, a negative ADF result implies that the given dataset is stationary. The well accepted threshold is 0.05 which is also used in this study to adjust MA filter length.

After the linear component is achieved with MA filter, a linear model is constructed as shown in Equation 7. The stationary component $l$ is modeled as a linear function of past values of the data series $l_{t-1},l_{t-2},...,l_{t-p}$ and random error series $\epsilon_{t-1},\epsilon_{t-2},...,\epsilon_{t-q}$ in Equation 1 using ARIMA model.

[TABLE]

where $g$ is a linear function of ARIMA.

Finally, nonlinear modeling ANN is used to implement functional relationship between components as indicated in Equation 4. The past observed data $y_{t-1},y_{t-2},...,y_{t-a}$, present ARIMA forecast result of the decomposed stationary data $\hat{L_{t}}$, and residuals of the data decomposition $r_{t-1},r_{t-2},...,r_{t-b}$ are fed to ANN as indicated in Equation 8:

[TABLE]

where $f$ is the nonlinear function of ANN, $a$ and $b$ are parameters of the model which show how much we will go back in time to use as features to ANN. Time series data determines how many of those features in the residual path and observed data path are going to be used in the nonlinear model. For example, if the given data does not show volatility, then the residual variable $b$ in Equation 8 might come out even as zero in tuning process. Likewise, $a$ variable in Equation 8 is also empirically determined in the tuning process.

The proposed model does not only exploit the unique strength of single models, but also eliminates the three strong assumptions performed by other hybrid methods. Therefore, risk of low forecasting performance in unexpected situations is highly avoided. The competitive performance of our proposed algorithm is shown in our experimental results by using various type of datasets.

4 Empirical Results

The performance results of the proposed hybrid method along with the other methods discussed in this paper are evaluated on four different datasets. Three of them are well-known benchmark datasets - the Wolf’s sunspot data, the Canadian lynx data, and the British pound/US dollar exchange rate data - which have been widely used in statistics and the neural network literature [15, 24, 25, 27]. The other dataset is publicly available electricity price of Turkey Intraday Market [31].

In the experiments, only one-step-ahead forecasting is considered. In order to compare accuracy performances, three evaluation metrics are used: Mean Absolute Error (MAE), Mean Squared Error (MSE) and Mean Absolute Scaled Error (MASE) whose formulations are indicated as follows respectively:

[TABLE]

where $e_{t}=y_{t}-\hat{y}_{t}$ and $y_{t}$ is the actual data value, $\hat{y}_{t}$ is the forecasted value at given time $t$. While MAE specifies the average of the absolute errors over the performed prediction, MSE measures the average of the squared error. Since both MAE and MSE results depend on the scale of the given data, when comparing time series which have different scales, they are not preferable. Therefore, scale-free error metric MASE can be used to compare forecast accuracy between series.

To assess the forecasting performance of the different methods, each dataset is divided into training and testing sets. While the training data is used for model development, the test data is used to evaluate the established model. In order to tune the hyper-parameters the methods, the last 20% of training dataset is used as validation set.

In addition, due to the fact that ANN performs random initialization and produces different results at each run, the methods which include ANN algorithm are executed 50 times and average results are reported. Table 1 gives the forecasting results of all examined methods on these all four datasets.

4.1 Forecasts for Sunspot dataset

The Wolf’s sunspot series, which contains annual activity of spots visible on the face of the sun, has been extensively used in numerous linear and nonlinear models [24]. The data includes the annual count of sunspots from 1700 to 1987 (see Figure 2) giving a total of 288 observations. ADF stationarity test result of the dataset is 0.083 which is greater than the threshold 0.05. This implies that there is a unit root on the dataset, thus the dataset can be regarded as non-stationary time series. 288 observations in the dataset is divided into two samples: 221 observations between 1700-1920 years are considered as training data to develop the model, the last 67 observations between 1921-1987 years are considered as test data and used to evaluate the model performance.

In Sunspot dataset, when ARIMA is individually used as a forecasting method, we also set the order of ARIMA to 9 (AR(9)) as same as the other many studies [15, 24, 25]. When ANN is individually used as a forecasting method, similar to these studies, three layered 4 x 4 x 1 ANN architecture is used which is composed of four input nodes, four hidden nodes, and one output node.

In the proposed method, the linear component comes out when the MA filter length is 15. After using the filter, the achieved linear component has 0.006 stationary test result which indicates its stationarity, since it is a value less than the threshold, 0.05. The best fitted neural network in the final step of the proposed hybrid method has 7 nodes in the input layer where 4 of them are observed values, 2 of them are residuals, and one node is assigned for the result of linear component forecast. According to our experiments, when the number of hidden nodes are adjusted to same number of nodes as in the input layer, the best fitted ANN model is achieved.

When numerical results of Sunspot dataset given in Table 1 are analyzed, individual methods such as ARIMA, ANN have apparently lowest performance as compared to hybrid methods. This suggests that either ARIMA or ANN, when individually used, do not capture all patterns in the data series. Therefore, combining two methods by taking advantage of each of them can be an effective way to overcome this limitation. We indeed observe in Table 1 that the hybrid methods Zhang, Khashei-Bijari, and Babu-Reddy methods produce better results as compared to individual ones. However, they produce lower forecasting performance than our proposed hybrid method. The assumptions those hybrid methods make can be restricting in many situations as mentioned in Section 3. Our proposed hybrid method eliminates those assumptions and yields better generalization performance.

4.2 Forecasts for Lynx dataset

The lynx dataset, which contains the number of lynx trapped per year in the Mackenzie River district of Northern Canada, is an another extensively analyzed time series data in the literature [15, 24]. The data shows a periodicity of approximately 10 years as seen in Figure 4. Moreover, ADF stationary test results of the dataset is 0.056 which implies that the dataset is almost stationary. There are 114 observations in the data, corresponding to the period of 1821-1934. The first 100 observations between 1821-1920 years are considered as training data to develop the model, the last 14 observations between 1921-1934 years are used as test data to evaluate the model performance. In addition, like in other studies [15, 24], the logarithms (to the base 10) of the data are used in the analysis.

When ARIMA is used as an individual model, we used the AR model of order 12 (AR(12)) for Lynx dataset which is also used by [15, 24]. Similar to these studies, three layered 7 X 5 X 1 ANN architecture is used when ANN is individually used as a forecasting method.

In the proposed method, the linear component is extracted from the Lynx dataset when the MA filter length is 5. The relatively short MA filter length was expected, since the ADF test result shows a certain level of stationarity in the data. As a result of MA filter, the achieved linear component has 0.006 stationary test result which indicates even more stationarity to be properly modeled by ARIMA. The best fitted neural network in the final step of the proposed hybrid method has 9 nodes in the input layer where 5 of them are observed values, 3 of them are residuals, and one node is assigned for the result of linear component forecast. According to our tuning experiments, when the number of hidden nodes are adjusted to the same number of nodes as in the input layer, the best fitted ANN model is achieved.

In this dataset, among the individually used methods, ARIMA gives better accuracy as compared to ANN in contrast to the Sunspot dataset (see in Table 1). This is most likely due to the fact that Lynx dataset is more stationary dataset compared to Sunspot dataset. Due to this relative stationarity, the effect of hybrid methods might not be easily observed. Although we have circumstances which do not necessarily favor the data decomposition and model combination, hybrid methods do not give a lower performance than the individual linear models, and they even provide better results. The highest performance is mostly achieved by our proposed hybrid method. Figure 5 compares the actual and forecast values for all examined methods.

4.3 Forecasts for Gbp/Usd dataset

The other benchmark dataset is the exchange rate between British pound and US dollar which contains weekly observation from 1980 to 1993, giving 731 data points in the time series. Predicting exchange rate is an important yet difficult task due to high volatility. ADF stationary test result of the dataset is 0.58 which is highly greater than the threshold 0.05. This implies that the dataset is highly volatile and non-stationary. This non-stationarity can be even seen in the plot, given in Figure 6, which shows numerous changing turning points in the series. Similar to other datasets, the experimental setup is same as in previous hybrid studies [15, 24] where data is transformed using natural logarithmic function and separated into two samples. The first 679 observations from 1980-1992 years are considered as training data to develop the model, the last 52 observations between 1992-1993 years are used as test data to evaluate the model performance.

In this dataset, when ARIMA is individually used, rather than using regression type of model in the ARIMA itself, random walk model is chosen as best-fitted ARIMA model. This approach has been used by Zhang [15] and also been suggested by many studies in the exchange rate literature [32]. In this model, the most recent observation is the best guide for the next forecast. When ANN is individually used as a forecasting method, the best fitted ANN is set as three layered 7 x 6 x 1 architecture.

In the proposed method, in order to decompose this highly volatile data, MA filter length comes out to be 40. As a result of this decomposition, ADF test result of the obtained component is 0.007 which indicates the stationarity of the component. To compute the best final forecast in the proposed model, the ANN is constructed as three layered 9 x 9 x 1 architecture. In this architecture, input layer is composed of the last 5 of observed values, the last 3 of residuals, and the result of linear component forecast.

Results of the Gbp/Usd dataset forecasts are compared in Figure 7. Both ANN and hybrid methods have much better performance than the individual ARIMA method for a highly fluctuating forecast horizon. The proposed hybrid method is able to capture this volatile pattern much better and outperforms the other methods in all error metrics.

4.4 Forecasts for Intraday dataset

The last analyzed dataset is Turkey intraday electricity market price data which is publicly available [31]. The dataset contains 581 observations which consist of daily averaged prices from July 2015 to December 2017. As compared to datasets such as Sunspot and Lynx, data pattern in this dataset is also highly fluctuating (see Figure 8). ADF stationary test result of the dataset shows 0.27 value which is highly greater than the threshold of 0.05. Natural logarithmic transformation is applied on the dataset for scaling purposes.

In this dataset, when ARIMA is individually used as a forecasting method, we found autoregressive model of order 9 (AR(9)) to be the most parsimonious one among all ARIMA models. The best fitted individual ANN model is achieved in three layered 3 X 6 X 1 architecture, after our tuning process.

In this dataset, the length of MA filter which we use in the proposed hybrid method comes out as 6. The obtained linear component after the MA filter has 0.004 ADF stationary test result; this implies that the component can be properly modeled by ARIMA method in our architecture. The best fitted ANN in the final step of the proposed method has 17 x 17 x 1 architecture. Input layer is composed of the last 8 of original data, the last 8 of residuals and the result of linear component forecast.

Similar to Gbp/Usd dataset, Intraday dataset is highly volatile and non-stationary which cannot be effectively modeled by using only a linear model. As can be seen in Table 1, ANN and hybrid methods significantly outperform the individual linear model of ARIMA. Furthermore, the proposed method gives remarkably superior accuracy as compared to other hybrid methods in all error metrics. The comparison of the actual and forecast values for all examined method are given in Figure 9.

5 Discussion and Improvement

There are several important results obtained in our experiments Firstly, when individual methods’ results are compared among themselves, we see that ARIMA outperforms ANN for the datasets which present more linearity and vice versa (see Table 1). Moreover, hybrid methods have better performance as compared to individual ones especially in more fluctuating datasets. Finally, the assumptions made by other hybrid methods degenerate the forecasting performance when unexpected situations occur in the data. Our proposed hybrid method which avoids these assumptions apparently creates more general models and outperforms the other examined methods.

Figure 10 compares the distribution of MASE results of datasets with ADF test results of the corresponding dataset. In this figure, boxplots are drawn by using MASE values of all examined methods for each dataset, presented in Table 1. When we compare error results among time series data by using a scale-invariant error metric MASE, it is observed that the more non-stationarity in a dataset leads to a higher error value. For example, Lynx dataset, which turns out to be the most linear among all datasets according to ADF test results, has the lowest MASE results. On the other hand, Gbp/Usd dataset, which shows the most non-linearity according to ADF test results, has the highest MASE results. As a result of this, we can conclude that having more regular data distribution in a time series leads to more accurate results in forecasting. This conclusion motivates us to propose an improvement on our already best performing proposed hybrid method. This improved method aims to produce more stationary subseries from given time series by using a multi-scale decomposition technique. Then, those achieved linear subseries can be modeled with a higher accuracy using the proposed hybrid method.

In the literature there are several multi-scale decomposition methods such as Empirical Mode Decomposition (EMD), Wavelet Packet Decomposition (WPD), Fourier Transform (FT) and etc. [33]. Since EMD does not make a priori assumption about the given time series and preserves time scale of the data throughout the decomposition, it is a more prefer technique than the others for decomposing time series [26].

The main principle of EMD is to decompose a given time series data into a sum of several subseries. Those subseries are called Intrinsic Mode Functions (IMFs) and the remaining component after subtracting the summation of IMFs from the original data is called residue. These subseries have two important properties which allow them to be easily modeled: Each subseries has its own local characteristic time scale and they are relatively stationary subseries. Let $y(t)$ be a given time series data, and then the EMD calculation can be described as follows:

[TABLE]

where $IMF_{i}(t)$ $(i=1,2,...,n)$ represents the different subseries, and $R_{n}(t)$ is the residue after summation of $n$ IMFs are subtracted from the original data.

The EMD-based methods includes three main steps, as seen in Figure 11: In the first step, the original time series data is decomposed into IMFs. In the second step, forecasting is performed by usinf our hybr’dd method on each IMF. In the last step, forecast results of each individual model are summed up to achieve the final forecast of the original time series. We use an additive function in the end to capture the additive relation of IMFs with the original data.

The EMD-based methods are evaluated on the same datasets by using the experimental setup in Figure 11. In order to evaluate the effect of EMD, all examined methods are executed in the second step of the algorithm (see Figure 11). Table 2 gives the forecasting results of all examined methods with EMD on four all datasets. When these results are compared with the previous ones showed in Table 1, the methods with EMD give significantly higher accuracies. The percentage improvement for all datasets are presented in Table 3. The improvements at each dataset varies between 23% and 89% for all error metrics. We also provide a bar chart (see Figure 12) that shows the MASE results for each method which are averaged over the results of all datasets. The chart indicates that methods with EMD achieve remarkably less error in their forecasts. Further analysis shows that EMD-based methods give even greater improvements in the accuracies for non-stationary datasets. This is most likely due to the fact that EMD is able to resolve high volatility problem in time series data. For example, while our hybrid method with EMD improves MASE results of relatively stationary datasets Sunspot and Lynx 30% and 28% respectively, this improvement jump over 50% in Gbp/Usd and Intraday datasets which are highly non-stationary. Another indication of EMD’s capability of solving volatility problem is that ARIMA is able to achieve better results than ANN for all datasets (see Table 3), although ANN without using EMD was better in non-stationary datasets (see Table 1).

In final, we want to point out that our proposed hybrid method with EMD gives the best results as compared to other methods (see Table 3). Subseries obtained from EMD are relatively stationary as compared to original data however, they still show fluctuations in their frequency range. Performing one more decomposition on these subseries using MA filter and constituting functional relation between the stationary part, residuals, and the original data values outperforms all other remaining methods.

6 Conclusions

Time series forecasting is an important yet often a challenging task used in many different application domains. The studies in the litreture mainly focus on either linear or nonlinear modelings individually or a combination of them. While linear models such as ARIMA gives better forecasting accuracy with stationary time series data, nonlinear methods such as ANN is more appropriate for non-stationary datasets. In order to take advantage of the unique strength of each different type of methods in a more general setting, hybrid methods are proposed. Hybrid methods basically use linear and nonlinear modeling, ARIMA and ANN respectively on the corresponding decomposed components and then combine the results. Hybrid ARIMA-ANN methods give better results in general as compared to the cases where they are individually used. However, they generally suffer from the assumptions they make while constructing their model. These assumptions lead to produce inconsistent results and give low accuracies in overall if unexpected situations occur.

In this study, a new hybrid ARIMA-ANN model based forecasting method is proposed to overcome three main assumptions made by traditional hybrid ARIMA-ANN models. Firstly, the proposed method removes the assumption that the linear component is the ARIMA model output of the given data. Rather, it extracts the linear component by using MA filter. It is known that data showing linear characteristics can be more accurately modeled by linear methods. Therefore, properly decomposed data yields more accurate linear forecasting and consequently more accurate final results in the hybrid methods. Secondly, the proposed method does not directly model residuals via a nonlinear method ANN, since the assumption of that residuals might comprise valid nonlinear patterns, does not always hold. Thirdly, the proposed method does not restrict linear and nonlinear component modeling and also combining the results of them. Rather, it can capture structures of the linear and nonlinear components in a better way, and produces more general models than those existing hybrid models.

In the light of our experimental results, we can conclude that forecasting performance gets better if more stationary time series data is provided. This result motives us to make original time series data more stationary in order to improve accuracy results. We show that when EMD multi-scale data decomposition is combined with all examined methods, accuracy results can be remarkably improved. Our experimental results indicate that our hybrid method with EMD gives remarkably superior accuracy as compared to all other examined methods.