Forecasting the carbon price sequence in the Hubei emissions exchange using a hybrid model based on ensemble empirical mode decomposition

The prediction of carbon price is exceedingly essential for the regulators, investors, and participants of the carbon trading market. It is the basis for formulating market policies and improving risk management capabilities. China's carbon price series are nonlinear and nonstationary, so it is difficult to predict accurately with traditional models. This paper proposes a multiscale ensemble prediction model based on ensemble empirical mode decomposition (EEMD‐ADD) to improve the prediction accuracy of carbon price. Firstly, EEMD is used to decompose the carbon price sequence into several intrinsic mode functions (IMFs), and these IMFs are divided into high‐frequency component, low‐frequency component and the trend component. Then, LSSVM, PSO‐LSSVM, and BA‐LSSVM are used to predict the three components respectively after comparative analysis. Finally, the results are combined to obtain the final prediction value. In the empirical analysis of the Hubei Emissions Exchange, the proposed model outperforms other comparative models. The RMSE, MAE, and MAPE values of the EEMD‐ADD model are 0.6180, 0.4726, and 1.6342, and the DS, CP, and CD values are 94.36, 92.16, and 96.48. In addition, the model performed best in other time periods. The results suggest that the proposed model is effective and could predict carbon prices more accurately.


| 2709
WU and LIU for 20% of primary energy consumption, carbon emissions will reach the peak around 2030 and strive for an early peak. Nevertheless, China's rapid economic growth has led to a sharp increase in energy demand and a gradual increase in carbon emissions. And the traditional methods are difficult to meet the needs of emission reduction. So, the Chinese government has proposed the carbon emission trading scheme (ETS) as an important measure to mitigate climate change. 4 In other words, the market regulates carbon emissions by itself, allowing participants to realize their interests and needs. This method is considered to be the most effective way to reduce emissions 5 and has become an important tool for the Chinese government to formulate policies. Since 2011, China has successively set up several carbon emission exchanges. With the development and improvement of carbon emission trading market, the frequency and liveness of carbon trading have gradually increased. The forecasting of carbon prices plays an increasingly significant role in understanding China's carbon market dynamics and in making decisions about carbon emissions reductions. 6 Therefore, it is intensely crucial to select appropriate methods to achieve higher prediction accuracy.
With the maturity of carbon trading market at home and abroad, the research on carbon price has increased dramatically. These studies are mainly divided into two categories, one of which is based on econometric and statistical models for prediction. Among them, the generalized autoregressive conditional heteroskedasticity (GARCH) model is widely used. Byun and Cho applied GARCH model to predict the carbon price of EU ETS. 7 Li and Lu used empirical mode decomposition (EMD) and generalized autoregressive conditional heteroskedasticity (GARCH) model to study the price data of five pilot projects in China and predicted the carbon price. 8 Paolella and Taschini used the mixed orthodox GARCH model to predict the carbon price of the first stage in EU ETS. 9 Benz and Truck adopted MS-AR-GARCH model to conduct a recursive rolling prediction of carbon price gains in the first stage of EU ETS. 10 However, these methods cannot meet the requirements of high accuracy because of the nonlinear and nonstationary of carbon price series. To overcome the shortcomings of traditional methods, intelligent algorithms with self-learning ability have been introduced by researchers. Zhu and Wei used LSSVM model to forecast the carbon price of EU ETS from April 2005 to February 2011. 11 In 2013, Zhu and Wei proposed ARIMA-LSSVM model to forecast DEC 10 and DEC 12 carbon futures contracts of the European Climate Exchange (ECX). 12 Tian and Hao established an improved carbon price analysis and prediction system composed of analysis module and prediction module, and the empirical results show that the HOA-based model has the best prediction performance compared with other traditional models. 13 In order to predict carbon price more accurately, Hao and Tian proposed a kernel-based extremum learning machine optimization model with good generalization and stability and based on this, established a mixed prediction framework considering multiple factors. 14 Lu et al 15 used six different artificial intelligence models to predict the daily carbon price and trading volume of eight carbon markets in China, by comparing the prediction results, a more suitable model for specific carbon market was obtained. In order to make up for the defects of ARMA model and accurately predict the carbon price, Wang used Boosting algorithm for optimal subset search. 16 Artificial intelligence algorithm effectively overcomes the limitations of traditional models, but its forecasting ability for nonlinear and nonstationary carbon price series is still insufficient. In recent years, the multiscale ensemble hybrid forecasting model represented by EMD has been widely used. Yao and Lv studied the formation mechanism and carbon price prediction of carbon market in Hubei province by virtue of EMD-SVM model. 17 Zhu et al 18 establish a carbon price forecasting model based on EMD and neural networks. Zhu and Han used EMD-LSSVM method to predict carbon price. 19 Gao used the EMD-PSO-SVM model to predict the carbon price. 20 Huang and Sun proposed a novel hybrid carbon price forecasting model based on EMD-VMD-GA-BP. 21 Although EMD can improve the model's fitting degree, it is easier to cause mode mixing, thus reducing the prediction accuracy. Wu and Huang 22 proposed the EEMD which adds the Gaussian white noise so as to avoid the mode mixing problem. For this reason, mounting number of scholars use EEMD-types 23-27 model to predict nonlinear and nonstationary sequences. It is also widely used in carbon price forecasting, for instance, Zhang and Hu used EEMD and BP to forecast the error sequence, and the results show that the method has higher precision accuracy on carbon price. 28 Feng et al 29 used variance ratio and EEMD to analyze the carbon price fluctuation factors in the EU ETS.
However, the research on carbon price prediction is in its initial stage, and there are still some shortcomings: (a) the existing prediction methods are not robust for accurately forecasting carbon prices that have the characteristics of the nonlinear and nonstationary. Krogh et al 30 pointed out that when the single model constituting the combination model is accurate and diversified, the combination model is bound to achieve better results than the single model. Therefore, this paper analyses the characteristics of different frequency components decomposed by EEMD and chooses a suitable prediction model for each frequency component. Finally, all classification prediction results are integrated. (b) Most of the existing studies focus on the international carbon trading market and seldom take China's carbon price as a sample. Through Augmented-Dickey-Fuller (ADF) and Brock-Dechert-Scheinkman (BDS) tests, it is known that the regional carbon price in China ETS has the characteristics of nonlinear and nonstationary, which is similar to the fluctuation trend of the international carbon market. So, this paper can use the method of studying the international carbon price to analyze and predict the China ETS. 31 In this paper, a multiscale ensemble forecasting model is constructed based on EEMD, which is called EEMD-ADD. The original price sequence of Hubei Emissions Exchange is taken as sample for empirical analysis. To verify the validity of the model, the EEMD-ADD model is compared with LSSVM, PSO, BA, PSO-LSSVM, and BA-LSSVM. The innovations and main contributions of this paper are as follows: (a) EEMD algorithm is applied to predict carbon price. The IMFs obtained from EEMD is adopted to study the carbon price series, which reduces the data irregularity and improves the prediction accuracy. (b) Combining EEMD, LSSVM, PSO-LSSVM, and BA-LSSVM, the hybrid model of carbon price forecasting is first proposed, which might provide a new idea for carbon price forecasting, and effectively make up for the shortcomings of the single model. (c) Data from the carbon markets in Hubei are used for empirical analysis, while five other relevant models and other time periods are used for performance comparison to prove the robustness and effectiveness of the models. In addition, different evaluation indicators are used to analysis the model performance, including RMSE, MAE, MAPE, and DS, CP, CD.
The rest of this article is arranged as follows: Section 2 introduces the algorithms mentioned in this article. Section 3 introduces the specific structure of the hybrid prediction model. Section 4 carries on the demonstration analysis based on the carbon price sequence of the Hubei Emission Exchange. Section 5 verifies the prediction performance of the EEMD-ADD. Finally, Section 6 provides some conclusions and future work.

| METHODOLOGY
This paper first introduces four basic algorithms EEMD, LSSVM, PSO, BA and then proposes two combined algorithms based on them. The details are as follows.

| Ensemble empirical mode decomposition
Empirical mode decomposition (EMD) is a method suitable for analyzing and processing of nonlinear and nonstationary signals. [32][33][34] The method filters arbitrary signals into deterministic IMFs of different scales and a residual component.
The core is to screen the signals layer by layer and decompose the signals of different frequencies at different levels. The IMFs of each order must meet the following two requirements: (a) in the whole data sequence, the number of extremum and Zero-crossings must be equal or at most one difference; (b) the local mean of the data sequence is zero in any time. However, EMD is relatively easy to cause mode mixing. In order to avoid this problem, Wu and Huang 22 proposed EEMD that increases Gaussian white noise. Gauss white noise has the characteristics of an independent random sequence and a zero mean. And different scales of the signal can be compared and identified in an eigenmode function, then eliminate the white noise sequence by averaging the whole signal, so as to reduce the appearance of mode mixing problem. 35 The basic steps of EEMD algorithm are as follows: Step 1: Add a set of Gaussian white noise into the original sequence Step 2: Decompose the integrated sequence using the EMD algorithm to obtain IMFs Step 3: Repeat steps 1-2 with different scales of Gaussian white noise each time to obtain the corresponding IMFs where x (t) is the original signal, each C i (t) represents the i-th IMF, and R (t) is the residual component.

| Least squares support vector machine
On the basis of support vector machine, Suykens et al 36 proposed the least square support vector machine (LSSVM), which is an improvement on the standard support vector machine algorithm, making the algorithm simpler, easy to implement and enhance the computing speed, thus reducing the problem of high complexity of support vector machine. In LSSVM model, the training sample set S = { x i ,y i |i = 1,2 … t},x i = R n ,y i = R. Then, the optimal function is constructed by a high dimensional feature space, it described as follows: where f (x) stands for prediction value, (x) represents the nonlinear mapping function from the input space to the high-dimensional feature space, is weight and b is bias.
Its structural risk minimization is shown as follows: where 2 is the complexity; c stands for error penalty factor;R emp denotes the error range function. The objective function of LSSVM is as follows: WU and LIU i indicates the relaxation factor,y i is a constraint condition of the objective function.
Lagrange is used to solve the optimization problem: where i is the Lagrange multipliers According to the optimization conditions: We can get: According to the above formula, the optimization problem can be transformed into a process of solving linear equations, as follows: The final form of the LSSVM model is as follows: The radial basis function (RBF) is:

| Particle swarm optimization
Particle swarm optimization (PSO) was proposed by Dr Eberhart and Dr Kennedy in 1995 to study the predatory behavior of birds and was extended to solve optimization problems. 37 When using PSO algorithm to solve the optimization problem, the solution of each problem is regarded as a particle with no mass and no volume in the space, and the position of each particle is regarded as a solution of the space. The flight speed can be adjusted at any time according to the analysis results of the flying experience of individuals and groups. In each iteration of D-dimensional space, the particle optimizes its position and velocity by searching for two solutions, the individual extreme pbest and the global extreme gbest. The velocity and position update formula of the K-th iteration of particle I can be expressed as: where V t+1 i and V t i are the velocities of particle i at the next moment and the current moment; X t+1 i and X t i are the positions of particle i at the next moment and the current moment. represents the inertia weight, the larger makes the algorithm more inclined to global search, while the smaller makes the algorithm more inclined to local search. C 1 and C 2 are two learning factors, nonnegative constant. R 1 and R 2 are two uniform random sequences sampled from U [0,1].

| Bat alogrithm
Bat algorithm (BA) is a new random search optimization algorithm 38 proposed by Yang, which is inspired by the process of searching for food and avoiding obstacles by miniature bats in nature. It searches the solution space of problem by simulating the process of bat foraging, so as to achieve the optimization purpose. In order to find the optimal solution conveniently, we assume that bats forage in D-dimensional space, set at t-time, the position is X t i , the speed is V t i , and the current global optimal position is X i . The updated formulas of bat's position and velocity at (t + 1)-time are as follows: The random number ∈ [0,1], and F i is the search frequency emitted by bat i, F i ∈ F min ,F max .
The formula for updating the intensity and frequency of sound wave (pulse) are:

| PSO-LSSVM
For the LSSVM model with radial basis function, the kernel parameter and the normalized parameter have a great influence on the final prediction accuracy. 39 Therefore, we adopt PSO algorithm to optimize it. The steps are as follows: Step 1: Initialize the parameters of PSO and LSSVM Step 2: Calculate the fitness of each particle. Compare the current fitness value of each particle with the optimum fitness value and select the optimal one Step 3: Comparing the fitness value of each particle's optimal position with that of the population's optimal position, select the best and keep the current state of the particle Step 4: Judgment of optimizing termination conditions. If the error requirements or the maximum number of evolutions are reached, the optimization will end. Otherwise, repeat Step 2 to continue the optimization

| BA-LSSVM
Because the forecasting ability of LSSVM depends on the parameters, the BA is used to obtain the optimal parameters. Based on the BA-LSSVM model; the basic steps of parameter optimization can be described as follows: Step 1: Parameters setting. They are bat number M, maximum iteration number N, maximum pulse intensity A, and pulse frequency Step 2: Initialization population. Initialize the position X and speed V of the bat population. Each bat position is composed of LSSVM ( , ). The formula is: Step 3: Update parameters. Calculate the population fitness, find the current optimal solution, and update the flight speed and position of bats according to formula (11)- (13) Step 4: Update loudness and pulse frequency according to formula (14) and (15) Step 5: Output global optimal parameters. Repeat Steps 2-4 until the maximum number of iterations is reached and output the global optimal parameters. Then, a carbon price prediction model is performed.

PROPOSED HYBRID MODEL
In this paper, we used a multiscale ensemble prediction model EEMD-ADD to predict the carbon price sequence of Hubei ETS, as shown in Figure 1. The steps are as follows: Step 1: The carbon price sequence is decomposed into several IMFs and a residual by EEMD Step 2: According to the fine-to-coarse reconstruction algorithm, the IMFs are divided into low-frequency component and high-frequency component, the residual as the trend component

| Data selection and test
To push forward low-carbon procedure, the National Development and Reform Commission of China (NDRC) has established seven ETS pilots including Beijing, Tianjin, Shenzhen, Shanghai, Chongqing, Guangdong Province, and Hubei Province. 40 Among them, Hubei has the largest trading volume and a comprehensive business structure. The total volume of trade is 52.879 million tons, accounting for 35.4% of China; and the total volume of trade is 111.445 million-yuan, accounting for 31.4% of the total volume of trade in China. This paper analyses the carbon price series of Hubei 2014 (HB2014) from 2 April 2014 to 5 December 2018. Excluding the daily data of zero trading volume, there are 1129 data (data from China Carbon Emission Trading Network). The carbon price sequence is shown in Figure 2. Statistical descriptions on carbon price series are shown in Table 1.
HB2014 carbon price sequence is clearly irregular. In the initial stage, the sequence was relatively stable, followed by rapid decline and sharp fluctuations, then goes up by big margins, and finally tended to be stable. ADF is used to test the stability of carbon price series and BDS is used to test its linear characteristics. Table 2 is the ADF stability test result. This statistic meets the critical values of the 1%, 5% and 10% level. And the P-value is .1916. Obviously, the price sequence is nonstationary. Table 3 is the result of BDS test. The BDS value increases with the increase of embedding dimension, and the corresponding Z-statistics value is larger than the critical value. And the relevant P-value is 0, less than .05. Consequently, the carbon price sequence is nonlinear. To sum up, the price sequence of Hubei carbon trading market is nonstationary and nonlinear, making the sequence fluctuation characteristics basically the same as that of EU ETS.

| The analysis of IMF statistics
In this paper, the price sequence of Hubei carbon market is decomposed into nine IMFs and one residual by EEMD. As shown in Figure 3, the IMFs are arranged in order of frequency from high to low. To further investigate the results of modal decomposition, the mean period of each IMF, the F I G U R E 1 The framework of the EEMD-ADD model correlation between each IMF and the original sequence, and the variance and variance percentage of each IMF were calculated, see Table 4. The average period of IMF1-IMF5 is shorter than 100 days, which indicates that these fluctuations are quite intense. And the average cycle of the IMF6-IMF9 is over 200 days, with less fluctuation. The linear correlation between two sequences coefficient is described by the Pearson correlation, and whether the two sequences have a trend of simultaneous change is tested via the Kendall correlation coefficient. The correlation between the residual and the original sequence was the highest, reaching 0.3403 and 0.1985, but there was little fluctuation. The correlation coefficient of IMF6-IMF9 is relatively high, while that of IMF9 is negatively correlated with the original signal. The higher the variance ratio, the greater the influence of correlation IMF on the original sequence. The variance ratio of the residuals was 4.11%, which greatly influenced the original sequence. The trend of residuals is similar to the original sequence. Meanwhile, the sum of variance ratios of IMF6-IMF9 is 0.2%, it has little resemblance to the original sequence. It should be emphasized that due to the rounding error and the nonlinearity of the original sequence, the total variance of the IMF and residuals is not equal to the variance of the original sequence.
Through the fine-to-coarse reconstruction algorithm, the IMFs are decomposed into high-frequency and low-frequency components, 41 and the residual as the trend component. Calculate the mean and the T test of the sum of superposition sequences S i (i = 1,2, ⋯ ,9) from C i (t) to C i (t): The results are as Table 5.
The T test of IMF 5 (mean > 0.01) was significantly deviated from zero. Therefore, IMF1-IMF5 was taken as the high-frequency component, IMF6-IMF9 was taken as the low-frequency component, and residual as the trend component. Figure 4 shows the three components and the original sequence. Table 6 shows the statistical indicators of these sequences, including the Pearson and Kendall correlation between each component and the original sequence, the variance and variance percentage.

| Low-frequency componentimpacts of major events
The impact of major events on carbon price is mainly reflected in the low-frequency component. Such heterogeneous environment has marvelous impact on carbon prices, lasting for a long time, with a larger range, but with a lower frequency. For example, in April 2016, China set a goal of Intended Nationally Determined Contribution in the Paris Agreement, which caused significant fluctuations in carbon F I G U R E 3 The decomposition of the HB2014 sequence based on EEMD prices, reduced to less than 20 yuan and lasted until mid-July 2018. During this period, China held the 19th National People's Congress of the Communist Party of China. The government attaches great importance to energy conservation and emission reduction, resulting in the continuous rise in carbon prices. IMF6-IMF9 belongs to low-frequency component, with an average period of more than 200 days, which indicates that the impact of major events is difficult to eliminate and has a great impact on carbon price fluctuation. The Pearson and Kendall coefficients of low-frequency component are 0.7639 and 0.5443, shows that there is an obvious correlation between low-frequency component and the original sequence. And its variance contribution rate is higher than that of high-frequency component, so the fluctuations could be explained by major events. All in all, market participants should pay more attention to the occurrence of major events and judge the trend of carbon price according to its natures.

| High-frequency component-shortterm market fluctuations
Although the products in the carbon trading market are special, they still have general market attributes, that is, they are affected by the supply and demand mechanism of the market. The price rises when demand increases; the price decreases when demand decreases. Under the influence of market mechanism, the supply and demand of carbon emission rights are affected by energy prices, weather, political situation, and so on. These factors are short-lived and high in frequency, so their effect on carbon prices is mainly reflected in the high-frequency component. However, its impact tends to be relatively small.
The short average period and small variance contribution of high-frequency components indicate that they have little effect on the long-term trend of carbon prices, but have a greater short-term impact. Consequently, the establishment of short-term price rules is of great significance.

| Results
As shown in Figure 5,   Table 9. The analysis shows that: (a) The EEMD-ADD model has the highest accuracy, and the accuracy of PSO-LSSVM, BA-LSSVM, LSSVM, PSO, and BA decreases in turn. (b) The prediction accuracy of compound model is higher than that of single model. As can be seen from that the prediction trend performance of the multiscale ensemble prediction model is the highest. The above results show that the EEMD-ADD model is superior to signal models and some hybrid models at predicting the carbon price sequence.

| STABILIT Y TEST
In order to further verify the predictive performance of EEMD-ADD model, the other three-year carbon price sequence of Hubei Province and five other models are given for stability detection. The data are shown in Table 10. As shown in Table 11, the analyses and their conclusions are similar to those in 4.3.3. For these sequences, EEMD-ADD has the highest prediction accuracy and the best trend prediction performance. Moreover, when the sample size gradually decreases, the prediction accuracy will gradually decline, which shows that the prediction accuracy of large samples is higher. The empirical analysis in this paper not only collected a large number of sample data, but also decomposed the carbon price series into low-frequency, high-frequency components and trend terms through EEMD, which greatly reduced the impact of price fluctuation, so as to achieve higher accuracy. Thus, it can be verified that the EEMD-ADD prediction model is more effective than signal models and some hybrid models.

WORK
In this paper, a multiscale ensemble forecasting model is proposed for carbon price predicting. First, EEMD is used to decompose the original carbon price sequence into nine IMFs and one residual. Second, reconstruct these IMFs into three different frequency components by a fine-to-coarse reconstruction algorithm. Third, after comparative analysis, LSSVM, PSO-LSSVM, and BA-LSSVM are used to predict the three components, respectively. Fourth, combine the results to get the final predicted value. Using the carbon price of HB2014 for samples, the prediction performance of EEMD-ADD model has been compared with LSSVM, PSO, BA, PSO-LSSVM, and BA-LSSVM. And three additional carbon price series can also prove the performance of EEMD-ADD.
In general, we can draw some conclusions from the empirical analysis results: (a) EEMD can decompose the carbon price sequence into components with different frequencies, so it can capture the characteristics of each component more accurately. (b) By making full use of the advantages of different models, the complex characteristics of carbon spot price can be well described, so as to obtain more reliable and accurate prediction results. Therefore, the proposed model could provide significant improvements for carbon price prediction and is competitive at predicting nonlinear and nonstationary carbon prices.
In the future, (a) we should establish more advanced optimization algorithms and prediction models to further improve prediction performances. (b) We should study the changing characteristics of carbon price from regional and stage perspectives so as to predict its changing trend. (c) We need to analyze carbon prices in accordance with national policies and market rules in order to maximize profits.