Forecasting China's energy intensity by using an improved DVCGM (1, N) model considering the hysteresis effect

Energy consumption is considered as an important driver of economic growth. Energy intensity (EI), which is the ratio of total energy consumption to gross domestic product (GDP) of a country or region in a certain period of time, measures the performance of energy utilization. According to the energy section of the National Bureau of Statistics (NBS), China's largest energy consumption is 4,640 million tons of standard coal, and China has become the world largest energy consumer. Entering a new normal period, the government focuses on a crucial rebalancing and diversifying economy, with higher requirements for sustainable development. In the 13th Five-Year Plan for economic and social development, a total reduction of 15% has been set as the energy performance target. It is imperative to estimate and predict the EI to evaluate energy conservation and emission reduction. However, due to the complexity and dynamics of social economy, the implementation of energy conservation and emission reduction policies will not immediately reduce EI, and there is often a certain time lag. This inevitable time-lag should be taken into account to estimate and forecast the EI in China accurately.

Numerous research studies have studied the influencing factors and driving forces of EI through econometrics methods such as cointegration analysis, metrology and decomposition analysis and scenario analysis. Zhu et al. (2015) applied cointegration analysis method on the large and state-owned enterprises and found that energy-saving regulations in China are one of the most important factors in reducing aggregate EI. Karimu et al. (2017) studied the EI and convergence of Swedish industry by combining metrology and decomposition analytical methods. Ma and Yu (2017) used panel data model to discuss the driving factors which lead to EI decline. It is revealed that industrial structure, energy conservation regulations and EI are closely related. Tan et al. (2018) used index decomposition analysis and production decomposition analysis methods to analyze the factors which are related to the decline in the EI and pointed out that technology improvement effect is the most significant factor. In EI forecasting, Pao et al. (2012) used improved grey models (GMs) to predict China's CO₂ emissions, energy consumption and economic growth. Dong et al. (2018) estimated the driving force of regional EI in China and forecasted the potential of regional energy conservation with scenario analysis. Wu et al. (2018) used a new multivariable GM to predict energy consumption in Shandong Province.

The grey system theory is an interdisciplinary theory proposed by Deng (1982) and the grey prediction method applies well to small size data forecasting. As an important part of grey prediction theory, the GM (1, N) model is the basic model of multivariable grey system modeling approach. In recent years, numerous scholars have thoroughly discussed the model parameters optimization (Tien, 2005, 2010), the model accuracy improvement (Tien, 2011; Wang et al., 2016) and the expansion of the GM (Guo et al., 2013; Kose and Tasci, 2015, 2019; Ding et al., 2017).

Based on grey multivariable model with time-lagged system, Zhai et al. (1996) introduced the lag term into the GM (1, 2) model and determined the delay parameters with the goal of minimizing the modeling error. Hao (2011) used grey correlation analysis to determine time-lag period between variables and then on this basis to establish forecast model of GM (1, N). Zhang et al. (2015) constructed a time-delay multivariable discrete GM, DDGM (1, N) model, by introducing a time-delay control factor and solved the time-delay parameters by using the grey dimensionally expanding identification method, which obtained a good application effect. Ma and Yu (2017) used a novel time-delay multivariable GM to predict the natural gas consumption in China. Dang et al. (2017) constructed the discrete delay grey multivariable DDGMD (1, N) model by introducing the driving information control adjustment coefficient T and the action coefficient λ and solved the coefficient respectively by grey dimension expansion method and particle swarm optimization (PSO) algorithm. Xiong (2019) built a multivariable time-delay discrete MGM (1, m, t) model and studied the mechanism of modeling and the process of modeling, and the calculation method of time delay is given. The hysteresis effect is discussed by an example to verify the validity of the model.

GM (1, N) model, in spite of successfully applying in various fields, sometimes ignores the influence of virtual variables such as policy on the main system in practical applications. Zhang (2016) considered the influence of dummy variables on system behavior variables and built a discrete multivariable prediction model based on dummy variables, which further expanded the application scope of the model. Ding et al. (2018) introduced the dummy variables into the GM (1, N) model, gave the concrete model construction method in mathematics and verified the effectiveness of the new model with cases.

Generally speaking, we have abundant literature discussing the influencing factors of EI and its forecasting. Furthermore, there are some practices targeting the time-lag phenomenon using grey theory. The existing literature has explored and studied the GM from multiple angles, but there is still room for the GM to expand in the combination of dummy variables and time-delay systems. As the common GM (1, N) model does not take into account the time delay between variables, and there are dummy variables in the system which are difficult to be measured by quantity, our optimized DVCGM (1, N) model, which is short for dummy variables control grey model of N variables, incorporates the above two features and improves prediction accuracy to a satisfying level. It is practical significant because the time delay of dummy variables, i.e. government policy, is often seen in real world but is difficult to measure. The optimized DVCGM (1, N) model takes these variables into consideration, solves the practical problem and expands the grey system theory and the grey prediction method system, which improves the accuracy of grey prediction model.

This paper's main interest is to estimate and forecast EI by considering the influence of government policies and to test the accuracy of the improved GM through a comparison study. The main contribution of the study to the literature is to consider hysteresis effect and increase the forecast accuracy. It is imperative for optimizing energy structure, improving energy utilization efficiency and ensuring energy security. The results show that the improved GM produces better results than the other three conventional models. The rest of the paper is organized as follows. Section 2 briefly introduces grey theory and multivariable grey prediction models. A nonlinear optimization method based on PSO algorithm is constructed to calculate the hysteresis parameter in the improved model. In Section 3, EI is estimated by considering the hysteresis effect of energy-saving policies by an improved DVCGM (1, N) model. We compare the results with the other two GMs and one econometric model. The improved DVCGM (1, N) model has the best performance in the estimation comparison and is applied to forecast future EI. Section 4 is the conclusion of the study with the limitation and future path.

The three grey prediction models used in this paper are interrelated and in a progressive order. The GM (1, N) model is a traditional multivariable grey prediction model. The DVCGM (1, N) model introduces virtual variables, taking into account the influence of policy and other factors on the basis of the GM (1, N) model. Time-delay parameter is introduced in the improved DVCGM (1, N) model, considering the hysteresis effects of historical variables, which further enriches the existing grey prediction theory.

Grey prediction model can be regarded as two levels of work. At a lower level, the original sequence produces the sequence of generation by one accumulative generation (1-AGO), and then it forms the sequence of mean generation of consecutive neighbors; similarly, the sequence of influencing factors generates the sequence of generation by 1-AGO; Constructing B and Y matrix and calculating system parameters through ordinary least squares (OLS) regression. Once the system parameters a and b are determined, we can obtain the time response function (TRF) by solving the whitenization equation. At a higher level, the continuous differential equation with initial values is used as the reflection equation, and the discrete data are mapped to a manageable function, which is further restored to the TRF as the simulation basis. Typical procedures can be described briefly by the program in Figure 1.

$X1(1)$ is the 1 - AGO sequence from $X1(0)$⁠, and $Xi(1)$ is the sequences of the relevant factors, where

$Z1(1)$ is the sequence of mean generation of consecutive neighbors from $X1(1)$⁠, where

where k = 2, 3, n, a is the development coefficient and $b1,b2,⋯,bn$ are the grey input coefficients obtained by the least squares method. To determine these coefficients, the matrix B and $YN$ are defined as follows:

The values of the coefficients a and $b1,b2,⋯,bn$ can be determined by the following equation:

After determining the coefficients of a and $b1,b2,⋯,bn$⁠, the differential equation of the GM can be determined by Eq (2). The solution of the above differential equation is as follows:

where $Xˆ1(1)(k)$ is the prediction of the AGO of the original sequence. By considering that the estimation of the first element of the first AGO of a sequence is equal to the first element of the sequence, the following relation is determined: $Xˆ1(1)(1)=X1(0)(1)$

Finally, in order to predict the elements of the original sequence, the inverse accumulated generating operation should be performed. Therefore, the predicted values can be determined as follows:

where $xˆ1(0)(n)$ is an estimation of the original sequence, which is simulation values, $xˆ1(0)(n+1),xˆ1(0)(n+2),⋯$ are predictive values.

Traditional GM (1, N) model ignores the influence of virtual variables, which will inevitably lead to significant errors in practical applications. Therefore, it is necessary to construct a new multivariable predictive model with virtual variable control, based on the traditional GM (1, N) model, i.e. the DVCGM (1, N) model. The modeling steps for DVCGM (1, N) can be illustrated in Figure 2.

Virtual variable sequence is

$Xi(1)$ and $Dj(1)$ are the 1-AGO sequence, $X1(0)$ is the behavior sequence of the system, $Xi(1)(i=2,⋯,M)$ and $Dj(1)(j=M+1,⋯,N)$ is the driving factor sequence. Then GM (1, N) model with dummy variable can be expressed as:

named DVCGM (1, N) model (dummy variables control grey model of N variables). $Z1(1)$ is the sequence of mean generation of consecutive neighbors from $X1(1)$ where,$Z1(1)(k)=0.5X1(1)(k−1)+0.5X1(1)(k),(k=2,3,…,n).$

where $∑i=2Mbixi(1)(k)$ is independent quantization variable driver, $∑j=M+1Nbjdj(1)(k)$ is the dummy variable driver. If virtual variables are not considered or b_j=0, the model can be transformed to a traditional GM (1, N) model.

as mentioned in Definition 4. The parameter column of the model is

The matrix B and $YN$ are defined as follows:

As mentioned in Definition 4, the parameter column of the model is

Let Eq (5) be defined as the differential equation of DVCGM (1, N) model:

After determining the coefficients of a and $b1,b2,⋯,bn$⁠, the differential equation of the GM can be determined by Eq (5). The solution of the above differential equation is as follows:

Finally, in order to predict the elements of the original sequence, the inverse accumulated generating operation should be performed. Therefore, the predicted values can be determined as follows:

where $xˆ1(0)(n)$ is an estimation of the original sequence, which is simulation values, $xˆ1(0)(n+1),xˆ1(0)(n+2),⋯$ are predictive values.

The classic multivariable grey prediction models, such as traditional GM (1, N) and DVCGM (1, N) models, can reflect the influences of current driving-variables on the present system behavior and innately ignore the hysteresis effect of historical variables. Therefore, an improved DVCGM (1, N) model is proposed, integrating these above prediction model.

In this section, the hysteresis parameter $λi$ is innovatively introduced into the DVCGM (1, N) model to improve the prediction accuracy. Supported by PSO algorithm, detailed process and algorithm can be described as following:

As mentioned in Definition 4. The parameter column of the model is

The matrix B and $YN$ are defined as follows:

The least square estimation of the parameter column satisfies the following requirements:

1. When n = N+1, $bˆ=B−1Y,|B|≠0;$

2. When n > N+1, $bˆ=(BTB)−1BTY,|BTB|≠0;$

3. When n < N+1, $bˆ=BT(BTB)−1Y,|BTB|≠0;$

Proof: Substitute k = 2, 3,…, n into the model, you can get

That is, by the least square method, $Y=Bbˆ$

1. When n = N+1, B has an inverse matrix, the equations have a unique solution, we can get $bˆ=B−1Y,|B|≠0;$

2. When n > N+1, B is column full rank, the full rank decomposition of B is B = DC. Then the generalized inverse matrix of B can be obtained:

Because B is a full rank matrix, C can be taken as a unit matrix, B = D, so

• (3)When n < N+1, B is a row full rank matrix, D can be taken as a unit matrix, B = C, so

Where k = 2, 3, n, a is the development coefficient and $b1,b2,⋯,bn$ are the grey input coefficients, $λi$ is the hysteresis parameter. PSO algorithm is used to determine the coefficients of a, $b1,b2,⋯,bn$⁠, and $λi$⁠. Then the differential equation of the GM can be determined by Eq (7). The solution of Eq (7) is as follows:

When the range of the driving factor sequence is small, the driver term can be viewed as a grey constant, and then the approximate TRF sequence of the grey differential equation of the model is

Finally, in order to predict the elements of the original sequence, the inverse accumulated generating operation should be performed. The predicted values can be determined as follows:

where $xˆ1(0)(n)$ is an estimation of the original sequence, which is simulation values, $xˆ1(0)(n+1),xˆ1(0)(n+2),⋯$ are predictive values.

The most important part of the improved DVCGM (1, N) model is estimating the time-lag parameter, which directly affects the accuracy of the model. However, the time-lag parameters must be determined in advance, followed by B and Y matrix construction and system parameters calculation through OLS. Once the system parameters a and b are determined, we can obtain the TRF and the simulation and prediction value of the model.

In this paper, a nonlinear optimization model is established by using the Least One Multiplication. Then, the time-lag parameter is determined. When the range of the driving factor sequence is small, Eq (9) is used as the TRF, $λi$ can be solved by the following nonlinear programming model:

The model takes the relationship between structural parameters as the constraint condition and minimizes the average simulation relative error of the system characteristic variables, which can improve accuracy to the greatest extent.

The above optimization problem can be solved by PSO (Kiran and Mustafa, 2017; Mason et al., 2018; Chen et al., 2018). According to Eq (10), a nonlinear optimization method based on PSO algorithm is constructed to obtain the hysteresis parameter. PSO sets a certain number of particles in feasible region to find the best location and can be used to seek optimal values of $λi$⁠. Denote $λi$ in Eq (10) and construct the fitness function of each particle, according to Eq (12).

Obviously, the average simulation relative error of the system characteristic variable sequence varies depending on the number of lag periods. The value of lag period should be selected to make the average simulation relative error as small as possible. Therefore, the improved DVCGM (1, N) model can well describe the hysteresis effect between the system characteristic variables. Once the hysteresis parameter is determined, the structural parameters of the model are set accordingly, and the simulation and prediction results can be obtained according to Eq (9).

Detailed procedure of the improved DVCGM (1, N) model is illustrated as follows.

• Step 1. Collect raw data and establish original sequence $X1(0)$⁠, $X1(1)$ is the 1-AGO sequence from $X1(0)$⁠, $Xi(1)$ is the sequence of the relevant factors. $Z1(1)$ is the sequence of mean generation of consecutive neighbors from $X1(1)$⁠. Then, determining virtual variable sequence

$Dj(1)$ is also the sequence of the relevant factors.

• Step 2. Solving delay parameters $λi$ by PSO according (10), then constructing vector Y and matrix B. Using the least square method, the values of the coefficients a and $b1,b2,⋯,bn$ can be determined.

• Step 3. After considering the lag effect, the differential equation of the grey model can be determined by Eq (7).

• Step 4. The time response function of the new model is established to generate prediction data according to Eq (9).

Forecasting EI can be considered as a grey problem, because EI is greatly affected by technological progress, population factor, industrial structure and so on. These factors influence EI through a dynamic and complicated mechanism. The uncertain impacts and limited number of data provides a good basis for grey theory application. There are four procedures in this part, including data collection, parameter estimation, result comparisons and future forecasts. Three competing models, namely GM (1, N), DVCGM (1, N) and ARIMA model, are employed to test the accuracy of the improved DVCGM (1, N) on EI forecasting.

The three GMs used in this paper are interrelated and are from basic to the advanced. The GM (1, N) model is a traditional grey multivariate prediction model. On its basis, DVCGM (1, N) model introduces dummy variables, taking into account the influence of policies and other factors. The improved DVCGM (1, N) model introduces time-lag parameters, which further enriches the existing grey prediction theory by considering the hysteresis effect of policies. The above GMs are used to estimate China's EI, and the optimal model is selected by comparing their performance, to predict China's EI in the next five years. In addition, we use ARIMA model as a comparison of grey methods to show that the improved DVCGM (1, N) model is not only better than the traditional GM, but also better than the non-grey econometric model.

The indicators selected in this paper are as follows. EI is the ratio of total energy consumption to GDP. Population factor is the employed population. The ratio of the added value of industrial production to the added value of energy consumed by the industry is used as the substitution variable of technological progress (Yan, 2011). Industrial structure is measured by the ratio of output value of each industry to GDP. For the consistency of the statistical scope, we choose a time scale of 2001–2017, and all data are collected from China Statistical Yearbook. With small sample size (17 periods' real measurement values) and insufficient information, this case fits well with the grey system.

First, we calculate the grey correlation between EI and population, technological progress and industrial structure. As shown in Table 1, technological progress has the highest correlation with EI. Therefore, we select technological progress as the driving variable, and EI is the system behavior variable. Then, we can build GM (1, N) model for EI forecast.

As shown in Figure 3, China's EI increased rapidly from 2001, rising sharply to a peak of 1.4542 tons standard coal per 10,000 Yuan in 2005 and declined steadily afterward. In 2017, the EI has decreased by 62.67%, leaving only 0.5428 tons of standard coal per 10,000 Yuan. The 11th Five-Year Plan in 2006 is a kind of watershed for Chinese EI. The government introduced a series of strict policies of energy conservation and emission reduction, and China's EI declined obviously. As a result, we need to put policy and the hysteresis effect of policy into consideration when estimate EI of China.

Four models are applied in the simulation of EI of China. First, we use the GM (1, N) model mentioned in Section 2.1 to predict EI. We select technological progress as the driving variable, EI as the system behavior variable and then establish the GM (1, N) model. According to 2.1, the TRF is obtained as follows:

The simulation results of GM (1, N) model are illustrated in Table 3, and the relative errors and average relative errors can be calculated, as shown in Table 3.

Secondly, we build DVCGM (1, N) model for EI estimate. We select technological progress as the driving variable, EI as the system behavior variable. Energy conservation policy (P) is introduced as a dummy variable. Before 2006, as the strict energy-saving policies had not been implemented, P value is 0; after 2006, P value is 1. According to Section 2.2, the TRF is obtained as follows:

The simulation results of DVCGM (1, N) model are illustrated in Table 3, and the relative errors and average relative errors can be calculated, as shown in Table 3.

Thirdly, we use the improved DVCGM (1, N) model to estimate EI. As discussed in Section 2.3, DVCGM (1, N) model takes consideration of the hysteresis effect of dummy variables. We select technological progress as the driving variable, EI as the system behavior variable and energy-saving policy as virtual variable. Different from DVCGM (1, N) model, we need to determine the delay parameters, along with the structure parameters and build the TRF to get the simulation results.

1. Determination of time-delay parameter.

As stated by Eq (10-12), a nonlinear optimization method based on PSO algorithm is constructed to obtain the hysteresis parameter. According to the optimization model shown in Eq (9), the average relative percentage errors (APE) of the model under different lag periods is calculated in Table 2. When the hysteresis parameter is 2, the average relative percentage errors (APE) of the model is the smallest (0.9754%). Therefore, the value of time lag is two years. In the view of the complexity of socio-economic ecology, adaptive adjustments are made according to the changes of policies. These energy-saving policies influence economic operation through a certain transmission mechanism and gradually affect the EI.

• (2)Determination of structural parameters.

The hysteresis parameter is substituted into the B matrix and Y matrix. According to Theorem 2, the matrix operation of least square regression is used to obtain the structural parameter values b₁, b₂. B^T= (1.251669, 0.938994).

• (3)Calculate the simulated and predicted values.

Putting the values of the estimated structural parameter and hysteresis parameter in Eq (7-9), we can obtain the optimal TRF as follows:

Then, the simulation results of improved DVCGM (1, N) model are illustrated in Table 3, and the relative percentage errors (PE) and average relative percentage errors (APE) can be calculated, as shown in Table 3.

Finally, we employ ARIMA (autoregressive composite moving average) model to test our results from the non-grey perspective. The estimation results of ARIMA model are shown in Table A1-4 in the  Appendix. All the coefficients are statistically significant and the model is well fitted (R-squared is 0.850452). The ARIMA (2,1,1) model of time series is determined as follows:

The performance of the ARIMA model and the calculated relative errors and average relative errors are presented in Table 3.

As illustrated in Table 3, there is great consistency between simulated values and real values for the improved DVCGM (1, N) model. For APE, which is the performance prediction index, its values of the improved DVCGM (1, N) is the smallest (1.76% in the in-sample periods and 1.92% in the out-sample periods) among all four models.

As shown in Figure 4a, the overall trends of the GM (1, N), DVCGM (1, N) and improved DVCGM (1, N) model and ARIMA model comply with the true curve to some extent. However, the performance evaluation is distinct (Figure 4b). The APE value of the improved DVCGM (1, N) model is the smallest with the minimum fluctuation, demonstrating the efficacy and reliability of the model. It also suggests that the improved DVCGM (1, N) model is better than the traditional grey multivariable models. The ARIMA model has the largest values of APE among these four competing models. GM (1, N) model obtains the second largest values of APE, suggesting DVCGM (1, N) is the suboptimal choice.

In view of the hysteresis effect of energy-saving policies, the improved DVCGM (1, N) model has a much lower error than the other three models. Therefore, we choose the improved DVCGM (1, N) model as the best model for predicting EI.

Based on the improved DVCGM (1, N) model, we can forecast the output value of China's EI from 2018 to 2022 using data of technological progress. However, the data of technological progress from 2018 to 2022 is unknown and needs to be predicted in advance. The one-step rolling GM (1, 1) model is considered for the data prediction.

We used the GM (1, 1) model for in-sample simulation and out-of-sample prediction of technological progress. The results (Table 4) show that the in-sample simulation error is 2.97%, and the out-of-sample hind cast error is 2.11%, which indicates that the GM (1, 1) model gives satisfying results and can be used to predict technological progress from 2018 to 2022.

It is also noted that energy conservation policy (P) is introduced as a dummy variable. Before 2020, P value is 0; after 2020, P value is 1. It is because that the 13th Five-Year Plan is from 2016 to 2020, new performance target and new policies for the next stage are formulating and will be effective after 2020. Therefore, the data we need to substitute into the improved DVCGM (1, N) model is shown in Table 5.

The predicted results of China's EI are shown in Table 6. In addition, we also draw a line graph to make the results more iconic in Figure 5.

As illustrated in Figure 5, there is a downward trend of the EI in the next five years. By 2020, the EI is expected to decrease by 20% or more than it was in 2016. That is to say, during the 13th Five-Year Plan period (2016–2020), EI will drop by more than 15%, meeting the country's energy performance target. Therefore, government policies have a profound influence on EI. When formulating energy conservation and emission reduction policies, we should consider the hysteresis effect of the policies and make adjustments accordingly to achieve the goal.

Over the past 20 years, China is gradually shifting from a resource-intensive and energy-driven economy to a more sustained economy. EI in China has fallen almost continuously while China focuses on the industrial upgrading and promoting transformation of the economic structure. Energy-saving policies and regulations are introduced to help China reach its energy performance target. However, few studies have been carried out to consider the hysteresis effects of policies on the estimation of EI. Therefore, to address such a challenge problem, an improved grey multivariable model is designed to forecast China's EI considering the hysteresis effect of government policies. To further improve its forecasting capability, a nonlinear optimization method based on PSO algorithm is constructed to calculate the hysteresis parameter. In addition, three conventional models, namely GM (1, N), DVCGM (1, N) and ARIMA models, are applied to test the accuracy of this improved DVCGM (1, N) model. The empirical results demonstrate that the proposed model considering the hysteresis effects of energy conservation policies performs best and matches well with the actual observations. Accordingly, this proposed model is used to forecast EI value from 2018 to 2022. The main conclusions are as follows:

1. The improved DVCGM (1, N) model can solve the modeling problem of small sample systems with time-delay causality. A nonlinear optimization method based on PSO algorithm is constructed to calculate the hysteresis parameter. It overcomes the defects of traditional GMs and econometric models.

2. GM (1, N), DVCGM (1, N) and ARIMA model are taken as comparative models. The accuracy of improved DVCGM (1, N) model was tested by the average relative percentage errors. The results show that the Improved DVCGM (1, N) model notes the hysteresis effect of government policies and significantly improves the prediction accuracy of China's EI than the other three models. As suggested by APEs, the overall fitting in descending order is improved DVCGM (1, N) model, DVCGM (1, N), GM (1, N) and ARIMA model.

3. China's EI is greatly influenced by technological progress and is much of policy-driven. When formulating energy conservation and emission reduction policies, we should fully consider the hysteresis effect of the policies, so as to make adjustment of the relative policies and better achieve the national energy performance target.

A few caveats are appropriate. It is an interesting further path to work out the hysteresis parameter directly from the nonlinear programming model. Besides, the sustainability of the hysteresis effect of policy is worth considering. Furthermore, population factors and industrial structure also have good correlation with EI. These will be investigated in our further studies.

This work was supported by the National Social Science Foundation of China (Grant No. 19BJY151).