A Decomposition-Integration Framework of Carbon Price Forecasting Based on Econometrics and Machine Learning Methods

Abstract

Carbon price forecasting and pricing are critical for stabilizing carbon markets, mitigating investment risks, and fostering economic development. This paper presents an advanced decomposition-integration framework which seamlessly integrates econometric models with machine learning techniques to enhance carbon price forecasting. First, the complete ensemble empirical mode decomposition with adaptive noise (CEEMDAN) method is employed to decompose carbon price data into distinct modal components, each defined by specific frequency characteristics. Then, Lempel–Ziv complexity and dispersion entropy algorithms are applied to analyze these components, facilitating the identification of their unique frequency attributes. The framework subsequently employs GARCH models for predicting high-frequency components and a gated recurrent unit (GRU) neural network optimized by the grey wolf algorithm for low-frequency components. Finally, the optimized GRU model is utilized to integrate these predictive outcomes nonlinearly, ensuring a comprehensive and precise forecast. Empirical evidence demonstrates that this framework not only accurately captures the diverse characteristics of different data components but also significantly outperforms traditional benchmark models in predictive accuracy. By optimizing the GRU model with the grey wolf optimizer (GWO) algorithm, the framework enhances both prediction stability and adaptability, while the nonlinear integration approach effectively mitigates error accumulation. This innovative framework offers a scientifically rigorous and efficient tool for carbon price forecasting, providing valuable insights for policymakers and market participants in carbon trading.

1. Introduction

As human society progresses rapidly, there is a consistent annual increase in global greenhouse gas emissions. In 2023, worldwide carbon emissions reached 37.4 billion tons, reflecting a 1.1% rise from 2022 (CO₂ Emissions in 2023-Analysis-IEA). In response, various countries and organizations have implemented measures to curb carbon dioxide emissions. One effective strategy for reducing carbon emissions is carbon pricing [1], which encompasses tools such as carbon trading markets and carbon taxes [2,3]. Carbon trading, as a market-based instrument, not only encourages energy conservation, but also accelerates advancements in low-carbon technologies and enhances corporate investments in clean energy [4,5]. When compared to carbon taxes, carbon trading tends to offer superior economic benefits in effectively mitigating carbon emissions [6], making it a particularly cost-effective approach.

According to the International Energy Agency and the International Carbon Action Partnership (Emissions Trading Worldwide: 2024 ICAP Status Report ∣ International Carbon Action Partnership (https://icapcarbonaction.com)), as of March 2024, there were 36 operational carbon markets worldwide, which collectively address nearly 20% of global greenhouse gas emissions and represent 60% of the world’s GDP. Furthermore, 22 countries and regions are in the process of contemplating the establishment of carbon markets. Notably, China, the largest emitter of carbon dioxide, accounted for 33% of global carbon emissions in 2023. To play a more significant role in global carbon reduction initiatives, China is actively investigating and implementing carbon trading markets. Given that carbon pricing is recognized as a highly effective approach to combating global warming, the accurate forecasting of carbon prices has emerged as a key focus in research concerning low-carbon transitions.

In September 2020, during the 75th session of the United Nations General Assembly, China unveiled its ambitious objectives of reaching peak carbon dioxide emissions by 2030 and attaining carbon neutrality by 2060, referred to as the “Dual Carbon Goals”. Since 2011, China has launched seven pilot carbon trading programs in regions such as Hubei, Guangdong, Beijing, Shanghai, Tianjin, Chongqing, and Shenzhen aimed at effectively managing greenhouse gas emissions through market-based approaches while fostering the growth of a green economy. Nonetheless, given its late initiation, the Chinese carbon market faces numerous challenges, including an underdeveloped regulatory framework, overly generous quota distributions, limited industrial participation and trading options, low trading volumes, and notable liquidity fluctuations [7]. Insights from these pilot markets reveal that China’s carbon market exhibits significant volatility in both trading volumes and prices. This dynamic behavior is complex and marked by frequent fluctuations, impacting various economic sectors such as industry, energy, and agriculture. Therefore, accurately forecasting carbon price trends is crucial for government policy formulation, enterprise strategic optimization, and investor decision-making [8,9].

Factors such as energy policy, economic growth, and unforeseen events significantly influence carbon price fluctuations [10,11,12,13]. These influences impart notable non-stationary and nonlinear characteristics to carbon pricing [14], complicating accurate forecasting efforts. Traditional econometric techniques, including the Autoregressive Integrated Moving Average (ARIMA) and Generalized Autoregressive Conditional Heteroskedasticity (GARCH), often fall short in effectively capturing these nonlinear aspects [15]. Consequently, many researchers have shifted towards machine learning approaches for carbon price prediction [16]. However, the reliability of these opaque predictive models raises concerns, especially when addressing the inherent nonlinearity and substantial volatility associated with carbon prices.

In recent years, advancements in signal processing technology have drawn significant attention to decomposition-integration frameworks [17]. The varying influences on carbon prices across different time scales [18] have led to a growing application of these frameworks in forecasting [19,20]. However, few studies have explored frameworks that effectively combine traditional econometric approaches with machine learning techniques. This gap highlights the necessity to identify an effective method for forecasting carbon prices in China, thereby assisting investors in navigating the current volatile market.

This paper introduces a hybrid forecasting framework, CEEMDAN-GWO-GRU/GARCH, which leverages complete ensemble empirical mode decomposition with adaptive noise (CEEMDAN) [21], the Grey Wolf Optimization (GWO) algorithm, the GARCH model, and a gated recurrent unit (GRU). The framework is designed to predict the carbon prices of the China Emissions Exchange Guangzhou (CEEX Guangzhou) from December 2013 to December 2022. Initially, the price series is decomposed into multiple subsequences with varying frequencies using CEEMDAN. Following this, high-frequency and low-frequency sequences are determined through Lempel–Ziv complexity and dispersion entropy (DispEn). For each group of subsequences, models based on GARCH and the GRU, optimized by GWO, are developed. Ultimately, the optimized GRU model aggregates all predictive values to produce the final forecast.

The contributions and innovations presented in this paper are summarized as follows:

(1) We propose a new carbon price forecasting decomposition-integration framework, CEEMDAN-GWO-GRU/GARCH-GRU, which is particularly effective for handling complex, non-stationary, and nonlinear carbon price time series data. The framework decomposes the carbon price series into several components with different frequency characteristics using CEEMDAN’s adaptive noise mechanism. The frequency characteristics of each component are then identified using the Lempel–Ziv algorithm and dispersion entropy. This decomposition and identification process ensures the model’s thorough utilization of different frequency components, thereby improving forecasting accuracy.

(2) Unlike studies that rely solely on econometric models or machine learning methods for forecasting, the new decomposition-integration framework proposed in this paper effectively combines the strengths of both econometric models and machine learning methods to optimize predictions of different frequency components. The GARCH model is used to forecast high-frequency components, which leverages its strength in capturing volatility clustering in high-frequency data. For low-frequency components, the framework employs a GWO-optimized GRU model, which effectively integrates the GRU’s ability to handle long-term dependencies with GWO’s global optimization capabilities. This effective fusion of econometric models and machine learning methods enables the model to accurately capture the different frequency characteristics of carbon price data, thereby enhancing overall forecasting performance.

(3) We introduce the GWO-GRU model to nonlinearly integrate the forecast results of different frequency components, significantly improving the issue of error accumulation commonly found in traditional linear integration methods and ensuring optimal parameter configuration during the integration process. The comprehensive integration technique not only reduces the errors of traditional single models or standard hybrid methods but also enhances the forecasting framework’s adaptability to carbon price fluctuations, thereby increasing the overall accuracy and stability of the predictions.

Additionally, a comprehensive series of empirical analyses is conducted to validate the effectiveness and robustness of the framework. This evaluation includes metrics such as the Root Mean Square Error (RMSE), Mean Absolute Error (MAE), and Mean Absolute Percentage Error (MAPE). Insights gained from linear regression and the Diebold–Mariano test (DM test) [22,23] comparing predicted and actual values underscore the framework’s strengths in carbon price forecasting. Moreover, robustness tests reveal that the framework maintains strong performance even when variations occur in model hyperparameters and training set lengths.

The subsequent sections are organized as follows: Section 2 reviews the existing literature on carbon price forecasting and identifies current research limitations. Section 3 presents the underlying methodologies and theories. Section 4 elaborates on the development of the proposed framework. Section 5 focuses on the empirical analysis, while Section 6 concludes the paper and suggests avenues for future research.

2. Literature Review

Carbon prices are influenced by various factors, including economic conditions and policy decisions, resulting in characteristics such as nonlinearity, volatility, and uncertainty [24]. Currently, two primary approaches dominate the research landscape regarding carbon price forecasting. The first approach focuses on predicting prices using trending factors, while the second relies on historical data [20].

Within the realm of forecasting based on trending factors, Ref. [15] argues that numerous external elements impact carbon prices. Relying solely on historical data may hinder the development of advanced algorithms capable of accurately representing these influencing factors. Fortunately, a wealth of literature has thoroughly examined these variables. For instance, Ref. [25] investigated the European carbon market from 2005 to 2007, revealing a connection between carbon prices, energy prices, and temperature variations. Ref. [26] established a positive correlation between carbon prices and industrial production using a Markov transformation model, noting that Brent crude oil prices exert the most significant influence. Additionally, Ref. [27] identified fuel prices and the dynamics of supply and demand in the energy market as primary drivers of carbon prices. Ref. [28] employed mixed-data sampling regression models (MIDAS) to forecast European carbon prices, considering economic and energy indicators, and discovered that predictions based on coal prices yield the best outcomes, while those based on crude oil prices produce comparatively higher results. Ref. [29] posits that incorporating energy price factors, including oil, coal, and natural gas, can enhance predictive accuracy. However, the selection of influential factors is often subjective and incomplete, which can lead to cumulative errors and inaccurate forecasts, particularly when addressing the nonlinear nature of carbon prices [30].

The second category, derived from historical data, encompasses three primary methodologies: traditional econometric models, individual machine learning models, and decomposition-integration frameworks [4]. In the realm of traditional econometric models, ARIMA and various autoregressive approaches are commonly employed for forecasting carbon prices [31,32]. Ref. [26] demonstrated the effectiveness of a non-parametric model in predicting the BlueNext spot price and future prices from the European Climate Exchange (ECX), revealing that its prediction error rate outperforms that of autoregressive models. Ref. [33] found that GARCH models yield superior predictive performance compared to implied volatility and the K-nearest neighbor model for carbon price fluctuations. Nonetheless, these traditional econometric models struggle to accommodate the evident nonlinear characteristics of carbon price volatility [15].

In the context of big data and machine learning advancements, researchers have increasingly turned to artificial intelligence models for carbon price predictions. Ref. [34] successfully employed a Back-Propagation Neural Network (BPNN) to forecast prices in the Fujian carbon exchange. Ref. [35] introduced an explainable machine learning (XML) framework, effectively predicting and interpreting carbon prices across three major markets in Shenzhen, Hubei, and Beijing, China. Ref. [36] ascertained that a Genetic Algorithm combined with Extreme Learning Machines (ELM) offers effective forecasting for the Beijing carbon exchange, utilizing a robust performance evaluation system. Ref. [37] focused on the Chinese carbon market, constructing a CNN-LSTM model to predict trading prices, achieving accuracy comparable to other leading models, including CEEMDAN, Boosting, and GRUs. Furthermore, LSTM, GRUs, and related models have gained considerable traction in carbon price forecasting [38,39,40].

Decomposition-integration frameworks typically encompass five key stages: decomposition, reconstruction, prediction, integration, and evaluation [41]. Breaking down carbon prices allows for the effective capture of sequence characteristics, revealing inherent patterns within various frequency data and ultimately enhancing predictive accuracy [42]. Ref. [43] suggests that a hybrid model, integrating Empirical Mode Decomposition (EMD), Least-Squares Support Vector Machines (LSSVM), and GARCH with the Lempel–Ziv complexity algorithm, is highly effective for forecasting carbon prices. Due to modal aliasing issues associated with the EMD method [44], Ref. [20] developed an EEMD-BPNN-ELM model using ensemble empirical mode decomposition (EEMD), while [44] introduced a multi-scale ensemble prediction model built on EEMD.

In earlier studies, Ref. [5] employed Complementary Ensemble Empirical Mode Decomposition (CEEMD) to break down and reconstruct carbon prices into three distinct elements: trend, periodic fluctuations, and random fluctuations, based on the attributes of each subsequence. Their framework combined a co-integration model (CIM), ant colony algorithm (ACA), grey neural network (GNN), and GARCH, achieving impressive predictive performance. Ref. [45] applied six machine learning frameworks to forecast the carbon price across eight Chinese exchanges, finding that the radial basis function neural network (RBFNN) and kernel-based nonlinear extension of the Arps decline model (KNEA), optimized using GWO with CEEMDAN, yielded the most accurate and consistent predictions. Ref. [46] proposed a refined decomposed ensemble prediction model that employs time convolutional networks with varying structures to forecast carbon prices, enhancing predictive performance while minimizing structural complexity and computational costs. Empirical findings highlight the proposed model’s superiority and resilience. Ref. [47] developed a hybrid forecasting model based on a decomposition-integration framework called VMD-CNN-BiLSTM-MLP to forecast EUA carbon futures prices for two different five-year periods before and after the implementation of emission reduction policies. The results indicate that their model outperformed others. Alongside these mainstream decomposition approaches, some researchers have applied secondary decomposition to time series data, extracting additional characteristic information [48,49].

Despite substantial advances in forecasting carbon prices using trend factors and historical data, existing studies still exhibit the following shortcomings in practical applications: (1) Carbon price forecasting employs traditional econometric models, AI models, and hybrids, each with distinct advantages but also limitations in addressing the non-stationary, multi-frequency, and nonlinear dynamics of carbon prices. The models’ sensitivity and adaptability often hinder forecast accuracy and stability. There is a clear need to refine decomposition ensemble frameworks to mitigate modal mixing and optimize ensemble strategies, ultimately enhancing forecast precision. (2) Econometric models are well suited for capturing relationships between economic variables, while machine learning excels at handling nonlinear, high-dimensional data. However, better integrating these approaches to improve the accuracy, reliability, and practical application of carbon price forecasts remains a key challenge. (3) Error accumulation is a pervasive issue in time series forecasting, especially in volatile markets like carbon. Advancing nonlinear ensemble methods to flexibly adjust forecast weights and combinations, thereby dispersing errors and improving stability, is a critical area for further research. (4) Manual hyperparameter adjustment is costly in terms of calculation and easily falls to the local optimum, making the optimal set difficult to find. Therefore, we introduce a new decomposition-integration framework for carbon price forecasting, CEEMDAN-GWO-GRU/GARCH-GWGRU, which integrates econometric models with machine learning to address these limitations.

3. Methodology

3.1. CEEMDAN

Adaptive noise-complete ensemble empirical modal decomposition (CEEMDAN) comprises a noise-assisted decomposition algorithm designed to overcome the limitations of empirical modal decomposition (EMD) [50] and ensemble empirical modal decomposition (EEMD) [21] regarding mode aliasing and the effects of white noise. Its primary objective is to more effectively address high-frequency noise and the non-stationary characteristics present in signals.

The EEMD [51] algorithm mitigates mode aliasing by adding auxiliary white noise with different frequencies and uniform distributions to the target signal multiple times, followed by separate EMDs. This process involves averaging the intrinsic modal functions (IMFs) and the trend components obtained from each decomposition. While EEMD represents a significant advancement over traditional EMD, the introduction of white noise can still interfere with the decomposition process, affecting the purity of the modal components and consequently diminishing the accuracy of subsequent analyses.

To address this limitation, the CEEMDAN algorithm incorporates a more refined method for adaptive noise processing. During each iteration of the EMD, it adds adaptive white noise that matches the characteristics of the current remaining signal. This targeted approach effectively smooths out the interference impulses introduced by the white noise, thereby reducing its negative impact on the decomposition results and further enhancing both the accuracy and the stability of the process.

The operational procedure of the CEEMDAN algorithm is as follows:

(1) The CEEMDAN algorithm decomposes the original sequence $x\left(t\right)+a_0{m}_j$ in N calculations, the parameter $a_0$ is used to control the signal-to-noise ratio between the additional noise and the original signal, and the first modal component can be expressed as follows:

$$\overline{IM{F}_1\left(t\right)}={\displaystyle \frac{1}{N}}\sum _{j=1}^{N}IM{F}_1^j\left(t\right).$$

(1)

(2) At this point, the residual signal can be expressed as $r_1\left(t\right)=x\left(t\right)-\overline{IM{F}_1\left(t\right)}$.

(3) Before the next decomposition step, the white noise signal is introduced and self-decomposed, combining the first stage’s components and the residual signal. The signal to be decomposed can be expressed as $r_1\left(t\right)+a_{1}em{d}_1\left(m_j\left(t\right)\right)$, where the $em{d}_1\left(\right)$ operator denotes the ith modal component generated using the EMD algorithm. Subsequently, the second modal component is integrated and averaged using the EEMD algorithm to obtain the following:

$$\overline{IM{F}_2\left(t\right)}={\displaystyle \frac{1}{N}}\sum _{j=1}^{N}em{d}_1(r_1\left(t\right)+a_{1}em{d}_1\left(m_j\left(t\right)\right)).$$

(2)

(4) Continuing to calculate the residual signal using the method of step (2), the $m+1$ modal components are obtained as follows:

$$\overline{IM{F}_{m+1}\left(t\right)}={\displaystyle \frac{1}{N}}\sum _{j=1}^{N}em{d}_1(r_m\left(t\right)+a_{m}em{d}_m\left(m_j\left(t\right)\right)).$$

(3)

(5) The decomposition algorithm terminates when the final residual signal obtained satisfies the iteration conditions. Assuming the final division into N modal components, the original sequence can be represented as follows:

$$x\left(t\right)=\sum _{j=1}^N\overline{IM{F}_j\left(t\right)}+r\left(t\right).$$

(4)

3.2. Lempel–Ziv Algorithm

The Lempel–Ziv complexity algorithm has a wide range of applications in time series analysis, mainly through replication and insertion operations to measure the complexity of the time series, which are then used for classification, identification, and other tasks. The main principle is that, if the data series contains many repetitive patterns, then it is easy to compress, so its complexity is low. On the contrary, if the data series contains many unique patterns, then it is more difficult to compress, so its complexity is high [52]. The specific procedure for calculating the Lempel–Ziv complexity is as follows:

Step 1 Data Preprocessing: Given a sequence $X_N=\{x_1,x_2,\dots ,x_N\}$, the original sequence $X_N$ is first converted to a binary sequence $S_N=\{s_1,s_2,\dots ,s_N\}$. This is achieved by assigning a value of 1 to the elements of the sequence that are greater than or equal to the median, and otherwise assigning a value of 0. The purpose of this step is to convert the time series into a format suitable for processing by the Lempel–Ziv algorithm.

Step 2 Complexity Computation: Next, the temporary character variables $S_{v,0}$ and $Q_0$ are initialized to be empty characters and the initial value of complexity is set to $C_N\left(0\right)=0$. Then, each character of the binary sequence $S_N$ is traversed. For each character $s_r$, we define the string $Q_r=\left\{Q_{r-1}{s}_r\right\}$, and determine whether $Q_r$ already appears in the previous string $S_{v,r-1}=\left\{S_{v,r-2}{s}_{r-1}\right\}$:

• If $Q_r$ already exists, the complexity is not increased and the process continues with the next character.
• If $Q_r$ does not exist, the complexity is increased by 1 ($C_N\left(r\right)=C_N(r-1)+1$), $Q_r=\left\{\right\}$, and then the process continues with the next character. This process continues until the last character in the sequence.

Step 3 Normalization: The computed complexity $C_N\left(N\right)$ is a non-negative integer value, but we usually normalize it in order to facilitate comparison between different sequences. The normalized complexity formula is as follows:

$$C_{nN}={\displaystyle \frac{C_N\left(N\right)}{{lim}_{N\to \infty }{C}_N\left(N\right)}}\approx {\displaystyle \frac{C_N\left(N\right){log}_{k}N}{N}},$$

(5)

where k is the number of different characters in the binary sequence $S_N$.

3.3. Dispersion Entropy (DispEn)

The concept of entropy is derived from system thermodynamics and is often used to quantify the complexity and irregularity of data and to reveal dynamic properties and underlying patterns in time series data. Dispersion entropy (DispEn) is a nonlinear dynamics method used to analyze the complexity of time series. Compared to traditional methods such as sample entropy, dispersion entropy is able to detect changes in signal amplitude and frequency at the same time, thus assessing the signal complexity more comprehensively. Its calculation steps are as follows:

(1) First, the normal distribution function is used as follows:

$$y_j={\displaystyle \frac{1}{\sigma \sqrt{2\pi }}}{\int }_{-\infty }^{x_j}exp\left(-{\displaystyle \frac{{(t-\mu )}^2}{2{\sigma }^2}}\right)dt.$$

(6)

Mapping the timing data x to $y=\{y_1,y_2,\dots ,y_N\}$, $y\in (0,1)$. $\sigma$ and $\mu$ denote the standard deviation and mean, respectively.

(2) Each $y_j$ is assigned to an integer of $[1,2,\dots ,c]$ by a linear transformation, i.e.,

$$u_j^c=round(c{y}_j+0.5),$$

(7)

where $u_j^c$ represents the jth element of the categorized signal obtained after linear transformation; c represents the number of categories; $round$ represents the rounding function.

(3) The embedding vector $u_i^{m,c}$ is computed. An embedding vector $u_i^{m,c}$ is created with embedding dimension m and time delay d, i.e.,

$$u_i^{m,c}=\{u_i^c,u_{i+d}^c,\dots ,u_{i+(m-1)d}^c\},i=1,2,\dots ,N-(m-1)d.$$

(8)

(4) Each sequence $u_i^{m,c}$ is mapped into ${\pi }_{v_0{v}_1\dots v_{m-1}}$, which is the diffuse mode $u_i^{m,c}$, i.e.,

$$u_i^c=v_0,u_{i+d}^c=v_1\dots u_{i+(m-1)d}^c=v_{m-1},$$

(9)

where the number of diffuse modes is $c^m$.

(5) For $c^m$ diffuse modes ${\pi }_{v_0{v}_1\dots v_{m-1}}$, the probability of each diffuse mode is as follows:

$$p\left({\pi }_{v_0{v}_1\dots v_{m-1}}\right)={\displaystyle \frac{Number\left\{{\pi }_{v_0{v}_1\dots v_{m-1}}\right\}}{N-(m-1)d}}.$$

(10)

where $Number\left\{{\pi }_{v_0{v}_1\dots v_{m-1}}\right\}$ denotes the number of diffuse patterns ${\pi }_{v_0{v}_1\dots v_{m-1}}$ assigned to the sequence $u_i^{m,c}$ and $p\left({\pi }_{v_0{v}_1\dots v_{m-1}}\right)$ is the number of diffuse patterns assigned to the sequence divided by the total number of embedding vectors.

(6) Finally, based on the definition of Shannon’s entropy, the value of Dispersal Entropy (DisPEn) is calculated as follows:

$$DisPEn(x,m,c,d)=-\sum _{\pi =1}^{c^m}p\left({\pi }_{v_0{v}_1\dots v_{m-1}}\right)ln\left(p\left({\pi }_{v_0{v}_1\dots v_{m-1}}\right)\right).$$

(11)

From Equation (11), it can be seen that the embedding dimension m, the number of classes c, and the delay d affect the results.

3.4. GARCH

The GARCH model is an important extension of the ARCH model [53]. Although its structure resembles ARCH’s, the lag order is often less. The GARCH model is defined by Formulas (12) and (13), where $p\ge 0$, $q\ge 0$, $a_0>0$, $a_i\ge 0(i=1,\dots ,q)$, ${\beta }_i\ge 0(i=1,\dots ,p)$.

$$y_t={\mu }_t+{ϵ}_t,{ϵ}_t=e_t{\sigma }_t,e_t\sim i.i.N(0,1),$$

(12)

$${\sigma }_t^2={\alpha }_0+\sum _{i=1}^q{\alpha }_i{ϵ}_{t-i}^2+\sum _{j=1}^p{\beta }_j{\sigma }_{t-j}^2={\alpha }_0+\alpha \left(L\right){ϵ}_t^2+\beta \left(L\right){\sigma }_t^2.$$

(13)

The model satisfying the above conditions is called GARCH(p,q). GARCH(1,1) is the simplest GARCH model, whose conditional variance function is defined as Formula (14), where ${\alpha }_0>0$, $\alpha \ge 0$, $\beta \ge 0$. When $\alpha +\beta <1$, GARCH(1,1) is smooth.

$${\sigma }_t^2={\alpha }_0+\alpha {ϵ}_{t-1}^2+\beta {\sigma }_{t-1}^2.$$

(14)

3.5. GRU

The GRU, an improved algorithm based on LSTM, has a simple structure and can also avoid the long dependence problem [54]. Figure 1 illustrates the internal architecture of the GRU, which comprises an update gate and a reset gate. The update gate regulates the number of historical data retained at the current time step, while the reset gate controls the extent to which previous information is disregarded.

The memory of the last time and the current state is expressed as $h_{t-1}$ and $h_t$, respectively, the input time series is $x_t$, and $y_t$ represents the output sequence. First, the reset gate and update gate states are calculated by Formulas (15) and (16), with the input vector $x_t$ and $h_{t-1}$ transferred by the sigmoid function. Second, Formula (17) obtains new memory contents ${\tilde{h}}_t$ using the tanh function with the linear transformation of the reset gate states $r_t$, $h_{t-1}$, and $x_t$. Third, Formula (18) calculates the current state $h_t$. Finally, the final result $y_t$ is output by the sigmoid function of Formula (19). Formulas (15)–(19) show the calculation process of the GRU, where $W_r$, $W_z$, $W_{{\tilde{h}}_t}$, and $W_o$ are weight matrices, ${\tilde{h}}_t$ represents new memory contents, $\sigma (·)$ is the sigmoid activation function, and $tanh(·)$ is the tanh activation function:

$$r_t=\sigma \left(W_t[h_{t-1},x_t]\right),$$

(15)

$$z_t=\sigma \left(W_z[h_{t-1},x_t]\right),$$

(16)

$${\tilde{h}}_t=tanh\left(W_{{\tilde{h}}_t}[h_{t-1}\times r_t,x_t]\right),$$

(17)

$$h_t=(1-z_t)\times h_{t-1}+z_t\times {\tilde{h}}_t,$$

(18)

$$y_t=\sigma \left(W_o{h}_t\right).$$

(19)

3.6. GWO

Grey Wolf Optimization (GWO) is an intelligent optimization algorithm which completes optimization tasks by simulating the hunting mechanism of grey wolves in nature, including searching, surrounding, and attacking [55]. The algorithm divides the data into four fitness levels: $\alpha$, $\beta$, $\delta$, and $\omega$, where $\alpha$ means the optimal solution, $\beta$ means the sub-optimal solution, and $\delta$ means the worst solution. The others are $\omega$ wolves, responsible for rounding up prey. In the process of hunting, the grey wolves gradually narrow the solution space by rounding up prey, and finally approach the optimal solution. The location of the prey corresponds to the solution of the optimization problem.

In this study, we apply the GWO algorithm to optimize the hyperparameters of the GRU carbon price forecasting model. Specifically, we optimize three key hyperparameters of the model through the GWO algorithm: learning rate, batch size, epochs. These hyperparameters are crucial to the training effectiveness and prediction accuracy of the model, so their optimization can significantly improve model performance.

In this scenario, the “location” in the GWO algorithm represents the value of the GRU hyperparameter. In the GRU model, the “position” of the grey wolves can be defined as the combined vector of hyperparameters, $X=\{learningrate,batchsize,epochs\}$. By iteratively updating the location of the grey wolves, we can continuously adjust these hyperparameters to achieve the optimal prediction effect. The core of the GWO algorithm consists of the following steps:

(1) The distance vector D between the wolves and the prey is determined:

$$D=|C·X_p\left(t\right)-X\left(t\right)|,C=2r_1,$$

(20)

where $X_p\left(t\right)$ represents the location of the prey at t generation; $X\left(t\right)$ represents the location of the grey wolves at t generation; C is the constant vector; $r_1$ is a random number in the range of $[0,1]$.

(2) The way grey wolves update as they move toward the prey can be described as follows:

$$\left\{\begin{array}{c}{D}_{\alpha }=\mid C_1·X_{\alpha }\left(t\right)-X_u\left(t\right)\mid ,\hfill \\ D_{\beta }=\mid C_2·X_{\beta }\left(t\right)-X_u\left(t\right)\mid ,\hfill \\ D_{\delta }=\mid C_3·X_{\delta }\left(t\right)-X_u\left(t\right)\mid ,\hfill \\ X_1=X_{\alpha }-A_1·D_{\alpha },\hfill \\ X_2=X_{\beta }-I_2·D_{\beta },\hfill \\ X_3=X_{\delta }-A_3·D_{\delta },\hfill \\ X_u(t+1)={\displaystyle \frac{X_1+X_2+X_3}{3}},\hfill \end{array}\right.$$

(21)

$$A=2\alpha r_2-\alpha ,$$

(22)

where A is the convergence factor, and $r_2$ is a random number in the range of $[0,1]$. a decreases from 2 to 0 as the iteration increases. In this way, GWO is able to efficiently explore the hyperparameter space and optimize the learning rate, batch size, and epochs in the GRU model, thereby improving the predictive performance of the model.

4. The Proposed Framework

This paper combines CEEMDAN, GARCH, and GWO for hyperparameter optimization to propose a nonlinear decomposition-integration framework named CEEMDAN-GWO-GRU/GARCH-GRU. The model construction process mainly includes three stages, as shown in Figure 2.

4.1. Stage 1: Original Series Decomposition and Frequency Identification

To begin with, the framework utilizes CEEMDAN to break down the original time series into multiple simple subsequences that vary in frequency and amplitude. These subsequences, referred to as Intrinsic Mode Functions (IMFs), are organized from the highest to the lowest frequency. The decomposed carbon price consists of three main types of component: high frequency (volatility), medium frequency (periodicity), and low frequency (tendency), each of which has a different significance for policy and trading decisions. High-frequency volatility reflects the short-term market sentiment and unexpected events. For traders, high-frequency volatility provides important information for capturing short-term arbitrage opportunities, but also an important basis for quantitative trading model optimization. However, these volatilities can lead to market instability, so policymakers need to mitigate risks by introducing price volatility limits or providing market stability instruments, such as carbon futures. The periodicity of the medium frequency is usually related to seasonal factors (such as changes in energy demand) and the economic cycle. This component can provide a reference for the formulation of medium-term carbon emission targets. For example, policymakers can use the mid-frequency cycle trend to predict the range of carbon price fluctuations over a period of time in the future, thereby developing flexible emission quota allocation strategies. The low-frequency trend reflects the long-term development of the carbon market, which is influenced by macroeconomic policies and carbon reduction targets. It provides an important basis for policymakers who are developing long-term carbon pricing mechanisms. Traders can also adjust their portfolios according to low-frequency trends. By combining the characteristics of each component, policymakers and market participants can more precisely understand the dynamics of the carbon market and make decisions accordingly, thus achieving both economic benefits and emission reduction targets.

This paper uses the Lempel–Ziv algorithm and DispEn to identify the sequence frequency. Through different types of complexity recognition algorithms, high-frequency and low-frequency subsequences can be identified more objectively and effectively. The specific identification process is as follows:

(1) The complexity of the ith subsequence is calculated by the Lempel–Ziv algorithm or DispEn, where $C_i$, $i=1,2,\dots ,N$.

(2) The critical value ${\lambda }_0$ (≥0.8) [56] is set, and the minimum value m is found, satisfying Formula (25).

$$\eta ={\displaystyle \frac{{\sum }_{i=1}^m{C}_i}{{\sum }_{i=1}^n{C}_i}}\ge {\lambda }_0,m\in k.$$

(23)

(3) Subsequences 1 to m are recognized as high-frequency subsequences and $m+1$ to N as low-frequency subsequences.

4.2. Stage 2: Different-Frequency Subsequence Forecasting

Due to the distinct nature of high-frequency and low-frequency subsequences, this study adopts tailored predictive techniques for each type. Specifically, the GARCH model is employed to forecast high-frequency subsequences, while the GRU is utilized for predicting low-frequency ones.

High-frequency subsequences are characterized by considerable randomness and unpredictability, which reflects the short-term volatility of carbon prices. The GARCH model is well suited for capturing the features of high volatility, making it an effective tool for forecasting these subsequences [57].

In contrast, low-frequency subsequences exhibit relatively lower volatility and possess simpler structures, primarily reflecting the cyclical and trend-based nature of carbon prices. According to prior studies, the GRU outperforms other econometric and machine learning methods in predicting cyclical and trend series [58,59]. As a result, the GRU is the preferred choice for forecasting low-frequency subsequences. To further enhance the predictive performance of the GRU, the GWO algorithm is integrated for hyperparameter optimization, improving the overall accuracy of the framework.

4.3. Stage 3: Nonlinear Integration of Subsequences

This study employs a nonlinear integration approach to combine the predictive values derived from high-frequency and low-frequency subsequences. The input consists of the summed predictive values of each subsequence from the training set, with the actual prices from that same set serving as the output. The GWO-GRU model is trained to understand the error relationships between the predicted values and the actual prices. Subsequently, the total predictive values from the test set are fed into the trained GWO-GRU model, which corrects the prediction errors to generate the final carbon price results.

4.4. Evaluation and Statistical Tests

This paper compares various forecasting models and frameworks to verify their accuracy and validity. Three common loss functions, RMSE, MAE, and MAPE, are used to evaluate the model. Compared with a single evaluation indicator, applying three different evaluation indicators can obtain a more objective and comprehensive model and framework. Their definitions are shown as Formulas (24)–(26), where $x\left(t\right)$ and $\widehat{x}\left(t\right)$ represent the real value and predictive value of the test set at the moment t, and n is the sample size of the test set. The three evaluation indicators reflect the errors between the predictive and real values from different perspectives. The smaller the value of the evaluation indicators, the better and more reliable the prediction effect of the models and frameworks.

$$RMSE=\sqrt{{\displaystyle \frac{1}{n}}\sum _{t=1}^n{(x\left(t\right)-\widehat{x}\left(t\right))}^2},$$

(24)

$$MAE={\displaystyle \frac{1}{n}}\sum _{t=1}^n\mid x\left(t\right)-\widehat{x}\left(t\right)\mid ,$$

(25)

$$MAPE={\displaystyle \frac{1}{n}}\sum _{t=1}^n{\displaystyle \frac{\mid x\left(t\right)-\widehat{x}\left(t\right)\mid }{x\left(t\right)}}.$$

(26)

Additionally, the Diebold–Mariano test [22] is also used to test the significance of the prediction results of the models and frameworks. The DM test, a classical statistical test method, is often used to test whether the prediction effect of two models is significantly different. Through the RMSE, MAE, MAPE, and DM tests, this paper can make an effective and objective evaluation of the proposed decomposition-integration framework, CEEMDAN-GWO-GRU/GARCH-GRU.

5. Empirical Research

5.1. Data Description and Normalization

This paper selects the daily carbon trading prices from the Guangzhou Carbon Emissions Exchange as the research sample. The sample dataset covers trading data from 19 December 2013 to 3 March 2023, totaling 1969 data points. The data are obtained from the official website. The China Emissions Exchange Guangzhou (CEEX Guangzhou), one of the earliest carbon markets in China, started trial operation in December 2013. The trading volume of CEEX Guangzhou continues to lead the whole of China. As of August 2023, CEEX Guangzhou continues to maintain a leading position in the national carbon market, with a total carbon allowance trading volume exceeding 298 million tons and a total transaction value of over CNY 6.789 billion, making it the world’s largest spot carbon allowance trading market. Figure 3 shows the carbon price series of CEEX Guangzhou.

In this study, the first 90% of the data are classified as the training set and the remaining samples are the test set. In

Table A1. CEEX Guangzhou carbon price data description.

CEEX Guangzhou	Sample		Size	Date Range
Data	Total sample		1969	2013.12.09∼2023.03.03
	Training sample		1772	2013.12.09∼2022.05.19
	Test sample		197	2022.05.20∼2023.03.03
Description	Mean	Median	Std.	Kurtosis	Skewness
	32.201	25.62	21.987	−0.051	1.145
S-W test	W-statistic			p-value
	0.804			<0.001
ADF test	T-statistic	p-value	1%	5%	10%
	−0.131	0.946	−3.434	−2.863	−2.568

, we give the relevant sample data along with the descriptive statistics. According to the Shapiro–Wilk (S-W) test, the dataset does not conform to normal distribution. Combined with the ADF test and the sample standard deviation, it indicates that the price of CEEX Guangzhou is unstable.

To reduce the impact of data noise and improve the prediction performance of the model, this paper normalizes the carbon price data and processes the data to the range of 0 to 1 using Formula (27), where $y_{max}$ and $y_{min}$ equal 1 and 0, respectively, representing the maximum and minimum values after data normalization, and $x_{max}$ and $x_{min}$ are the maximum and minimum values of the original data. The predictive value is finally de-normalized by Formula (28), where $\widehat{y}$ is the predictive value of the original data after normalization processing, and $\widehat{x}$ is the final predictive result.

$$y={\displaystyle \frac{x-x_{min}}{x_{max}-x_{min}}},$$

(27)

$$\widehat{x}=\widehat{y}(x_{max}-x_{min})+x_{min}.$$

(28)

5.2. Carbon Price Sequence Decomposition and Complexity Identification

The descriptive statistics, the results from the ADF test presented in Table A1, and the price trend illustrated in Figure 3 indicate that the carbon pricing at CEEX Guangzhou exhibits both nonlinearity and non-stationarity. Consequently, it is necessary to decompose the carbon price time series to obtain stable subsequences. To achieve this, the study utilizes the CEEMDAN algorithm, which effectively mitigates mode aliasing issues and decreases computational requirements. The detailed outcomes of this decomposition are illustrated in Figure 4.

The subsequences are arranged according to frequency from high to low, and the last item (IMF10) is the residual. Compared with the original series, we find that the decomposed subsequence is characterized by clearer features, stronger regularity, and a smoother trend, which helps to improve the accuracy of carbon price forecasting.

Following the decomposition of the carbon price for CEEX Guangzhou, the subsequent phase involves applying both the Lempel–Ziv complexity algorithm and DispEn to discern high-frequency and low-frequency subsequences. The results of this complexity analysis are displayed in Figure 5. Both methods indicate that IMF5 attains the critical threshold, with values of 87.33% and 80.47%, respectively. Thus, the high-frequency components are identified as IMF1 through IMF5, while the low-frequency components consist of IMF6 to IMF10.

5.3. One-Step-Ahead Forecasting

5.3.1. Model Selection and Hyperparameter Setting

We select GARCH, the GRU, LSTM, Bidirectional LSTM (BiLSTM), Support Vector Regression (SVR), an Elman neural network, and their combined models as comparison models to validate the performance and reliability of our proposed framework.

Figure 6 illustrates that, during the model selection process, we have the option to either apply the CEEMDAN method for decomposing carbon prices in the initial phase or move directly to the prediction stage. In the second phase, GARCH or ELMAN is chosen for high-frequency sequences, while the GRU, LSTM, and BiLSTM are employed for low-frequency sequences. For the final integration stage, one can opt for either a straightforward summation method or a nonlinear integration approach using the GRU. Furthermore, the use of the GWO algorithm for model optimization is also an option. Consequently, we arrive at a total of 29 models, ultimately selecting the GARCH and GWO-GRU models.

This study evaluates the proposed framework against 29 models through one-step-ahead forecasting to demonstrate its effectiveness. One-step-ahead forecasting entails predicting the next day’s price based solely on preceding actual prices. We select a variety of established carbon price forecasting models to assess the performance and reliability of both the decomposition-integration approach and our proposed framework. These models include GARCH, the GRU, LSTM, Bidirectional LSTM (BiLSTM), Support Vector Regression (SVR), and Elman Recurrent Neural Networks (ELMAN). It is noteworthy that GARCH is deemed unsuitable for the decomposition-integration method due to the absence of the ARCH effect in the decomposed low-frequency sequences.

As depicted in Figure 6, the initial stage of sequence decomposition allows the use of CEEMDAN to decompose carbon prices or allows one to proceed directly to forecasting without decomposition. In the subsequent stage of model selection, GARCH or ELMAN is applied to high-frequency sequences, whereas the GRU, LSTM, and BiLSTM are utilized for low-frequency sequences. In the final selection of integration methods, one can choose between linear approaches, such as simple summation, or nonlinear methods like the GRU. Additionally, the decision to incorporate the GWO algorithm for model optimization is available. Thus, we conclude with 29 models and select the GARCH and GWO-GRU configurations. Further details regarding the comparisons are provided in Section 5.3.2, Section 5.3.3, Section 5.3.4, Section 5.3.5 and Section 5.4.

Establishing an appropriate GARCH model needs to test the ARCH effect of the sequence, and it is found that the ARCH effect exists in IMF1 to IMF5. According to the Akaike Information Criterion (AIC) and Bayesian Information Criterion (BIC), the optimal number of autoregressive terms and moving average terms $(p,q)$ is selected within a range measuring no more than 3, as shown in

Table A2. Hyperparameter settings of the framework.

Model	Subsequence	p		q
GARCH	IMF1	2		2
	IMF2	2		1
	IMF3	2		1
	IMF4	1		1
	IMF5	2		2
	Subsequence	Learning rate	Batch size	Epochs
GRU	IMF6	0.008	116	161
	IMF7	0.006	124	146
	IMF8	0.006	96	265
	IMF9	0.012	78	235
	IMF10	0.005	126	366
	Method	Learning rate	Batch size	Epochs
	Nonlinear	0.006	64	188
	integration

In Grid Search, the scope of learning rate is $[0.001,0.005,0.01,0.02,0.03,0.04,0.05]$, the scope of batch size is $[16,32,64,128,256]$, the scope of epochs is $[100,200,300,400,500]$.

.

Since hyperparameters significantly affect the predictive accuracy of machine learning models [60], this paper focuses on optimizing iterations epochs, learning rate, and batch size. When forecasting different data, the optimal hyperparameters of a model are also different. Adjusting hyperparameters for each subsequence to achieve optimal prediction effects is necessary.

Although GWO has been widely applied and has achieved excellent results in machine learning models [17,61,62], GWO is prone to fall into local optimization [63]. Therefore, we combine GWO and Grid Search [64] to fully leverage their advantages and improve predictive effectiveness. Grid Search is used to search for the approximate range of hyperparameter combinations, and GWO is used to reduce the range further and find an optimal hyperparameter combination. The GWO parameters are set as follows: the population of grey wolves is 20 and the maximum number of iterations is 30. The hyperparameter settings of each subsequence in the proposed framework, CEEMDAN-GWO-GRU/GARCH-GRU, are shown in Table A2. For the remaining hyperparameters not listed in the table, all models are set to the same values. For example, the number of responses is set to 1, the number of hidden units to 100, and the number of features to 3. In addition, we use the Adam optimizer to update the weights and multiply the learning rate by 0.15 after 200 iterations automatically.

5.3.2. Forecasting Results and Evaluation

After model selection and parameter setting, we obtain the proposed framework and 29 models and frameworks, including CEEMDAN-GWO-GRU/GARCH-GRU (the proposed model, framework P for short), CEEMDAN-GWO-ELMAN/LSTM-GRU (CGELG), CEEMDAN-GARCH/GRU-GRU (CGGG), CEEMDAN-ELMAN/BILSTM (CEB), VMD-GARCH/LSTM-LSTM(VGLL), CEEMDAN-GARCH/GRU (CGG), CEEMDAN-GARCH/BILSTM (CGB), CEEMDAN-GRU (CG), CEEMDAN-SVR (CS), the GRU, LSTM, and ELMAN.

The prediction results of all the above models and frameworks on the test set are shown in Figure 7, where the red bold line and the blue bold line represent the prediction results of Framework P and the real value of the carbon price, respectively. Obviously, the prediction and trend of Framework P are closer to the real carbon price and it performs better than the other compared models. In addition, the two small figures show the local prediction results and residuals, respectively, further proving the effectiveness of the model.

To provide a more comprehensive assessment of the superiority of Framework P’s predictions, the results of Framework P and the comparison models under the three assessment criteria of RMSE, MAE, and MAPE are shown in

Table 1. Model and framework evaluation by RMSE, MAE, and MAPE.

Models	RMSE	Rank	MAE	Rank	MAPE	Rank
Framework P	0.2493	1	0.1895	1	0.0024	1
CEEMDAN-GARCH/GRU-GRU	0.3865	2	0.3026	2	0.0039	2
CEEMDAN-GWO-ELMAN/GRU-GRU	0.3997	3	0.3220	3	0.0042	3
CEEMDAN-GWO-GARCH/LSTM-GRU	0.4178	4	0.3423	4	0.0044	4
CEEMDAN-ELMAN/GRU-GRU	0.4564	5	0.3582	5	0.0046	5
CEEMDAN-GARCH/GRU	0.4978	6	0.3726	6	0.0048	6
CEEMDAN-GARCH/LSTM-GRU	0.5062	7	0.3780	7	0.0049	7
CEEMDAN-GWO-ELMAN/LSTM-GRU	0.5313	8	0.3968	8	0.0051	8
CEEMDAN-ELMAN/GRU	0.5335	9	0.4086	9	0.0053	9
VMD-GARCH/LSTM-LSTM	0.5644	10	0.4366	11	0.0056	11
CEEMDAN-ELMAN/LSTM-GRU	0.5749	11	0.4274	10	0.0055	10
CEEMDAN-GARCH/LSTM	0.6794	12	0.5461	12	0.0071	12
CEEMDAN-ELMAN/LSTM	0.6904	13	0.5565	13	0.0072	13
GARCH	0.9856	14	0.6755	14	0.0087	14
LSTM	0.9957	15	0.7508	15	0.0097	15
CEEMDAN-GWO-GARCH/BiLSTM-GRU	1.0061	16	0.8752	18	0.0113	18
CEEMDAN-GRU	1.0092	17	0.9003	20	0.0116	20
CEEMDAN-ELMAN	1.0147	18	0.8804	19	0.0113	19
CEEMDAN-LSTM	1.0253	19	0.9234	21	0.0120	21
CEEMDAN-GARCH/BiLSTM-GRU	1.0895	20	0.9376	22	0.0122	22
SVR	1.0986	21	0.7655	16	0.0099	16
GRU	1.1063	22	0.7869	17	0.0102	17
CEEMDAN-GWO-ELMAN/BiLSTM-GRU	1.1298	23	0.9461	23	0.0122	23
CEEMDAN-SVR	1.1323	24	1.0038	25	0.0130	25
ELMAN	1.2350	25	0.9778	24	0.0125	24
CEEMDAN-ELMAN/BiLSTM-GRU	1.2445	26	1.0936	26	0.0142	26
CEEMDAN-GARCH/BiLSTM	1.5552	27	1.3432	27	0.0175	27
CEEMDAN-ELMAN/BiLSTM	1.5787	28	1.3466	28	0.0176	28
CEEMDAN-BILSTM	1.9799	29	1.6186	29	0.0218	29
BILSTM	2.8130	30	2.3724	30	0.0305	30

. Framework P ranks first and performs best, with an RMSE, MAE, and MAPE of 0.2493, 0.1895, and 0.24%. Compared with the framework proposed by [56,65], Framework P can better combine the decomposition-integration method, econometrics, machine learning, and the nonlinear integration method, giving full play to their respective advantages to achieve a greater effect. Additionally, we also find that most models have the same ranking, but some models change under the three evaluation criteria. Among them, the single SVR and GRU models have the biggest change in ranking, both changing by five places. We find that all single forecasting models have a large change in ranking, indicating that the prediction results of single models are unreliable.

The comparison results for all models and frameworks based on the three evaluation metrics, RMSE, MAE, and MAPE, are presented in Figure 8. The bars corresponding to Framework P are consistently the shortest, demonstrating that this proposed framework achieves the highest accuracy in forecasting carbon prices with the lowest prediction error. Furthermore, it is observed that models utilizing the nonlinear integration method generally exhibit lower prediction errors compared to their counterparts, suggesting that this method effectively enhances predictive accuracy. Additionally, it is noted that models exhibiting larger errors tend to be simpler, highlighting the challenges single models face in accurately predicting carbon prices.

We establish 30 linear regression models with the actual carbon price as the dependent variable and predictive carbon prices as the independent variable. The results of linear regression parameters are shown in Figure 9. The R² of Framework P is 0.9860, close to 1, which indicates that the prediction result of Framework P is closer to the real value than the other 29 comparison models and frameworks, proving the accuracy and reliability of Framework P’s prediction.

Based on the three loss functions of the MSE, MAE, and MAPE, we conduct the Diebold–Mariano test on Framework P and the other comparable models and frameworks shown in

Table 2. Model and framework performance by RMSE, MAE, and MAPE.

Models	MSE	MAE	MAPE
CEEMDAN-GWO-ELMAN/GRU-GRU	−6.1771 ***	−7.0221 ***	−7.0538 ***
CEEMDAN-GWO-GARCH/LSTM-GRU	−7.5423 ***	−11.5173 ***	−11.5183 ***
CEEMDAN-GWO-ELMAN/LSTM-GRU	−6.5809 ***	−9.1506 ***	−9.1324 ***
CEEMDAN-GWO-GARCH/BILSTM-GRU	−13.4917 ***	−21.2599 ***	−21.3997 ***
CEEMDAN-GWO-ELMAN/BILSTM-GRU	−12.0551 ***	−18.0871 ***	−18.1672 ***
VMD-GARCH/LSTM-LSTM	−6.7997 ***	−9.2261 ***	−9.2359 ***
CEEMDAN-GARCH/GRU-GRU	−5.1761 ***	−5.7379 ***	−5.7630 ***
CEEMDAN-ELMAN/GRU-GRU	−7.3982 ***	−8.7560 ***	−8.7531 ***
CEEMDAN-GARCH/LSTM-GRU	−6.2647 ***	−8.2370 ***	−8.2471 ***
CEEMDAN-ELMAN/LSTM-GRU	−7.1920 ***	−10.2632 ***	−10.2495 ***
CEEMDAN-GARCH/BILSTM-GRU	−12.6294 ***	−17.0233 ***	−16.9984 ***
CEEMDAN-ELMAN/BILSTM-GRU	−14.7298 ***	−19.4080 ***	−19.3736 ***
CEEMDAN-GARCH/GRU	−6.0847 ***	−7.8530 ***	−7.8559 ***
CEEMDAN-ELMAN/GRU	−6.6917 ***	−9.1777 ***	−9.2004 ***
CEEMDAN-GARCH/LSTM	−8.9931 ***	−13.2971 ***	−13.3296 ***
CEEMDAN-ELMAN/LSTM	−8.9677 ***	−13.6142 ***	−13.6215 ***
CEEMDAN-GARCH/BILSTM	−14.1627 ***	−19.8066 ***	−19.5371 ***
CEEMDAN-ELMAN/BILSTM	−13.8912 ***	−19.1060 ***	−18.8426 ***
CEEMDAN-GRU	−13.5660 ***	−22.6894 ***	−22.7061 ***
CEEMDAN-LSTM	−14.2671 ***	−20.0551 ***	−20.0200 ***
CEEMDAN-ELMAN	−12.6400 ***	−18.9933 ***	−19.2058 ***
CEEMDAN-SVR	−14.2333 ***	−22.6436 ***	−22.6531 ***
CEEMDAN-BILSTM	−11.7816 ***	−17.3966 ***	−17.0744 ***
GRU	−6.2386 ***	−12.2136 ***	−12.1221 ***
LSTM	−6.7956 ***	−13.3869 ***	−13.2468 ***
ELMAN	−9.4926 ***	−16.0688 ***	−16.1014 ***
SVR	−5.9837 ***	−11.2275 ***	−11.1336 ***
GARCH	−5.9129 ***	−10.6305 ***	−10.5727 ***
BILSTM	−12.9678 ***	−20.3269 ***	−20.5668 ***

*** means that the null hypothesis is rejected at the 1% significance level.

. In 87 comparison examples with 29 comparison models and frameworks, Framework P rejects the null hypothesis at the significance level of 1%. The error of the loss function is statistically significant, indicating that Framework P is superior to the other comparable models and frameworks.

5.3.3. Comparison of Single Models and Decomposition-Integration Frameworks

We further analyze and compare Framework P with other models and frameworks through various comparative analyses to obtain deep conclusions. Since the models’ rankings under the three evaluation indexes of RMSE, MAE, and MAPE are the same, only the RMSE is selected as the evaluation index for comparison in this section. Figure 10 compares the prediction results of single models and decomposition-integration frameworks. The radar map shows that, although CEEMDAN-BILSTM has a great improvement compared with the single BILSTM model, the improvement effect of the other four models is not obvious, and even the effect of the LSTM and SVR models becomes worse after applying the decomposition-integration method.

To further explore the reasons, we draw partial prediction results in Figure 10. We find that BILSTM cannot learn the characteristics of carbon prices well, with the predictive result randomly fluctuating, and other single models have obvious lag. The single framework prediction results in Figure 10 indicate that the five models, except for the BILSTM model, only take last-day carbon prices as the predictive values and cannot learn the fundamental characteristics of carbon prices. Although their results are good in terms of the evaluation indicators, the five models have no practical significance. Accordingly, we believe single models are ineffective in carbon price forecasting.

In contrast, the decomposition-integration framework prediction results in Figure 10 show that decomposition-integration frameworks overcome the lag problem of single models and can predict carbon prices more accurately. Compared with the BILSTM model, the CEEMDAN-BILSTM model can learn carbon price characteristics more effectively, and the forecasting result fluctuates less and is more accurate. According to the above analysis, it is proved that decomposition-integration frameworks can predict carbon prices more accurately than single models, which is consistent with the views of [4,40,43]. Therefore, we can obtain a series of simple and stable subsequences with features easily learned through carbon price decomposition to improve predictive accuracy.

5.3.4. Comparison of Single and Hybrid Methods in Decomposition-Integration Frameworks

Figure 11 presents the Root Mean Square Error (RMSE) for single methods and hybrid methods within the decomposition-integration framework. The red squares represent single models, while the blue circles represent hybrid models. Clearly, the hybrid models are more accurate than the single models, with an average RMSE reduction of 35%, which significantly enhances predictive performance. Compared with the studies of [17,65], our study reveals the advantages of using hybrid methods for decomposition-integration frameworks. Hybrid methods apply more appropriate models to subsequences based on data characteristics, while single methods blindly use models for all predictions. For example, the GRU is an excellent forecasting model for low-frequency sequences but is not ideal for high-frequency sequences. GARCH can capture the complex randomness, time variability, and aggregation of high-frequency sequences. Combining GARCH and the GRU can give full play to their strengths and obtain more accurate forecasting results than using only the GRU to predict all subsequences.

5.3.5. Model and Method Selections in Decomposition-Integration Frameworks

We further analyze the hybrid methods within the decomposition-integration framework. As shown in Figure 12, the framework using GARCH tends to be more accurate than the one using ELMAN, particularly in the nonlinear integration framework. The RMSE of the nonlinear integration GARCH framework is nearly 0.07 lower than that of the ELMAN framework, while the RMSE of the linear integration GARCH framework is reduced by less than 0.04. This finding suggests that GARCH is better suited for predicting high-frequency series, consistent with the earlier conclusion that the GARCH model can effectively capture the characteristics of high-frequency series.

As shown in Figure 13, the prediction performance of the decomposition-integration framework based on the GRU outperforms the LSTM and BILSTM frameworks, indicating that the GRU is more suitable for predicting low-frequency series. Additionally, the prediction performance of the BILSTM framework is significantly lower than that of the GRU and LSTM, demonstrating that BILSTM is not suitable for forecasting low-frequency carbon price series.

In Figure 14, a comparison is made between the nonlinear integration method and the linear integration method. The findings reveal that the RMSE for the nonlinear integration framework is, on average, 23.1% lower than that of the linear integration framework. When contrasted with the straightforward linear integration methods applied by [65,66,67], our results highlight the enhanced effectiveness and accuracy of the nonlinear integration approach. The main reason for this difference is that simple linear summation of the subsequences aggregates the errors of each subsequence. The nonlinear integration method, on the other hand, learns the relationship between the simple linear sum and the real values, and thus is able to correct the errors and improve the prediction accuracy.

We also investigated the influence of intelligent optimization algorithms on prediction outcomes. Our results indicate that, compared to the findings of [4,56], these algorithms can significantly enhance prediction accuracy. As shown in Figure 15, the framework using the GWO algorithm performs better than the framework without optimization algorithms. Intelligent optimization algorithms maximize the performance of the framework by finding the optimal combination of hyperparameters.

5.4. Multi-Step-Ahead Forecasting

The previous discussion mainly focuses on out-of-sample one-step-ahead forecasting. However, in the context of the actual operation of the carbon trading market, multi-step-ahead forecasting of carbon prices is equally important. Since multi-step forecasting cannot incorporate the latest known information immediately, it often faces the issue of error accumulation as the forecast horizon extends, a problem that is particularly evident in complex decomposition-integration forecasting frameworks. Therefore, this section primarily compares and analyzes the effectiveness and robustness of the decomposition-integration framework in multi-step-ahead forecasting.

For simplicity, this section focuses on the nine models that performed best in the one-step forecasting scenario as a comparison. Multi-step-ahead forecasts of the CEEX Guangzhou carbon price are conducted with forecast horizons of 3, 5, and 7 steps. Figure 16 shows the RMSE of the models at different forecast horizons. As the forecast horizon increases, the difficulty of capturing the time dependency of the series also increases, resulting in a general decline in the predictive accuracy of all models. As shown in Figure 16, the RMSE of Framework P, proposed in this paper, is significantly lower than that of other comparison models across all forecast horizons, indicating that Framework P outperforms the other models and significantly improves predictive accuracy. Compared to the next-best model, Framework P achieves an average RMSE improvement of 23.2%. Additionally, as the forecast horizon increases, the error curves of models using linear integration methods (such as CEEMDAN-GARCH/LSTM, CEEMDAN-ELMAN/GRU, and CEEMDAN-GARCH/GRU) are steeper compared to those using nonlinear integration methods. This phenomenon suggests that the disadvantage of linear integration methods in multi-step-ahead forecasting is further amplified due to the error accumulation effect. In contrast, nonlinear integration methods, by learning the differences between predicted and actual values, demonstrate greater stability in multi-step-ahead forecasting.

In summary, compared to other models, the proposed Framework P has a lower and more stable error curve, further demonstrating its robustness in carbon price forecasting.

5.5. Robustness Test

5.5.1. Robustness Test of Hyperparameters

The framework proposed in this paper and most comparison frameworks use neural networks whose performance is easily affected by hyperparameters. If a model is less affected by hyperparameter changes, it has strong robustness. Since the learning rate, iteration epochs, and batch size have been optimized in Section 5.5, this section only adjusts and researches the number of hidden units significantly impacting the models’ performance. We test the parameter values of 10, 50, 100, 200, and 300 and record the RMSE of each experiment, as shown in Figure 17.

The box diagram for BILSTM and its derivative frameworks is much longer and higher than that for other comparable models, again indicating that BILSTM frameworks have a poor effect on carbon price forecasting. Although the minimum RMSE of CEEMDAN-ELMAN is lower than that of the single ELMAN model, the maximum RMSE is significantly higher. On the one hand, since the single ELMAN model only learns the simple features of the last-day carbon price, it has no practical prediction significance even though the RMSE is low. On the other hand, ELMAN is good at extracting data features of high-frequency sequences, while its prediction effect on low-frequency sequences is not stable [4,68]. With changes in hyperparameters, this defect is further amplified, making the RMSE of CEEMDAN-ELMAN high overall and it has a large degree of dispersion because low-frequency sequences account for a large proportion of the original sequence.

We exclude BILSTM, its derivative frameworks, and CEEMDAN-ELMAN to further analyze the hyperparameter sensitivity of the other models and frameworks shown in Figure 18. Obviously, the fluctuation range of LSTM and its derivative frameworks is relatively large, which means that they are more sensitive to the change of hyperparameters. In contrast, the GRU and its derivative frameworks are more stable with the change, and their median and average are significantly lower than that of LSTM, proving that the GRU is more suitable for predicting carbon prices. Additionally, although GARCH and ELMAN have similar box lengths, the overall performance of GARCH frameworks is significantly lower, indicating the effectiveness of GARCH. The box of Framework P, proposed in this paper, is significantly lower than other comparable models and frameworks. The RMSE value is relatively stable in response to changes in hyperparameters, indicating that Framework P has a relatively good performance in carbon price forecasting.

5.5.2. Robustness Test of the Sample Ratio

This paper further analyzes the robustness of Framework P under different training sets and test set lengths. Instead of 10%, we select the first 80% of data as the training set, with 1575 samples from 9 December 2013 to 22 July 2021, and the last 20% of data as the test set, with 394 samples from 22 July 2022 to 3 March 2023. As BILSTM and its derivative frameworks are unsuitable for carbon price forecasting, we exclude them and obtain results for the other models and frameworks shown in

Table 3. Linear regression analysis of the carbon price.

Models	RMSE	MAE	MAPE	MSE	MAE	MAPE
Framework P	0.7200	0.5195	0.0077
CEEMDAN-GARCH/GRU-GRU	0.9346	0.7879	0.0134	−4.7897 ***	−9.9407 ***	−10.7727 ***
CEEMDAN-GWO-ELMAN/GRU-GRU	0.9504	0.7358	0.0118	−3.7273 ***	−6.3850 ***	−7.8867 ***
VMD-GARCH/LSTM-LSTM	1.0528	0.8529	0.0132	−5.1316 ***	−5.5379 ***	−5.5628 ***
CEEMDAN-ELMAN/GRU-GRU	1.1651	0.9869	0.0148	−7.4339 ***	−13.9543 ***	−14.9470 ***
GARCH	1.4256	0.8598	0.0130	−4.4203 ***	−7.0562 ***	−7.8696 ***
CEEMDAN-GWO-GARCH/LSTM-GRU	1.6283	1.3237	0.0185	−14.4173 ***	−17.4572 ***	−17.8334 ***
CEEMDAN-GWO-ELMAN/LSTM-GRU	1.7670	1.4217	0.0195	−13.1729 ***	−17.4485 ***	−17.3448 ***
CEEMDAN-GARCH/LSTM-GRU	1.9300	1.6418	0.0231	−17.1490 ***	−22.2719 ***	−23.6570 ***
GRU	1.9390	1.4422	0.0205	−7.1105 ***	−16.0831 ***	−16.9759 ***
CEEMDAN-GARCH/GRU	1.9473	1.7061	0.0238	−18.3099 ***	−26.0535 ***	−27.5366 ***
CEEMDAN-ELMAN/GRU	1.9864	1.7045	0.0237	−19.2076 ***	−26.4907 ***	−28.2950 ***
CEEMDAN-ELMAN/LSTM-GRU	1.9908	1.6543	0.0231	−15.8921 ***	−20.1434 ***	−20.9796 ***
CEEMDAN-GRU	2.0528	1.7982	0.0252	−17.6675 ***	−26.6100 ***	−28.5906 ***
SVR	2.5305	2.0495	0.0282	−9.9736 ***	−22.2362 ***	−23.9457 ***
CEEMDAN-GARCH/LSTM	2.8715	2.4311	0.0327	−21.3015 ***	−25.4924 ***	−25.8639 ***
CEEMDAN-ELMAN/LSTM	2.9090	2.4694	0.0333	−21.5075 ***	−25.7659 ***	−26.3088 ***
CEEMDAN-LSTM	2.9189	2.4706	0.0334	−21.3924 ***	−25.3399 ***	−26.0190 ***
CEEMDAN-ELMAN	3.1121	2.6654	0.0366	−20.1226 ***	−27.6298 ***	−29.3190 ***
LSTM	3.2580	2.6762	0.0365	−12.9878 ***	−24.1499 ***	−26.2976 ***
CEEMDAN-SVR	3.2695	2.7476	0.0373	−19.7956 ***	−25.7724 ***	−26.9241 ***
ELMAN	3.3062	2.6721	0.0366	−10.4115 ***	−23.7687 ***	−25.5273 ***

*** means that the null hypothesis is rejected at the 1% significance level.

. The RMSE, MAE, and MAPE of Framework P are, respectively, 0.7200, 0.5195, and 0.77%, ranking first among the 21 models and frameworks. The Diebold–Mariano test results show that Framework P rejects the null hypothesis at the significance level of 1% in 63 comparison examples, indicating that the loss is statistically significant. Therefore, the prediction performance of Framework P under different sample lengths is significantly better than that of other comparable models and frameworks, which strongly proves the validity and robustness of the proposed Framework P in carbon price forecasting.

6. Conclusions

Accurate carbon price prediction is crucial for businesses to mitigate operational risks and enhance decision-making. For policymakers, it plays a vital role in maintaining market stability and promoting environmental sustainability. However, predicting carbon prices remains a significant challenge due to their non-stationary behavior, nonlinear complexities, and high volatility. This study introduces a novel decomposition-integration framework, CEEMDAN-GWO-GRU/GARCH-GRU, for carbon price forecasting in China, leveraging the strengths of both econometric models and machine learning techniques from a multi-scale perspective.

Initially, carbon trading prices are decomposed using the CEEMDAN model, and the distinct characteristics of modal components are identified through the Lempel–Ziv complexity algorithm and dispersion entropy. Subsequently, a divide-and-conquer approach is applied, where the GARCH model and the GWO-GRU model are used for high-frequency and low-frequency components, respectively, ensuring tailored predictions for each. This method circumvents the potential pitfalls of applying the same model uniformly to all modal components. In the final step, carbon price forecasts are obtained by nonlinearly integrating the predictions of different components using the GWO-GRU model, yielding results that outperform single and traditional hybrid models. The effectiveness of this new framework is demonstrated using carbon price data from the Guangzhou Carbon Exchange. Through statistical analyses, including the RMSE, MAE, MAPE, linear regression, DM tests, and sample ratio robustness tests, the study confirms the improved accuracy and robustness of the proposed forecasting approach for non-stationary and nonlinear carbon price series.

The main conclusions of this paper are as follows:

(1) The decomposition-integration method can effectively solve the lag problem of single models and significantly improve the prediction accuracy and tendency.

(2) Hybrid methods in decomposition-integration frameworks are better than single methods, in which the GRU and GARCH perform better than LSTM, BILSTM, and ELMAN, and are not susceptible to hyperparameters.

(3) The nonlinear integration approach has the capability to understand the relationship between predicted and actual values, allowing it to rectify prediction errors. By using the GWO algorithm, we can identify the optimal hyperparameter combinations, thereby enhancing the accuracy of the predictions.

(4) We compare CEEMDAN-GWO-GRU/GARCH-GRU with 29 models and frameworks, proving its superiority and effectiveness, which may play a role in reducing risks, stabilizing the market, and promoting economic development.

Despite these successes, the study has some limitations. The current model relies solely on historical carbon prices, without incorporating external factors such as policy changes or economic conditions, even though considering these external variables could further improve the forecasting performance. Additionally, the proposed model could be adapted for other time series forecasting applications, including those in fields such as environmental science, finance, and energy, broadening its relevance to various complex scenarios.