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

Purpose – The purpose of this paper is to examine the effectiveness of an improved dummy variables control grey model (DVCGM) considering the hysteresis effect of government policies in China’s energy intensity (EI) forecasting. Design/methodology/approach – Energy consumption is considered as an important driver of economic development. China has introduced policies those aim at the optimization of energy structure and EI. In this study, EI is forecasted by an improved DVCGM, considering the hysteresis effect of energy-saving policies of the government. A nonlinear optimization method based on particle swarm optimization (PSO) algorithm is constructed to calculate the hysteresis parameter. A one-step rolling mechanism is applied to provide input data of the prediction model. Grey model (GM) (1, N), DVCGM (1, N) and ARIMAmodel are applied to test the accuracy of the improved DVCGM (1, N) model prediction. Findings –The results show that the improvedDVCGMprovides reliable results andworkswell in simulation and predictions using multivariable data in small sample size and time-lag virtual variable. Accordingly, the improvedDVCGM notes the hysteresis effect of government policies and significantly improves the prediction accuracy of China’s EI than the other three models. Originality/value – This study estimates the EI considering the hysteresis effect of energy-saving policies in China by using an improved DVCGM. The main contribution of this paper is to propose a model to estimate EI, considering the hysteresis effect of energy-saving policies and improve forecasting accuracy.


Introduction
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 2 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(Tien, , 2010, the model accuracy improvement (Tien, 2011;Wang et al., 2016) and the expansion of the GM (Guo et al., 2013;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 GS (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.

Methodology
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.

GM (1, N) model
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 Hysteresis effect of China's EI forecasting 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.
where k 5 2, 3, n, a is the development coefficient and b 1 ; b 2 ; Á Á Á ; b n are the grey input coefficients obtained by the least squares method. To determine these coefficients, the matrix B and Y N are defined as follows: The values of the coefficients a and b 1 ; b 2 ; Á Á Á ; b n can be determined by the following equation: Definition 3. Let Eq (2) be defined as the differential equation (or) called grey reflection equation:

Hysteresis effect of China's EI forecasting
After determining the coefficients of a and b 1 ; b 2 ; Á Á Á ; b n , the differential equation of the GM can be determined by Eq (2). The solution of the above differential equation is as follows: 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: 1 ð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:

DVCGM (1, N) model
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.
is the behavior sequence of the system, j ðj ¼ M þ 1; Á Á Á ; N Þ is the driving factor sequence. Then GM (1, N) model with dummy variable can be expressed as: i ðkÞ is independent quantization variable driver, The matrix B and Y N 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 b 1 ; b 2 ; Á Á Á ; b n , 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: 1 ðnÞ is an estimation of the original sequence, which is simulation values, b x ð0Þ 1 ðn þ 1Þ; b x ð0Þ 1 ðn þ 2Þ; Á Á Á are predictive values.

Improved DVCGM (1, N) model
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 GS behavior and innately ignore the hysteresis effect of historical variables. Therefore, an improved DVCGM (1, N) model is proposed, integrating these above prediction model.
2.3.1 Construction of the improved DVCGM (1, N) 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: The matrix B and Y N are defined as follows: . . .  The least square estimation of the parameter column satisfies the following requirements: (1) When Proof: Substitute k 5 2, 3,. . ., n into the model, you can get x x 1 ðnÞ ¼ −ax That is, by the least square method, Hysteresis effect of China's EI forecasting (1) When n 5 Nþ1, B has an inverse matrix, the equations have a unique solution, we can (2) When n > Nþ1, B is column full rank, the full rank decomposition of B is B 5 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 5 D, so (3) When n < Nþ1, B is a row full rank matrix, D can be taken as a unit matrix, B 5 C, so Definition 6. Let Eq (7) be defined as the differential equation of improved DVCGM (1, N) model: Where k 5 2, 3, n, a is the development coefficient and b 1 ; b 2 ; Á Á Á ; b n are the grey input coefficients, λ i is the hysteresis parameter. PSO algorithm is used to determine the coefficients of a, b 1 ; b 2 ; Á Á Á ; b n , and λ i . Then the differential equation of the GM can be determined by Eq (7). The solution of Eq (7) is as follows: x 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 b x  (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;. 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).

Hysteresis effect of China's EI forecasting
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).

Modeling procedure.
Detailed procedure of the improved DVCGM (1, N) model is illustrated as follows.
Step 1. Collect raw data and establish original sequence X j 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 b 1 ; b 2 ; Á Á Á ; b n 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).

Application
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 GS 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.

Variables selection and data collection
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  Energy intensity Technological progress tons standard coal per 10,000 Yuan Table 1.

Hysteresis effect of China's EI forecasting
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.

Simulation of energy intensity in 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.
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.
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 1 , b 2 . B T 5 (1.251669, 0.938994). Putting the values of the estimated structural parameter and hysteresis parameter in Eq (7-9), we can obtain the optimal TRF as follows:  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: ΔEI ¼ À 0:0311 þ 1:52ΔEI tÀ1 À 0:94ΔEI tÀ2 þ ε tÀ0:99 ε tÀ1 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.  Table 3.  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.

Simulation of energy intensity in China by
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.

Conclusions
Over the past 20 years, China is gradually shifting from a resource-intensive and energydriven 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 policydriven. 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 GS 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. Test Table A4. Estimation of ARIMA model