# A New Approach for Reconstruction of IMFs of Decomposition and Ensemble Model for Forecasting Crude Oil Prices

Authors: Peng Xu; Muhammad Aamir; Ani Shabri; Muhammad Ishaq; Adnan Aslam; Li Li
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

Accurate forecasting for the crude oil price is important for government agencies, investors, and researchers. To cope with this issue, in this paper, a new paradigm is designed for the reconstruction of intrinsic mode functions (IMFs) of decomposition and ensemble models to reduce the complexity in computation and to enhance the forecasting accuracy. Decomposition and ensemble methodologies significantly enhance the forecasting accuracy under the framework of “divide and conquer” with the proposed reconstruction of IMFs method. The proposed approach used the autocorrelation at lag 1 of all IMFs for the reconstruction. The ensemble empirical mode decomposition (EEMD) technique is employed to decompose the data into different IMFs. Models that utilized the decomposed data relatively perform well, as compared to its application to the undecomposed data. However, sometimes, the decomposition may produce poor results due to the error accumulation at the end. Thus, in this study, the reconstruction of IMFs is proposed for minimizing the aforementioned error, thereby increasing the forecasting accuracy. The Brent and West Texas Intermediate (WTI) datasets (daily and weekly) are exploited to compare the forecasting performance of autoregressive integrated moving average (ARIMA) along with artificial neural network (ANN) models with the decomposed data. The results have proven that the new paradigm of reconstruction of IMFs through autocorrelation was a better and simple strategy that significantly improved the performance of single models including ARIMA and ANN. Hence, it is concluded that the proposed model takes less computational time and achieved higher forecasting accuracy with the reconstruction of IMFs as opposed to using all IMFs.


---

## Body

## 1. Introduction

Crude oil is a very important commodity in the world because of its unique nature, as it affects the life of every individual in many ways. According to the IEA, in early 2018, the world currently consumes 99.3 million barrels of oil and liquid fuels daily. In global market, it is the most active and heavily traded commodity. Due to the high demand of crude oil in every field of life, it needs more attention as compared to other commodities. Oil is a nonrenewable commodity, but the world consumes it in different ways; thus, it is a challenge for mathematicians, statisticians, and econometricians to develop a better strategy for understanding the price changing aspect of crude oil. Due to the irregular and stochastic nature of oil prices, it is a very complex and challenging task for researchers to develop an appropriate model for COPs forecasting. The compound and complex nature of the COPs make this area widely opened for researchers to develop many different procedures to forecast feasibly.In the existing literature, several models have been proposed for forecasting the time series data and applied for forecasting the COPs, which is generally described in [1, 2], yt+l=fYt+ϵt+l, where yt represents the price of crude oil at time t, Yt=yt−1,yt−2,…,yt−h is the previous value of time t with lagging of h, l is the horizon of the prediction, and ϵt+l denotes the prediction error which is independent and identically distributed random variable. The function f∗ is distinguishably designed with the respective parameter evaluation techniques, while the current COPs forecasting models can be divided broadly into three categories. In the first category, traditional econometric models fall with strict assumptions such as stationarity and linearity in parameter estimation but with relatively simple fixed functions. The second category is called AI, which consists of flexible functions having the capability of self-learning at the stage of model training. The last category is very popular these days as the combination of several single models into a hybrid model systematically. This study focuses on the third category of hybridizing the popular single models into one for achieving the high forecasting accuracy with a simple technique.The first category consists of ARIMA, RW, GARCH, VAR, and ECM models which are more frequently used to forecast the COPs. For the prediction of monthly COPs, Y. Xiang and X. H. Zhuang [3] used the ARIMA model for Brent data covering a period from November 2012 to April 2013. For estimation of the volatility and the conditional mean [4], GARCH model is used for daily COPs data of WTI covering a period from 9 December 2000 to 2 January 2010. The VAR model is utilized in [5] for US monthly COPs covering from January 1980 to November 2002. The RW model is also used as a benchmark to forecast COPs moments [6]. The ECM model was used in [7] for the prediction of COPs of Brent and WTI covering the period from 1994 to 2002.The second category concerned about the AI methods, the predominant techniques, which are frequently used to forecast the COPs, and empirical studies which proved that their superiority over the traditional econometrics models are SVR, ANN, and LSSVR. For SVR models, W. Xie et al. [8] used this technique for the prediction of WTI monthly COPs covering a period from January 1970 to December 2003, and the empirical study showed that the SVR model provides more accurate predictions as compared to the ARIMA model and BPNN models. Similarly, the SVR model was used in [9] for the prediction of COPs of WTI covering a period from January 1986 to December 2009. Movagharnejad et al. [10] used the ANN model for forecasting the COPs data covering a period from January 2000 to April 2010. The evolutionary neural network was presented by Chiroma et al. [11] for forecasting the WTI monthly COPs covering a period from May 1987 to December 2011. Li et al. [12] utilized the LSSVR model for forecasting WTI COPs from 4 January 2008 to 18 October 2018, with a conclusion that BPNN, ARIMA, and SVR models were outperformed by LSSVR model. The AI models performed better than the traditional models; however, overfitting and best parameter selection are the weakness of these models [13].The third category is about the hybridization of the models into a single model. In this category, one of the most famous concepts “divide and conquer” [14] is used to decompose the original series into meaningful components and then ensemble them after applying different techniques on each component. Nowadays, the decomposition and ensemble strategy was used frequently to forecast COPs. This methodology has three main steps: in the first step, the time series is divided into small independent components; in the second step, a suitable technique is applied to each component and predict; and the last step, ensembles the individual predictions for the original series prediction [13, 15]. Using the “divide and conquer” concept in the decomposition and ensemble modelling techniques, some recent studies have demonstrated the usefulness of these techniques which identifies the enhancement of the forecasting accuracies. Yu et al. [15] decomposed the original series by using the EMD and predicted WTI and Brent daily COPs covering periods from 20 May 1987 to 30 September 2006 and 1 January 1986 to 30 September 2006, respectively. Subsequently, a decision is reached, whereby the “divide and conquer” approach performed well and enhanced the forecasting accuracy for both WTI and Brent COPs series. The complementary ensemble empirical mode decomposition technique is introduced in [16] for forecasting the COPs, the empirical results of the decomposition, and ensemble strategy which is supported by the theory of improving the forecasting accuracy of the model. The same outcomes were achieved in [17] which suggested a new AI learning-paradigm with compressed sensing for forecasting the daily COPs covering a period from 3 January 2011 to 17 July 2013. In the EEMD and EELM paradigm which is introduced in [18] for the prediction of WTI COPs covering a period from 2 January 1986 to 21 October 2103, the results supported the decomposition and ensemble methodologies. In conclusion, the decomposition and ensemble methodologies worked well and significantly improved the performance of the models used for forecasting the COPs.Moreover, Lin and Sun [19] also proposed a hybrid method that is a combination of the CEEMDAN and MLGRU neural network to forecast the COPs and concluded that the new model enhanced the forecasting accuracy. Xian et al. [20] suggested a new model integrating the two methods called EEMD-ICA which has been proposed for financial time-series signal processing. Wu et al. [21] also proposed a novel model based on EEMD and LSTM. The proposed model EEMD-LSTM is applied to the COPs of WTI, and the numerical results have proved the superiority of the new technique.Furthermore, the complex data of COPs will be effectively handled by the decomposition ensemble models as compared to the ARIMA, GARCH, and ANN single models. However, an essential issue could arise regarding the model computational time, cost, and complexity. Due to the break down of decomposition ensemble models, the time series data are divided into several independent components called IMFs. After the decomposition, the modelling of all IMFs could be a time-consuming method due to the fact of individual IMF modelling and forecasting. Moreover, at the end, sometimes, produced poor results due to the estimation errors of all models are accumulated in the last ensemble step of forecasting [2, 22]. Addressing the problem of reconstruction of IMFs is incorporated after the decomposition step to forecast the individual IMF. During the IMFs reconstruction step, all the IMFs, after the decomposition of the original series, are reconstructed into some meaningful components. To overcome the above problem in decomposition and ensemble models, B. Zhu et al. [23] proposed a new multiscale paradigm through a kernel function incorporating EEMD, PSO, and LSSVM. Using the EEMD method, the original series was decomposed into IMFs and reconstructed into groups such as HFs, LFs, and trend. Later, the HFs were forecasted using the ARIMA model, while LFs and trend were forecasted by the PSO method, and for the final prediction, the whole forecasted results were simply combined. Aamir et al. [24] also reconstructed the IMFs of the EEMD model using the order of the ARIMA model for forecasting the weekly COPs of Brent and WTI. The authors decided that the reconstruction of IMFs was an effective technique that reduces time and cost and enhances the forecasting accuracy. Rios and De Mello [25] empirically proved that the acquired IMFs present various levels of stochasticity. This was the motivation to develop formal research to prove that the decomposed IMFs might be separated into two components, namely, stochastic and deterministic. The problem arises as how to divide these IMFs into SD components. The IMFs were divided in [25] into two components based on the RP and AMI and concluded that dividing the IMFs into two components enhanced the forecasting accuracy. Furthermore, the authors in [26] claimed that the IMFs could be divided into two additive components called stochastic and deterministic, and the researchers can work to better model these components. Aamir and Shabri [27] also reconstructed the IMFs in two components called stochastic and deterministic through AMI. The authors concluded that the reconstructed components enhanced the forecasting accuracy taken WTI and Brent daily COPs as a sample. Furthermore, Du et al. [28] also proposed the BPE index method to solve the problem based on PEEMD. The authors also used the reconstruction of IMFs procedure by implementing the T-test.Most of the above studies have used one or two criteria for reconstructing the IMFs into components, such as HFs, LFs, trend, RP, AMI, the order of ARIMA model,T-test, and average, while avoiding some other data characteristics. In order to capture some important components from the data dynamics efficiently, a comprehensive data analysis or some effective technique is strongly recommended for all the decomposed modes as well as a better procedure for forecasting reconstructed modes. In this study, the EEMD technique is used to decompose the time series for analyzing a wider class of time series.To overcome the above problem and further improve the forecasting efficiency, in this study, we aim at proposing a novel method of the use of decomposed IMFs obtained from EEMD. The innovation and contributions of this study are outlined as (i) a new forecasting approach for COPs of utilizing the reconstructed IMFs was proposed, following the well-known “decomposition and ensemble” framework. (ii) To the best of our knowledge, this is the first time that IMFs are reconstructed in this format and used for time-series forecasting. (iii) A threshold value of autocorrelation is fixed based on simulations and statistical tests for separating the IMFs into SD components. (iv) In publicly accessible Brent and WTI (daily and weekly) COPs, extensive experiments were conducted, and it was shown that the proposed approach outperformed several state-of-the-art methods for forecasting COPs. (v) We further analyzed the characteristics of the SC which needs more consideration and can help in improving the forecasting accuracy of the COPs significantly.In the following section, the framework of the proposed study, methods, and materials are outlined in detail.


## 2. Materials and Methods

### 2.1. Data Description

The data used in this study were obtained from DataStream and can be found in the supplementary files. A total of four datasets are used to evaluate the performance of the proposed procedure. All datasets consist of the daily and weekly COPs of Brent and WTI. There are a different number of observations in each dataset to check the correctness, robustness, and generalizability of the suggested procedure. Every time series was divided into training and testing sets of observations. The first 80% of the total observations in every time series were used as training set, whereas the rest 20% were used as the testing set. The details of every time series are discussed in brief in subsequent sections.


#### 2.1.1. Daily

The Brent COPs time series contain a total of 12436 observations (spanning 2 February 1970 to 1 October 2017), the first 9,949 observations belong to the training series, and the rest 2487 part of the testing series. The next is the WTI COPs time series containing a total of 9060 observations (spanning 12 January 1983 to 1 October 2017), the first 7248 observations belong to the training series, while the rest of 1812 part of the testing series.


#### 2.1.2. Weekly

The Brent COPs time series contain a total of 2488 observations (spanning 2 February 1970 to 1 October 2017), the first 1991 observations belong to the training series, and the rest 497 part of the testing series. The next is the WTI COPs time series containing a total of 1812 observations (spanning 12 January 1983 to 1 October 2017), the first 1450 observations belong to the training series, while the rest 362 part of the testing series.


### 2.2. Framework of the Proposed Study

Reconstruction of IMFs obtained from EEMD is proposed in this study. Firstly, the obtained IMFs from EEMD are divided into two components, namely, stochastic and deterministic. In this study, the proposed procedure for dividing IMFs into two components is autocorrelation (correlation (yt,yt−1)) of every IMF. If the autocorrelation is less than a threshold value, then that IMF will be considered stochastic, otherwise deterministic. To obtained two components add all stochastic IMFs for SC and deterministic for DC. The idea of reconstructing the IMFs into SD components should be utilized in more detail for further improvement because the SC requires more attention which can help in improving the forecasting accuracy significantly. The IMFs being part of the SC handled separately because most of the variation in total shared by the SC [29]. Thus, in this study, the stochastic IMFs modelled separately, for every stochastic IMF a different model is chosen and for final forecasting ensemble all the models’ forecast, whereas all deterministic IMFs would be treated as one component. The ARIMA and ANN models will be fitted for stochastic IMFs and DC and for the final output, add all these components, so the proposed models are denoted by EEMD-R-ARIMA and EEMD-R-ANN. For ANN, the AR terms of the ARIMA model would be used as inputs. The complete framework of the proposed method is demonstrated in Figure 1.Figure 1
The framework of the proposed method.


### 2.3. Experimental Study

It is discussed that the main aim of this study is to use the minimum number of IMFs with higher forecasting accuracy. For a minimum number of IMFs use, the reconstruction of IMFs is proposed in this study. Thus, first obtain all IMFs from EEMD and divide them into two components, namely, stochastic and deterministic. However, a procedure or a threshold value is required which decides that the given IMF is stochastic or deterministic. The proposed procedure is the autocorrelation of IMFs from which a decision will be taken that the given IMF has stochastic or deterministic influences. Based on four different threshold values, i.e.,ρyt,yt−1=0.500,0.900,0.950,0.990 of autocorrelation, the obtained IMFs would be divided into two components. The SC consist of all those IMFs whose autocorrelation is less than the threshold value; otherwise, the IMF belongs to the DC. To obtain the best threshold value for the division of IMFs into SD components, synthetic time series would be used which is characterized through the additive noise and organized into four different scenarios as follows:(1)
In a time series composed of a sine function and normal distribution, the sine function represents the DC, while the normal distribution represents SC:(1)Yt=sin2πt+ε0,1.(2)
In a time series composed of a sine function and ARMA model, the sine function represents the DC, while the ARMA model represents the SC with an error of 0.25:(2)Yt=sin2πt+ARMA1,1.(3)
In a time series composed of a sine function, ARMA model, and normal distribution, the sine function represents the DC, while the ARMA model and normal distribution represent SC:(3)Yt=sin2πt+ARMA1,1+ε0,1.(4)
In a time series composed of a Lorenz system and normal distribution, the Lorenz system represents the DC, while the normal distribution represents SC:(4)Yt=Lorenz system10,28,83+ε0,1.Therefore, to divide the IMFs into SD components, in this study, autocorrelation is used. IMFs have stochastic influences being part of the SC while having deterministic influences being part of the DC. To obtain the original series, add both SD components. The following formula is used to compute the autocorrelation of all decomposed IMFs:(5)rYt,Yt−1=∑t=1T−1Yt−Y¯Yt−1−Y¯∑t=1TYt−Y¯2.The problem is how to decide that the given IMF has a stochastic or deterministic influence. To fix this problem, an experimental study is performed to obtain a threshold value rather a confidence interval from which we will be able to detach the IMFs regarding their influences. When these IMFs are detached according to their influences, the next step is to obtain SD components by summing up all IMFs, respectively. A correlation test is performed to check whether the detached components are independent. For testing, the following test statistic is used with the alternative hypothesisHA:ρSD≠0:(6)t=rSDn−21−rSD2∼tn−2,where n denotes the total number of observations and rSD represents the correlation among the SD components and is defined as follows:(7)rSD=∑t=1TYS,t−Y¯SYD,t−Y¯D∑t=1TYS,t−Y¯S2∑t=1TYD,t−Y¯D2.To double-check the distribution of these two components, another assessment is also steered which can test the difference of the two correlations produced by the synthetic data SD components and the reconstructing SD components obtained from IMFs. For testing, the following test statistic is used with an alternative hypothesis that the difference between the two correlations is not equal to zero, i.e.,H0:ρSDs≠ρSDdorρSDs−ρSDd≠0:(8)z=z1−z2δz1−z2∼N0,1,where z1=1/2log1+rSDs/1−rSDs, z2=1/2log1+rSDd/1−rSDd, δz1−z2=2/T−3, “s” represents the synthetic data, and “d” is used for decomposed data. The test is performed for different threshold values of the autocorrelation which divided the IMFs into two components (SD). The best threshold value will be chosen for the real-life datasets and should implement.The above four different scenarios are used to check how well the proposed method extracted the SD components. The purpose of using different series is to check the generalizability and robustness of the proposed method. One of the advantages of the use of synthetic data is that it provided a controlled scenario and better supports the conclusions.


### 2.4. Methods Used in This Study

#### 2.4.1. The ARIMA Box–Jenkins Approach

The inclusion of both AR and MA terms makes the Box-Jenkins models exceptionally flexible. The interdependencies in time series (Yt) are measured by the AR terms, while the dependencies on preceding error terms are measured by the MA terms [30, 31]. An ARMA model of order (p, q) for a univariate series has the following form:(9)Yt=γ0+∑i=1pγiYt−i+∑j=1qβjεt−j+εt,where p and q represent the number of AR and MA terms, respectively, γ represents the constant and AR coefficients, whereas β is used for MA terms. The white noise process is denoted by ϵt∼N0,σ2. Using the Box–Jenkins procedure which needed stationary data, this can be achieved by differencing the time series many times (but at most two times usually). The ACF and PACF plots can be used to select the appropriate order of the polynomials and compare them with the theoretical bounds of [32]. Rival ARIMA models can also be selected based on AIC [33] and BIC [34]. The next is the diagnostic checking that the fitted model is adequate by using the LB test [35]. After obtaining the appropriate ARIMA model, the last step is the forecasting of the series.


#### 2.4.2. Artificial Neural Network

ANNs have many distinguishing characteristics over the ARIMA models and very successful alternative to linear models to handle the nonlinearity factors when forecasting the time series. One of the best characteristics of an ANN is the universal approximation of any nonlinear function up to the best degree of accuracy of your desired [36, 37]. A most commonly used ANN is the FNN with a single hidden layer having one output node for forecasting the time series applications [38, 39]. An m×n×1 ANN model has the following form:(10)Yt=μ0+∑j=1nμjgϕ0j+∑i=1mϕijYt−i+εt,where μjj=0,1,2,…,n, ϕiji=1,2,…,m,j=1,2,…,n denotes the weights, μ0andϕ0j represent the bias terms, ϵt is used for the white noise, Yt−ii=1,2,…,m represents the input patterns, Yt represents the single output of the model, m is used for the number of input variables, n is used for the number of hidden nodes, and g⋅ is used for the transfer function of the hidden layer. Equation (11) performs nonlinear mapping from past observations Yt−1,Yt−2,…,Yt−m, for the future value Yt. i.e.,(11)Yt=ϑYt−1,Yt−2,…,Yt−m,ω+εt,where ϑ is a function determined by the neural network training and ω is a vector representing all parameters. Therefore, from some viewpoint, the FNN model looks like or equivalent to an NAR model [40, 41]. The ANN effectively modelled the nonlinear time series, but for linear problems, it provided mixed results [38].


#### 2.4.3. EEMD

The new opportunities are provided in [42] to extract the time opposing elements from data after the advancement of transient local and versatile strategy for EMD. The transient region is a standout amongst the most critical qualities of EMD, which is procured from the spline fitting through maxima/minima of inputted data. The spline fitting has a high transient area and can be easily affirmed using numerical programming. The transient domain of EMD is all around spared if the quantity of filtering is smaller and fixed. The filtering process of EMD is worse with the excessive use of both transient region and oscillatory components. One of the benefits of the EMD procedure is that it automatically bypasses the stationarity assumption of the data. For the physical interpretability of the results, EMD territory gives a significant condition but not adequate. An EMD approach involves continuing with the levels of removing oscillatory sections of lesser replication, and the noise-wrapped results at one level can provoke continuing with the twisting of subsequent oscillatory components, making EMD happening hardly physically interpretable. Hence, due to the absence of this robustness problem, the EMD procedure becomes inadequate.Several attempts have been made to overcome the problem of robustness. For example, Huang et al. [42] proposed the discontinuity test in which variation is controlled by one element. Due to the controlled scenarios, it reduced the adaptiveness properties of EMD, demanding for more fascinating methodologies. To fulfill this robustness problem, a noise-assisted data analysis procedure, namely, EEMD was introduced in [43] which exactly solved the highlighted issue raised in EMD. The following are the steps involved in EEMD:(1)
A series of unique white noise summed with true seriesYt(2)
In the next step, the new series obtained in (9) can be decomposed accordingly into IMFs(3)
Repeat (9) and (10) by summing the different white noise series to the true series Yt(4)
As the outcome attained the ensemble means of the respective IMFsThe EEMD utilized the properties of white noise in the above steps and are as follows: the minima and maxima distribution should be uniform transiently and EMD being adequately a dual-channel bank for white-noise [44–46]. The minima and maxima distribution on all-time spans are relatively uniform with the added white noise. Another property of EEMD is the control on oscillatory elements after adding the white noise to the series which works like a dual filter bank which significantly reduces the mode mixing problem. Thus, as an outcome, the decomposition turns out to be steadier and physically significant.


#### 2.4.4. EEMD-ARIMA and ANN Model

Wang et al. [47] summarized the hybrid model of EEMD-ARIMA by incorporating the steps as follows:(1)
Firstly, the true seriesYt be decomposed into different IMFs(2)
Secondly, for everykth IMF the best ARIMA model will be selected and forecast accordingly(3)
The last step is to obtain the targeted series summed all thekth IMFs outputThe EEMD-ANN model should be fitted in the same way above. Furthermore, for every ANN model, the delaying parameter or starting input parameter be the AR terms of every ARIMA model ofkth IMF.


### 2.5. Forecast Evaluation Criterions

Forecasting accuracy is used to measure the degree of exactness among the proposed and state-of-the-art methods. When there are competing models, the forecasting accuracy measures are the most important criteria. In this study, four-typical forecasting accuracy measures are used consisting of MAE, RMSE, MAPE, and DS. Typically, these methods are defined as follows:(1)
MAE(12)MAE=1F∑t=1FY^t−Yt.(2)
RMSE(13)RMSE=1F∑t=1FY^t−Yt2.(3)
MAPE(14)MAPE=1F∑t=1FY^t−YtYt×100,whereF represents the total number of forecasts of the testing set, the original value of the given time series at time t is denoted by Yt and forecasted value is denoted by Y^t. Four different datasets are used in this study for validation and justification purposes. The next is the directional prediction measure which is more appropriate for business purposes since the investors are concerned with the market trend. Thus, DS is used for competing models to measure the directional forecasting accuracy and defined as follows.(4)
DS(15)DS=1F∑t=1FDt×100,whereDt=1, if Yt+1−Yt×Y^t+1−Yt≥0, else Dt=0, and F corresponds to the number of total forecasts of testing set. Apparently, the lower the values of MAPE, RMSE, and MAE, the higher the forecast accuracy in contrast with higher the values of DS produced more accurate forecasts. Next, to compare the two models prediction errors, we employed the DM test which is defined as follows.(5)
DM TEST(16)DM=d¯1/nγ0+2∑k=1nγk.It is used to test the statistical significance of the forecast accuracy of two forecast methods. Variabledt is the subtraction of absolute errors of the two methods. d¯=1/n∑t=1ndt is the mean of dt. γk=Covdt,dt−k is the autocovariance at lag k. DM∼N(0, 1), if the p value > α, we conclude there is no significant difference between the two forecasts.


## 2.1. Data Description

The data used in this study were obtained from DataStream and can be found in the supplementary files. A total of four datasets are used to evaluate the performance of the proposed procedure. All datasets consist of the daily and weekly COPs of Brent and WTI. There are a different number of observations in each dataset to check the correctness, robustness, and generalizability of the suggested procedure. Every time series was divided into training and testing sets of observations. The first 80% of the total observations in every time series were used as training set, whereas the rest 20% were used as the testing set. The details of every time series are discussed in brief in subsequent sections.


### 2.1.1. Daily

The Brent COPs time series contain a total of 12436 observations (spanning 2 February 1970 to 1 October 2017), the first 9,949 observations belong to the training series, and the rest 2487 part of the testing series. The next is the WTI COPs time series containing a total of 9060 observations (spanning 12 January 1983 to 1 October 2017), the first 7248 observations belong to the training series, while the rest of 1812 part of the testing series.


### 2.1.2. Weekly

The Brent COPs time series contain a total of 2488 observations (spanning 2 February 1970 to 1 October 2017), the first 1991 observations belong to the training series, and the rest 497 part of the testing series. The next is the WTI COPs time series containing a total of 1812 observations (spanning 12 January 1983 to 1 October 2017), the first 1450 observations belong to the training series, while the rest 362 part of the testing series.


## 2.1.1. Daily

The Brent COPs time series contain a total of 12436 observations (spanning 2 February 1970 to 1 October 2017), the first 9,949 observations belong to the training series, and the rest 2487 part of the testing series. The next is the WTI COPs time series containing a total of 9060 observations (spanning 12 January 1983 to 1 October 2017), the first 7248 observations belong to the training series, while the rest of 1812 part of the testing series.


## 2.1.2. Weekly

The Brent COPs time series contain a total of 2488 observations (spanning 2 February 1970 to 1 October 2017), the first 1991 observations belong to the training series, and the rest 497 part of the testing series. The next is the WTI COPs time series containing a total of 1812 observations (spanning 12 January 1983 to 1 October 2017), the first 1450 observations belong to the training series, while the rest 362 part of the testing series.


## 2.2. Framework of the Proposed Study

Reconstruction of IMFs obtained from EEMD is proposed in this study. Firstly, the obtained IMFs from EEMD are divided into two components, namely, stochastic and deterministic. In this study, the proposed procedure for dividing IMFs into two components is autocorrelation (correlation (yt,yt−1)) of every IMF. If the autocorrelation is less than a threshold value, then that IMF will be considered stochastic, otherwise deterministic. To obtained two components add all stochastic IMFs for SC and deterministic for DC. The idea of reconstructing the IMFs into SD components should be utilized in more detail for further improvement because the SC requires more attention which can help in improving the forecasting accuracy significantly. The IMFs being part of the SC handled separately because most of the variation in total shared by the SC [29]. Thus, in this study, the stochastic IMFs modelled separately, for every stochastic IMF a different model is chosen and for final forecasting ensemble all the models’ forecast, whereas all deterministic IMFs would be treated as one component. The ARIMA and ANN models will be fitted for stochastic IMFs and DC and for the final output, add all these components, so the proposed models are denoted by EEMD-R-ARIMA and EEMD-R-ANN. For ANN, the AR terms of the ARIMA model would be used as inputs. The complete framework of the proposed method is demonstrated in Figure 1.Figure 1
The framework of the proposed method.


## 2.3. Experimental Study

It is discussed that the main aim of this study is to use the minimum number of IMFs with higher forecasting accuracy. For a minimum number of IMFs use, the reconstruction of IMFs is proposed in this study. Thus, first obtain all IMFs from EEMD and divide them into two components, namely, stochastic and deterministic. However, a procedure or a threshold value is required which decides that the given IMF is stochastic or deterministic. The proposed procedure is the autocorrelation of IMFs from which a decision will be taken that the given IMF has stochastic or deterministic influences. Based on four different threshold values, i.e.,ρyt,yt−1=0.500,0.900,0.950,0.990 of autocorrelation, the obtained IMFs would be divided into two components. The SC consist of all those IMFs whose autocorrelation is less than the threshold value; otherwise, the IMF belongs to the DC. To obtain the best threshold value for the division of IMFs into SD components, synthetic time series would be used which is characterized through the additive noise and organized into four different scenarios as follows:(1)
In a time series composed of a sine function and normal distribution, the sine function represents the DC, while the normal distribution represents SC:(1)Yt=sin2πt+ε0,1.(2)
In a time series composed of a sine function and ARMA model, the sine function represents the DC, while the ARMA model represents the SC with an error of 0.25:(2)Yt=sin2πt+ARMA1,1.(3)
In a time series composed of a sine function, ARMA model, and normal distribution, the sine function represents the DC, while the ARMA model and normal distribution represent SC:(3)Yt=sin2πt+ARMA1,1+ε0,1.(4)
In a time series composed of a Lorenz system and normal distribution, the Lorenz system represents the DC, while the normal distribution represents SC:(4)Yt=Lorenz system10,28,83+ε0,1.Therefore, to divide the IMFs into SD components, in this study, autocorrelation is used. IMFs have stochastic influences being part of the SC while having deterministic influences being part of the DC. To obtain the original series, add both SD components. The following formula is used to compute the autocorrelation of all decomposed IMFs:(5)rYt,Yt−1=∑t=1T−1Yt−Y¯Yt−1−Y¯∑t=1TYt−Y¯2.The problem is how to decide that the given IMF has a stochastic or deterministic influence. To fix this problem, an experimental study is performed to obtain a threshold value rather a confidence interval from which we will be able to detach the IMFs regarding their influences. When these IMFs are detached according to their influences, the next step is to obtain SD components by summing up all IMFs, respectively. A correlation test is performed to check whether the detached components are independent. For testing, the following test statistic is used with the alternative hypothesisHA:ρSD≠0:(6)t=rSDn−21−rSD2∼tn−2,where n denotes the total number of observations and rSD represents the correlation among the SD components and is defined as follows:(7)rSD=∑t=1TYS,t−Y¯SYD,t−Y¯D∑t=1TYS,t−Y¯S2∑t=1TYD,t−Y¯D2.To double-check the distribution of these two components, another assessment is also steered which can test the difference of the two correlations produced by the synthetic data SD components and the reconstructing SD components obtained from IMFs. For testing, the following test statistic is used with an alternative hypothesis that the difference between the two correlations is not equal to zero, i.e.,H0:ρSDs≠ρSDdorρSDs−ρSDd≠0:(8)z=z1−z2δz1−z2∼N0,1,where z1=1/2log1+rSDs/1−rSDs, z2=1/2log1+rSDd/1−rSDd, δz1−z2=2/T−3, “s” represents the synthetic data, and “d” is used for decomposed data. The test is performed for different threshold values of the autocorrelation which divided the IMFs into two components (SD). The best threshold value will be chosen for the real-life datasets and should implement.The above four different scenarios are used to check how well the proposed method extracted the SD components. The purpose of using different series is to check the generalizability and robustness of the proposed method. One of the advantages of the use of synthetic data is that it provided a controlled scenario and better supports the conclusions.


## 2.4. Methods Used in This Study

### 2.4.1. The ARIMA Box–Jenkins Approach

The inclusion of both AR and MA terms makes the Box-Jenkins models exceptionally flexible. The interdependencies in time series (Yt) are measured by the AR terms, while the dependencies on preceding error terms are measured by the MA terms [30, 31]. An ARMA model of order (p, q) for a univariate series has the following form:(9)Yt=γ0+∑i=1pγiYt−i+∑j=1qβjεt−j+εt,where p and q represent the number of AR and MA terms, respectively, γ represents the constant and AR coefficients, whereas β is used for MA terms. The white noise process is denoted by ϵt∼N0,σ2. Using the Box–Jenkins procedure which needed stationary data, this can be achieved by differencing the time series many times (but at most two times usually). The ACF and PACF plots can be used to select the appropriate order of the polynomials and compare them with the theoretical bounds of [32]. Rival ARIMA models can also be selected based on AIC [33] and BIC [34]. The next is the diagnostic checking that the fitted model is adequate by using the LB test [35]. After obtaining the appropriate ARIMA model, the last step is the forecasting of the series.


### 2.4.2. Artificial Neural Network

ANNs have many distinguishing characteristics over the ARIMA models and very successful alternative to linear models to handle the nonlinearity factors when forecasting the time series. One of the best characteristics of an ANN is the universal approximation of any nonlinear function up to the best degree of accuracy of your desired [36, 37]. A most commonly used ANN is the FNN with a single hidden layer having one output node for forecasting the time series applications [38, 39]. An m×n×1 ANN model has the following form:(10)Yt=μ0+∑j=1nμjgϕ0j+∑i=1mϕijYt−i+εt,where μjj=0,1,2,…,n, ϕiji=1,2,…,m,j=1,2,…,n denotes the weights, μ0andϕ0j represent the bias terms, ϵt is used for the white noise, Yt−ii=1,2,…,m represents the input patterns, Yt represents the single output of the model, m is used for the number of input variables, n is used for the number of hidden nodes, and g⋅ is used for the transfer function of the hidden layer. Equation (11) performs nonlinear mapping from past observations Yt−1,Yt−2,…,Yt−m, for the future value Yt. i.e.,(11)Yt=ϑYt−1,Yt−2,…,Yt−m,ω+εt,where ϑ is a function determined by the neural network training and ω is a vector representing all parameters. Therefore, from some viewpoint, the FNN model looks like or equivalent to an NAR model [40, 41]. The ANN effectively modelled the nonlinear time series, but for linear problems, it provided mixed results [38].


### 2.4.3. EEMD

The new opportunities are provided in [42] to extract the time opposing elements from data after the advancement of transient local and versatile strategy for EMD. The transient region is a standout amongst the most critical qualities of EMD, which is procured from the spline fitting through maxima/minima of inputted data. The spline fitting has a high transient area and can be easily affirmed using numerical programming. The transient domain of EMD is all around spared if the quantity of filtering is smaller and fixed. The filtering process of EMD is worse with the excessive use of both transient region and oscillatory components. One of the benefits of the EMD procedure is that it automatically bypasses the stationarity assumption of the data. For the physical interpretability of the results, EMD territory gives a significant condition but not adequate. An EMD approach involves continuing with the levels of removing oscillatory sections of lesser replication, and the noise-wrapped results at one level can provoke continuing with the twisting of subsequent oscillatory components, making EMD happening hardly physically interpretable. Hence, due to the absence of this robustness problem, the EMD procedure becomes inadequate.Several attempts have been made to overcome the problem of robustness. For example, Huang et al. [42] proposed the discontinuity test in which variation is controlled by one element. Due to the controlled scenarios, it reduced the adaptiveness properties of EMD, demanding for more fascinating methodologies. To fulfill this robustness problem, a noise-assisted data analysis procedure, namely, EEMD was introduced in [43] which exactly solved the highlighted issue raised in EMD. The following are the steps involved in EEMD:(1)
A series of unique white noise summed with true seriesYt(2)
In the next step, the new series obtained in (9) can be decomposed accordingly into IMFs(3)
Repeat (9) and (10) by summing the different white noise series to the true series Yt(4)
As the outcome attained the ensemble means of the respective IMFsThe EEMD utilized the properties of white noise in the above steps and are as follows: the minima and maxima distribution should be uniform transiently and EMD being adequately a dual-channel bank for white-noise [44–46]. The minima and maxima distribution on all-time spans are relatively uniform with the added white noise. Another property of EEMD is the control on oscillatory elements after adding the white noise to the series which works like a dual filter bank which significantly reduces the mode mixing problem. Thus, as an outcome, the decomposition turns out to be steadier and physically significant.


### 2.4.4. EEMD-ARIMA and ANN Model

Wang et al. [47] summarized the hybrid model of EEMD-ARIMA by incorporating the steps as follows:(1)
Firstly, the true seriesYt be decomposed into different IMFs(2)
Secondly, for everykth IMF the best ARIMA model will be selected and forecast accordingly(3)
The last step is to obtain the targeted series summed all thekth IMFs outputThe EEMD-ANN model should be fitted in the same way above. Furthermore, for every ANN model, the delaying parameter or starting input parameter be the AR terms of every ARIMA model ofkth IMF.


## 2.4.1. The ARIMA Box–Jenkins Approach

The inclusion of both AR and MA terms makes the Box-Jenkins models exceptionally flexible. The interdependencies in time series (Yt) are measured by the AR terms, while the dependencies on preceding error terms are measured by the MA terms [30, 31]. An ARMA model of order (p, q) for a univariate series has the following form:(9)Yt=γ0+∑i=1pγiYt−i+∑j=1qβjεt−j+εt,where p and q represent the number of AR and MA terms, respectively, γ represents the constant and AR coefficients, whereas β is used for MA terms. The white noise process is denoted by ϵt∼N0,σ2. Using the Box–Jenkins procedure which needed stationary data, this can be achieved by differencing the time series many times (but at most two times usually). The ACF and PACF plots can be used to select the appropriate order of the polynomials and compare them with the theoretical bounds of [32]. Rival ARIMA models can also be selected based on AIC [33] and BIC [34]. The next is the diagnostic checking that the fitted model is adequate by using the LB test [35]. After obtaining the appropriate ARIMA model, the last step is the forecasting of the series.


## 2.4.2. Artificial Neural Network

ANNs have many distinguishing characteristics over the ARIMA models and very successful alternative to linear models to handle the nonlinearity factors when forecasting the time series. One of the best characteristics of an ANN is the universal approximation of any nonlinear function up to the best degree of accuracy of your desired [36, 37]. A most commonly used ANN is the FNN with a single hidden layer having one output node for forecasting the time series applications [38, 39]. An m×n×1 ANN model has the following form:(10)Yt=μ0+∑j=1nμjgϕ0j+∑i=1mϕijYt−i+εt,where μjj=0,1,2,…,n, ϕiji=1,2,…,m,j=1,2,…,n denotes the weights, μ0andϕ0j represent the bias terms, ϵt is used for the white noise, Yt−ii=1,2,…,m represents the input patterns, Yt represents the single output of the model, m is used for the number of input variables, n is used for the number of hidden nodes, and g⋅ is used for the transfer function of the hidden layer. Equation (11) performs nonlinear mapping from past observations Yt−1,Yt−2,…,Yt−m, for the future value Yt. i.e.,(11)Yt=ϑYt−1,Yt−2,…,Yt−m,ω+εt,where ϑ is a function determined by the neural network training and ω is a vector representing all parameters. Therefore, from some viewpoint, the FNN model looks like or equivalent to an NAR model [40, 41]. The ANN effectively modelled the nonlinear time series, but for linear problems, it provided mixed results [38].


## 2.4.3. EEMD

The new opportunities are provided in [42] to extract the time opposing elements from data after the advancement of transient local and versatile strategy for EMD. The transient region is a standout amongst the most critical qualities of EMD, which is procured from the spline fitting through maxima/minima of inputted data. The spline fitting has a high transient area and can be easily affirmed using numerical programming. The transient domain of EMD is all around spared if the quantity of filtering is smaller and fixed. The filtering process of EMD is worse with the excessive use of both transient region and oscillatory components. One of the benefits of the EMD procedure is that it automatically bypasses the stationarity assumption of the data. For the physical interpretability of the results, EMD territory gives a significant condition but not adequate. An EMD approach involves continuing with the levels of removing oscillatory sections of lesser replication, and the noise-wrapped results at one level can provoke continuing with the twisting of subsequent oscillatory components, making EMD happening hardly physically interpretable. Hence, due to the absence of this robustness problem, the EMD procedure becomes inadequate.Several attempts have been made to overcome the problem of robustness. For example, Huang et al. [42] proposed the discontinuity test in which variation is controlled by one element. Due to the controlled scenarios, it reduced the adaptiveness properties of EMD, demanding for more fascinating methodologies. To fulfill this robustness problem, a noise-assisted data analysis procedure, namely, EEMD was introduced in [43] which exactly solved the highlighted issue raised in EMD. The following are the steps involved in EEMD:(1)
A series of unique white noise summed with true seriesYt(2)
In the next step, the new series obtained in (9) can be decomposed accordingly into IMFs(3)
Repeat (9) and (10) by summing the different white noise series to the true series Yt(4)
As the outcome attained the ensemble means of the respective IMFsThe EEMD utilized the properties of white noise in the above steps and are as follows: the minima and maxima distribution should be uniform transiently and EMD being adequately a dual-channel bank for white-noise [44–46]. The minima and maxima distribution on all-time spans are relatively uniform with the added white noise. Another property of EEMD is the control on oscillatory elements after adding the white noise to the series which works like a dual filter bank which significantly reduces the mode mixing problem. Thus, as an outcome, the decomposition turns out to be steadier and physically significant.


## 2.4.4. EEMD-ARIMA and ANN Model

Wang et al. [47] summarized the hybrid model of EEMD-ARIMA by incorporating the steps as follows:(1)
Firstly, the true seriesYt be decomposed into different IMFs(2)
Secondly, for everykth IMF the best ARIMA model will be selected and forecast accordingly(3)
The last step is to obtain the targeted series summed all thekth IMFs outputThe EEMD-ANN model should be fitted in the same way above. Furthermore, for every ANN model, the delaying parameter or starting input parameter be the AR terms of every ARIMA model ofkth IMF.


## 2.5. Forecast Evaluation Criterions

Forecasting accuracy is used to measure the degree of exactness among the proposed and state-of-the-art methods. When there are competing models, the forecasting accuracy measures are the most important criteria. In this study, four-typical forecasting accuracy measures are used consisting of MAE, RMSE, MAPE, and DS. Typically, these methods are defined as follows:(1)
MAE(12)MAE=1F∑t=1FY^t−Yt.(2)
RMSE(13)RMSE=1F∑t=1FY^t−Yt2.(3)
MAPE(14)MAPE=1F∑t=1FY^t−YtYt×100,whereF represents the total number of forecasts of the testing set, the original value of the given time series at time t is denoted by Yt and forecasted value is denoted by Y^t. Four different datasets are used in this study for validation and justification purposes. The next is the directional prediction measure which is more appropriate for business purposes since the investors are concerned with the market trend. Thus, DS is used for competing models to measure the directional forecasting accuracy and defined as follows.(4)
DS(15)DS=1F∑t=1FDt×100,whereDt=1, if Yt+1−Yt×Y^t+1−Yt≥0, else Dt=0, and F corresponds to the number of total forecasts of testing set. Apparently, the lower the values of MAPE, RMSE, and MAE, the higher the forecast accuracy in contrast with higher the values of DS produced more accurate forecasts. Next, to compare the two models prediction errors, we employed the DM test which is defined as follows.(5)
DM TEST(16)DM=d¯1/nγ0+2∑k=1nγk.It is used to test the statistical significance of the forecast accuracy of two forecast methods. Variabledt is the subtraction of absolute errors of the two methods. d¯=1/n∑t=1ndt is the mean of dt. γk=Covdt,dt−k is the autocovariance at lag k. DM∼N(0, 1), if the p value > α, we conclude there is no significant difference between the two forecasts.


## 3. Results

### 3.1. ARIMA Model

Data stationarity is one of the assumptions which is required when using the ARMA model. In most of the cases, the time series for real data is not stationary because of the trend factor, so the ARIMA model applicability comes out. To get the stationary time series, we take the successive differences of the original series [48, 49]. ADF test is used to check the stationarity of a time series [50]. Using the ARIMA model, first identify the appropriate model by selecting the order of AR and MA terms after obtaining the stationary time series. ACF and PACF plots are used for the selection of best orders of AR and MA terms. Other approaches are also available for obtaining the orders of AR and MA terms called theoretical approaches consisting of AIC and BIC. For model adequacy checking, the LB test is used. In the last step, we forecast the appropriate time series. In this section, the ARIMA model is fitted for all training periods and forecast the testing periods. Table 1 represents the order of best ARIMA models for all datasets.Table 1
p values of ADF unit root test and order of selected ARIMA models for all data.
DataSeriesADF testp valueADF test 1st diff p valueOrder of ARIMA modelBrentDaily0.990<0.010(1, 1, 0)Weekly0.990<0.010(2, 1, 2)WTIDaily0.349<0.010(3, 1, 2)Weekly0.054<0.010(5, 1, 5)The fitted ARIMA models are found in Table2, where the estimated values of their respective parameters are presented with the standard error in parenthesis.Table 2
Estimated parameters and their standard errors of the ARIMA models.
DataSeriesAR coefficientsMA coefficientsBrentDaily−0.0291 (0.010)Weekly0.9111 (0.181)−0.7531 (0.115)−0.9081 (0.168)0.7921 (0.102)WTIDaily−1.1741 (0.053)−0.9071 (0.047)−0.0761 (0.013)1.1301 (0.052)0.8091 (0.050)Weekly0.6031 (0.035)0.2421 (0.026)0.2241 (0.022)−0.6751 (0.033)−0.2001 (0.034)−0.1611 (0.031)0.6631 (0.027)−0.8751 (0.032)−0.6311 (0.034)0.8571 (0.031)1p value less than 0.01.The next step after estimating the parameters is the extraction of the fitted model residuals and their testing to check whether the residuals are uncorrelated. The model diagnostic tests are also performed for all four datasets that the fitted models are adequate and can be used for forecasting.From Figure2, it is observed that the LB test reveals that all the fitted models are adequate and could be used for forecasting. Figure 2 shows on the third graph of LB statistic that all p values for all datasets are greater than 0.05, which means that the hypothesis of no serial correlation among fitted model residuals is not significant. So, all these models are adequate and provide the best forecast for the future. The last is the forecasting accuracy of these models which are presented in Tables 3 and 4 for daily and weekly datasets, respectively.Figure 2
Diagnostic test results: (a) Brent daily, (b) WTI daily, (c) Brent weekly, and (d) WTI weekly.
(a)(b)(c)(d)Table 3
Measurements of testing forecasts accuracies of daily COPs.
ModelBrentWTIMAERMSEMAPEDSMAERMSEMAPEDSEEMD-R-ANN0.3910.5020.54289.780.3720.4760.57291.44EEMD-R-ARIMA0.4080.5340.55988.250.3750.4790.57790.07EEMD-SD-ANN0.7510.9991.02276.030.6910.9131.01877.42EEMD-SD-ARIMA0.7571.0081.02775.380.7200.9571.05575.99EEMD-AMI-ANN1.0311.3691.40562.470.9911.2631.40462.07EEMD-AMI-ARIMA1.0361.3821.40461.970.9991.2711.41060.92EEMD-ANN0.5040.6440.68181.130.4100.5280.62083.38EEMD-ARIMA0.5190.6670.69779.830.4190.5330.62981.84ANN1.0501.4601.42752.291.0021.3821.47855.55ARIMA1.0631.4871.43351.651.0091.3891.48953.37Table 4
Measurements of testing forecasts accuracies of weekly COPs.
ModelBrentWTIMAERMSEMAPEDSMAERMSEMAPEDSEEMD-R-ANN0.6820.8620.93393.160.6370.8240.93694.46EEMD-R-ARIMA0.6840.8650.93692.080.6460.8380.95292.81EEMD-SD-ANN1.7732.2832.35975.651.7772.2492.57978.39EEMD-SD-ARIMA1.7892.3452.36575.011.7862.3022.58277.71EEMD-AMI-ANN2.4323.2943.29556.742.2792.9443.37158.17EEMD-AMI-ARIMA2.5343.3943.43955.812.3122.9763.40657.49EEMD-ANN0.7470.9711.02181.890.7280.9201.06983.66EEMD-ARIMA0.7670.9951.05280.170.7860.9971.16281.93ANN2.4503.3323.32253.322.3183.0153.41154.14ARIMA2.5723.4673.47151.292.3573.0643.46552.94


### 3.2. ANN Model

The input structure of ANN is subject to the order of AR terms of the best chosen ARIMA model. To fix the input parameters of ANN, the order of AR termsp is used provided in Table 1. After fulfilling the following requirements, the ANN structure is constructed:(1)
Feedforward backpropagation used as a network type.(2)
ANN inputs for WTI weekly data areYt−1,Yt−2,Yt−3,Yt−4,Yt−5 and the same for Brent are Yt−1,Yt−2. These inputs are automatically created by MATLAB ntstool once the training and target time series are imported to ANN structure.(3)
Levenberg–Marquardt training algorithms are used as a training function in this study.(4)
The weights are determined randomly depending on the strategies of the ANN toolboxes.(5)
The maximum neurons in the hidden layer are equal to2×k+1, where k is the number of inputs [51]. The maximum number of hidden neurons for WTI weekly data is 11 and for Brent is 5.(6)
The nonlinearity of crude oil prices is obligated, so to filter the nonlinearity in the hidden layer, a nonlinear transfer function is chosen such as tan-sigmoid or log-sigmoid. In this study, the NAR problem is selected from the ntstool.The complete structure of ANN is presented in Figure3 for both daily and weekly datasets. After completing the ANN construction, the next step is the training and testing processes to obtain the fitted and forecasting series. Training and testing processes only need to import the input and target variables for training and testing periods, respectively.Figure 3
ANN structure of all datasets. (a) Brent daily. (b) WTI daily. (c) Brent weekly. (d) WTI weekly.
(a)(b)(c)(d)


### 3.3. Data Decomposition

For using EEMD, the two parameters should be fixed in advance which are the white noise amplitude and the number of ensembles because it affects the decomposition process, and Wu and Huang [43] described the rule for these parameters in the form of equation, i.e., εh=e/h, where h is the ensemble number, εh representing the error of standard deviation, and the amplitude of added white noise represented by e. One important thing should be kept in mind when added white noise to the original time series is that its amplitude should not be too large or small because if large, the EEMD will return some redundant IMFs; on the contrary, if small, the EEMD IMFs remain the same as EMD. The authors [43, 52] suggested that the amplitude of the added white noise is equal to 20% of the standard deviation. About the white noise amplitude and size of the ensemble number, more details can be found in the study of [43]. In this study the EEMD technique is used. The two parameters of EEMD, the size of the ensemble number, and the white noise amplitude are equal to 100 and 0.2 times the standard deviation, respectively.EEMD technique is employed to all datasets using the above description. In this study, the r package “Rlibeemd” [53] is used for data decomposition. All the decomposed IMFs are presented in Figure 4 for all four datasets.Figure 4
Decomposition of all datasets using EEMD. (a) Brent daily. (b) WTI daily. (c) Brent weekly. (d) WTI weekly.
(a)(b)(c)(d)The total number of IMFs including residual for daily Brent and WTI is 13, while for weekly, Brent and WTI are 11 and 10, respectively. From Figures4(a)–4(d), IMF.1, IMF.2, IMF.3, and IMF.4 have no specific trend, and almost the observations distributed randomly. But after the IMF.4, the IMFs follow a trend and move slowly around the long-term mean, whereas the last IMF called residual represents the overall trend of the series. The overall decomposition provided by EEMD is physically meaningful, all IMFs are independent, and each IMF is consistent with the previous IMF [54]. Therefore, the decomposition obtained from EEMD is very helpful because it transforms the nonstationary and nonlinear series into stationary series which improve the forecasting accuracy of the models.


### 3.4. Simulations for Threshold Value Determination

This study aims to determine a threshold value for the division of IMFs into SD components. After obtaining the decomposed IMFs for all synthetic datasets, the autocorrelation is computed for every IMF. From autocorrelation, it is easy to determine that the current observation depends how much on its past observations. Thus, closing value to zero means there is no dependency on the current value on its past values, whereas closing value to one means the current value strongly influenced by its past values. The problem is how to fix a point that one can say that, below this value, the considered IMF is stochastic. To solve this issue different values of autocorrelation are considered for dividing the IMFs into two components. First, the IMFs are separated into two different components through autocorrelation using the selected four values, i.e.,ρyt,yt−1=0.500,0.900,0.950,0.990. The SC consists of all IMFs whose autocorrelation is less than the given threshold value otherwise included in the DC. Initially different values of the autocorrelation were selected and tested, i.e., ρyt,yt−1=0.50,0.60,0.70,0.80,0.90,0.95,0.99 for the reconstruction of IMFs, but after the computation, the values ρyt,yt−1=0.60,0.70,0.80,0.90 almost provide the same results; thus, for avoiding some redundant analysis, only ρyt,yt−1=0.900 is chosen from the proceeding group. Therefore, in this study, the following four values of autocorrelation are used for reconstruction of IMFs, i.e., ρyt,yt−1=0.500,0.900,0.950,0.990. Thus, in this section, only the synthetic data are used because it provided a controlled scenario and better support the conclusion. The number of IMFs being part of the SD components is presented in Table 5 for all scenarios.Table 5
The number of IMFs being part of the SD components for all threshold values and all scenarios.
Dataρyt,yt−1Scenario IScenario IIScenario IIIScenario IVSCDCSCDCSCDCSCDC5000.50012–812–812–8——0.9001-23–81-23–81–34–81-23–80.9501–34–81-23–81–34–81–34–80.9901–34–81–34–81–45–81–34–810000.50012–912–912–9——0.9001–34–91-23–91-23–91-23–90.9501–34–91–34–91–34–91–34–90.9901–45–91–45–91–45–91–34–920000.5001-23–1012–1012–1012–100.9001–34–101–34–101–34–101-23–100.9501–34–101–34–101–34–101–34–100.9901–45–101–45–101–45–101–34–1050000.50012–12——12–1212–120.9001–34–121–34–121–34–121–34–120.9501–45–121–34–121–34–121–34–120.9901–56–121–45–121–45–121–45–1210,0000.50012–13——12–1312–130.9001-23–131–34–131–34–131–34–130.9501–34–131–45–131–34–131–34–130.9901–45–131–56–131–56–121–45–13The next step after dividing the IMFs into two components using the above threshold values is the testing of the obtained two components that they are independently distributed. TheT-test is used to test the hypothesis of zero correlation among the SD components. The best threshold value is selected for which the null hypothesis is accepted for all scenarios and all sample sizes. Table 6 presented the p values that both SD components are independently distributed.Table 6
p values of Pearson correlation test for all samples and threshold values.
Dataρyt,yt−1Scenario IScenario IIScenario IIIScenario IV5000.5000.7030.3200.935—0.9000.8530.5180.4540.8900.9500.2520.5180.4540.3600.9900.252<0.0010.1530.36010000.5000.6730.2400.619—0.9000.7830.3000.3300.8580.9500.7830.3290.2170.9380.990<0.001<0.0010.0440.93820000.5000.6240.1920.4520.5090.9000.7270.5140.1930.7200.9500.7270.5140.1930.5280.9900.8380.5630.0010.52850000.500<0.001—0.0910.3570.9000.1310.1630.1010.6250.9500.6630.1630.1010.6250.9900.7470.2980.0010.00310,0000.5000.039—0.1730.1750.9000.1380.0190.1400.5240.9500.3690.0830.1400.5240.9900.6050.283<0.0010.006The boldp values in Table 6 are statistically significant at 5% and communicated that the two components are not independently distributed.To double-check the threshold values another test is also performed which is the difference of correlation test. This test will confirm that the correlation among the original components and the components obtained from the reconstruction of IMFs is equal to zero. So, only those threshold values of autocorrelation are selected, for which the null hypothesis of no difference is not rejected for all scenarios and all samples.From Tables6 and 7, it is observed that, from scenarios I and IV, the recommended threshold values for dividing the IMFs into two components are ρyt,yt−1=0.90,0.95, and from scenario III, the recommended values are ρyt,yt−1=0.50,0.90,0.95, whereas from scenario II, the recommended value is ρyt,yt−1=0.95. So, it is observed from all scenarios that the only one value which significantly divided the IMFs into two components is ρyt,yt−1=0.95. Thus, the recommended value of autocorrelation for dividing the IMFs into two components is ρyt,yt−1=0.95. The summary of the threshold value determination is presented in Table 8.Table 7
p values of correlation difference test of true correlation and among the components obtained from the reconstruction of IMFs.
Dataρyt,yt−1Scenario IScenario IIScenario IIIScenario IV5000.5000.7600.4900.550—0.9000.8500.3500.2400.7300.9500.4000.3500.2400.8400.9900.400<0.0010.1000.84010000.5000.5200.6900.360—0.9000.5900.7500.2100.6400.9500.5900.8100.1500.5200.990<0.001<0.0010.0480.52020000.5001.0000.5200.4100.4600.9000.7000.8500.2300.6000.9500.7000.8500.2300.4700.9900.6600.9000.0080.47050000.500<0.001—0.7200.9500.9000.0400.1800.7500.8100.9500.1200.1800.7500.8100.9900.1300.0500.1500.13010,0000.5000.140—0.2000.0200.9000.2900.0300.1600.1000.9500.5200.0900.1600.1000.9900.7100.190<0.0010.001The bold values indicate that they are statistically significant at 5% and communicated that the true correlation and among the components obtained from the reconstruction of IMFs are not independent.Table 8
Summary of the threshold values for all four scenarios.
Level of significance (%)ρyt,yt−1Scenario IScenario IIScenario IIIScenario IV50.500✗✗✓✗0.900✗✗✓✓0.950✓✓✓✓0.990✗✗✗✗10.500✗✗✓✗0.900✓✓✓✓0.950✓✓✓✓0.990✗✗✗✗✓The threshold value is recommended. ✗The threshold value is not recommended.The only recommended value at both levels of significance is0.95. Thus, for division of IMFs into SD components, the threshold value of 0.95 is recommended and could be used. In Section 3.5, the real-world applications are tested by utilizing the recommended threshold value of 0.95.


### 3.5. IMFs Reconstruction

In this study, the reconstruction of IMFs is proposed based on the autocorrelation of everykth IMF obtained from EEMD. Each IMF autocorrelation is compared with a threshold value of ρyt,yt−1=0.95. An IMF is considered stochastic if its autocorrelation is less than the threshold value; otherwise, it is considered as deterministic. For comparison, the IMFs are also reconstructed through AMI. The AMI is simply the visual inspection of the plots which are obtained for all IMFs of each dataset. The detailed analysis of every dataset is followed in the subsequent sections.From Section3.3, it is detected that the total number of IMFs for daily datasets is 12 plus residue. For simplicity, the residue term is replaced with the IMF, so the total number of IMFs for both datasets is 13. For weekly data, the total number of IMFs is 11 for Brent, while for WTI, it is 10. Autocorrelation of all IMFs is computed and presented in Table 9.Table 9
Autocorrelation of Brent and WTI daily and weekly data.
IMFsDailyWeeklyBrentWTIBrentWTIIMF.1−0.467−0.484−0.479−0.448IMF.20.4730.4880.5780.592IMF.30.8360.8500.8890.892IMF.40.9530.9540.9830.987IMF.50.9870.9880.9950.997IMF.60.9970.9980.9990.999IMF.70.9990.9990.9990.999IMF.80.9991.0001.0001.000IMF.91.0001.0001.0001.000IMF.101.0001.0001.0001.000IMF.111.0001.0001.000—IMF.121.0001.000——IMF.131.0001.000——It is observed from Table9 that the autocorrelation of the IMF.1, IMF.2, and IMF.3 for all datasets is less than 0.95. So, the first three IMFs are considered as stochastic in all datasets, while the rest is considered as deterministic. So, to make a SC, all stochastic IMFs are summed, and for DC, all deterministic IMFs are summed. Next, to obtain the SD components from the AMI method for Brent daily data, their plots are presented in Figure 5.Figure 5
AMI plots of Brent daily data for all IMFs.From Figure5, it is observed that the AMI plots have the same pattern after the 6th IMF, so the first six IMFs are considered as stochastic, and the remainder, as deterministic. To make two components, add the first six IMFs for SC and from seven to thirteen for DC. Next, to obtain the SD components of WTI daily data, the plots of AMI for all IMFs are presented in Figure 6.Figure 6
AMI plots of WTI daily data for all IMFs.From Figure6, it is observed that the AMI plots have the same pattern after the 6th IMF, so the first six IMFs are considered as stochastic, and the remainder, as deterministic. To make two components, add the first six IMFs for SC and from seven to thirteen for DC. Next, to obtain the SD components of Brent weekly data, the plots of AMI for all IMFs are presented in Figure 7.Figure 7
AMI plots of Brent weekly data for all IMFs.From Figure7, it is observed that the AMI plots have the same pattern after the 6th IMF, so the first six IMFs are considered as stochastic, and the remainder, as deterministic. To make two components, add the first six IMFs for SC and from seven to eleven for DC. Next, to obtain the SD components of WTI weekly data, the plots of AMI for all IMFs are presented in Figure 8.Figure 8
AMI plots of WTI weekly data for all IMFs.From Figure8, it is observed that the AMI plots have the same pattern after the 6th IMF, so the first six IMFs are considered as stochastic, and the remainder, as deterministic. To make two components, add the first six IMFs for SC and from seven to ten for DC. The next is the correlation test that the SD components are independently obtained from the reconstruction of IMFs.From Table10, it is observed that, for all reconstructed components, the hypothesis of zero correlation is not rejected, meaning that the division of SD components is independent at a 5% level of significance. Moreover, there is a theorem, put forward by Takens [55], known as the embedding theorem, stating that a time series can be reconstructed in vectors with m values (described as the embedding dimension). Each of these values corresponds to an observation that is spaced out in intervals, following a time delay (or separation dimension) called τ. As numerous studies have shown, these vectors represent the interdependent relationships between observations, and since they increase the accuracy of the modelling, they also enhance the accuracy of the prediction [56]. The plots of the AMI of the SD components are presented in Figure 9, showing the interdependency among the SD components.Table 10
Test of correlation among the two components for daily and weekly data.
DataMethodCorrelationt statisticd.fp valueBrentDailyR0.0040.534124340.592AMI0.0101.112124340.266WeeklyR0.0160.79124860.429AMI0.0261.29124860.197WTIDailyR0.0090.87390580.383AMI0.0100.92290580.356WeeklyR0.0170.76018100.447AMI0.0350.50118100.133Alternative hypothesis: true correlation is not equal to zero.Figure 9
AMI plots of SD components obtained from the proposed approach. (a) Daily Brent. (b) Daily WTI. (c) Weekly Brent. (d) Weekly WTI.
(a)(b)(c)(d)From Figures9(a)–9(d), it is observed that the division of IMFs into SD components is independently distributed. After obtaining the two components, the next step is the appropriate model selection for every component using the same steps discussed above. Every model is selected based on the LB test for which the p values are greater than 0.05 for both daily and weekly datasets, all IMFs, and components so that all the fitted models are adequate and could be used for forecasting. In the next step, all the selected models are used for forecasting, and the last 20% data of every series were retained as the testing set forecasted. For forecasting of ARIMA models, r software is used. For ANN, the data are transferred to MATLAB, and ntstool is used for forecasting every component. After forecasting the series, the accuracy measures are calculated and presented in Tables 3 and 4 for both daily and weekly series.


## 3.1. ARIMA Model

Data stationarity is one of the assumptions which is required when using the ARMA model. In most of the cases, the time series for real data is not stationary because of the trend factor, so the ARIMA model applicability comes out. To get the stationary time series, we take the successive differences of the original series [48, 49]. ADF test is used to check the stationarity of a time series [50]. Using the ARIMA model, first identify the appropriate model by selecting the order of AR and MA terms after obtaining the stationary time series. ACF and PACF plots are used for the selection of best orders of AR and MA terms. Other approaches are also available for obtaining the orders of AR and MA terms called theoretical approaches consisting of AIC and BIC. For model adequacy checking, the LB test is used. In the last step, we forecast the appropriate time series. In this section, the ARIMA model is fitted for all training periods and forecast the testing periods. Table 1 represents the order of best ARIMA models for all datasets.Table 1
p values of ADF unit root test and order of selected ARIMA models for all data.
DataSeriesADF testp valueADF test 1st diff p valueOrder of ARIMA modelBrentDaily0.990<0.010(1, 1, 0)Weekly0.990<0.010(2, 1, 2)WTIDaily0.349<0.010(3, 1, 2)Weekly0.054<0.010(5, 1, 5)The fitted ARIMA models are found in Table2, where the estimated values of their respective parameters are presented with the standard error in parenthesis.Table 2
Estimated parameters and their standard errors of the ARIMA models.
DataSeriesAR coefficientsMA coefficientsBrentDaily−0.0291 (0.010)Weekly0.9111 (0.181)−0.7531 (0.115)−0.9081 (0.168)0.7921 (0.102)WTIDaily−1.1741 (0.053)−0.9071 (0.047)−0.0761 (0.013)1.1301 (0.052)0.8091 (0.050)Weekly0.6031 (0.035)0.2421 (0.026)0.2241 (0.022)−0.6751 (0.033)−0.2001 (0.034)−0.1611 (0.031)0.6631 (0.027)−0.8751 (0.032)−0.6311 (0.034)0.8571 (0.031)1p value less than 0.01.The next step after estimating the parameters is the extraction of the fitted model residuals and their testing to check whether the residuals are uncorrelated. The model diagnostic tests are also performed for all four datasets that the fitted models are adequate and can be used for forecasting.From Figure2, it is observed that the LB test reveals that all the fitted models are adequate and could be used for forecasting. Figure 2 shows on the third graph of LB statistic that all p values for all datasets are greater than 0.05, which means that the hypothesis of no serial correlation among fitted model residuals is not significant. So, all these models are adequate and provide the best forecast for the future. The last is the forecasting accuracy of these models which are presented in Tables 3 and 4 for daily and weekly datasets, respectively.Figure 2
Diagnostic test results: (a) Brent daily, (b) WTI daily, (c) Brent weekly, and (d) WTI weekly.
(a)(b)(c)(d)Table 3
Measurements of testing forecasts accuracies of daily COPs.
ModelBrentWTIMAERMSEMAPEDSMAERMSEMAPEDSEEMD-R-ANN0.3910.5020.54289.780.3720.4760.57291.44EEMD-R-ARIMA0.4080.5340.55988.250.3750.4790.57790.07EEMD-SD-ANN0.7510.9991.02276.030.6910.9131.01877.42EEMD-SD-ARIMA0.7571.0081.02775.380.7200.9571.05575.99EEMD-AMI-ANN1.0311.3691.40562.470.9911.2631.40462.07EEMD-AMI-ARIMA1.0361.3821.40461.970.9991.2711.41060.92EEMD-ANN0.5040.6440.68181.130.4100.5280.62083.38EEMD-ARIMA0.5190.6670.69779.830.4190.5330.62981.84ANN1.0501.4601.42752.291.0021.3821.47855.55ARIMA1.0631.4871.43351.651.0091.3891.48953.37Table 4
Measurements of testing forecasts accuracies of weekly COPs.
ModelBrentWTIMAERMSEMAPEDSMAERMSEMAPEDSEEMD-R-ANN0.6820.8620.93393.160.6370.8240.93694.46EEMD-R-ARIMA0.6840.8650.93692.080.6460.8380.95292.81EEMD-SD-ANN1.7732.2832.35975.651.7772.2492.57978.39EEMD-SD-ARIMA1.7892.3452.36575.011.7862.3022.58277.71EEMD-AMI-ANN2.4323.2943.29556.742.2792.9443.37158.17EEMD-AMI-ARIMA2.5343.3943.43955.812.3122.9763.40657.49EEMD-ANN0.7470.9711.02181.890.7280.9201.06983.66EEMD-ARIMA0.7670.9951.05280.170.7860.9971.16281.93ANN2.4503.3323.32253.322.3183.0153.41154.14ARIMA2.5723.4673.47151.292.3573.0643.46552.94


## 3.2. ANN Model

The input structure of ANN is subject to the order of AR terms of the best chosen ARIMA model. To fix the input parameters of ANN, the order of AR termsp is used provided in Table 1. After fulfilling the following requirements, the ANN structure is constructed:(1)
Feedforward backpropagation used as a network type.(2)
ANN inputs for WTI weekly data areYt−1,Yt−2,Yt−3,Yt−4,Yt−5 and the same for Brent are Yt−1,Yt−2. These inputs are automatically created by MATLAB ntstool once the training and target time series are imported to ANN structure.(3)
Levenberg–Marquardt training algorithms are used as a training function in this study.(4)
The weights are determined randomly depending on the strategies of the ANN toolboxes.(5)
The maximum neurons in the hidden layer are equal to2×k+1, where k is the number of inputs [51]. The maximum number of hidden neurons for WTI weekly data is 11 and for Brent is 5.(6)
The nonlinearity of crude oil prices is obligated, so to filter the nonlinearity in the hidden layer, a nonlinear transfer function is chosen such as tan-sigmoid or log-sigmoid. In this study, the NAR problem is selected from the ntstool.The complete structure of ANN is presented in Figure3 for both daily and weekly datasets. After completing the ANN construction, the next step is the training and testing processes to obtain the fitted and forecasting series. Training and testing processes only need to import the input and target variables for training and testing periods, respectively.Figure 3
ANN structure of all datasets. (a) Brent daily. (b) WTI daily. (c) Brent weekly. (d) WTI weekly.
(a)(b)(c)(d)


## 3.3. Data Decomposition

For using EEMD, the two parameters should be fixed in advance which are the white noise amplitude and the number of ensembles because it affects the decomposition process, and Wu and Huang [43] described the rule for these parameters in the form of equation, i.e., εh=e/h, where h is the ensemble number, εh representing the error of standard deviation, and the amplitude of added white noise represented by e. One important thing should be kept in mind when added white noise to the original time series is that its amplitude should not be too large or small because if large, the EEMD will return some redundant IMFs; on the contrary, if small, the EEMD IMFs remain the same as EMD. The authors [43, 52] suggested that the amplitude of the added white noise is equal to 20% of the standard deviation. About the white noise amplitude and size of the ensemble number, more details can be found in the study of [43]. In this study the EEMD technique is used. The two parameters of EEMD, the size of the ensemble number, and the white noise amplitude are equal to 100 and 0.2 times the standard deviation, respectively.EEMD technique is employed to all datasets using the above description. In this study, the r package “Rlibeemd” [53] is used for data decomposition. All the decomposed IMFs are presented in Figure 4 for all four datasets.Figure 4
Decomposition of all datasets using EEMD. (a) Brent daily. (b) WTI daily. (c) Brent weekly. (d) WTI weekly.
(a)(b)(c)(d)The total number of IMFs including residual for daily Brent and WTI is 13, while for weekly, Brent and WTI are 11 and 10, respectively. From Figures4(a)–4(d), IMF.1, IMF.2, IMF.3, and IMF.4 have no specific trend, and almost the observations distributed randomly. But after the IMF.4, the IMFs follow a trend and move slowly around the long-term mean, whereas the last IMF called residual represents the overall trend of the series. The overall decomposition provided by EEMD is physically meaningful, all IMFs are independent, and each IMF is consistent with the previous IMF [54]. Therefore, the decomposition obtained from EEMD is very helpful because it transforms the nonstationary and nonlinear series into stationary series which improve the forecasting accuracy of the models.


## 3.4. Simulations for Threshold Value Determination

This study aims to determine a threshold value for the division of IMFs into SD components. After obtaining the decomposed IMFs for all synthetic datasets, the autocorrelation is computed for every IMF. From autocorrelation, it is easy to determine that the current observation depends how much on its past observations. Thus, closing value to zero means there is no dependency on the current value on its past values, whereas closing value to one means the current value strongly influenced by its past values. The problem is how to fix a point that one can say that, below this value, the considered IMF is stochastic. To solve this issue different values of autocorrelation are considered for dividing the IMFs into two components. First, the IMFs are separated into two different components through autocorrelation using the selected four values, i.e.,ρyt,yt−1=0.500,0.900,0.950,0.990. The SC consists of all IMFs whose autocorrelation is less than the given threshold value otherwise included in the DC. Initially different values of the autocorrelation were selected and tested, i.e., ρyt,yt−1=0.50,0.60,0.70,0.80,0.90,0.95,0.99 for the reconstruction of IMFs, but after the computation, the values ρyt,yt−1=0.60,0.70,0.80,0.90 almost provide the same results; thus, for avoiding some redundant analysis, only ρyt,yt−1=0.900 is chosen from the proceeding group. Therefore, in this study, the following four values of autocorrelation are used for reconstruction of IMFs, i.e., ρyt,yt−1=0.500,0.900,0.950,0.990. Thus, in this section, only the synthetic data are used because it provided a controlled scenario and better support the conclusion. The number of IMFs being part of the SD components is presented in Table 5 for all scenarios.Table 5
The number of IMFs being part of the SD components for all threshold values and all scenarios.
Dataρyt,yt−1Scenario IScenario IIScenario IIIScenario IVSCDCSCDCSCDCSCDC5000.50012–812–812–8——0.9001-23–81-23–81–34–81-23–80.9501–34–81-23–81–34–81–34–80.9901–34–81–34–81–45–81–34–810000.50012–912–912–9——0.9001–34–91-23–91-23–91-23–90.9501–34–91–34–91–34–91–34–90.9901–45–91–45–91–45–91–34–920000.5001-23–1012–1012–1012–100.9001–34–101–34–101–34–101-23–100.9501–34–101–34–101–34–101–34–100.9901–45–101–45–101–45–101–34–1050000.50012–12——12–1212–120.9001–34–121–34–121–34–121–34–120.9501–45–121–34–121–34–121–34–120.9901–56–121–45–121–45–121–45–1210,0000.50012–13——12–1312–130.9001-23–131–34–131–34–131–34–130.9501–34–131–45–131–34–131–34–130.9901–45–131–56–131–56–121–45–13The next step after dividing the IMFs into two components using the above threshold values is the testing of the obtained two components that they are independently distributed. TheT-test is used to test the hypothesis of zero correlation among the SD components. The best threshold value is selected for which the null hypothesis is accepted for all scenarios and all sample sizes. Table 6 presented the p values that both SD components are independently distributed.Table 6
p values of Pearson correlation test for all samples and threshold values.
Dataρyt,yt−1Scenario IScenario IIScenario IIIScenario IV5000.5000.7030.3200.935—0.9000.8530.5180.4540.8900.9500.2520.5180.4540.3600.9900.252<0.0010.1530.36010000.5000.6730.2400.619—0.9000.7830.3000.3300.8580.9500.7830.3290.2170.9380.990<0.001<0.0010.0440.93820000.5000.6240.1920.4520.5090.9000.7270.5140.1930.7200.9500.7270.5140.1930.5280.9900.8380.5630.0010.52850000.500<0.001—0.0910.3570.9000.1310.1630.1010.6250.9500.6630.1630.1010.6250.9900.7470.2980.0010.00310,0000.5000.039—0.1730.1750.9000.1380.0190.1400.5240.9500.3690.0830.1400.5240.9900.6050.283<0.0010.006The boldp values in Table 6 are statistically significant at 5% and communicated that the two components are not independently distributed.To double-check the threshold values another test is also performed which is the difference of correlation test. This test will confirm that the correlation among the original components and the components obtained from the reconstruction of IMFs is equal to zero. So, only those threshold values of autocorrelation are selected, for which the null hypothesis of no difference is not rejected for all scenarios and all samples.From Tables6 and 7, it is observed that, from scenarios I and IV, the recommended threshold values for dividing the IMFs into two components are ρyt,yt−1=0.90,0.95, and from scenario III, the recommended values are ρyt,yt−1=0.50,0.90,0.95, whereas from scenario II, the recommended value is ρyt,yt−1=0.95. So, it is observed from all scenarios that the only one value which significantly divided the IMFs into two components is ρyt,yt−1=0.95. Thus, the recommended value of autocorrelation for dividing the IMFs into two components is ρyt,yt−1=0.95. The summary of the threshold value determination is presented in Table 8.Table 7
p values of correlation difference test of true correlation and among the components obtained from the reconstruction of IMFs.
Dataρyt,yt−1Scenario IScenario IIScenario IIIScenario IV5000.5000.7600.4900.550—0.9000.8500.3500.2400.7300.9500.4000.3500.2400.8400.9900.400<0.0010.1000.84010000.5000.5200.6900.360—0.9000.5900.7500.2100.6400.9500.5900.8100.1500.5200.990<0.001<0.0010.0480.52020000.5001.0000.5200.4100.4600.9000.7000.8500.2300.6000.9500.7000.8500.2300.4700.9900.6600.9000.0080.47050000.500<0.001—0.7200.9500.9000.0400.1800.7500.8100.9500.1200.1800.7500.8100.9900.1300.0500.1500.13010,0000.5000.140—0.2000.0200.9000.2900.0300.1600.1000.9500.5200.0900.1600.1000.9900.7100.190<0.0010.001The bold values indicate that they are statistically significant at 5% and communicated that the true correlation and among the components obtained from the reconstruction of IMFs are not independent.Table 8
Summary of the threshold values for all four scenarios.
Level of significance (%)ρyt,yt−1Scenario IScenario IIScenario IIIScenario IV50.500✗✗✓✗0.900✗✗✓✓0.950✓✓✓✓0.990✗✗✗✗10.500✗✗✓✗0.900✓✓✓✓0.950✓✓✓✓0.990✗✗✗✗✓The threshold value is recommended. ✗The threshold value is not recommended.The only recommended value at both levels of significance is0.95. Thus, for division of IMFs into SD components, the threshold value of 0.95 is recommended and could be used. In Section 3.5, the real-world applications are tested by utilizing the recommended threshold value of 0.95.


## 3.5. IMFs Reconstruction

In this study, the reconstruction of IMFs is proposed based on the autocorrelation of everykth IMF obtained from EEMD. Each IMF autocorrelation is compared with a threshold value of ρyt,yt−1=0.95. An IMF is considered stochastic if its autocorrelation is less than the threshold value; otherwise, it is considered as deterministic. For comparison, the IMFs are also reconstructed through AMI. The AMI is simply the visual inspection of the plots which are obtained for all IMFs of each dataset. The detailed analysis of every dataset is followed in the subsequent sections.From Section3.3, it is detected that the total number of IMFs for daily datasets is 12 plus residue. For simplicity, the residue term is replaced with the IMF, so the total number of IMFs for both datasets is 13. For weekly data, the total number of IMFs is 11 for Brent, while for WTI, it is 10. Autocorrelation of all IMFs is computed and presented in Table 9.Table 9
Autocorrelation of Brent and WTI daily and weekly data.
IMFsDailyWeeklyBrentWTIBrentWTIIMF.1−0.467−0.484−0.479−0.448IMF.20.4730.4880.5780.592IMF.30.8360.8500.8890.892IMF.40.9530.9540.9830.987IMF.50.9870.9880.9950.997IMF.60.9970.9980.9990.999IMF.70.9990.9990.9990.999IMF.80.9991.0001.0001.000IMF.91.0001.0001.0001.000IMF.101.0001.0001.0001.000IMF.111.0001.0001.000—IMF.121.0001.000——IMF.131.0001.000——It is observed from Table9 that the autocorrelation of the IMF.1, IMF.2, and IMF.3 for all datasets is less than 0.95. So, the first three IMFs are considered as stochastic in all datasets, while the rest is considered as deterministic. So, to make a SC, all stochastic IMFs are summed, and for DC, all deterministic IMFs are summed. Next, to obtain the SD components from the AMI method for Brent daily data, their plots are presented in Figure 5.Figure 5
AMI plots of Brent daily data for all IMFs.From Figure5, it is observed that the AMI plots have the same pattern after the 6th IMF, so the first six IMFs are considered as stochastic, and the remainder, as deterministic. To make two components, add the first six IMFs for SC and from seven to thirteen for DC. Next, to obtain the SD components of WTI daily data, the plots of AMI for all IMFs are presented in Figure 6.Figure 6
AMI plots of WTI daily data for all IMFs.From Figure6, it is observed that the AMI plots have the same pattern after the 6th IMF, so the first six IMFs are considered as stochastic, and the remainder, as deterministic. To make two components, add the first six IMFs for SC and from seven to thirteen for DC. Next, to obtain the SD components of Brent weekly data, the plots of AMI for all IMFs are presented in Figure 7.Figure 7
AMI plots of Brent weekly data for all IMFs.From Figure7, it is observed that the AMI plots have the same pattern after the 6th IMF, so the first six IMFs are considered as stochastic, and the remainder, as deterministic. To make two components, add the first six IMFs for SC and from seven to eleven for DC. Next, to obtain the SD components of WTI weekly data, the plots of AMI for all IMFs are presented in Figure 8.Figure 8
AMI plots of WTI weekly data for all IMFs.From Figure8, it is observed that the AMI plots have the same pattern after the 6th IMF, so the first six IMFs are considered as stochastic, and the remainder, as deterministic. To make two components, add the first six IMFs for SC and from seven to ten for DC. The next is the correlation test that the SD components are independently obtained from the reconstruction of IMFs.From Table10, it is observed that, for all reconstructed components, the hypothesis of zero correlation is not rejected, meaning that the division of SD components is independent at a 5% level of significance. Moreover, there is a theorem, put forward by Takens [55], known as the embedding theorem, stating that a time series can be reconstructed in vectors with m values (described as the embedding dimension). Each of these values corresponds to an observation that is spaced out in intervals, following a time delay (or separation dimension) called τ. As numerous studies have shown, these vectors represent the interdependent relationships between observations, and since they increase the accuracy of the modelling, they also enhance the accuracy of the prediction [56]. The plots of the AMI of the SD components are presented in Figure 9, showing the interdependency among the SD components.Table 10
Test of correlation among the two components for daily and weekly data.
DataMethodCorrelationt statisticd.fp valueBrentDailyR0.0040.534124340.592AMI0.0101.112124340.266WeeklyR0.0160.79124860.429AMI0.0261.29124860.197WTIDailyR0.0090.87390580.383AMI0.0100.92290580.356WeeklyR0.0170.76018100.447AMI0.0350.50118100.133Alternative hypothesis: true correlation is not equal to zero.Figure 9
AMI plots of SD components obtained from the proposed approach. (a) Daily Brent. (b) Daily WTI. (c) Weekly Brent. (d) Weekly WTI.
(a)(b)(c)(d)From Figures9(a)–9(d), it is observed that the division of IMFs into SD components is independently distributed. After obtaining the two components, the next step is the appropriate model selection for every component using the same steps discussed above. Every model is selected based on the LB test for which the p values are greater than 0.05 for both daily and weekly datasets, all IMFs, and components so that all the fitted models are adequate and could be used for forecasting. In the next step, all the selected models are used for forecasting, and the last 20% data of every series were retained as the testing set forecasted. For forecasting of ARIMA models, r software is used. For ANN, the data are transferred to MATLAB, and ntstool is used for forecasting every component. After forecasting the series, the accuracy measures are calculated and presented in Tables 3 and 4 for both daily and weekly series.


## 4. Discussion

The performance of the proposed model EEMD-R-ANN is evaluated using four quantitative measures including RMSE, MAPE, MAE, and DS. The results of all these methods for all models are obtained and presented in Table3 for both markets of daily COPs. The benchmark models include the well-known single ARIMA and ANN models, the hybrid models including EEMD-ARIMA and EEMD-ANN, the hybrid models using only the SD components obtained from the proposed method of reconstruction including EEMD-SD-ANN and EEMD-SD-ARIMA, models using AMI for reconstruction of IMFs including EEMD-AMI-ARIMA, and EEMD-AMI-ANN, and the last is the proposed method of R of IMFs based on autocorrelation which are denoted by EEMD-R-ARIMA and EEMD-R-ANN.The RMSE is computed for all models included in this study and shown in Tables3 and 4 for both markets and all datasets of COPs and plotted in Figures 10 and 11 for a clearer picture. In all models, EEMD-R-ANN and EEMD-R-ARIMA almost achieved the same and the best (lowermost) values for both the COPs series. The reconstructed hybrid models EEMD-AMI-ARIMA and EEMD-AMI-ANN achieved nearly the similar values but high; EEMD-SD-ANN, EEMD-SD-ARIMA, and the hybrid models EEMD-ANN and EEMD-ARIMA achieved the lowermost values but higher than the EEMD-R-ANN. The single models ARIMA and ANN achieved the worst (highest) values among all considered models. Therefore, from the RMSE perspective, the proposed model EEMD-R-ANN confirmed their usefulness.Figure 10
Forecasting accuracy measures of daily Brent and WTI data.Figure 11
Forecasting accuracy measures of weekly Brent and WTI data.As another level performance of forecasting, the MAE values of all the models on WTI and Brent (daily and weekly) COPs are shown in Tables3 and 4 and plotted in Figures 10 and 11 for clear picture. It can be seen that EEMD-R-ANN outperformed all other models, and their interpretations are not different from RMSE for both markets and all datasets. Hence, from the MAE viewpoint, the proposed model EEMD-R-ANN confirmed their usefulness by attaining the lowest values.The MAPE is another level performance of forecasting. The MAPE values of all models on both markets WTI and Brent (daily and weekly) are presented in Tables3 and 4 and plotted in Figures 10 and 11 for clear picture. From Tables 3 and 4, the model EEMD-R-ANN attained the lowermost values for both markets and all datasets (daily and weekly). The three models EEMD-ANN, EEMD-ARIMA, and EEMD-R-ARIMA attained nearly similar values, whereas the other two models EEMD-SD-ANN and EEMD-SD-ARIMA attained nearly the similar values but higher than the preceding three models. The model’s ANN, ARIMA, EEMD-AMI-ANN, and EEMD-AMI-ARIMA achieved nearly similar values but highest and ranked high among all models. Therefore, from the MAPE perspective, the proposed model EEMD-R-ANN confirmed their usefulness by attaining the lowermost values.The model EEMD-R-ANN produced the lowest values of MAPE for Brent and WTI daily COPs which are 0.542% and 0.572%, respectively, and attained the category of highly accurate forecasts [57]. The values of the MAPE for Brent and WTI daily COPs of three models which are EEMD-R-ARIMA, EEMD-ANN, and EEMD-ARIMA in the range of 0.559% through 0.697% and achieved the same category as attained by the model EEMD-R-ANN but with high values. The values of the MAPE for models EEMD-SD-ARIMA, EEMD-SD-ANN, EEMD-AMI-ARIMA, EEMD-AMI-ANN, ARIMA, and ANN were in the range of 1.018% to 1.489% of both COPs and fall under the category of good accurate forecasts. Thus, from the MAPE perspective, the best model for obtaining highly accurate forecasts is EEMD-R-ANN for daily COPs of both markets, i.e., Brent and WTI.For weekly data, the model EEMD-R-ANN produced the lowest values of MAPE for Brent and WTI COPs, which are 0.933% and 0.936%, respectively, and attained the class of highly accurate forecasts [57]. The values of the MAPE for Brent and WTI weekly COPs of the model EEMD-R-ARIMA were 0.936% and 0.952% and achieved the same category as attained by the EEMD-R-ANN model but with high values. The values of the MAPE for models ANN, ARIMA, EEMD-SD-ARIMA, EEMD-SD-ANN, EEMD-AMI-ARIMA, EEMD-AMI-ANN, EEMD-ARIMA, and EEMD-ANN were in the range of 1.018% to 1.489% of both COPs and fall under the class of good accurate forecasts. Thus, from the MAPE perspective, the best model for obtaining highly accurate forecasts is EEMD-R-ANN for weekly COPs of both markets, i.e., Brent and WTI.The next forecast accuracy measure is the directional forecasts which are much more important for policymakers and investors since they look to the market trend either the COPs are going up or down. The DS values are presented in Tables3 and 4 for daily and weekly data, respectively, and plotted in Figure 12 for both markets. The DS values of the ANN and ARIMA models of both daily and weekly COPs for both markets are in the range of 51.29%–55.55% and almost equal to the random guess. Hence, it is very challenging to achieve the directional forecasts at a reasonable level for single models, and the reason may be the complex and compound nature of COPs. The next is the decomposition and ensemble model that used all IMFs obtained from EEMD of both COPs daily and weekly which are in the range of 80.17%–83.38%. Thus, decomposition and ensemble models achieved a higher forecasting accuracy and having a good capability of producing more accurate directional forecasts in comparison with single models. But the reason with these models is that they are too time-consuming. To cope with this issue, the idea of reconstruction of IMFs is introduced. So, the next is the directional forecasting performance of the models which used the reconstructed IMFs. The forecasting performance of the models which used the reconstructed IMFs obtained from AMI is in the range of 55.81%–62.47%. These results are better than the single models but not far from the random guess. The next is the forecasting performance of the models which used the proposed method of reconstruction of IMFs, and their DS values are in the range of 75.01%–78.39% and showed some satisfactory directional forecasts. In the last, the directional forecasts of the proposed models EEMD-R-ARIMA and EEMD-R-ANN are in the range of 88.25%–94.46% for both COPs as well as for both markets which are the highest among all models. Thus, the proposed reconstruction of the IMFs method attained the highest values for directional forecasts and outperformed all other benchmark models including ARIMA, ANN, EEMD-ARIMA, EEMD-ANN, EEMD-AMI-ARIMA, EEMD-AMI-ANN, EEMD-SD-ARIMA, and EEMD-SD-ANN. The model EEMD-R-ANN relatively performed better than the EEMD-R-ARIMA, so the highly recommended model regarding the DS values is EEMD-R-ANN for both COPs and as well as for both markets.Figure 12
DS plot of different models.Thus, the proposed method of IMFs reconstruction through autocorrelation significantly improved the performance of the well-known ARIMA and ANN models. Thus, from this study, the suggested model is EEMD-R-ANN which produced the highest values for DS and lowermost values for RMSE, MAPE, and MAE for both Brent and WTI (daily and weekly) datasets.Furthermore, for visual inspection, the forecasting accuracy measures are also plotted and presented in Figures10–12 for both Brent and WTI daily and weekly datasets.To check the superiority of the proposed model EEMD-R-ANN, the DM test is also conducted. The test results of DM statistics with theirp values are shown in Table 11 for both Brent and WTI COPs. The output of the DM test statistically validated the conclusion drawn from the above forecasting accuracy measures. The proposed model EEMD-R-ANN statistically outperformed ARIMA, ANN, EEMD-ANN, EEMD-ARIMA, EEMD-AMI-ANN, EEMD-AMI-ARIMA, EEMD-SD-ANN, and EEMD-SD-ARIMA models with their respective p values which were far below than 0.01 for both daily and weekly data as well as for both Brent and WTI COPs. However, the DM test statistic values were not statistically significant for EEMD-R-ANN and EEMD-R-ARIMA models for both markets and all datasets which confirm the usefulness of the proposed reconstruction strategy. The above analysis and discussions are summarized in Section 4.1.Table 11
DM test results of Brent and WTI daily and weekly data.
ModelEEMD-R-ARIMAEEMD-SD-ANNEEMD-SD-ARIMAEEMD-AMI-ANNEEMD-AMI-ARIMAEEMD-ANNEEMD-ARIMAANNARIMABrent dailyEEMD-R-ANN−1.38 (0.17)−21.14∗ (<0.01)−19.12∗ (<0.01)−19.20∗ (<0.01)−17.39∗ (<0.01)−14.21∗ (<0.01)−14.01∗ (<0.01)−19.1∗ (<0.01)−18.60∗ (<0.01)WTI dailyEEMD-R-ANN−1.28 (0.20)−17.61∗ (<0.01)−15.98∗ (<0.01)−19.10∗ (<0.01)−18.01∗ (<0.01)−7.39∗ (<0.01)−6.81∗ (<0.01)−16.2∗ (<0.01)−14.91∗ (<0.01)Brent weeklyEEMD-R-ANN−1.05 (0.29)−12.39∗ (<0.01)−11.89∗ (<0.01)−10.86∗ (<0.01)−9.91∗ (<0.01)−5.44∗ (<0.01)−4.98∗ (<0.01)−9.14∗ (<0.01)−8.90∗ (<0.01)WTI weeklyEEMD-R-ANN−1.01 (0.31)−11.44∗ (<0.01)−10.01∗ (<0.01)−9.33∗ (<0.01)−8.71∗ (<0.01)−4.67∗ (<0.01)−4.01∗ (<0.01)−8.24∗ (<0.01)−7.88∗ (<0.01)∗Significant at 1%.


### 4.1. Summary of the Work

From the above discussions and analysis, some important findings are summarized as follows:(1)
To accurately forecast the COPs, it is hard to attain satisfactory results when using single models (i.e., ARIMA and ANN) because of the nonstationary and nonlinearity structure of the data.(2)
The models which used all IMFs relatively performed well than the models which used the reconstructed IMFs into two components, namely, stochastic and deterministic.(3)
The recommended threshold value of the autocorrelationρyt,yt−1=0.95 based on synthetic data was used for the reconstruction of IMFs on the four real datasets. The analysis has proved that the reconstruction of IMFs through autocorrelation is a better and simple strategy that enhanced the performance of all models.(4)
With the benefits of EEMD, reconstruction, ARIMA, and ANN, the proposed ensemble models EEMD-R-ANN and EEMD-R-ARIMA significantly outperform all other models listed in this paper in terms of MAE, RMSE, MAPE, DS, and DM test.(5)
The performance of the proposed models EEMD-R-ANN and EEMD-R-ARIMA is statistically insignificant in terms of the DM test. However, in terms of MAE, RMSE, MAPE, and DS, the model EEMD-R-ANN relatively performed well compared to EEMD-R-ARIMA. Thus, the suggested model for forecasting COPs is EEMD-R-ANN.(6)
The experimental results demonstrated that the reconstruction of IMFs into SD components using the proposed method was effective but modelling the IMFs being part of the stochastic component separately was more effective and powerful approach for forecasting COPs.(7)
As one of the most important commodities in the world, we analyze the COPs for illustration and verification of the proposed method. The empirical results showed that our proposed reconstruction of decomposition and ensemble model-based analysis approach is vital and effective and could be tested for more complex tasks.


## 4.1. Summary of the Work

From the above discussions and analysis, some important findings are summarized as follows:(1)
To accurately forecast the COPs, it is hard to attain satisfactory results when using single models (i.e., ARIMA and ANN) because of the nonstationary and nonlinearity structure of the data.(2)
The models which used all IMFs relatively performed well than the models which used the reconstructed IMFs into two components, namely, stochastic and deterministic.(3)
The recommended threshold value of the autocorrelationρyt,yt−1=0.95 based on synthetic data was used for the reconstruction of IMFs on the four real datasets. The analysis has proved that the reconstruction of IMFs through autocorrelation is a better and simple strategy that enhanced the performance of all models.(4)
With the benefits of EEMD, reconstruction, ARIMA, and ANN, the proposed ensemble models EEMD-R-ANN and EEMD-R-ARIMA significantly outperform all other models listed in this paper in terms of MAE, RMSE, MAPE, DS, and DM test.(5)
The performance of the proposed models EEMD-R-ANN and EEMD-R-ARIMA is statistically insignificant in terms of the DM test. However, in terms of MAE, RMSE, MAPE, and DS, the model EEMD-R-ANN relatively performed well compared to EEMD-R-ARIMA. Thus, the suggested model for forecasting COPs is EEMD-R-ANN.(6)
The experimental results demonstrated that the reconstruction of IMFs into SD components using the proposed method was effective but modelling the IMFs being part of the stochastic component separately was more effective and powerful approach for forecasting COPs.(7)
As one of the most important commodities in the world, we analyze the COPs for illustration and verification of the proposed method. The empirical results showed that our proposed reconstruction of decomposition and ensemble model-based analysis approach is vital and effective and could be tested for more complex tasks.


## 5. Conclusion

This paper proposed a new approach for the reconstruction of IMFs of decomposition and ensemble model according to their influences, i.e., stochastic and deterministic. We considered four different synthetic series to fix a value of autocorrelation for reconstruction of IMFs with minimum human interference. All these analysis confirm that the new approach supports the fast selection of reconstruction of IMFs and helps researchers to make predictions in situations where there are time constraints. In that sense, we believe that our approach can improve the most complex stage of reconstruction of IMFs according to their influences. Hence, the reconstruction of IMFs based on autocorrelation is one of the best methods to reconstruct the IMFs according to their influences (Section3.5). Thus, from this study, the threshold value of 0.95 of autocorrelation is recommended for separating the IMFs into SD components. To evaluate the performance of the proposed approach, four real-world COPs series were considered. The experimental findings demonstrated that all methodologies including two single benchmark models and eight ensemble models were effective. However, the forecasting accuracy measures in terms of MAE, RMSE, MAPE, and DS highlighted that the models EEMD-R-ANN and EEMD-R-ARIMA based on reconstruction of IMFs through autocorrelation were the most efficient methods for forecasting COPs. Moreover, the performance of the model EEMD-R-ANN is slightly better than the EEMD-R-ARIMA model. Thus, finally, the recommended model for forecasting the world COPs is EEMD-R-ANN. Hence, the reconstruction of IMFs based on autocorrelation is the best method to reconstruct the IMFs according to their influences.The advantages of the new proposed reconstruction approach were the results of handling the stochastics IMFs separately in which the linear and nonlinear parts of the ARIMA and ANN models were combined to handle the stochastic uncertainty. Furthermore, the proposed method of reconstruction of IMFs has also overcome the problem of supervising the method in which the reconstruction is carried out through visual inspection using the AMI plots, while in the proposed method, the reconstruction of IMFs is done through the autocorrelation of each IMF from which every researcher can obtain the same number of reconstructed IMFs using the same data with high forecasting accuracy and less computational time.In future, the work could be extended in two aspects: (1) studying the stochastic IMFs in more detail to improve the performance on forecasting COPs and (2) applying the EEMD-R-ANN to forecast other time series of energy, such as wind speed and electricity price.



---

*Source: 1325071-2020-10-06.xml*