# International Gold Price Forecast Based on CEEMDAN and Support Vector Regression with Grey Wolf Algorithm

Authors: Wanbo Lu; Tingting Qiu; Wenhui Shi; Xiaojun Sun
Journal: Complexity
 (2022)

---

## Abstract

Considering the complexity pattern of the gold price, this paper adopts the decomposition-reconstruction-forecast-mergence scheme to perform the international gold price forecast. The original gold price data are decomposed into 12 intrinsic mode functions and a residual by the complete ensemble empirical mode decomposition with adaptive noise (CEEMDAN) method, and then the 13 sequences are reconstructed into a high-frequency subsequence (IMFH), a low-frequency subsequence (IMFL), and the residual (Res). According to the different characteristics of the subsequences, the IMFL and Res are forecasted by the support vector regression (SVR) model. Besides, in order to further improve the prediction accuracy of IMFH, we have developed a novel hybrid method based on the support vector regression (SVR) model and the grey wolf optimizer (GWO) algorithm with SVR for predicting the IMFH of gold prices, i.e., the CEEMDAN-GWO-SVR model. This hybrid model combines the methodology of complex systems with machine learning techniques, making it more appropriate for analyzing relationships such as high-frequency dependences and solving complex nonlinear problems. Finally, the final result is obtained by combining the forecasting results of the three subsequences. The empirical results show that the proposed model demonstrates the highest prediction ability among all of the investigated models in a comparison of prediction errors with other individual models.


---

## Body

## 1. Introduction

As a precious metal, gold also acts as a reserve. The central banks of many countries hold an amount of gold, for the purpose of mitigating potential risks, resisting inflation, and increasing economic growth. Gold plays an important role in reserve assets. Especially in the last few years, the international situation has undergone complex and deep changes due to the continuous escalation of the China-US trade war. In this context, it is of great significance to mitigate international financial risks and maintain economic stability. Central banks’ demand for gold has been increasing. In 2021, central banks added a total of 463 tons of gold, 82% more than in 2020, bringing global central bank gold reserves to the highest level in nearly 30 years.The trend of the gold price is different from other minerals [1], which is impacted by various factors, including demand and supply, inflation, and political issues, among others. In the financial crisis of 2008, most commodities and equities dropped about 40%, while the gold price remained at a high point and increased about 6% [2]. With the continuous adjustment of international patterns in recent years, gold is facing a more complex market environment [3], and its price has been subject to severe fluctuations. Figure 1 shows the movement of the gold price for 21 years from 2000 to 2021. Obviously, the price of gold fluctuates up and down, showing a high degree of nonlinearity and uncertainty. Therefore, an accurate gold price forecast is generally considered to be challenging [4, 5].Figure 1
Daily international gold price curve (January 3, 2000–October 21, 2022).The main purpose of this research is to contribute to the accurate prediction of gold price and movement. To this end, we adopt the decomposition-reconstruction-forecast-mergence scheme to perform the international gold price forecast, and propose the hybrid CEEMDAD-GWO-SVR model. Our experiments demonstrate that this scheme and the hybrid model can significantly improve the forecasting performance. The rest of this paper is organized as follows: in Section2, a series of recent works on gold price forecasting are presented, in Section 3, the CEEMDAN method, the GWO algorithm, and the SVR model are recalled for completeness, and the main steps of the hybrid CEEMDAD-GWO-SVR model are introduced, in Section 4, the experimental results are presented and forecasting accuracies are compared across different models, and in Section 5,the conclusions are discussed.


## 2. Literature Review

### 2.1. Forecasting Models

Research on gold price forecasting was performed for decades, and numerous methods have been proposed. The commonly used methods mainly include multiple regression, time-series analysis, neural network methods, and hybridization or combination models.The multiple regression model makes prediction of gold prices with multiple influencing factors. The time-series methods forecast the future gold price based on the past time series, including ARIMA (auto-regressive integrated moving average), the ARCH (auto-regressive Conditional Heteroscedasticity) class models, and the GARCH (Generalized Auto-Regressive Conditional Heteroscedasticity) class models among others. For example, Wang and Xia [6] analyzed the marginal and joint impact of economic variables on the gold price, and provided a short-term forecast of the gold price by multiple regression, achieving a satisfactory result. Yang [7] used the ARIMA model, and the results showed that the model can reflect the true gold price trend to a certain extent. However, the above two methods have certain limitations: the explanatory variables of the multiple regression are difficult to choose in empirical research, the time-series methods have too many parameters to be estimated, and the optimal lag order is hard to choose. In practice, these two methods may fail to characterize the nonlinear and complex behavior of the gold price time series and decrease the forecasting accuracy.The neural network methods are well suited to address nonlinear and nonstationary problem, and have been widely applied in time-series forecast. Parisi et al. [8] analyzed recursive and rolling neural network models to forecast one-step-ahead sign variations in gold price, and they found that the rolling ward networks perform well in forecasting gold price sign variation. Zhang and Ci [9] proposed a deep belief network (DBN) model for gold price forecast, and the empirical results showed that the proposed DBN model outperforms the BP (back propagation) neural network model, genetic algorithm optimization for BP model (GA-BP), and the ARIMA model. However, the neural network models have transparency problem and it is difficult to understand the economic implications and impact mechanisms [10].In recent years, research on gold price forecast has focused on two types of hybrid models and combined models, such as (1) the hybridization or combination of different models and (2) the hybridization or combination of models with metaheuristic algorithms or advanced intelligent approaches. Since not every model is suitable for hybridization or combination with others, the second type of hybridization has more applications than the first one. Most importantly, the hybridization of machine learning approaches achieved many successes and provided more useful findings about the gold price forecast. Weng et al. [11] proposed a novel online extreme learning machine (OS-ELM) model based on genetic algorithm (GA) and regularization to forecast gold price. The proposed hybrid model solves the singularity problem of OS-ELM model and outperforms the ARIMA, BP, ELM, and OS-ELM models. Alameer et al. [12] presented a hybrid model for accurately forecasting long-term monthly gold price fluctuations. The model uses whale optimization algorithm (WOA) as a trainer to learn the multilayer perceptron neural network (NN), and their empirical results indicate the superiority of hybrid WOA-NN model over the classic NN, particle swarm optimization for NN (PSO-NN), genetic algorithm for NN (GA-NN), and ARIMA models.Owing to its outstanding forecasting performance and interpretability, the support vector regression (SVR) model, especially as part of a hybrid model, has been extensively used in landslide displacement prediction [13–15], daily traffic peak flow management [16, 17], financial investment [18, 19], air quality prediction [20, 21], and gold price forecast [22, 23] as well. Compared with other machine learning methodologies, SVR shows rapid convergence [24] and efficiently solves multidimensional function estimation problems [25], but the selection of the SVR penalty parameter and the kernel function parameter has an important influence on its own performance [26]. Especially, SVR would have a large forecasting error in dealing with nonlinear high-frequency time series. Results in [26] showed that the hybrid models of SVR and metaheuristic algorithms can significantly improve the performance of the single SVR model.Metaheuristic algorithms are widely used by researchers for their remarkable ability to solve optimization problems. Since the commencement of the classical metaheuristics, namely, genetic [27], simulated annealing [28], tabu search [29], and particle swarm optimization [30] in the early 1970s to late 1990s, a new generation of metaheuristics has grown exponentially. In the past three years, nearly 20 metaheuristics has been proposed, including Monarch butterfly optimization [31], butterfly optimization algorithm [32], shuffled shepherd optimization algorithm [33], and water strider algorithm [34]. Currently, there are more than 200 metaheuristic techniques in the literature and several surveys [35–37] provided comprehensive overviews on them. Grey Wolf Optimizer (GWO) [38] is one of metaheuristic algorithms, which was proposed by Mirjalili in 2014. It is a swarm intelligent algorithm that mimics the leadership hierarchy and hunting mechanism of grey wolves in nature. It has the advantages of fast convergence speed and strong search ability and can simultaneously optimize the parameters and weights of the involved kernels [39]. Recent studies [38, 40, 41] have demonstrated that GWO provides competitive results compared to some well-known metaheuristics. Therefore, this paper hybridizes SVR with GWO algorithm (i.e., PSO-SVR) to forecast the nonlinear high-frequency part of the gold price time series and the subsequent experimental results verify the satisfactory learning precision and generalization ability with PSO-SVR.


### 2.2. Data Decomposition

Studies [42–45] have shown that data decomposition before forecasting can effectively reduce the nonstationary characteristics of the time series, and it can improve the forecasting performance significantly. Current decomposition algorithms include empirical mode decomposition (EMD) [42], ensemble empirical mode decomposition (EEMD) [43], complete ensemble empirical mode decomposition with adaptive noise (CEEMDAN) [44], and variational mode decomposition (VMD) [45]. EMD is an adaptive decomposition method for nonlinear and nonstationary processes. This method can decompose a complicated dataset into a finite and often small number of intrinsic mode functions (IMFs). However, EMD may encounter the mode-mixing problem caused by intermittent signals. EEMD effectively resolves this problem through an ensemble of the signal plus the white Gaussian noise. However, the added white Gaussian noise also creates some of the following new problems: the reconstructed signal includes residual noise and different realizations of signal plus noise may produce different number of IMFs [44]. To overcome these problems, CEEMDAN was proposed. It performs the EEMD over white Gaussian noise added at each stage of the decomposition and makes the decomposition complete. VMD is an entirely nonrecursive variational mode decomposition model, and the experimental results in the study [45] showed that VMD is more robust to sampling and noise than EMD and EEMD. However, the influence parameters of VMD need to be set in advance [46–48]. In summary, CEEMDAN not only provides a better spectral separation of the IMFs than EMD and EEMD but also overcomes the parameter selection limitation of VMD. In addition, CEEMDAN reduces the computational cost due to its fewer sifting iterations. Considering the importance of data decomposition and the ease of implementation of the CEEMDAN method, this paper uses CEEMDAN to preprocess the data before forecasting.


### 2.3. Main Contributions

In this study, we first adopt CEEMDAN to conduct the data preprocessing process and then use GWO-SVR model to forecast the price of gold. The proposed model, named CEEMDAN-GWO-SVR, shows better forecasting ability and performance than the state-of-the-art benchmark models in the numerical experiments. All in all, a list of abbreviations and variables has been added in Table1 in the appendix. The principal contributions of this paper are as follows:(1)
Although numerous methods have been proposed to forecast the gold price, especially the hybrid or combined models were performed for decades, there are few studies based on “decomposing first and then merging.” This paper adopts a decomposition-reconstruction-forecast-mergence scheme to forecast the international gold price, using the CEEMDAN method to decompose the original data into several IMFs and the residual and then usingT-test to reconstruct them into three subsequences, namely, high-frequency subsequence (IMFH), low-frequency subsequence (IMFL), and the residual (Res), which are, respectively, forecasted by SVR model and finally the three forecasting results are combined to the final result.(2)
Considering the different characteristics of the three subsequences and the large error of SVR when dealing with nonlinear high-frequency time series, different SVR models are used here. Specifically, for IMFL and Res, this paper uses SVR with a linear kernel function to forecast the price of gold, and for IMFH, the GWO algorithm and the SVR model with the Gaussian kernel function are hybridized to generate the GWO-SVR, which overcomes the problem of subjective selection of parameters to improve the accuracy of forecast.(3)
A comparison demonstrates that the CEEMDAN method can effectively capture effective information from the original gold price data, and the proposed CEEMDAN-GWO-SVR model provides significant improvements over the benchmark models.Table 1
List of abbreviations.
AbbreviationFull nameARCHAuto-regressive conditional heteroscedasticityARIMAAuto-regressive integrated moving averageBPBack propagationBPNNBack propagation neural networkCEEMDANComplete ensemble empirical mode decomposition with adaptive noiseDADirectional accuracyDBNDeep belief networkEEMDEnsemble empirical mode decompositionELMExtreme learning machineEMDEmpirical mode decompositionGAGenetic algorithmGARCHGeneralized auto-regressive conditional heteroscedasticityGWOGrey wolf optimizerIMFIntrinsic mode functionIMFHHigh-frequency subsequenceIMFLLow-frequency subsequenceMAPEMean absolute percentage errorNNNeural networkOS-ELMOnline extreme learning machinePSOParticle swarm optimizationResResidualRMSERoot mean square errorSVRSupport vector regressionVCRVariance contribution rateVMDVariational mode decompositionWOAWhale optimization algorithm


## 2.1. Forecasting Models

Research on gold price forecasting was performed for decades, and numerous methods have been proposed. The commonly used methods mainly include multiple regression, time-series analysis, neural network methods, and hybridization or combination models.The multiple regression model makes prediction of gold prices with multiple influencing factors. The time-series methods forecast the future gold price based on the past time series, including ARIMA (auto-regressive integrated moving average), the ARCH (auto-regressive Conditional Heteroscedasticity) class models, and the GARCH (Generalized Auto-Regressive Conditional Heteroscedasticity) class models among others. For example, Wang and Xia [6] analyzed the marginal and joint impact of economic variables on the gold price, and provided a short-term forecast of the gold price by multiple regression, achieving a satisfactory result. Yang [7] used the ARIMA model, and the results showed that the model can reflect the true gold price trend to a certain extent. However, the above two methods have certain limitations: the explanatory variables of the multiple regression are difficult to choose in empirical research, the time-series methods have too many parameters to be estimated, and the optimal lag order is hard to choose. In practice, these two methods may fail to characterize the nonlinear and complex behavior of the gold price time series and decrease the forecasting accuracy.The neural network methods are well suited to address nonlinear and nonstationary problem, and have been widely applied in time-series forecast. Parisi et al. [8] analyzed recursive and rolling neural network models to forecast one-step-ahead sign variations in gold price, and they found that the rolling ward networks perform well in forecasting gold price sign variation. Zhang and Ci [9] proposed a deep belief network (DBN) model for gold price forecast, and the empirical results showed that the proposed DBN model outperforms the BP (back propagation) neural network model, genetic algorithm optimization for BP model (GA-BP), and the ARIMA model. However, the neural network models have transparency problem and it is difficult to understand the economic implications and impact mechanisms [10].In recent years, research on gold price forecast has focused on two types of hybrid models and combined models, such as (1) the hybridization or combination of different models and (2) the hybridization or combination of models with metaheuristic algorithms or advanced intelligent approaches. Since not every model is suitable for hybridization or combination with others, the second type of hybridization has more applications than the first one. Most importantly, the hybridization of machine learning approaches achieved many successes and provided more useful findings about the gold price forecast. Weng et al. [11] proposed a novel online extreme learning machine (OS-ELM) model based on genetic algorithm (GA) and regularization to forecast gold price. The proposed hybrid model solves the singularity problem of OS-ELM model and outperforms the ARIMA, BP, ELM, and OS-ELM models. Alameer et al. [12] presented a hybrid model for accurately forecasting long-term monthly gold price fluctuations. The model uses whale optimization algorithm (WOA) as a trainer to learn the multilayer perceptron neural network (NN), and their empirical results indicate the superiority of hybrid WOA-NN model over the classic NN, particle swarm optimization for NN (PSO-NN), genetic algorithm for NN (GA-NN), and ARIMA models.Owing to its outstanding forecasting performance and interpretability, the support vector regression (SVR) model, especially as part of a hybrid model, has been extensively used in landslide displacement prediction [13–15], daily traffic peak flow management [16, 17], financial investment [18, 19], air quality prediction [20, 21], and gold price forecast [22, 23] as well. Compared with other machine learning methodologies, SVR shows rapid convergence [24] and efficiently solves multidimensional function estimation problems [25], but the selection of the SVR penalty parameter and the kernel function parameter has an important influence on its own performance [26]. Especially, SVR would have a large forecasting error in dealing with nonlinear high-frequency time series. Results in [26] showed that the hybrid models of SVR and metaheuristic algorithms can significantly improve the performance of the single SVR model.Metaheuristic algorithms are widely used by researchers for their remarkable ability to solve optimization problems. Since the commencement of the classical metaheuristics, namely, genetic [27], simulated annealing [28], tabu search [29], and particle swarm optimization [30] in the early 1970s to late 1990s, a new generation of metaheuristics has grown exponentially. In the past three years, nearly 20 metaheuristics has been proposed, including Monarch butterfly optimization [31], butterfly optimization algorithm [32], shuffled shepherd optimization algorithm [33], and water strider algorithm [34]. Currently, there are more than 200 metaheuristic techniques in the literature and several surveys [35–37] provided comprehensive overviews on them. Grey Wolf Optimizer (GWO) [38] is one of metaheuristic algorithms, which was proposed by Mirjalili in 2014. It is a swarm intelligent algorithm that mimics the leadership hierarchy and hunting mechanism of grey wolves in nature. It has the advantages of fast convergence speed and strong search ability and can simultaneously optimize the parameters and weights of the involved kernels [39]. Recent studies [38, 40, 41] have demonstrated that GWO provides competitive results compared to some well-known metaheuristics. Therefore, this paper hybridizes SVR with GWO algorithm (i.e., PSO-SVR) to forecast the nonlinear high-frequency part of the gold price time series and the subsequent experimental results verify the satisfactory learning precision and generalization ability with PSO-SVR.


## 2.2. Data Decomposition

Studies [42–45] have shown that data decomposition before forecasting can effectively reduce the nonstationary characteristics of the time series, and it can improve the forecasting performance significantly. Current decomposition algorithms include empirical mode decomposition (EMD) [42], ensemble empirical mode decomposition (EEMD) [43], complete ensemble empirical mode decomposition with adaptive noise (CEEMDAN) [44], and variational mode decomposition (VMD) [45]. EMD is an adaptive decomposition method for nonlinear and nonstationary processes. This method can decompose a complicated dataset into a finite and often small number of intrinsic mode functions (IMFs). However, EMD may encounter the mode-mixing problem caused by intermittent signals. EEMD effectively resolves this problem through an ensemble of the signal plus the white Gaussian noise. However, the added white Gaussian noise also creates some of the following new problems: the reconstructed signal includes residual noise and different realizations of signal plus noise may produce different number of IMFs [44]. To overcome these problems, CEEMDAN was proposed. It performs the EEMD over white Gaussian noise added at each stage of the decomposition and makes the decomposition complete. VMD is an entirely nonrecursive variational mode decomposition model, and the experimental results in the study [45] showed that VMD is more robust to sampling and noise than EMD and EEMD. However, the influence parameters of VMD need to be set in advance [46–48]. In summary, CEEMDAN not only provides a better spectral separation of the IMFs than EMD and EEMD but also overcomes the parameter selection limitation of VMD. In addition, CEEMDAN reduces the computational cost due to its fewer sifting iterations. Considering the importance of data decomposition and the ease of implementation of the CEEMDAN method, this paper uses CEEMDAN to preprocess the data before forecasting.


## 2.3. Main Contributions

In this study, we first adopt CEEMDAN to conduct the data preprocessing process and then use GWO-SVR model to forecast the price of gold. The proposed model, named CEEMDAN-GWO-SVR, shows better forecasting ability and performance than the state-of-the-art benchmark models in the numerical experiments. All in all, a list of abbreviations and variables has been added in Table1 in the appendix. The principal contributions of this paper are as follows:(1)
Although numerous methods have been proposed to forecast the gold price, especially the hybrid or combined models were performed for decades, there are few studies based on “decomposing first and then merging.” This paper adopts a decomposition-reconstruction-forecast-mergence scheme to forecast the international gold price, using the CEEMDAN method to decompose the original data into several IMFs and the residual and then usingT-test to reconstruct them into three subsequences, namely, high-frequency subsequence (IMFH), low-frequency subsequence (IMFL), and the residual (Res), which are, respectively, forecasted by SVR model and finally the three forecasting results are combined to the final result.(2)
Considering the different characteristics of the three subsequences and the large error of SVR when dealing with nonlinear high-frequency time series, different SVR models are used here. Specifically, for IMFL and Res, this paper uses SVR with a linear kernel function to forecast the price of gold, and for IMFH, the GWO algorithm and the SVR model with the Gaussian kernel function are hybridized to generate the GWO-SVR, which overcomes the problem of subjective selection of parameters to improve the accuracy of forecast.(3)
A comparison demonstrates that the CEEMDAN method can effectively capture effective information from the original gold price data, and the proposed CEEMDAN-GWO-SVR model provides significant improvements over the benchmark models.Table 1
List of abbreviations.
AbbreviationFull nameARCHAuto-regressive conditional heteroscedasticityARIMAAuto-regressive integrated moving averageBPBack propagationBPNNBack propagation neural networkCEEMDANComplete ensemble empirical mode decomposition with adaptive noiseDADirectional accuracyDBNDeep belief networkEEMDEnsemble empirical mode decompositionELMExtreme learning machineEMDEmpirical mode decompositionGAGenetic algorithmGARCHGeneralized auto-regressive conditional heteroscedasticityGWOGrey wolf optimizerIMFIntrinsic mode functionIMFHHigh-frequency subsequenceIMFLLow-frequency subsequenceMAPEMean absolute percentage errorNNNeural networkOS-ELMOnline extreme learning machinePSOParticle swarm optimizationResResidualRMSERoot mean square errorSVRSupport vector regressionVCRVariance contribution rateVMDVariational mode decompositionWOAWhale optimization algorithm


## 3. Materials and Methods

### 3.1. Complete Ensemble Empirical Mode Decomposition with Adaptive Noise (CEEMDAN)

Based on EMD and EEMD, the CEEMDAN method adds white Gaussian noise at each stage of the decomposition and then computes the unique residue to obtain each IMF of the original dataset. The original signal can be exactly reconstructed by summing the IMFs owing to the complete decomposition of this method. The detailed algorithm is as follows:Step 1. We add white Gaussian noiseωit (i=1,⋯,I) to the original dataset xt and then generate the new dataset(1)xit=xt+ε0ωit,where i is the number of the trials and ε0 is the signal-to-noise ratio.Step 2. Eachxit (i=1,⋯,I) is decomposed by EMD by getting their first modes IMF1it and computing the first IMF of xt(2)IMF1˜t=1I.∑iIMF1it.Step 3. We calculate the first residual,r1t, as(3)r1t=xt−IMF1˜t.Step 4. We define the operatorEj., which produces the j-th mode by EMD and then decomposes r1t+ε1E1ωit (i=1,⋯,I) until we get the second IMF,(4)IMF2˜t=1I.∑iE1r1t+ε1E1ωit.Step 5. We calculate thek-th (k=2,⋯,K) residue and rkt, in equation (5), where K is the total number of IMF(5)rkt=rk−1t−IMFk˜t.Step 6. We decomposerkt+εkEkωit (i=1,⋯,I) by EMD and compute the k+1-th IMF as(6)IMFk+1˜t=1I.∑iE1rkt+εkEkωit.Step 7. We repeated Step 4 to 6 until the residue has less than two extrema and cannot be decomposed any longer, so the final residue can be expressed as(7)Rt=xt−∑kIMFk˜t,where the original dataset xt is decomposed(8)xt=∑kIMFk˜t+Rt.


### 3.2. Grey Wolf Optimizer (GWO)

GWO is metaheuristic algorithm, which mimicks the leadership hierarchy and hunting mechanism of grey wolves in nature. Four types of grey wolves, includingα, β, δ, and ω, are employed for simulating the social dominant hierarchy. α wolves are the leaders and the decision makers in the pack. β wolves are the second group in the hierarchy that help α wolves in decision-making and establish a discipline for the pack. δ wolves obey α and β wolves, but dominate ω wolves. They play the roles of scouts, sentinels, elders, hunters, and caretakers. ω wolves, the lowest in hierarchy, always have to obey all the other types of wolves. The GWO algorithm mathematically models this social hierarchy and hunting process to collect randomly preferred solutions. Specifically, α, β, δ, and ω are considered as the fittest, the second preferred, the third preferred, and the rest solutions, respectively. The GWO algorithm contains the following three steps: tracking, encircling, and attacking the prey.In the encircling, grey wolves encircle prey during the hunt, and the mathematical formulation is given by(9)X⟶t+1=X⟶pt−A⟶.C⟶.X⟶pt−X⟶t,where t is the number of the current iteration, X⟶ and X⟶p are the positions of a grey wolf and the prey, respectively, A⟶=2a⟶.r⟶1−a⟶ and C⟶=2r⟶2 are coefficient vectors, components of a⟶ are linearly decreased from 2 to 0, and r⟶1 and r⟶2 are random vectors in [0, 1].In hunting, it is supposed that theα, β, and δ have better knowledge about the potential location of prey. Therefore, the current first three best solutions are saved, and the other search agents (including the ω) update their positions according to the position of the best search agent. The formulas are given by equations (10)–(13):(10)X⟶1=X⟶α−A⟶1.C⟶1.X⟶α−X⟶,(11)X⟶2=X⟶β−A⟶2.C⟶2.X⟶β−X⟶,(12)X⟶3=X⟶δ−A⟶3.C⟶3.X⟶δ−X⟶,(13)X⟶t+1=X⟶1+X⟶2+X⟶33.The grey wolf agents diverge in search of the prey and converge for attacking the prey, which can be mathematical modelled by the value ofA⟶. The grey wolf agents tend to diverge from the prey when A⟶>1 and converge towards the prey when A⟶<1. Finally, the GWO algorithm is terminated by minimizing of the fitness function, which would be discussed in Section 3.4.


### 3.3. Support Vector Regression (SVR)

SVR is a supervised machine learning algorithm. Based on structural risk minimization, SVR can solve nonlinear problems with high dimensionality and small samples. For a set of training data withN observations x1,y1,x2,y2,⋯,xN,yN, where xi∈ℝn are input vectors, yi∈ℝ are output vectors, and the SVR function can be described as the following linear regression function:(14)fx=wTΦx+b,where Φx is the feature mapping function and the parameters w and b are the weight and bias, respectively. Based on the structural risk minimization, w and b can be calculated as follows:(15)Min:12w2+C∑i=1Nξi+ξi∗,s.t.yi−wTΦx+b≤ε+ξi∗,wTΦx+b−yi≤ε+ξi∗,ξi,ξi∗≥0,where C>0 is a cost function, ξi and ξi∗ are slack variables, and ε is an insensitive loss function. By introducing the Lagrange multipliers αi and αi∗, the SVR function is expressed as(16)fx=∑i=1Nαi−αi∗kx,xi+b,where kx,xi=Φx.Φxi is the kernel function, which is the inner product of the feature mapping function.SVR would have a large forecasting error when dealing with nonlinear high-frequency time series. Therefore, unlike the linear kernel function used in forecasting the low-frequency time series, the Gaussian kernel,kx,xi=exp−γx−xi2, is used to forecast the high-frequency time series, owing to its advanced mapping capability and stable forecasting performance. Meanwhile, the selection of the SVR penalty parameter C and the kernel function parameter γ has an important influence on its forecasting accuracy [26]. Thus, this paper tunes these SVR parameters using the GWO algorithm.


### 3.4. The Hybrid CEEMDAD-GWO-SVR Model

The details of the hybrid CEEMDAD-GWO-SVR model are as follows, and the flowchart is presented in Figure2. It should be noted that the “decomposition-reconstruction-forecast-mergence” scheme is also shown in Figure 2.Step 1. We conduct CEEMDAN to decompose the original gold price data into several relatively simple IMFs and the residual by using equations (1)–(8). This step corresponds to the “decomposition” of the scheme.Step 2. We reconstruct the IMFs and the residual into three subsequences, namely, IMFH, IMFL, and Res. Specifically, since the high-frequency IMFs are affected by some short-term random factors, they can be assumed to fluctuate around the zero mean. Therefore, each IMF is reconstructed into high-frequency subsequence (IMFH) and low-frequency subsequence (IMFL) byT-test with zero mean. The residual is classified as the third subsequence, i.e., Res. This step corresponds to the “reconstruction” of the scheme.Step 3. Forecasting IMFH by GWO-SVR and forecasting IMFL and Res by SVR. The optimal GWO-SVR model is determined by minimizing the fitness function. In this paper, the mean square error, given by equation (17), is used to calculate the forecasting accuracy. This step is decisive for the prediction performance of the proposed model, corresponding to the “forecast” of the scheme.(17)MSE=1n∑i=1Nyi−y^i2.Step 4. This is the last part of the scheme. By summing up the forecasting results of these three subsequences, the final result is obtained.Figure 2
The flowchart of CEEMDAD-GWO-SVR model.


## 3.1. Complete Ensemble Empirical Mode Decomposition with Adaptive Noise (CEEMDAN)

Based on EMD and EEMD, the CEEMDAN method adds white Gaussian noise at each stage of the decomposition and then computes the unique residue to obtain each IMF of the original dataset. The original signal can be exactly reconstructed by summing the IMFs owing to the complete decomposition of this method. The detailed algorithm is as follows:Step 1. We add white Gaussian noiseωit (i=1,⋯,I) to the original dataset xt and then generate the new dataset(1)xit=xt+ε0ωit,where i is the number of the trials and ε0 is the signal-to-noise ratio.Step 2. Eachxit (i=1,⋯,I) is decomposed by EMD by getting their first modes IMF1it and computing the first IMF of xt(2)IMF1˜t=1I.∑iIMF1it.Step 3. We calculate the first residual,r1t, as(3)r1t=xt−IMF1˜t.Step 4. We define the operatorEj., which produces the j-th mode by EMD and then decomposes r1t+ε1E1ωit (i=1,⋯,I) until we get the second IMF,(4)IMF2˜t=1I.∑iE1r1t+ε1E1ωit.Step 5. We calculate thek-th (k=2,⋯,K) residue and rkt, in equation (5), where K is the total number of IMF(5)rkt=rk−1t−IMFk˜t.Step 6. We decomposerkt+εkEkωit (i=1,⋯,I) by EMD and compute the k+1-th IMF as(6)IMFk+1˜t=1I.∑iE1rkt+εkEkωit.Step 7. We repeated Step 4 to 6 until the residue has less than two extrema and cannot be decomposed any longer, so the final residue can be expressed as(7)Rt=xt−∑kIMFk˜t,where the original dataset xt is decomposed(8)xt=∑kIMFk˜t+Rt.


## 3.2. Grey Wolf Optimizer (GWO)

GWO is metaheuristic algorithm, which mimicks the leadership hierarchy and hunting mechanism of grey wolves in nature. Four types of grey wolves, includingα, β, δ, and ω, are employed for simulating the social dominant hierarchy. α wolves are the leaders and the decision makers in the pack. β wolves are the second group in the hierarchy that help α wolves in decision-making and establish a discipline for the pack. δ wolves obey α and β wolves, but dominate ω wolves. They play the roles of scouts, sentinels, elders, hunters, and caretakers. ω wolves, the lowest in hierarchy, always have to obey all the other types of wolves. The GWO algorithm mathematically models this social hierarchy and hunting process to collect randomly preferred solutions. Specifically, α, β, δ, and ω are considered as the fittest, the second preferred, the third preferred, and the rest solutions, respectively. The GWO algorithm contains the following three steps: tracking, encircling, and attacking the prey.In the encircling, grey wolves encircle prey during the hunt, and the mathematical formulation is given by(9)X⟶t+1=X⟶pt−A⟶.C⟶.X⟶pt−X⟶t,where t is the number of the current iteration, X⟶ and X⟶p are the positions of a grey wolf and the prey, respectively, A⟶=2a⟶.r⟶1−a⟶ and C⟶=2r⟶2 are coefficient vectors, components of a⟶ are linearly decreased from 2 to 0, and r⟶1 and r⟶2 are random vectors in [0, 1].In hunting, it is supposed that theα, β, and δ have better knowledge about the potential location of prey. Therefore, the current first three best solutions are saved, and the other search agents (including the ω) update their positions according to the position of the best search agent. The formulas are given by equations (10)–(13):(10)X⟶1=X⟶α−A⟶1.C⟶1.X⟶α−X⟶,(11)X⟶2=X⟶β−A⟶2.C⟶2.X⟶β−X⟶,(12)X⟶3=X⟶δ−A⟶3.C⟶3.X⟶δ−X⟶,(13)X⟶t+1=X⟶1+X⟶2+X⟶33.The grey wolf agents diverge in search of the prey and converge for attacking the prey, which can be mathematical modelled by the value ofA⟶. The grey wolf agents tend to diverge from the prey when A⟶>1 and converge towards the prey when A⟶<1. Finally, the GWO algorithm is terminated by minimizing of the fitness function, which would be discussed in Section 3.4.


## 3.3. Support Vector Regression (SVR)

SVR is a supervised machine learning algorithm. Based on structural risk minimization, SVR can solve nonlinear problems with high dimensionality and small samples. For a set of training data withN observations x1,y1,x2,y2,⋯,xN,yN, where xi∈ℝn are input vectors, yi∈ℝ are output vectors, and the SVR function can be described as the following linear regression function:(14)fx=wTΦx+b,where Φx is the feature mapping function and the parameters w and b are the weight and bias, respectively. Based on the structural risk minimization, w and b can be calculated as follows:(15)Min:12w2+C∑i=1Nξi+ξi∗,s.t.yi−wTΦx+b≤ε+ξi∗,wTΦx+b−yi≤ε+ξi∗,ξi,ξi∗≥0,where C>0 is a cost function, ξi and ξi∗ are slack variables, and ε is an insensitive loss function. By introducing the Lagrange multipliers αi and αi∗, the SVR function is expressed as(16)fx=∑i=1Nαi−αi∗kx,xi+b,where kx,xi=Φx.Φxi is the kernel function, which is the inner product of the feature mapping function.SVR would have a large forecasting error when dealing with nonlinear high-frequency time series. Therefore, unlike the linear kernel function used in forecasting the low-frequency time series, the Gaussian kernel,kx,xi=exp−γx−xi2, is used to forecast the high-frequency time series, owing to its advanced mapping capability and stable forecasting performance. Meanwhile, the selection of the SVR penalty parameter C and the kernel function parameter γ has an important influence on its forecasting accuracy [26]. Thus, this paper tunes these SVR parameters using the GWO algorithm.


## 3.4. The Hybrid CEEMDAD-GWO-SVR Model

The details of the hybrid CEEMDAD-GWO-SVR model are as follows, and the flowchart is presented in Figure2. It should be noted that the “decomposition-reconstruction-forecast-mergence” scheme is also shown in Figure 2.Step 1. We conduct CEEMDAN to decompose the original gold price data into several relatively simple IMFs and the residual by using equations (1)–(8). This step corresponds to the “decomposition” of the scheme.Step 2. We reconstruct the IMFs and the residual into three subsequences, namely, IMFH, IMFL, and Res. Specifically, since the high-frequency IMFs are affected by some short-term random factors, they can be assumed to fluctuate around the zero mean. Therefore, each IMF is reconstructed into high-frequency subsequence (IMFH) and low-frequency subsequence (IMFL) byT-test with zero mean. The residual is classified as the third subsequence, i.e., Res. This step corresponds to the “reconstruction” of the scheme.Step 3. Forecasting IMFH by GWO-SVR and forecasting IMFL and Res by SVR. The optimal GWO-SVR model is determined by minimizing the fitness function. In this paper, the mean square error, given by equation (17), is used to calculate the forecasting accuracy. This step is decisive for the prediction performance of the proposed model, corresponding to the “forecast” of the scheme.(17)MSE=1n∑i=1Nyi−y^i2.Step 4. This is the last part of the scheme. By summing up the forecasting results of these three subsequences, the final result is obtained.Figure 2
The flowchart of CEEMDAD-GWO-SVR model.


## 4. Empirical Study

This study selects the daily closing prices of international gold from January 3, 2000 to October 21, 2022, getting from the World Gold Council website (https://www.gold.org), with a total of 5,950 data points. Figure 3 displays the fluctuation diagram of the gold price. We take the first 5,655 data points (from January 3, 2000 to September 3, 2021) as the training set, and the last 295 data points (from September 5, 2021 to October 21, 2022) as the test set. Table 2 presents details of each set.Figure 3
Daily international gold price curve (January 3, 2000–October 21, 2022).Table 2
The summary of the daily international gold price data set.
Data setData sizeData periodStatistical descriptionsMeanVarianceRangeSkewnessKurtosisCoefficient of variationFull set595001/03/2000–09/03/20211040.18269829.341811.20−0.09331.73890.4994Training set565501/03/2000–09/03/20211000.17251272.091811.20−0.06251.74610.5012Test set29509/06/2021–21/10/20221807.146612.26407.350.18112.76720.0450


### 4.1. Data Decomposition and Reconstruction

#### 4.1.1. Data Decomposition

In the CEEMDAN, the parameters of ensemble sizeI and the signal-to-noise ratio ε0 are set as 500 and 0.2, respectively, and the maximal number of iterations is 2,000, as proposed in [44]. Decomposition results is shown in Figure 4. The gold price is decomposed into 13 sequences with 12 IMFs and one residual. The frequency of sequences is ranked from high to low. The high-frequency IMFs fluctuate relatively randomly but basically around zero mean, which can be regarded as the interference of short-term random factors (e.g., psychological factors, speculative factors, and short-term policies). The low-frequency IMFs have roughly similar trend to that of the original gold price data, reflecting the influence of important factors such as market policies and macroeconomics. The residual shows the overall upward trend, determining the long-term trend of gold price.Figure 4
Decomposition of the gold price by CEEMDAN (January 3, 2000–October 21, 2022).Meanwhile, Pearson correlation coefficient, Kendall rank correlation coefficient, and the variance contribution rate (VCR) are selected to measure the relationship between each sequence and the original gold price data, as depicted in Table3. From the two correlation coefficients, the high-frequency IMFs have a low association with the original gold price data, while the low-frequency IMFs and the residual have a high correlation with the original. The similar pattern is observed for VCR and the residual covers 66.6% of all variance, indicating that it has the greatest influence on the original data and reflects the internal trajectory of the gold price data. To sum up, the trend of gold price is mainly described by the residual. The low-frequency IMFs have a great influence on gold price in a certain period, while the high-frequency IMFs only have a small influence in a short period.Table 3
Correlation coefficient, VCR, andT test results for each IMF and the residual.
IMFs/residualPearson correlationKendall rank correlationVariance contribution rate (%)T testt statisticp valueIMF10.01290.01310.0144−0.31070.7561IMF20.01880.01460.0020−0.86140.3891IMF30.01210.00880.02220.13870.8897IMF40.01290.00500.0279−0.19740.8435IMF50.02320.01720.03700.20000.8415IMF60.04300.03120.0634−1.01890.3083IMF70.03700.02050.0874−0.28240.7776IMF80.04750.04010.1458−1.33090.1833IMF90.07660.04090.29760.09820.9218IMF100.10670.04660.3195−1.91830.0551IMF110.37270.13912.48372.59790.0094IMF120.64440.34546.6239−10.45650.0001Residual0.92100.739968.1219190.14370.0001


#### 4.1.2. Data Reconstruction

The high-frequency IMFs have a small influence on gold price in a short period, so they can be combined into one sequence to improve the efficiency of modeling, but the problem is how to determine which sequences are high-frequency IMFs. Since the high-frequency IMFs are affected by some short-term random factors, they can be assumed to fluctuate around zero mean. Therefore, IMFs can be reconstructed byT test with zero mean for each IMF. The results are shown in Table 3. The p values of IMF1∼IMF9 are greater than 0.01.Since the nine IMFs fail to reject the null hypothesis fort test, so they are all divided into the high-frequency subsequence (IMFH) and IMFH=∑i=19IMFi. The remaining sequences, IMF10, IMF11, and IMF12, reject the null hypothesis for t test at 1% significance level, so we define the low-frequency subsequence (IMFL) as IMFL=∑i=1012IMFi. The residual is classified as the third subsequence of the gold price data, i.e., Res. Table 4 presents the reconstructed results.Table 4
Reconstructed results of the original gold price data.
SubsequencesSequences/residualPearson correlationKendall rank correlationVariance contribution rate (%)IMFHIMF1, IMF2, IMF3, IMF4, IMF5, IMF6, IMF7, IMF8, IMF90.09330.05370.9831IMFLIMF10, IMF11, IMF120.59690.268014.9269ResResidual0.92100.739968.1219Figure5 shows the reconstructed subsequence trend, and it can be found that: (1) IMFH fluctuates relatively randomly but basically around zero mean, which is mainly interfered by some short-term random factors; (2) IMFL has roughly similar trend to that of the original gold price, which is mainly affected by market policies and macroeconomics; and (3) Res shows the overall upward trend, conducting the trend of gold price.Figure 5
Reconstructed subsequence trend of the gold price (january 3, 2000–october 21, 2022).


### 4.2. Forecasting Results and Discussion

#### 4.2.1. Evaluation Measures

The forecasting accuracies of the proposed model and the benchmark models are evaluated with the following three measures: the mean absolute percentage error (MAPE), the root mean square error (RMSE), and the directional accuracy (DA). Their formulas are given by:(18)MAPE=1n.∑i=1nyi−fiyi,RMSE=1/n.∑i=1nyi−fi2,DA=1n.∑i=1ndi⋅100%,di=1,yi+1−yifi+1−yi≥0,0,yi+1−yifi+1−yi<0,where n is the number of forecasting points, yi is the actual gold price at point i, and fi is the gold price forecast at point i. Obviously, the smaller the MAPE and RMSE values, the higher the forecasting accuracy of the model, and the larger the DA values, the higher the direction accuracy of the model.


#### 4.2.2. In-Sample Forecast (January 1, 2016–September 3, 2021)

The original gold price data is decomposed and reconstructed into three subsequences: IMFH, IMFL, and Res, which can be expressed as(19)y=IMFH+IMFL+Res,where y is the actual gold price. The three subsequences are forecasted separately, and then combined to obtain the final forecast of the gold price. Considering the low forecasting accuracy of SVR when dealing with nonlinear high-frequency time series, IMFH is forecasted by GWO-SVR to improve the parameter selection, and IMFL and Res forecasts are conducted by SVR.In GWO algorithm, the optimization range is set at (0.001 and 200), and the population size and the maximum number of iterations are set to 20 and 50, respectively. For IMFH, the optimal penalty parameterC and the kernel function parameter γ of SVR are −132.7510 and −1.8545, selected by GWO. All calculations are done on the Ryzen R7-4800u CPU, with 16.0 GB DDR4 3200Mhz RAM, and windows 11 system.Figure6 shows the forecasting results of the three subsequences and the original gold price data. It can be observed that the blue lines (forecasting results) are very close to the red lines (the reconstructed subsequences), which indicate that the proposed CEEMDAN-GWO-SVR model achieves a relatively good forecasting result. Table 5 presents the measures of the forecasting accuracy. IMFL and Res have relatively small errors and high directional accuracy, achieving good forecasting performance. IMFH has lower accuracy than the other two, and it also reduces the forecasting accuracy of the original gold price data. This is because IMFH is affected by random factors, which increases the difficulty of forecasting.Figure 6
In-sample forecasting results of the original gold price data.Table 5
Forecasting accuracy indexes of the three subsequences and the original gold price data.
SubsequenceModelRMSEMAPEDAIMFHCEEMDAN-GWO-SVR11.75561.031776.0000IMFLCEEMDAN-SVR0.43810.005676.9067IMFRCEEMDAN-SVR1.01740.001160.8711Table6 presents the accuracy measures of the forecasting results of the proposed CEEMDAN-GWO-SVR model, SVR-based models, and non-SVR-based models. Among them, the SVR-based models include CEEMDAN-SVR, GWO-SVR, and SVR, and the non-SVR-based models include ARIMA, GARCH, ANN, CART, and BPNN (back propagation neural network). Obviously, the proposed CEEMDAN-GWO-SVR model is superior to the alternatives. For the SVR-based models, the CEEMDAN-SVR, CEEMDAN-PSO-SVR, and CEEMDAN-GA-SVR models can certainly improve the performance of SVR. While the hybrid CEEMDAN-GWO-SVR model achieves the best forecasting result, which shows that the proposed model can critically capture effective information from the original data, and it also demonstrates the superiority of the decomposition-reconstruction-forecast-mergence scheme. Figure 7 shows the convergence plots of GWO, which reflected that the objective can be minimized by the algorithm before the maximum number of iterations.Table 6
In-sample forecasting accuracy measures of the CEEMDAN-GWO-SVR and benchmark models.
ModelRMSEMAPEDACEEMDAN-GWO-SVR11.81840.007452.0178CEEMDAN-SVR11.89580.007451.6800CEEMDAN-GA-SVR11.80490.007451.7511CEEMDAN-PSO-SVR11.88990.007452.0356SVR12.41880.007849.4756ANN11.88670.007449.8667BPNN11.98880.007549.6356CART7.08800.005766.9511ARIMA (1, 1, 1)3.68810.001987.0933ARIMA (1, 1, 1)-GARCH (1, 1)1.09550.001787.0400Figure 7
The convergence plots of GWO.


#### 4.2.3. Out-of-Sample Forecast (September 6, 2021–October 21, 2022)

The above results prove that the CEEMDAN-GWO-SVR model outperforms the benchmark models for in-sample forecast. In order to further test the proposed model for out-of-sample forecasts, this paper estimates the gold price for the next 295 consecutive days (September 6, 2021–October 2021). Parameters of the CEEMDAN-GWO-SVR model here are set in the same way as above.Out-of-sample forecasting accuracy indexes comparison are respectively shown in Table7. The superiority of the proposed CEEMDAN-GWO-SVR model can also be demonstrated. Among the three SVR-based models, the CEEMDAN-SVR model obtains the best performance, and at the same time, it also has superior forecast performance than the non-SVR-based models. This shows that the CEEMDAN method could improve the forecasting performance. Moreover, when CEEMDAN-SVR and GWO are hybridized, the forecasting accuracy of gold price is improved further, which shows the superiority of the CEEMDAN-GWO-SVR model proposed in this paper. It should be noted that in short-term forecasts, although the CEEMDAN-GWO-SVR model has relatively lower directional accuracy, it still has smaller forecasting errors (MAPE and RMSE) than alternatives. Therefore, the proposed CEEMDAN-GWO-SVR model can be considered as the best model.Table 7
Out-of-sample forecasting accuracy measures.
ModelRMSEMAPEDACEEMDAN-GWO-SVR16.16480.006746.0377CEEMDAN-SVR16.25550.006846.4151CEEMDAN-GA-SVR16.17430.006748.3019CEEMDAN-PSO-SVR16.18850.006746.4151SVR17.08560.007049.0566ANN16.32750.006748.6792BPNN16.60900.006849.8113CART24.48740.010846.7925ARIMA (1, 1, 1)87.39970.038848.6792ARIMA (1, 1, 1)-GARCH (1, 1)87.13620.034449.4340


## 4.1. Data Decomposition and Reconstruction

### 4.1.1. Data Decomposition

In the CEEMDAN, the parameters of ensemble sizeI and the signal-to-noise ratio ε0 are set as 500 and 0.2, respectively, and the maximal number of iterations is 2,000, as proposed in [44]. Decomposition results is shown in Figure 4. The gold price is decomposed into 13 sequences with 12 IMFs and one residual. The frequency of sequences is ranked from high to low. The high-frequency IMFs fluctuate relatively randomly but basically around zero mean, which can be regarded as the interference of short-term random factors (e.g., psychological factors, speculative factors, and short-term policies). The low-frequency IMFs have roughly similar trend to that of the original gold price data, reflecting the influence of important factors such as market policies and macroeconomics. The residual shows the overall upward trend, determining the long-term trend of gold price.Figure 4
Decomposition of the gold price by CEEMDAN (January 3, 2000–October 21, 2022).Meanwhile, Pearson correlation coefficient, Kendall rank correlation coefficient, and the variance contribution rate (VCR) are selected to measure the relationship between each sequence and the original gold price data, as depicted in Table3. From the two correlation coefficients, the high-frequency IMFs have a low association with the original gold price data, while the low-frequency IMFs and the residual have a high correlation with the original. The similar pattern is observed for VCR and the residual covers 66.6% of all variance, indicating that it has the greatest influence on the original data and reflects the internal trajectory of the gold price data. To sum up, the trend of gold price is mainly described by the residual. The low-frequency IMFs have a great influence on gold price in a certain period, while the high-frequency IMFs only have a small influence in a short period.Table 3
Correlation coefficient, VCR, andT test results for each IMF and the residual.
IMFs/residualPearson correlationKendall rank correlationVariance contribution rate (%)T testt statisticp valueIMF10.01290.01310.0144−0.31070.7561IMF20.01880.01460.0020−0.86140.3891IMF30.01210.00880.02220.13870.8897IMF40.01290.00500.0279−0.19740.8435IMF50.02320.01720.03700.20000.8415IMF60.04300.03120.0634−1.01890.3083IMF70.03700.02050.0874−0.28240.7776IMF80.04750.04010.1458−1.33090.1833IMF90.07660.04090.29760.09820.9218IMF100.10670.04660.3195−1.91830.0551IMF110.37270.13912.48372.59790.0094IMF120.64440.34546.6239−10.45650.0001Residual0.92100.739968.1219190.14370.0001


### 4.1.2. Data Reconstruction

The high-frequency IMFs have a small influence on gold price in a short period, so they can be combined into one sequence to improve the efficiency of modeling, but the problem is how to determine which sequences are high-frequency IMFs. Since the high-frequency IMFs are affected by some short-term random factors, they can be assumed to fluctuate around zero mean. Therefore, IMFs can be reconstructed byT test with zero mean for each IMF. The results are shown in Table 3. The p values of IMF1∼IMF9 are greater than 0.01.Since the nine IMFs fail to reject the null hypothesis fort test, so they are all divided into the high-frequency subsequence (IMFH) and IMFH=∑i=19IMFi. The remaining sequences, IMF10, IMF11, and IMF12, reject the null hypothesis for t test at 1% significance level, so we define the low-frequency subsequence (IMFL) as IMFL=∑i=1012IMFi. The residual is classified as the third subsequence of the gold price data, i.e., Res. Table 4 presents the reconstructed results.Table 4
Reconstructed results of the original gold price data.
SubsequencesSequences/residualPearson correlationKendall rank correlationVariance contribution rate (%)IMFHIMF1, IMF2, IMF3, IMF4, IMF5, IMF6, IMF7, IMF8, IMF90.09330.05370.9831IMFLIMF10, IMF11, IMF120.59690.268014.9269ResResidual0.92100.739968.1219Figure5 shows the reconstructed subsequence trend, and it can be found that: (1) IMFH fluctuates relatively randomly but basically around zero mean, which is mainly interfered by some short-term random factors; (2) IMFL has roughly similar trend to that of the original gold price, which is mainly affected by market policies and macroeconomics; and (3) Res shows the overall upward trend, conducting the trend of gold price.Figure 5
Reconstructed subsequence trend of the gold price (january 3, 2000–october 21, 2022).


## 4.1.1. Data Decomposition

In the CEEMDAN, the parameters of ensemble sizeI and the signal-to-noise ratio ε0 are set as 500 and 0.2, respectively, and the maximal number of iterations is 2,000, as proposed in [44]. Decomposition results is shown in Figure 4. The gold price is decomposed into 13 sequences with 12 IMFs and one residual. The frequency of sequences is ranked from high to low. The high-frequency IMFs fluctuate relatively randomly but basically around zero mean, which can be regarded as the interference of short-term random factors (e.g., psychological factors, speculative factors, and short-term policies). The low-frequency IMFs have roughly similar trend to that of the original gold price data, reflecting the influence of important factors such as market policies and macroeconomics. The residual shows the overall upward trend, determining the long-term trend of gold price.Figure 4
Decomposition of the gold price by CEEMDAN (January 3, 2000–October 21, 2022).Meanwhile, Pearson correlation coefficient, Kendall rank correlation coefficient, and the variance contribution rate (VCR) are selected to measure the relationship between each sequence and the original gold price data, as depicted in Table3. From the two correlation coefficients, the high-frequency IMFs have a low association with the original gold price data, while the low-frequency IMFs and the residual have a high correlation with the original. The similar pattern is observed for VCR and the residual covers 66.6% of all variance, indicating that it has the greatest influence on the original data and reflects the internal trajectory of the gold price data. To sum up, the trend of gold price is mainly described by the residual. The low-frequency IMFs have a great influence on gold price in a certain period, while the high-frequency IMFs only have a small influence in a short period.Table 3
Correlation coefficient, VCR, andT test results for each IMF and the residual.
IMFs/residualPearson correlationKendall rank correlationVariance contribution rate (%)T testt statisticp valueIMF10.01290.01310.0144−0.31070.7561IMF20.01880.01460.0020−0.86140.3891IMF30.01210.00880.02220.13870.8897IMF40.01290.00500.0279−0.19740.8435IMF50.02320.01720.03700.20000.8415IMF60.04300.03120.0634−1.01890.3083IMF70.03700.02050.0874−0.28240.7776IMF80.04750.04010.1458−1.33090.1833IMF90.07660.04090.29760.09820.9218IMF100.10670.04660.3195−1.91830.0551IMF110.37270.13912.48372.59790.0094IMF120.64440.34546.6239−10.45650.0001Residual0.92100.739968.1219190.14370.0001


## 4.1.2. Data Reconstruction

The high-frequency IMFs have a small influence on gold price in a short period, so they can be combined into one sequence to improve the efficiency of modeling, but the problem is how to determine which sequences are high-frequency IMFs. Since the high-frequency IMFs are affected by some short-term random factors, they can be assumed to fluctuate around zero mean. Therefore, IMFs can be reconstructed byT test with zero mean for each IMF. The results are shown in Table 3. The p values of IMF1∼IMF9 are greater than 0.01.Since the nine IMFs fail to reject the null hypothesis fort test, so they are all divided into the high-frequency subsequence (IMFH) and IMFH=∑i=19IMFi. The remaining sequences, IMF10, IMF11, and IMF12, reject the null hypothesis for t test at 1% significance level, so we define the low-frequency subsequence (IMFL) as IMFL=∑i=1012IMFi. The residual is classified as the third subsequence of the gold price data, i.e., Res. Table 4 presents the reconstructed results.Table 4
Reconstructed results of the original gold price data.
SubsequencesSequences/residualPearson correlationKendall rank correlationVariance contribution rate (%)IMFHIMF1, IMF2, IMF3, IMF4, IMF5, IMF6, IMF7, IMF8, IMF90.09330.05370.9831IMFLIMF10, IMF11, IMF120.59690.268014.9269ResResidual0.92100.739968.1219Figure5 shows the reconstructed subsequence trend, and it can be found that: (1) IMFH fluctuates relatively randomly but basically around zero mean, which is mainly interfered by some short-term random factors; (2) IMFL has roughly similar trend to that of the original gold price, which is mainly affected by market policies and macroeconomics; and (3) Res shows the overall upward trend, conducting the trend of gold price.Figure 5
Reconstructed subsequence trend of the gold price (january 3, 2000–october 21, 2022).


## 4.2. Forecasting Results and Discussion

### 4.2.1. Evaluation Measures

The forecasting accuracies of the proposed model and the benchmark models are evaluated with the following three measures: the mean absolute percentage error (MAPE), the root mean square error (RMSE), and the directional accuracy (DA). Their formulas are given by:(18)MAPE=1n.∑i=1nyi−fiyi,RMSE=1/n.∑i=1nyi−fi2,DA=1n.∑i=1ndi⋅100%,di=1,yi+1−yifi+1−yi≥0,0,yi+1−yifi+1−yi<0,where n is the number of forecasting points, yi is the actual gold price at point i, and fi is the gold price forecast at point i. Obviously, the smaller the MAPE and RMSE values, the higher the forecasting accuracy of the model, and the larger the DA values, the higher the direction accuracy of the model.


### 4.2.2. In-Sample Forecast (January 1, 2016–September 3, 2021)

The original gold price data is decomposed and reconstructed into three subsequences: IMFH, IMFL, and Res, which can be expressed as(19)y=IMFH+IMFL+Res,where y is the actual gold price. The three subsequences are forecasted separately, and then combined to obtain the final forecast of the gold price. Considering the low forecasting accuracy of SVR when dealing with nonlinear high-frequency time series, IMFH is forecasted by GWO-SVR to improve the parameter selection, and IMFL and Res forecasts are conducted by SVR.In GWO algorithm, the optimization range is set at (0.001 and 200), and the population size and the maximum number of iterations are set to 20 and 50, respectively. For IMFH, the optimal penalty parameterC and the kernel function parameter γ of SVR are −132.7510 and −1.8545, selected by GWO. All calculations are done on the Ryzen R7-4800u CPU, with 16.0 GB DDR4 3200Mhz RAM, and windows 11 system.Figure6 shows the forecasting results of the three subsequences and the original gold price data. It can be observed that the blue lines (forecasting results) are very close to the red lines (the reconstructed subsequences), which indicate that the proposed CEEMDAN-GWO-SVR model achieves a relatively good forecasting result. Table 5 presents the measures of the forecasting accuracy. IMFL and Res have relatively small errors and high directional accuracy, achieving good forecasting performance. IMFH has lower accuracy than the other two, and it also reduces the forecasting accuracy of the original gold price data. This is because IMFH is affected by random factors, which increases the difficulty of forecasting.Figure 6
In-sample forecasting results of the original gold price data.Table 5
Forecasting accuracy indexes of the three subsequences and the original gold price data.
SubsequenceModelRMSEMAPEDAIMFHCEEMDAN-GWO-SVR11.75561.031776.0000IMFLCEEMDAN-SVR0.43810.005676.9067IMFRCEEMDAN-SVR1.01740.001160.8711Table6 presents the accuracy measures of the forecasting results of the proposed CEEMDAN-GWO-SVR model, SVR-based models, and non-SVR-based models. Among them, the SVR-based models include CEEMDAN-SVR, GWO-SVR, and SVR, and the non-SVR-based models include ARIMA, GARCH, ANN, CART, and BPNN (back propagation neural network). Obviously, the proposed CEEMDAN-GWO-SVR model is superior to the alternatives. For the SVR-based models, the CEEMDAN-SVR, CEEMDAN-PSO-SVR, and CEEMDAN-GA-SVR models can certainly improve the performance of SVR. While the hybrid CEEMDAN-GWO-SVR model achieves the best forecasting result, which shows that the proposed model can critically capture effective information from the original data, and it also demonstrates the superiority of the decomposition-reconstruction-forecast-mergence scheme. Figure 7 shows the convergence plots of GWO, which reflected that the objective can be minimized by the algorithm before the maximum number of iterations.Table 6
In-sample forecasting accuracy measures of the CEEMDAN-GWO-SVR and benchmark models.
ModelRMSEMAPEDACEEMDAN-GWO-SVR11.81840.007452.0178CEEMDAN-SVR11.89580.007451.6800CEEMDAN-GA-SVR11.80490.007451.7511CEEMDAN-PSO-SVR11.88990.007452.0356SVR12.41880.007849.4756ANN11.88670.007449.8667BPNN11.98880.007549.6356CART7.08800.005766.9511ARIMA (1, 1, 1)3.68810.001987.0933ARIMA (1, 1, 1)-GARCH (1, 1)1.09550.001787.0400Figure 7
The convergence plots of GWO.


### 4.2.3. Out-of-Sample Forecast (September 6, 2021–October 21, 2022)

The above results prove that the CEEMDAN-GWO-SVR model outperforms the benchmark models for in-sample forecast. In order to further test the proposed model for out-of-sample forecasts, this paper estimates the gold price for the next 295 consecutive days (September 6, 2021–October 2021). Parameters of the CEEMDAN-GWO-SVR model here are set in the same way as above.Out-of-sample forecasting accuracy indexes comparison are respectively shown in Table7. The superiority of the proposed CEEMDAN-GWO-SVR model can also be demonstrated. Among the three SVR-based models, the CEEMDAN-SVR model obtains the best performance, and at the same time, it also has superior forecast performance than the non-SVR-based models. This shows that the CEEMDAN method could improve the forecasting performance. Moreover, when CEEMDAN-SVR and GWO are hybridized, the forecasting accuracy of gold price is improved further, which shows the superiority of the CEEMDAN-GWO-SVR model proposed in this paper. It should be noted that in short-term forecasts, although the CEEMDAN-GWO-SVR model has relatively lower directional accuracy, it still has smaller forecasting errors (MAPE and RMSE) than alternatives. Therefore, the proposed CEEMDAN-GWO-SVR model can be considered as the best model.Table 7
Out-of-sample forecasting accuracy measures.
ModelRMSEMAPEDACEEMDAN-GWO-SVR16.16480.006746.0377CEEMDAN-SVR16.25550.006846.4151CEEMDAN-GA-SVR16.17430.006748.3019CEEMDAN-PSO-SVR16.18850.006746.4151SVR17.08560.007049.0566ANN16.32750.006748.6792BPNN16.60900.006849.8113CART24.48740.010846.7925ARIMA (1, 1, 1)87.39970.038848.6792ARIMA (1, 1, 1)-GARCH (1, 1)87.13620.034449.4340


## 4.2.1. Evaluation Measures

The forecasting accuracies of the proposed model and the benchmark models are evaluated with the following three measures: the mean absolute percentage error (MAPE), the root mean square error (RMSE), and the directional accuracy (DA). Their formulas are given by:(18)MAPE=1n.∑i=1nyi−fiyi,RMSE=1/n.∑i=1nyi−fi2,DA=1n.∑i=1ndi⋅100%,di=1,yi+1−yifi+1−yi≥0,0,yi+1−yifi+1−yi<0,where n is the number of forecasting points, yi is the actual gold price at point i, and fi is the gold price forecast at point i. Obviously, the smaller the MAPE and RMSE values, the higher the forecasting accuracy of the model, and the larger the DA values, the higher the direction accuracy of the model.


## 4.2.2. In-Sample Forecast (January 1, 2016–September 3, 2021)

The original gold price data is decomposed and reconstructed into three subsequences: IMFH, IMFL, and Res, which can be expressed as(19)y=IMFH+IMFL+Res,where y is the actual gold price. The three subsequences are forecasted separately, and then combined to obtain the final forecast of the gold price. Considering the low forecasting accuracy of SVR when dealing with nonlinear high-frequency time series, IMFH is forecasted by GWO-SVR to improve the parameter selection, and IMFL and Res forecasts are conducted by SVR.In GWO algorithm, the optimization range is set at (0.001 and 200), and the population size and the maximum number of iterations are set to 20 and 50, respectively. For IMFH, the optimal penalty parameterC and the kernel function parameter γ of SVR are −132.7510 and −1.8545, selected by GWO. All calculations are done on the Ryzen R7-4800u CPU, with 16.0 GB DDR4 3200Mhz RAM, and windows 11 system.Figure6 shows the forecasting results of the three subsequences and the original gold price data. It can be observed that the blue lines (forecasting results) are very close to the red lines (the reconstructed subsequences), which indicate that the proposed CEEMDAN-GWO-SVR model achieves a relatively good forecasting result. Table 5 presents the measures of the forecasting accuracy. IMFL and Res have relatively small errors and high directional accuracy, achieving good forecasting performance. IMFH has lower accuracy than the other two, and it also reduces the forecasting accuracy of the original gold price data. This is because IMFH is affected by random factors, which increases the difficulty of forecasting.Figure 6
In-sample forecasting results of the original gold price data.Table 5
Forecasting accuracy indexes of the three subsequences and the original gold price data.
SubsequenceModelRMSEMAPEDAIMFHCEEMDAN-GWO-SVR11.75561.031776.0000IMFLCEEMDAN-SVR0.43810.005676.9067IMFRCEEMDAN-SVR1.01740.001160.8711Table6 presents the accuracy measures of the forecasting results of the proposed CEEMDAN-GWO-SVR model, SVR-based models, and non-SVR-based models. Among them, the SVR-based models include CEEMDAN-SVR, GWO-SVR, and SVR, and the non-SVR-based models include ARIMA, GARCH, ANN, CART, and BPNN (back propagation neural network). Obviously, the proposed CEEMDAN-GWO-SVR model is superior to the alternatives. For the SVR-based models, the CEEMDAN-SVR, CEEMDAN-PSO-SVR, and CEEMDAN-GA-SVR models can certainly improve the performance of SVR. While the hybrid CEEMDAN-GWO-SVR model achieves the best forecasting result, which shows that the proposed model can critically capture effective information from the original data, and it also demonstrates the superiority of the decomposition-reconstruction-forecast-mergence scheme. Figure 7 shows the convergence plots of GWO, which reflected that the objective can be minimized by the algorithm before the maximum number of iterations.Table 6
In-sample forecasting accuracy measures of the CEEMDAN-GWO-SVR and benchmark models.
ModelRMSEMAPEDACEEMDAN-GWO-SVR11.81840.007452.0178CEEMDAN-SVR11.89580.007451.6800CEEMDAN-GA-SVR11.80490.007451.7511CEEMDAN-PSO-SVR11.88990.007452.0356SVR12.41880.007849.4756ANN11.88670.007449.8667BPNN11.98880.007549.6356CART7.08800.005766.9511ARIMA (1, 1, 1)3.68810.001987.0933ARIMA (1, 1, 1)-GARCH (1, 1)1.09550.001787.0400Figure 7
The convergence plots of GWO.


## 4.2.3. Out-of-Sample Forecast (September 6, 2021–October 21, 2022)

The above results prove that the CEEMDAN-GWO-SVR model outperforms the benchmark models for in-sample forecast. In order to further test the proposed model for out-of-sample forecasts, this paper estimates the gold price for the next 295 consecutive days (September 6, 2021–October 2021). Parameters of the CEEMDAN-GWO-SVR model here are set in the same way as above.Out-of-sample forecasting accuracy indexes comparison are respectively shown in Table7. The superiority of the proposed CEEMDAN-GWO-SVR model can also be demonstrated. Among the three SVR-based models, the CEEMDAN-SVR model obtains the best performance, and at the same time, it also has superior forecast performance than the non-SVR-based models. This shows that the CEEMDAN method could improve the forecasting performance. Moreover, when CEEMDAN-SVR and GWO are hybridized, the forecasting accuracy of gold price is improved further, which shows the superiority of the CEEMDAN-GWO-SVR model proposed in this paper. It should be noted that in short-term forecasts, although the CEEMDAN-GWO-SVR model has relatively lower directional accuracy, it still has smaller forecasting errors (MAPE and RMSE) than alternatives. Therefore, the proposed CEEMDAN-GWO-SVR model can be considered as the best model.Table 7
Out-of-sample forecasting accuracy measures.
ModelRMSEMAPEDACEEMDAN-GWO-SVR16.16480.006746.0377CEEMDAN-SVR16.25550.006846.4151CEEMDAN-GA-SVR16.17430.006748.3019CEEMDAN-PSO-SVR16.18850.006746.4151SVR17.08560.007049.0566ANN16.32750.006748.6792BPNN16.60900.006849.8113CART24.48740.010846.7925ARIMA (1, 1, 1)87.39970.038848.6792ARIMA (1, 1, 1)-GARCH (1, 1)87.13620.034449.4340


## 5. Conclusions

This paper proposes a decomposition-reconstruction-forecast-mergence scheme to forecast the international gold price. In this scheme, a useful gold price forecasting model is proposed by combining the CEEMDAN method with GWO-SVR. The experimental results verify the validity of the CEEMDAN-GWO-SVR model for gold price forecast. The following conclusions are drawn:(1)
Due to the complex pattern of the gold price, data preprocessing in the early stage of forecasting is necessary. The CEEMDAN method adds white Gaussian noise at each stage of the decomposition and sufficiently reduces the noise of the original data. Also, experimental results show that the CEEMDAN method can improves the forecasting performance of SVR.(2)
Although the CEEMDAN-SVR model can certainly improve the performance of SVR, the forecasting accuracy can be further improved by hybridizing CEEMDAN-SVR and GWO. This is because GWO optimizes the selection of SVR parameters for high-frequency time series decomposed by CEEMDAN, thus ultimately improving the forecasting performance.(3)
The CEEMDAN-GWO-SVR model significantly outperforms the other six benchmark models, which also demonstrates the superiority of the decomposition-reconstruction-forecast-mergence scheme.However, this paper still has certain limitations. This paper only uses time series to forecast gold prices. Although the proposed model has achieved good forecasting results, the research on the main influencing factors of gold price is still an important research issue.Further research can be carried out from the following aspects in the future. The GWO algorithm has certain problems in population diversity and the balance between the exploitation and exploration [49]. Therefore, the algorithm could be improved to further improve the prediction performance of the SVR model. Moreover, except GWO, finding a suitable approach to hybridize more recent metaheuristic algorithms with SVR is also an important research issue.



---

*Source: 1511479-2022-12-17.xml*