# Use of BP Neural Networks to Determine China’s Regional CO2 Emission Quota

Authors: Yawei Qi; Wenxiang Peng; Ran Yan; Guangping Rao
Journal: Complexity
 (2021)

---

## Abstract

China declared a long-term commitment at the United Nations General Assembly (UNGA) in 2020 to reduce CO2 emissions. This announcement has been described by Reuters as “the most important climate change commitment in years.” The allocation of China’s provincial CO2 emission quotas (hereafter referred to as quotas) is crucial for building a unified national carbon market, which is an important policy tool necessary to achieve carbon emissions reduction. In the present research, we used historical quota data of China’s carbon emission trading policy pilot areas from 2014 to 2017 to identify alternative features of corporate CO2 emissions and build a backpropagation neural network model (BP) to train the benchmark model. Later, we used the model to calculate the quotas for other regions, provided they implement the carbon emission trading policy. Finally, we added up the quotas to obtain the total national quota. Additionally, considering the perspective of carbon emission terminal, a new characteristic system of quota allocation was proposed in order to retrain BP including the following three aspects: enterprise production, household consumption, and regional environment. The results of the benchmark model and the new models were compared. This feature system not only builds a reasonable quota-related indicator framework but also perfectly matches China’s existing “bottom-up” total control quota approach. Compared with the previous literature, the present report proposes a quota allocation feature system closer to China’s policy and trains BP to obtain reasonable feature weights. The model is very important for the establishment of a unified national carbon emission trading market and the determination of regional quotas in China.


---

## Body

## 1. Introduction

Because of the 2020 COVID-19 epidemic, human beings have become more aware of their relationship with nature and of the importance of sustaining a harmonious coexistence of man and nature. In a time of significant crises, including the COVID-19 epidemic and climate change, the international community agrees that only through the development and implementation of green and low-carbon technologies, society can achieve high-quality economic recovery [1–3]. On Sept. 22, during the General Debate of the 75th Session of the UNGA, Chinese President Xi Jinping declared that China aims to reach CO2 emissions peak before 2030 and achieve carbon neutrality before 2060. According to this, China plans to restore its economy by promoting low-carbon technologies and lifestyle. Reducing CO2 emissions has become an important goal of China’s 14th Five-Year Plan.During the “13th Five-Year Plan” period, the Chinese government learned from the successful experience of the European Union, Japan, and other economies in reducing CO2 emissions and began to explore the use of market-based methods: CO2 emission trading systems [4, 5]. From 2014 to 2019, the central government implemented a pilot CO2 emission trading policy in 7 provinces and cities, and each regional government formulated relevant trading standards and rules. The government in the pilot regions implemented the “allocation + trading” quota management rules for emission-control enterprises, that is, the emission-control enterprises received free quotas issued by the government at the beginning of the performance period. These quotas are determined after an enterprise’s self-inspection and CO2 emissions report is issued and a third-party verification is performed. If an enterprise’s CO2 emission is exceeded/remained during the period, it can be bought/sold in the carbon trading market. After six years, China’s pilot carbon market has increased and has become the world’s second largest carbon market in terms of quota trading volume. Preliminary statistics have shown that a total of 2,837 emission-control agencies, 1082 nonemission-control agencies, and 11,169 individuals have participated in the pilot market. The cumulative transaction volume of the 7 pilot market quotas has reached 406 million tons, and the cumulative transaction volume is about 9.28 billion yuan. By the end of 2019, China’s carbon intensity was reduced by about 48.1% as compared with data from 2005. In addition, nonfossil energy accounted for 15.3% of the primary energy consumption, which means that China has achieved its 2020 emission reduction target ahead of schedule. Thus, China is currently applying market mechanisms to control and reduce greenhouse gas emissions and to promote the green and low-carbon transformation of economic development. Moreover, the implementation of the CO2 emissions trading market represents not only an important institutional innovation for China but also an important policy tool to implement international agreements for emissions reduction. Given the significant emissions reduction results due to the implementation of the CO2 emission trading market policies in the pilot regions, the Chinese government has initiated the creation of a unified national CO2 emission trading market to help all regions in the country to reduce carbon emissions. It is expected that with these new rules, the CO2 emission peak target by 2030 will be achieved as soon as possible. Therefore, the design and implementation of a unified national CO2 emission market is an issue that needs to be studied urgently. The foundation for a proper design and implementation of a CO2 emission trading market program to achieve the intended emission reductions resides in the correct determination of the national quota and quotas for each region (province and city).After finding alternative features of corporate carbon emissions, we used quota data on China’s carbon emission trading policy pilot areas from 2014 to 2017 and the BP model to calculate the quotas of other regions in the sample interval, provided that such regions implemented the carbon emissions trading policy. We obtained the national total quota. With respect to carbon emission terminal, we divided the quota allocation system into three aspects: enterprise production, household consumption, and regional environment and then retrained the BP to obtain new results which were compared to those of the benchmark model. We found out that (1) from 2014 to 2016, China’s total quota displayed a yearly increase and a sudden decrease in 2017. During the initial stage of the national CO2 emission trading market program, the national quota is expected to maintain a relatively high level. Later, during the following 3 to 4 years, through adjustment and adaptation, quotas in each region are expected to show a downward trend and increase in the change rate. This forces enterprises to either participate in the CO2 emission trading market or improve their technology to reduce emissions. (2) Considering the feature system built by adding household consumption and regional environment, the training model displays a smaller loss rate, and the test results (other regional quotas) describe the actual situation in a more accurate way. Thus, when building a unified national CO2 emission trading market and determining quotas for various regions in addition to enterprise production, it is more reasonable to consider a feature system that takes into account household consumption and regional environment. At the same time, this feature system can be used in combination with China’s “bottom-up” total control and postadjustment method. It not only allows regional quota decision makers to predict CO2 emissions in advance through the existing data in the feature system and the trained model before obtaining the final real summary of the CO2 emissions but also allows enterprises considering CO2 price when making investment decisions and trying to make profits in the market.


## 2. Literature Review

Previous research has mainly studied the following two aspects: (a) initial distribution of quotas for CO2 emissions according to different principles and methods and (b) the allocation efficiency according to the initial allocation of carbon emission quotas.The publications in the first area reported the study of distribution subjects and distribution methods. CO2 emission quota allocation can be divided according to the two following perspectives: region and industry. Region includes the initial quota allocation among countries and the initial quota allocation among different provinces within the same country. In this context, the earliest research studied CO2 emission quota allocation among countries. In 1992, the United Nations Conference on Environment and Development established “common but differentiated responsibilities” as the principle of international environmental cooperation. Later, in 1998, Rose et al. introduced the principle of CO2 emission rights allocation, which should include equity and efficiency. However, what kind of allocation method involves “equity” and how to balance the importance of “equity” and “efficiency” are controversial topics [6]. Some scholars in developed countries believe that allocating CO2 emission rights based on population size is in line with the principle of equity [7]. However, other experts in developing countries have proposed that “accumulated emissions per capita” is more related to equity. According to this concept, the allowed CO2 emissions per capita in developing countries during the development stage should be higher than those in developed countries [8, 9]. Later, more researchers used different features and methods to quantitatively analyze the initial quota [10–12]. Based on population, GDP, and CO2 emission data for 132 countries, Wang et al. proposed the Gini coefficient optimization model that optimizes the historical CO2 emissions quota for various countries and is able to project future quotas [13].Chinese scholars have focused on CO2 quota allocation among provinces and industries in China [14–19]. Song et al. considered the comprehensive distribution principle of three indicators: hereditary system, egalitarianism, and payment, to create the provincial environmental fixed cost allocation optimization model (FCAM), which will be able to determine the allocation of provincial CO2 emission rights for 2020, with a more balanced equity and efficiency [20]. From the perspective of equity and efficiency, Yu and Wu used a master-slave hierarchical interactive iterative algorithm based on satisfaction to build a two-level planning model (the upper-level planning model based on equity and the lower-level planning model based on efficiency) to optimize the allocation of CO2 emission rights between provinces [21]. Qian et al. used Chinese enterprise carbon patent data, and from a consumption and production perspective created the stochastic frontier model to measure regional CO2 emission efficiency. Later, they used the estimated efficiency value to numerically simulate the regional allocation of CO2 emission rights [22]. Wu et al. proposed a coupling model of China’s multiregional CGE and CO2 trade (CGE-3MS). The model showed the decision-making and optimization process of trading CO2 units, and they analyzed the impact of the carbon market on the economy and emission-control industries in China under different initial quotas [23].Compared with the first piece of literature, the second one focused on whether the existing allocation of quotas between regions and industries is reasonable and effective. In general, allocation efficiency is used to measure whether the allocation is reasonable. Some Chinese scholars used the original DEA and the zero-sum DEA models to measure quota allocation for Chinese provinces. They found out that the results of the zero-sum DEA model were better than those obtained with the original DEA model [24]. In addition, other scholars revised the original DEA model and proposed a new allocation method that included the evaluation of DEA efficiency and historical CO2 emissions. Later, they used China’s emission commitment as the decision-making guidance and selected the maximization of the average efficiency as the final goal [25, 26]. Besides, after proving that the original DEA model was inefficient to determine CO2 emissions quotas, some scholars used other models to redistribute the quotas of various industries [27–29]. Huang and Zhang used the SBM and RE/CE models to get a more comprehensive efficiency that reveals Chinese energy use and the CO2 emission situation. From the empirical study of 30 regions in China, they found that the southern region of China has the most efficient score, while northeastern China has poor performance. Price factor has a significant influence on energy use and CO2 emission efficient score of some provinces [30].In summary, we have identified two points that should be taken into account when creating regional quota allocation systems and building the corresponding models. First, feature selection should be performed with caution. Factors related to the regional CO2 emission accounting as well as those related to CO2 emissions of the accounting entities (enterprises participating in the CO2 emission trading market) should be selected. Second, when using the feature system data in the allocation model, it is necessary to carefully determine the weight of each feature to ensure the scientific, rigorous, and accurate quota allocation.Therefore, in the present research, we first selected features that are closely related to the accounting entities (emission-control enterprises), such as the number of industrial enterprises exceeding designated size in the region, the energy consumption per unit of industrial added value, and the proportion of coal-fired energy included in the total energy consumption. We also considered features that are related to regional accounting CO2 emissions, such as per capita energy consumption, per capita carbon emissions, forest coverage, and green areas. Then, we built a feature system that can be used in conjunction with China’s “bottom-up” total control and postadjustment method to provide a predictive model. In addition, this model allows enterprises considering CO2 prices when making investment decisions, thereby stimulating corporate green innovation. In addition, for the first time, we considered the use of the BP in order to determine the weight of each feature based on historical quota data. This is different from the previous literature, which used subjective weight determination methods. In addition, our methods are more in line with China’s actual CO2 emission situation and more accurate during calculations. Finally, the feature system and the corresponding BP neural network model proposed in the present research can be used to calculate the quotas of other regions during the same period and also predict future quotas for the same region.


## 3. Empirical Analysis

### 3.1. Data and Feature Statistics

#### 3.1.1. Data

On October 29, 2011, China´s Development and Reform Commission indicated that China should start implementing the pilot carbon emission trading policy. Specifically, Beijing, Tianjin, Shanghai, Chongqing, Guangdong, Hubei, and Shenzhen were labeled as the pilot regions. This allowed each regional government to determine the transaction start time, corresponding transaction subject enterprises, and quotas according to local conditions. Although the central government promulgated the pilot policy in 2011, each pilot region actually started the carbon emission trading system between 2013 and 2014. Therefore, considering the complete CO2 quota and the availability of additional data, we selected 31 regions (5 cities and 26 provinces) as the research sample for the period 2014–2017. With respect to sample pretreatment, the following aspects were considered: (1) given that Shenzhen belongs to Guangdong Province and that both of them are pilot regions, Shenzhen was analyzed separately; (2) regions with serious missing values in historical feature data were eliminated (Tibet); (3) taking into account that the different dimensions of the data may affect the prediction, all features were normalized using(1)bijt=aijt−minaijtmaxaijt−minaijt,where aijt indicates the jth original feature of the i region in the t year, minaijt and maxaijt represent the minimum and maximum values of the original feature j in 31 regions in the t year, respectively, and bjit indicates the j feature of i region in the t year after normalization.Regional features (except that for Shenzhen) were obtained from the National Bureau of Statistics and China Statistical Yearbook. Features of Shenzhen were obtained from the Shenzhen Statistical Yearbook, and the historical quotas data came from the 2014–2018 Beijing Carbon Market Annual Report.China and its regions have not released official data on CO2 emissions. Thus, we used the CO2 emission calculation method given in the 2006 IPPC National Greenhouse Inventory Guidelines and reported by Qi [31]:(2)CO2=∑i=1nEFi×Ei×10−6=∑i=1nCEFi×NCVi×ORi×4412×Ei×10−6,where CO2 indicates carbon dioxide emissions t and EFi represents the carbon dioxide emissions factor for a specific fuel kg−CO2/kg,m3. We considered 11 types of energy consumption including raw coal, coal, coke, crude oil, fuel oil, gasoline, kerosene, diesel, liquefied petroleum gas, refinery dry gas, and natural gas n=11; Ei is the fuel consumption for fuel ikg,m3; CEFi represents the carbon content of the fuel itC/TJ; NCVi indicates the average low calorific value of fuel ikJ/kg,m3; ORi denotes the carbon oxidation rate of fuel i%. The carbon dioxide emission coefficients of various energy sources are shown in Table 1.Table 1
The CO2 emission coefficients of various energy sources.
EnergyNCVCEFOREFRaw coal20 90826.370.941.900 3Coal2090825.80.921.8300Coke2843529.50.932.860 4Crude oil4181620.10.983.020 2Fuel oil4181621.10.983.170 5Gasoline4307018.90.982.925 1Kerosene4307019.50.983.017 9Diesel4265220.20.983.095 9Liquefied petroleum gas5017917.20.983.101 3Refinery dry gas4605518.20.983.011 9Natural gas3893115.30.992.162 2Notes: (1) NCV comes from China “General Principles of Comprehensive Energy Consumption Calculation” (GB/T 2589–2008); (2) CEF and OR come from China’s “Provincial Greenhouse Gas Inventory Compilation Guide” (NDRC Office [2011] No. 1041).


#### 3.1.2. Feature Statistics

Considering the calculation order, there are two carbon allocation methods:(1)
The “up-bottom” allocation method is applied from a macro (provincial) perspective, according to the general principle of efficiency and equity. This methodology considers population, provincial economic level (GDP), degree of industrialization (industrial structure), historical factors (accumulated carbon emissions per capita), and natural endowments to determine quotas for the different regions and for enterprises. In addition, the method allows determining the amount of quotas in advance, in such a way that participating entities are able to consider CO2 prices in decision-making processes. However, due to the lack of data regarding actual emissions of the participating entities, when the quotas correspond to the actual situation and, in consequence, entities are motivated to participate in market transactions emissions reduction, these processes display great uncertainty.(2)
In the case of the “bottom-up” allocation method, the level-by-level summary determines the quotas for each region through the calculation of emissions at the microlevel (emission terminal). Considering carbon dioxide production terminal in human society (mainly consumers and enterprises), the enterprise component considers the number of equipment that emits CO2, the scale of the enterprise, the scientific research level of enterprise emissions reduction (patents and R&D investment), and so on. The consumer component involves population size, per capita carbon emissions, and so on. Taking into account the emission data of the participating entities, this method is suitable for adjusting quotas afterwards (without specifying quotas in advance). However, it may eventually result in oversupply of quotas due to the excessively high emission limits. This situation may discourage participants to consider CO2 prices when they are making investment decisions. In consequence, participants are not motivated to reduce emissions.At the same time, given the problem of excessive CO2 emissions caused by humans, people are also planting trees (forest carbon sinks) and using photosynthesis to reduce emissions. Therefore, resident life should also be considered when calculating quotas.In 2020, the central government will formally begin to build a national carbon emission trading market. First, it will issue a national carbon emission trading quota setting and allocation implementation plan for the power industry (2019–2020). The plan involves a “bottom-up” quota determination method, that is, relevant provincial departments will be required to determine the list of key emitters and their actual output. Later, they will identify key emitters’ quotas based on the benchmark method (the free quotas in each region for 2019–2020 were preallocated according to the 70% of the power (heat) supply of each key emitter in 2018). Then, after the quotas of all the key emitters in each region are verified, they will be added up to form the total quota of the region, and the regional quotas will be further added to obtain the total quota of the country.According to China’s “bottom-up” quota allocation method, it is most accurate to use the CO2 emission equipment data of emission-control enterprises in the regions. However, since China’s national carbon emission trading market is still in its infancy, relevant data (corporate power supply (heat) units, actual output) of emission-control enterprises in regions have not been released. Also, “bottom-up” quota allocation methodologies present several limitations, which may result in oversupply of quotas, insufficient demand, low CO2 price, and inactive market due to excessively high emission limits. Because of this, participants usually do not consider CO2 prices when they are in the process of investment decision-making and cannot effectively motivate participants to reduce emissions.Therefore, based on the operating experience of The European Union Emission Trading Scheme (EU ETS) and Regional Greenhouse Gas Initiative (RGGI), we proposed a more reasonable “top-down” allocation method in combination with China’s existing “bottom-up” quota allocation method [32]. This method was divided into two steps. The first step involved the development of a feature system that is more suitable for China’s provincial carbon quota allocation. This was performed by selecting features that are related to participating entities. For example, (a) the current emission-control enterprises in China are mainly industrial enterprises; thus, the region that contains more industrial enterprises will have more quotas; (b) energy consumption, including coal-fired energy that is used by emission-control enterprises in production activities ranks first in regional energy consumption; thus, energy structure and energy consumption per unit GDP are both factors affecting regional quotas. The second step involves training suitable BP to calculate nonpilot regional quotas based on historical pilot regional quotas and feature system data.Most ETS are based on total CO2 emission control. Determination of total CO2 emissions should not only consider the overall emission target, but also the regional differences (level of economic development, technological differences, and forest carbon sink).Based on previous reports and considering different factors included in the “top-down” quota control approach adopted by the EU ETS after phase II [33–35], we divided the features into three categories: (a) enterprise production, (b) household consumption, and (c) regional environment.(a)
Household consumption: in 2019, 30% of China’s CO2 emissions corresponded to consumer and man-made emissions. Household consumption, as the main body of society, is one of the main causes of CO2 emissions. Regions with large populations display frequent economic activities. Therefore, resident consumption factors depend on per capita GDP, energy consumption, carbon emissions, cumulative carbon emissions, and disposable income.(b)
Enterprise production: in 2019, 70% of China’s CO2 emissions were the result of industrial production or generative emissions. Among them, the carbon dioxide emissions related to the power industry accounted for more than 40%. In addition, those from the steel industry, which are part of the manufacturing industry, accounted for about 15%. Industrial enterprises exceeding permitted size are not only the main body of the industry but also the main entities participating in the carbon emissions trading market. Therefore, factors related to CO2 emissions from enterprises include the number of industrial enterprises exceeding permitted size, full-time equivalent of R&D personnel in industrial enterprises exceeding designated size, industrial structure, and energy structure, energy consumption per GDP unit, energy consumption per unit of industrial added value, and electricity consumption per GDP unit.(c)
Regional environment: in addition to factors related to resident and enterprises, regional environment also affects CO2 emissions and decomposition. We considered three aspects, regional economic level, technological level, and green resources, mainly including total freight volume, total passenger volume, degree of openness to the outside world, total gas supply, total supply of liquefied petroleum gas, urbanization rate, technological level, green area, and forest carbon sink. Specific explanations of features are given in Table 2.Table 2
Specific explanation of features.
FeatureExplanationGDP per capitaGDP/populationEnergy consumption per capitaEnergy consumption/populationDegree of opening to the outside worldTotal import and export/GDPCarbon emissions per capitaTotalCO2emissions/populationCumulative carbon emissions per capitaCumulativeCO2emissions/Cumulative populationTotal passenger volumeThe number of passengers actually carried by the means of transport in a certain period of timeDisposable income per capitaDisposable income/populationIndustrial structureAdded value of tertiary industry/GDPEnergy consumption per unit of GDPEnergy consumption/GDPEnergy consumption per unit of industrial added valueEnergy consumption/industrial added valueElectricity consumption per unit of GDPElectricity consumption/GDPEnergy structureCoal-fired energy/energy consumptionTotal freight volumeIn a certain period of time, the actual weight of the cargo carried by the means of transportThe number of industrial enterprises above designated sizeIncluding independent auditor industrial enterprises and affiliated industrial production unitsFull-time equivalent of R&D personnel in industrial enterprises above designated sizeThe number of full-time staff plus part-time staff is converted to the total number of full-time staff based on workloadTechnological levelNumber of patents granted/10,000 peopleForest carbon sinkForest area/total land areaGreen areaIncluding park area, production green area, protective green area, subsidiary green area, and other green areaTotal gas supplyRefers to the amount of natural gas supplied by gas units during the reporting periodTotal supply of liquefied petroleum gasRefers to the amount of liquefied petroleum gas supplied by gas units during the reporting periodUrbanization rateUrban population/population


### 3.2. Empirical Model

#### 3.2.1. Specifications of the Backpropagation Neural Network Model (BP)

Compared with the subjective weight assignment methods (AHP, Expert Evaluation Method, TOPSIS, etc.) used in the formation of the “top-down” model and reported in previous literature, BP involves a multilayer feedforward network and error direction propagation-learning algorithm. Because of its unique adaptability, learning ability, and strong generalization ability, it is widely used in the fields of automatic identification, predictive estimation, engineering, biology, and medicine, among others. For the purpose of the present research, BP can more objectively and accurately quantify the impact of features on quotas and dynamically reflect the nonlinear impact of features on quotas at different stages [36]. Therefore, after training the BP based on pilot regional quotas and feature system data, we calculated the nonpilot regional quotas.BP is composed of an input layer, a hidden layer, and an output layer. These three basic elements are fully connected during the whole network training. For example, for a neural network model with only one hidden layer, the process of BP neural network is mainly divided into two stages.The first stage is the forward propagation of the signal, from the input layer to the hidden layer and finally to the output layer. Assuming that the number of samples isA, input layer has m nodes, output layer has n nodes, and hidden layer has p nodes; xai is the input/output of the input layer, a = 1, 2, …, A, i = 1, 2, …, m; Baj, baj is the input/output of the hidden layer, respectively, j = 1, 2, …, p; Yak, and yak are the input/output of the output layer, respectively, k = 1, 2, …, n; yak∗ is the expected label of the output layer (historical quotas); wij and wjk are the weights from input to hidden and hidden to output, respectively; ，θj and γk are the biases from input to hidden and hidden to output, respectively.Input layer to hidden layer: determine input function,Baj=xaiwij+θj, and then transform Baj into baj through activation function baj=fBaj.Hidden layer to output layer: determine input function,Yak=bajwjk+γk, and then transform Yak into yak through activation function, yak=fYak, where yak is the final result.Later, the loss function is determined and loss is calculated according toyak and yak∗. When the loss is either smaller than the set range or reaches the upper limit of the number of iterations, the model ends the training; otherwise, it enters a second stage.The second stage is the backpropagation of the loss. The loss information is returned along the original propagation path through the learning signal. Starting from the last layer, the weight and bias are corrected layer by layer, and finally the loss is within the set range.Weight and bias update formula:(3)Wl∗=Wl+ΔWl,θl∗=θl+Δθl,ΔWl=−η∂Loss∂Wl,Δθl=−η∂Loss∂θl,where Wl∗ and θl∗ indicate the updated weight and bias of the l layer, Wl and θl indicate the original weight and bias of the l layer, ΔWl and Δθl represent the correction part of weight and bias of the l layer, and η is a fixed value that indicates the learning rate.The fundamental part of backpropagation is to minimize the loss through the update of weights and biases, using the gradient descent method (actually using the chain partial derivative). The specific derivation process is given below.The following assumptions are considered:The activation function from the input layer to the hidden layer isfx=1/1+e−x (sigmoid).The activation function from the input layer to the hidden layer isfx=x.Loss function isLoss=1/2×∑k=1nyak∗−yak2=1/2×∑k=1neak2, eak=yak∗−yak.The weights update process of the hidden-output layer:(4)∂Loss∂wjk=∂Loss∂eak×∂eak∂yak×∂yak∂Yak×∂Yak∂wjk=eak×−1×1×baj=−bajeak,wjk∗=wjk−η∂Loss∂wjk=wjk+ηbajeak.The weight update process of the input-hidden layer:(5)∂Loss∂wij=∂Loss∂eak×∂eak∂yak×∂yak∂Yak×∂Yak∂baj×∂baj∂Bij×∂Bij∂wij=eak×−1×wjk×baj1−baj×xai=−baj1−bajwjkxaieak,wij∗=wij−η∂Loss∂wij=wij+ηbaj1−bajeakwjkxai.In the same way, the updated value of the bias can be obtained:Hidden-output layer:γk∗=γk+ηeakInput-hidden layer:θj∗=θj+ηbaj1−bajwjkeak


#### 3.2.2. Key Parameters

According to the model principle previously mentioned, there are 5 key parameters that determine the learning effect of the BP: activation function, loss function, learning rate, the number of hidden layer nodes, and gradient descent algorithm.(1)
Activation function:For the BP both, the hidden layer and the output layer need to use an activation function.For the hidden layer, the activation function is generally a nonlinear function. The reason for this is that, if the activation function is a linear function, the output is a linear combination of the input, which is equivalent to the effect of the no hidden layer (the hidden layer is invalid). The introduction of a nonlinear function as the hidden layer activation function makes the network more powerful, increases its ability to learn complex data, and reflects the nonlinear relationship between input and output. Therefore, we introduced four nonlinear activation functions that are widely used (Tables3 and 4).Although ReLU has two problems, it is currently the most commonly used activation function for BP. In addition, it is the default activation function used by most feedforward neural networks.For the output layer, the choice of its activation function depends on whether the problem is a regression problem or a classification problem. In the event, it is a classification problem, the sigmoid activation function represents a good choice; for regression problems, a linear activation function is more appropriate.(2)
Loss function:With regards to the problem and the output layer activation function to match different loss functions,(1)
Cross-entropy function: it is suitable for binary classification problems, and the output layer activation function is sigmoid(2)
Log-likelihood cost: it is suitable for multiclassification problems, and the output layer activation function is softmax(3)
Mean square error (MSE): it is suitable for regression problems, and the output layer activation function is a linear function(3)
Learning rate:The learning rate value is an important part of the BP, which represents the speed of information accumulation in the neural network over time, and its value is between [0, 1]. Under ideal circumstances, we would start with a large learning rate and gradually reduce the speed until the loss value no longer diverges (if the learning rate is set too low, the training progress will be very slow because only very few adjustments to the weight of the network are made. However, if the learning rate is set too high, it may bring undesirable consequences on the loss function (Figure1)).(4)
The number of hidden layer nodes:The number of hidden layer nodes has a great influence on the prediction accuracy of BP; if the number of nodes is too small, the network cannot perform a proper learning process, and it will need more times to train. In addition, the training accuracy is also affected. When the number of nodes is too large, the training time increases and the network will result in overfitting. However, there is no conventional formula for determining the number of nodes. Some empirical formulas are given below for reference:l<n−1l<m+n+a, l=log2n,l≥k×n/n+m, where n indicates the number of input layer nodes, l indicates the number of hidden layer nodes, m is the number of output layer nodes, a represents a constant between 0–10, and k corresponds to the number of samples.In fact, the number of hidden layer nodes can be roughly calculated according to the reference formula. Later, trial and error is used to find the optimal number of nodes. Generally speaking, the BP error shows a trend where it first decreases and later increases with the increase of hidden layer nodes.(5)
Gradient descent algorithm:We introduce six well-known gradient descent algorithms (Tables5 and 6).According to this analysis, there are no perfect key parameters that can suit all conditions. The appropriate selection of key parameters depends on the specific problem of study. In the present research, we studied a regression problem. Thus, we chose ReLU andfx=x as the activation function and MSE as the loss function. Also, Adam may be appropriate as the gradient descent algorithm; however, the learning rate and the number of hidden layer nodes cannot be determined in advance. In summary, the final determination of all key parameters needs BP training. For this reason, the determination of key parameters is provided in Section 4.Table 3
Activation functions and figures.
NameFunctionDerivative functionSigmoidSx=1/1+e−xS′x=Sx1−SxSoftmaxSoftmaxZi=eZi/∑c=1CeZc=pi∂softmaxZi/∂Zj=pi1−pi,j=i,−pj⋅pi,j≠i.tanhtanhx=ex−e−x/ex+e−xT′x=1−Tx2ReLUReLUx=max0,xReLU′x=0,x≤0,1,x>0.Table 4
Advantages and disadvantages of activation functions.
NameAdvantagesDisadvantagesSigmoid(1) It can smoothly map the real-number field to [0, 1](1) Gradient vanishing, that is, in the process of backpropagation, the derivative will gradually become 0; thus, the parameters cannot be updated and the neural network cannot be optimized(2) Monotonically increasing, continuous derivable, and its derivative form is very simple(2) Nonzero-centered: the output value of the function is always greater than 0, which will slow down the convergence speed of the model training(3) Suitable for handling binary classification problems(3) Exponentiation is relatively time-consumingSoftmax(1) It maps the output value to (0, 1), and the sum of the mapped output value is 1(1) The operation of Softmax involves the calculation of exponential function; in consequence, an “overflow problem” for computers occurs(2) It divides the entire hyperspace according to the number of classifications(3) Suitable for multiclassification problems(2) Not suitable for face recognition taskstanh(1) It can smoothly map the real-number field to [−1, 1](1) Gradient vanishing, that is, in the process of backpropagation, the derivative will gradually become 0; thus, the parameters cannot be updated and the neural network cannot be optimized(2) Solve nonzero-centered problem(3) Suitable for handling binary classification problemsReLU(1) Solve the gradient vanishing in the positive interval(1) Nonzero-centered(2) The calculation is simple, no exponential calculation is required, and the activation value can be obtained with only one value(2) Dead ReLU problem: it is “vulnerable” during training; whenx < 0, the gradient is 0; the gradient of these nodes and subsequent nodes are always 0, and no longer responds to any data, causing the corresponding parameters to never be updated(3) The convergence speed is much faster than sigmoid and tanhFigure 1
The loss function results of different learning rates.
(a)(b)(c)Table 5
Gradient descent functions and parameter description.
NameGradient descent functionParameter descriptionBGDθ∗=θ−η⋅g⌢=θ−η⋅∇θJθg⌢ denotes gradient estimateθ denotes weightSGDθ∗=θ−η⋅g⌢=θ−η⋅∇θJθ;xi;yig⌢ denotes gradient estimateθ denotes weightMomentumvt=γvt−1+η⋅∇θJθvt denotes gradient estimateθ∗=θ−vtθ denotes weightAdaGradg⌢t,i=∇θJθig⌢t,i denotes gradient estimateθt+1,i∗=θt,i−η/g⌢t,i⋅g⌢t,i+ε⋅g⌢t,iθt,i denotes weightRMSpropgt=∇θJθt−1γ=0.9Eg2t denotes Exponential weighted averageEg2t=γEg2t−1+1−γgt2gt denotes gradient estimateθt+1∗=θt−η/Eg2t+εgtθt denotes weightAdamgt=∇θJθt−1β1=0.9β2=0.999ε=10e−8mt=β1mt−1+1−β1gtgt denotes gradient estimate θt denotes weightvt=β2vt−1+1−β2gt2m⌢t=mt/1−β1t, v⌢t=vt/1−β2tθt+1∗=θt−η/v⌢t+εm⌢tTable 6
Advantages and disadvantages of gradient descent algorithms.
NameAdvantagesDisadvantagesBGDThe principle of gradient descent is simple(1) Calculation is very slow(2) Difficult to handle a large dataset(3) Cannot add new data to update the modelSGD(1) Compared with BGD, SGD training speed is faster(1) Frequent updates may cause severe oscillations in the loss function(2) New data can be added to update the modelMomentum(1) Consider the speed of the previous step and the new gradient(2) Can speed up the convergence and suppress the shockAdaGrad(1) Compared with SGD, it adds a denominator(1) If gradient update is frequent, it may cause the subsequent gradient updates be slow or disappear(2) Handle the case where the number of gradient updates is smallRMSprop(1) Similar to momentum, it can reduce fluctuations(2) Overcome the problem of the sharp decrease or disappearance of gradient in AdaGrad(3) It performs better than SGD, momentum, and AdaGrad, based on the nonstationary objective functionAdam(1) Combine momentum and RMSProp(2) Integrate the contents of gradient descent, momentum, Adagrad, and RMSprop with certain improvement(3) Easy to use, insensitive to gradient scaling, can be used for large data, processing sparse data, easy to adjust hyperparameters, etc.


## 3.1. Data and Feature Statistics

### 3.1.1. Data

On October 29, 2011, China´s Development and Reform Commission indicated that China should start implementing the pilot carbon emission trading policy. Specifically, Beijing, Tianjin, Shanghai, Chongqing, Guangdong, Hubei, and Shenzhen were labeled as the pilot regions. This allowed each regional government to determine the transaction start time, corresponding transaction subject enterprises, and quotas according to local conditions. Although the central government promulgated the pilot policy in 2011, each pilot region actually started the carbon emission trading system between 2013 and 2014. Therefore, considering the complete CO2 quota and the availability of additional data, we selected 31 regions (5 cities and 26 provinces) as the research sample for the period 2014–2017. With respect to sample pretreatment, the following aspects were considered: (1) given that Shenzhen belongs to Guangdong Province and that both of them are pilot regions, Shenzhen was analyzed separately; (2) regions with serious missing values in historical feature data were eliminated (Tibet); (3) taking into account that the different dimensions of the data may affect the prediction, all features were normalized using(1)bijt=aijt−minaijtmaxaijt−minaijt,where aijt indicates the jth original feature of the i region in the t year, minaijt and maxaijt represent the minimum and maximum values of the original feature j in 31 regions in the t year, respectively, and bjit indicates the j feature of i region in the t year after normalization.Regional features (except that for Shenzhen) were obtained from the National Bureau of Statistics and China Statistical Yearbook. Features of Shenzhen were obtained from the Shenzhen Statistical Yearbook, and the historical quotas data came from the 2014–2018 Beijing Carbon Market Annual Report.China and its regions have not released official data on CO2 emissions. Thus, we used the CO2 emission calculation method given in the 2006 IPPC National Greenhouse Inventory Guidelines and reported by Qi [31]:(2)CO2=∑i=1nEFi×Ei×10−6=∑i=1nCEFi×NCVi×ORi×4412×Ei×10−6,where CO2 indicates carbon dioxide emissions t and EFi represents the carbon dioxide emissions factor for a specific fuel kg−CO2/kg,m3. We considered 11 types of energy consumption including raw coal, coal, coke, crude oil, fuel oil, gasoline, kerosene, diesel, liquefied petroleum gas, refinery dry gas, and natural gas n=11; Ei is the fuel consumption for fuel ikg,m3; CEFi represents the carbon content of the fuel itC/TJ; NCVi indicates the average low calorific value of fuel ikJ/kg,m3; ORi denotes the carbon oxidation rate of fuel i%. The carbon dioxide emission coefficients of various energy sources are shown in Table 1.Table 1
The CO2 emission coefficients of various energy sources.
EnergyNCVCEFOREFRaw coal20 90826.370.941.900 3Coal2090825.80.921.8300Coke2843529.50.932.860 4Crude oil4181620.10.983.020 2Fuel oil4181621.10.983.170 5Gasoline4307018.90.982.925 1Kerosene4307019.50.983.017 9Diesel4265220.20.983.095 9Liquefied petroleum gas5017917.20.983.101 3Refinery dry gas4605518.20.983.011 9Natural gas3893115.30.992.162 2Notes: (1) NCV comes from China “General Principles of Comprehensive Energy Consumption Calculation” (GB/T 2589–2008); (2) CEF and OR come from China’s “Provincial Greenhouse Gas Inventory Compilation Guide” (NDRC Office [2011] No. 1041).


### 3.1.2. Feature Statistics

Considering the calculation order, there are two carbon allocation methods:(1)
The “up-bottom” allocation method is applied from a macro (provincial) perspective, according to the general principle of efficiency and equity. This methodology considers population, provincial economic level (GDP), degree of industrialization (industrial structure), historical factors (accumulated carbon emissions per capita), and natural endowments to determine quotas for the different regions and for enterprises. In addition, the method allows determining the amount of quotas in advance, in such a way that participating entities are able to consider CO2 prices in decision-making processes. However, due to the lack of data regarding actual emissions of the participating entities, when the quotas correspond to the actual situation and, in consequence, entities are motivated to participate in market transactions emissions reduction, these processes display great uncertainty.(2)
In the case of the “bottom-up” allocation method, the level-by-level summary determines the quotas for each region through the calculation of emissions at the microlevel (emission terminal). Considering carbon dioxide production terminal in human society (mainly consumers and enterprises), the enterprise component considers the number of equipment that emits CO2, the scale of the enterprise, the scientific research level of enterprise emissions reduction (patents and R&D investment), and so on. The consumer component involves population size, per capita carbon emissions, and so on. Taking into account the emission data of the participating entities, this method is suitable for adjusting quotas afterwards (without specifying quotas in advance). However, it may eventually result in oversupply of quotas due to the excessively high emission limits. This situation may discourage participants to consider CO2 prices when they are making investment decisions. In consequence, participants are not motivated to reduce emissions.At the same time, given the problem of excessive CO2 emissions caused by humans, people are also planting trees (forest carbon sinks) and using photosynthesis to reduce emissions. Therefore, resident life should also be considered when calculating quotas.In 2020, the central government will formally begin to build a national carbon emission trading market. First, it will issue a national carbon emission trading quota setting and allocation implementation plan for the power industry (2019–2020). The plan involves a “bottom-up” quota determination method, that is, relevant provincial departments will be required to determine the list of key emitters and their actual output. Later, they will identify key emitters’ quotas based on the benchmark method (the free quotas in each region for 2019–2020 were preallocated according to the 70% of the power (heat) supply of each key emitter in 2018). Then, after the quotas of all the key emitters in each region are verified, they will be added up to form the total quota of the region, and the regional quotas will be further added to obtain the total quota of the country.According to China’s “bottom-up” quota allocation method, it is most accurate to use the CO2 emission equipment data of emission-control enterprises in the regions. However, since China’s national carbon emission trading market is still in its infancy, relevant data (corporate power supply (heat) units, actual output) of emission-control enterprises in regions have not been released. Also, “bottom-up” quota allocation methodologies present several limitations, which may result in oversupply of quotas, insufficient demand, low CO2 price, and inactive market due to excessively high emission limits. Because of this, participants usually do not consider CO2 prices when they are in the process of investment decision-making and cannot effectively motivate participants to reduce emissions.Therefore, based on the operating experience of The European Union Emission Trading Scheme (EU ETS) and Regional Greenhouse Gas Initiative (RGGI), we proposed a more reasonable “top-down” allocation method in combination with China’s existing “bottom-up” quota allocation method [32]. This method was divided into two steps. The first step involved the development of a feature system that is more suitable for China’s provincial carbon quota allocation. This was performed by selecting features that are related to participating entities. For example, (a) the current emission-control enterprises in China are mainly industrial enterprises; thus, the region that contains more industrial enterprises will have more quotas; (b) energy consumption, including coal-fired energy that is used by emission-control enterprises in production activities ranks first in regional energy consumption; thus, energy structure and energy consumption per unit GDP are both factors affecting regional quotas. The second step involves training suitable BP to calculate nonpilot regional quotas based on historical pilot regional quotas and feature system data.Most ETS are based on total CO2 emission control. Determination of total CO2 emissions should not only consider the overall emission target, but also the regional differences (level of economic development, technological differences, and forest carbon sink).Based on previous reports and considering different factors included in the “top-down” quota control approach adopted by the EU ETS after phase II [33–35], we divided the features into three categories: (a) enterprise production, (b) household consumption, and (c) regional environment.(a)
Household consumption: in 2019, 30% of China’s CO2 emissions corresponded to consumer and man-made emissions. Household consumption, as the main body of society, is one of the main causes of CO2 emissions. Regions with large populations display frequent economic activities. Therefore, resident consumption factors depend on per capita GDP, energy consumption, carbon emissions, cumulative carbon emissions, and disposable income.(b)
Enterprise production: in 2019, 70% of China’s CO2 emissions were the result of industrial production or generative emissions. Among them, the carbon dioxide emissions related to the power industry accounted for more than 40%. In addition, those from the steel industry, which are part of the manufacturing industry, accounted for about 15%. Industrial enterprises exceeding permitted size are not only the main body of the industry but also the main entities participating in the carbon emissions trading market. Therefore, factors related to CO2 emissions from enterprises include the number of industrial enterprises exceeding permitted size, full-time equivalent of R&D personnel in industrial enterprises exceeding designated size, industrial structure, and energy structure, energy consumption per GDP unit, energy consumption per unit of industrial added value, and electricity consumption per GDP unit.(c)
Regional environment: in addition to factors related to resident and enterprises, regional environment also affects CO2 emissions and decomposition. We considered three aspects, regional economic level, technological level, and green resources, mainly including total freight volume, total passenger volume, degree of openness to the outside world, total gas supply, total supply of liquefied petroleum gas, urbanization rate, technological level, green area, and forest carbon sink. Specific explanations of features are given in Table 2.Table 2
Specific explanation of features.
FeatureExplanationGDP per capitaGDP/populationEnergy consumption per capitaEnergy consumption/populationDegree of opening to the outside worldTotal import and export/GDPCarbon emissions per capitaTotalCO2emissions/populationCumulative carbon emissions per capitaCumulativeCO2emissions/Cumulative populationTotal passenger volumeThe number of passengers actually carried by the means of transport in a certain period of timeDisposable income per capitaDisposable income/populationIndustrial structureAdded value of tertiary industry/GDPEnergy consumption per unit of GDPEnergy consumption/GDPEnergy consumption per unit of industrial added valueEnergy consumption/industrial added valueElectricity consumption per unit of GDPElectricity consumption/GDPEnergy structureCoal-fired energy/energy consumptionTotal freight volumeIn a certain period of time, the actual weight of the cargo carried by the means of transportThe number of industrial enterprises above designated sizeIncluding independent auditor industrial enterprises and affiliated industrial production unitsFull-time equivalent of R&D personnel in industrial enterprises above designated sizeThe number of full-time staff plus part-time staff is converted to the total number of full-time staff based on workloadTechnological levelNumber of patents granted/10,000 peopleForest carbon sinkForest area/total land areaGreen areaIncluding park area, production green area, protective green area, subsidiary green area, and other green areaTotal gas supplyRefers to the amount of natural gas supplied by gas units during the reporting periodTotal supply of liquefied petroleum gasRefers to the amount of liquefied petroleum gas supplied by gas units during the reporting periodUrbanization rateUrban population/population


## 3.1.1. Data

On October 29, 2011, China´s Development and Reform Commission indicated that China should start implementing the pilot carbon emission trading policy. Specifically, Beijing, Tianjin, Shanghai, Chongqing, Guangdong, Hubei, and Shenzhen were labeled as the pilot regions. This allowed each regional government to determine the transaction start time, corresponding transaction subject enterprises, and quotas according to local conditions. Although the central government promulgated the pilot policy in 2011, each pilot region actually started the carbon emission trading system between 2013 and 2014. Therefore, considering the complete CO2 quota and the availability of additional data, we selected 31 regions (5 cities and 26 provinces) as the research sample for the period 2014–2017. With respect to sample pretreatment, the following aspects were considered: (1) given that Shenzhen belongs to Guangdong Province and that both of them are pilot regions, Shenzhen was analyzed separately; (2) regions with serious missing values in historical feature data were eliminated (Tibet); (3) taking into account that the different dimensions of the data may affect the prediction, all features were normalized using(1)bijt=aijt−minaijtmaxaijt−minaijt,where aijt indicates the jth original feature of the i region in the t year, minaijt and maxaijt represent the minimum and maximum values of the original feature j in 31 regions in the t year, respectively, and bjit indicates the j feature of i region in the t year after normalization.Regional features (except that for Shenzhen) were obtained from the National Bureau of Statistics and China Statistical Yearbook. Features of Shenzhen were obtained from the Shenzhen Statistical Yearbook, and the historical quotas data came from the 2014–2018 Beijing Carbon Market Annual Report.China and its regions have not released official data on CO2 emissions. Thus, we used the CO2 emission calculation method given in the 2006 IPPC National Greenhouse Inventory Guidelines and reported by Qi [31]:(2)CO2=∑i=1nEFi×Ei×10−6=∑i=1nCEFi×NCVi×ORi×4412×Ei×10−6,where CO2 indicates carbon dioxide emissions t and EFi represents the carbon dioxide emissions factor for a specific fuel kg−CO2/kg,m3. We considered 11 types of energy consumption including raw coal, coal, coke, crude oil, fuel oil, gasoline, kerosene, diesel, liquefied petroleum gas, refinery dry gas, and natural gas n=11; Ei is the fuel consumption for fuel ikg,m3; CEFi represents the carbon content of the fuel itC/TJ; NCVi indicates the average low calorific value of fuel ikJ/kg,m3; ORi denotes the carbon oxidation rate of fuel i%. The carbon dioxide emission coefficients of various energy sources are shown in Table 1.Table 1
The CO2 emission coefficients of various energy sources.
EnergyNCVCEFOREFRaw coal20 90826.370.941.900 3Coal2090825.80.921.8300Coke2843529.50.932.860 4Crude oil4181620.10.983.020 2Fuel oil4181621.10.983.170 5Gasoline4307018.90.982.925 1Kerosene4307019.50.983.017 9Diesel4265220.20.983.095 9Liquefied petroleum gas5017917.20.983.101 3Refinery dry gas4605518.20.983.011 9Natural gas3893115.30.992.162 2Notes: (1) NCV comes from China “General Principles of Comprehensive Energy Consumption Calculation” (GB/T 2589–2008); (2) CEF and OR come from China’s “Provincial Greenhouse Gas Inventory Compilation Guide” (NDRC Office [2011] No. 1041).


## 3.1.2. Feature Statistics

Considering the calculation order, there are two carbon allocation methods:(1)
The “up-bottom” allocation method is applied from a macro (provincial) perspective, according to the general principle of efficiency and equity. This methodology considers population, provincial economic level (GDP), degree of industrialization (industrial structure), historical factors (accumulated carbon emissions per capita), and natural endowments to determine quotas for the different regions and for enterprises. In addition, the method allows determining the amount of quotas in advance, in such a way that participating entities are able to consider CO2 prices in decision-making processes. However, due to the lack of data regarding actual emissions of the participating entities, when the quotas correspond to the actual situation and, in consequence, entities are motivated to participate in market transactions emissions reduction, these processes display great uncertainty.(2)
In the case of the “bottom-up” allocation method, the level-by-level summary determines the quotas for each region through the calculation of emissions at the microlevel (emission terminal). Considering carbon dioxide production terminal in human society (mainly consumers and enterprises), the enterprise component considers the number of equipment that emits CO2, the scale of the enterprise, the scientific research level of enterprise emissions reduction (patents and R&D investment), and so on. The consumer component involves population size, per capita carbon emissions, and so on. Taking into account the emission data of the participating entities, this method is suitable for adjusting quotas afterwards (without specifying quotas in advance). However, it may eventually result in oversupply of quotas due to the excessively high emission limits. This situation may discourage participants to consider CO2 prices when they are making investment decisions. In consequence, participants are not motivated to reduce emissions.At the same time, given the problem of excessive CO2 emissions caused by humans, people are also planting trees (forest carbon sinks) and using photosynthesis to reduce emissions. Therefore, resident life should also be considered when calculating quotas.In 2020, the central government will formally begin to build a national carbon emission trading market. First, it will issue a national carbon emission trading quota setting and allocation implementation plan for the power industry (2019–2020). The plan involves a “bottom-up” quota determination method, that is, relevant provincial departments will be required to determine the list of key emitters and their actual output. Later, they will identify key emitters’ quotas based on the benchmark method (the free quotas in each region for 2019–2020 were preallocated according to the 70% of the power (heat) supply of each key emitter in 2018). Then, after the quotas of all the key emitters in each region are verified, they will be added up to form the total quota of the region, and the regional quotas will be further added to obtain the total quota of the country.According to China’s “bottom-up” quota allocation method, it is most accurate to use the CO2 emission equipment data of emission-control enterprises in the regions. However, since China’s national carbon emission trading market is still in its infancy, relevant data (corporate power supply (heat) units, actual output) of emission-control enterprises in regions have not been released. Also, “bottom-up” quota allocation methodologies present several limitations, which may result in oversupply of quotas, insufficient demand, low CO2 price, and inactive market due to excessively high emission limits. Because of this, participants usually do not consider CO2 prices when they are in the process of investment decision-making and cannot effectively motivate participants to reduce emissions.Therefore, based on the operating experience of The European Union Emission Trading Scheme (EU ETS) and Regional Greenhouse Gas Initiative (RGGI), we proposed a more reasonable “top-down” allocation method in combination with China’s existing “bottom-up” quota allocation method [32]. This method was divided into two steps. The first step involved the development of a feature system that is more suitable for China’s provincial carbon quota allocation. This was performed by selecting features that are related to participating entities. For example, (a) the current emission-control enterprises in China are mainly industrial enterprises; thus, the region that contains more industrial enterprises will have more quotas; (b) energy consumption, including coal-fired energy that is used by emission-control enterprises in production activities ranks first in regional energy consumption; thus, energy structure and energy consumption per unit GDP are both factors affecting regional quotas. The second step involves training suitable BP to calculate nonpilot regional quotas based on historical pilot regional quotas and feature system data.Most ETS are based on total CO2 emission control. Determination of total CO2 emissions should not only consider the overall emission target, but also the regional differences (level of economic development, technological differences, and forest carbon sink).Based on previous reports and considering different factors included in the “top-down” quota control approach adopted by the EU ETS after phase II [33–35], we divided the features into three categories: (a) enterprise production, (b) household consumption, and (c) regional environment.(a)
Household consumption: in 2019, 30% of China’s CO2 emissions corresponded to consumer and man-made emissions. Household consumption, as the main body of society, is one of the main causes of CO2 emissions. Regions with large populations display frequent economic activities. Therefore, resident consumption factors depend on per capita GDP, energy consumption, carbon emissions, cumulative carbon emissions, and disposable income.(b)
Enterprise production: in 2019, 70% of China’s CO2 emissions were the result of industrial production or generative emissions. Among them, the carbon dioxide emissions related to the power industry accounted for more than 40%. In addition, those from the steel industry, which are part of the manufacturing industry, accounted for about 15%. Industrial enterprises exceeding permitted size are not only the main body of the industry but also the main entities participating in the carbon emissions trading market. Therefore, factors related to CO2 emissions from enterprises include the number of industrial enterprises exceeding permitted size, full-time equivalent of R&D personnel in industrial enterprises exceeding designated size, industrial structure, and energy structure, energy consumption per GDP unit, energy consumption per unit of industrial added value, and electricity consumption per GDP unit.(c)
Regional environment: in addition to factors related to resident and enterprises, regional environment also affects CO2 emissions and decomposition. We considered three aspects, regional economic level, technological level, and green resources, mainly including total freight volume, total passenger volume, degree of openness to the outside world, total gas supply, total supply of liquefied petroleum gas, urbanization rate, technological level, green area, and forest carbon sink. Specific explanations of features are given in Table 2.Table 2
Specific explanation of features.
FeatureExplanationGDP per capitaGDP/populationEnergy consumption per capitaEnergy consumption/populationDegree of opening to the outside worldTotal import and export/GDPCarbon emissions per capitaTotalCO2emissions/populationCumulative carbon emissions per capitaCumulativeCO2emissions/Cumulative populationTotal passenger volumeThe number of passengers actually carried by the means of transport in a certain period of timeDisposable income per capitaDisposable income/populationIndustrial structureAdded value of tertiary industry/GDPEnergy consumption per unit of GDPEnergy consumption/GDPEnergy consumption per unit of industrial added valueEnergy consumption/industrial added valueElectricity consumption per unit of GDPElectricity consumption/GDPEnergy structureCoal-fired energy/energy consumptionTotal freight volumeIn a certain period of time, the actual weight of the cargo carried by the means of transportThe number of industrial enterprises above designated sizeIncluding independent auditor industrial enterprises and affiliated industrial production unitsFull-time equivalent of R&D personnel in industrial enterprises above designated sizeThe number of full-time staff plus part-time staff is converted to the total number of full-time staff based on workloadTechnological levelNumber of patents granted/10,000 peopleForest carbon sinkForest area/total land areaGreen areaIncluding park area, production green area, protective green area, subsidiary green area, and other green areaTotal gas supplyRefers to the amount of natural gas supplied by gas units during the reporting periodTotal supply of liquefied petroleum gasRefers to the amount of liquefied petroleum gas supplied by gas units during the reporting periodUrbanization rateUrban population/population


## 3.2. Empirical Model

### 3.2.1. Specifications of the Backpropagation Neural Network Model (BP)

Compared with the subjective weight assignment methods (AHP, Expert Evaluation Method, TOPSIS, etc.) used in the formation of the “top-down” model and reported in previous literature, BP involves a multilayer feedforward network and error direction propagation-learning algorithm. Because of its unique adaptability, learning ability, and strong generalization ability, it is widely used in the fields of automatic identification, predictive estimation, engineering, biology, and medicine, among others. For the purpose of the present research, BP can more objectively and accurately quantify the impact of features on quotas and dynamically reflect the nonlinear impact of features on quotas at different stages [36]. Therefore, after training the BP based on pilot regional quotas and feature system data, we calculated the nonpilot regional quotas.BP is composed of an input layer, a hidden layer, and an output layer. These three basic elements are fully connected during the whole network training. For example, for a neural network model with only one hidden layer, the process of BP neural network is mainly divided into two stages.The first stage is the forward propagation of the signal, from the input layer to the hidden layer and finally to the output layer. Assuming that the number of samples isA, input layer has m nodes, output layer has n nodes, and hidden layer has p nodes; xai is the input/output of the input layer, a = 1, 2, …, A, i = 1, 2, …, m; Baj, baj is the input/output of the hidden layer, respectively, j = 1, 2, …, p; Yak, and yak are the input/output of the output layer, respectively, k = 1, 2, …, n; yak∗ is the expected label of the output layer (historical quotas); wij and wjk are the weights from input to hidden and hidden to output, respectively; ，θj and γk are the biases from input to hidden and hidden to output, respectively.Input layer to hidden layer: determine input function,Baj=xaiwij+θj, and then transform Baj into baj through activation function baj=fBaj.Hidden layer to output layer: determine input function,Yak=bajwjk+γk, and then transform Yak into yak through activation function, yak=fYak, where yak is the final result.Later, the loss function is determined and loss is calculated according toyak and yak∗. When the loss is either smaller than the set range or reaches the upper limit of the number of iterations, the model ends the training; otherwise, it enters a second stage.The second stage is the backpropagation of the loss. The loss information is returned along the original propagation path through the learning signal. Starting from the last layer, the weight and bias are corrected layer by layer, and finally the loss is within the set range.Weight and bias update formula:(3)Wl∗=Wl+ΔWl,θl∗=θl+Δθl,ΔWl=−η∂Loss∂Wl,Δθl=−η∂Loss∂θl,where Wl∗ and θl∗ indicate the updated weight and bias of the l layer, Wl and θl indicate the original weight and bias of the l layer, ΔWl and Δθl represent the correction part of weight and bias of the l layer, and η is a fixed value that indicates the learning rate.The fundamental part of backpropagation is to minimize the loss through the update of weights and biases, using the gradient descent method (actually using the chain partial derivative). The specific derivation process is given below.The following assumptions are considered:The activation function from the input layer to the hidden layer isfx=1/1+e−x (sigmoid).The activation function from the input layer to the hidden layer isfx=x.Loss function isLoss=1/2×∑k=1nyak∗−yak2=1/2×∑k=1neak2, eak=yak∗−yak.The weights update process of the hidden-output layer:(4)∂Loss∂wjk=∂Loss∂eak×∂eak∂yak×∂yak∂Yak×∂Yak∂wjk=eak×−1×1×baj=−bajeak,wjk∗=wjk−η∂Loss∂wjk=wjk+ηbajeak.The weight update process of the input-hidden layer:(5)∂Loss∂wij=∂Loss∂eak×∂eak∂yak×∂yak∂Yak×∂Yak∂baj×∂baj∂Bij×∂Bij∂wij=eak×−1×wjk×baj1−baj×xai=−baj1−bajwjkxaieak,wij∗=wij−η∂Loss∂wij=wij+ηbaj1−bajeakwjkxai.In the same way, the updated value of the bias can be obtained:Hidden-output layer:γk∗=γk+ηeakInput-hidden layer:θj∗=θj+ηbaj1−bajwjkeak


### 3.2.2. Key Parameters

According to the model principle previously mentioned, there are 5 key parameters that determine the learning effect of the BP: activation function, loss function, learning rate, the number of hidden layer nodes, and gradient descent algorithm.(1)
Activation function:For the BP both, the hidden layer and the output layer need to use an activation function.For the hidden layer, the activation function is generally a nonlinear function. The reason for this is that, if the activation function is a linear function, the output is a linear combination of the input, which is equivalent to the effect of the no hidden layer (the hidden layer is invalid). The introduction of a nonlinear function as the hidden layer activation function makes the network more powerful, increases its ability to learn complex data, and reflects the nonlinear relationship between input and output. Therefore, we introduced four nonlinear activation functions that are widely used (Tables3 and 4).Although ReLU has two problems, it is currently the most commonly used activation function for BP. In addition, it is the default activation function used by most feedforward neural networks.For the output layer, the choice of its activation function depends on whether the problem is a regression problem or a classification problem. In the event, it is a classification problem, the sigmoid activation function represents a good choice; for regression problems, a linear activation function is more appropriate.(2)
Loss function:With regards to the problem and the output layer activation function to match different loss functions,(1)
Cross-entropy function: it is suitable for binary classification problems, and the output layer activation function is sigmoid(2)
Log-likelihood cost: it is suitable for multiclassification problems, and the output layer activation function is softmax(3)
Mean square error (MSE): it is suitable for regression problems, and the output layer activation function is a linear function(3)
Learning rate:The learning rate value is an important part of the BP, which represents the speed of information accumulation in the neural network over time, and its value is between [0, 1]. Under ideal circumstances, we would start with a large learning rate and gradually reduce the speed until the loss value no longer diverges (if the learning rate is set too low, the training progress will be very slow because only very few adjustments to the weight of the network are made. However, if the learning rate is set too high, it may bring undesirable consequences on the loss function (Figure1)).(4)
The number of hidden layer nodes:The number of hidden layer nodes has a great influence on the prediction accuracy of BP; if the number of nodes is too small, the network cannot perform a proper learning process, and it will need more times to train. In addition, the training accuracy is also affected. When the number of nodes is too large, the training time increases and the network will result in overfitting. However, there is no conventional formula for determining the number of nodes. Some empirical formulas are given below for reference:l<n−1l<m+n+a, l=log2n,l≥k×n/n+m, where n indicates the number of input layer nodes, l indicates the number of hidden layer nodes, m is the number of output layer nodes, a represents a constant between 0–10, and k corresponds to the number of samples.In fact, the number of hidden layer nodes can be roughly calculated according to the reference formula. Later, trial and error is used to find the optimal number of nodes. Generally speaking, the BP error shows a trend where it first decreases and later increases with the increase of hidden layer nodes.(5)
Gradient descent algorithm:We introduce six well-known gradient descent algorithms (Tables5 and 6).According to this analysis, there are no perfect key parameters that can suit all conditions. The appropriate selection of key parameters depends on the specific problem of study. In the present research, we studied a regression problem. Thus, we chose ReLU andfx=x as the activation function and MSE as the loss function. Also, Adam may be appropriate as the gradient descent algorithm; however, the learning rate and the number of hidden layer nodes cannot be determined in advance. In summary, the final determination of all key parameters needs BP training. For this reason, the determination of key parameters is provided in Section 4.Table 3
Activation functions and figures.
NameFunctionDerivative functionSigmoidSx=1/1+e−xS′x=Sx1−SxSoftmaxSoftmaxZi=eZi/∑c=1CeZc=pi∂softmaxZi/∂Zj=pi1−pi,j=i,−pj⋅pi,j≠i.tanhtanhx=ex−e−x/ex+e−xT′x=1−Tx2ReLUReLUx=max0,xReLU′x=0,x≤0,1,x>0.Table 4
Advantages and disadvantages of activation functions.
NameAdvantagesDisadvantagesSigmoid(1) It can smoothly map the real-number field to [0, 1](1) Gradient vanishing, that is, in the process of backpropagation, the derivative will gradually become 0; thus, the parameters cannot be updated and the neural network cannot be optimized(2) Monotonically increasing, continuous derivable, and its derivative form is very simple(2) Nonzero-centered: the output value of the function is always greater than 0, which will slow down the convergence speed of the model training(3) Suitable for handling binary classification problems(3) Exponentiation is relatively time-consumingSoftmax(1) It maps the output value to (0, 1), and the sum of the mapped output value is 1(1) The operation of Softmax involves the calculation of exponential function; in consequence, an “overflow problem” for computers occurs(2) It divides the entire hyperspace according to the number of classifications(3) Suitable for multiclassification problems(2) Not suitable for face recognition taskstanh(1) It can smoothly map the real-number field to [−1, 1](1) Gradient vanishing, that is, in the process of backpropagation, the derivative will gradually become 0; thus, the parameters cannot be updated and the neural network cannot be optimized(2) Solve nonzero-centered problem(3) Suitable for handling binary classification problemsReLU(1) Solve the gradient vanishing in the positive interval(1) Nonzero-centered(2) The calculation is simple, no exponential calculation is required, and the activation value can be obtained with only one value(2) Dead ReLU problem: it is “vulnerable” during training; whenx < 0, the gradient is 0; the gradient of these nodes and subsequent nodes are always 0, and no longer responds to any data, causing the corresponding parameters to never be updated(3) The convergence speed is much faster than sigmoid and tanhFigure 1
The loss function results of different learning rates.
(a)(b)(c)Table 5
Gradient descent functions and parameter description.
NameGradient descent functionParameter descriptionBGDθ∗=θ−η⋅g⌢=θ−η⋅∇θJθg⌢ denotes gradient estimateθ denotes weightSGDθ∗=θ−η⋅g⌢=θ−η⋅∇θJθ;xi;yig⌢ denotes gradient estimateθ denotes weightMomentumvt=γvt−1+η⋅∇θJθvt denotes gradient estimateθ∗=θ−vtθ denotes weightAdaGradg⌢t,i=∇θJθig⌢t,i denotes gradient estimateθt+1,i∗=θt,i−η/g⌢t,i⋅g⌢t,i+ε⋅g⌢t,iθt,i denotes weightRMSpropgt=∇θJθt−1γ=0.9Eg2t denotes Exponential weighted averageEg2t=γEg2t−1+1−γgt2gt denotes gradient estimateθt+1∗=θt−η/Eg2t+εgtθt denotes weightAdamgt=∇θJθt−1β1=0.9β2=0.999ε=10e−8mt=β1mt−1+1−β1gtgt denotes gradient estimate θt denotes weightvt=β2vt−1+1−β2gt2m⌢t=mt/1−β1t, v⌢t=vt/1−β2tθt+1∗=θt−η/v⌢t+εm⌢tTable 6
Advantages and disadvantages of gradient descent algorithms.
NameAdvantagesDisadvantagesBGDThe principle of gradient descent is simple(1) Calculation is very slow(2) Difficult to handle a large dataset(3) Cannot add new data to update the modelSGD(1) Compared with BGD, SGD training speed is faster(1) Frequent updates may cause severe oscillations in the loss function(2) New data can be added to update the modelMomentum(1) Consider the speed of the previous step and the new gradient(2) Can speed up the convergence and suppress the shockAdaGrad(1) Compared with SGD, it adds a denominator(1) If gradient update is frequent, it may cause the subsequent gradient updates be slow or disappear(2) Handle the case where the number of gradient updates is smallRMSprop(1) Similar to momentum, it can reduce fluctuations(2) Overcome the problem of the sharp decrease or disappearance of gradient in AdaGrad(3) It performs better than SGD, momentum, and AdaGrad, based on the nonstationary objective functionAdam(1) Combine momentum and RMSProp(2) Integrate the contents of gradient descent, momentum, Adagrad, and RMSprop with certain improvement(3) Easy to use, insensitive to gradient scaling, can be used for large data, processing sparse data, easy to adjust hyperparameters, etc.


## 3.2.1. Specifications of the Backpropagation Neural Network Model (BP)

Compared with the subjective weight assignment methods (AHP, Expert Evaluation Method, TOPSIS, etc.) used in the formation of the “top-down” model and reported in previous literature, BP involves a multilayer feedforward network and error direction propagation-learning algorithm. Because of its unique adaptability, learning ability, and strong generalization ability, it is widely used in the fields of automatic identification, predictive estimation, engineering, biology, and medicine, among others. For the purpose of the present research, BP can more objectively and accurately quantify the impact of features on quotas and dynamically reflect the nonlinear impact of features on quotas at different stages [36]. Therefore, after training the BP based on pilot regional quotas and feature system data, we calculated the nonpilot regional quotas.BP is composed of an input layer, a hidden layer, and an output layer. These three basic elements are fully connected during the whole network training. For example, for a neural network model with only one hidden layer, the process of BP neural network is mainly divided into two stages.The first stage is the forward propagation of the signal, from the input layer to the hidden layer and finally to the output layer. Assuming that the number of samples isA, input layer has m nodes, output layer has n nodes, and hidden layer has p nodes; xai is the input/output of the input layer, a = 1, 2, …, A, i = 1, 2, …, m; Baj, baj is the input/output of the hidden layer, respectively, j = 1, 2, …, p; Yak, and yak are the input/output of the output layer, respectively, k = 1, 2, …, n; yak∗ is the expected label of the output layer (historical quotas); wij and wjk are the weights from input to hidden and hidden to output, respectively; ，θj and γk are the biases from input to hidden and hidden to output, respectively.Input layer to hidden layer: determine input function,Baj=xaiwij+θj, and then transform Baj into baj through activation function baj=fBaj.Hidden layer to output layer: determine input function,Yak=bajwjk+γk, and then transform Yak into yak through activation function, yak=fYak, where yak is the final result.Later, the loss function is determined and loss is calculated according toyak and yak∗. When the loss is either smaller than the set range or reaches the upper limit of the number of iterations, the model ends the training; otherwise, it enters a second stage.The second stage is the backpropagation of the loss. The loss information is returned along the original propagation path through the learning signal. Starting from the last layer, the weight and bias are corrected layer by layer, and finally the loss is within the set range.Weight and bias update formula:(3)Wl∗=Wl+ΔWl,θl∗=θl+Δθl,ΔWl=−η∂Loss∂Wl,Δθl=−η∂Loss∂θl,where Wl∗ and θl∗ indicate the updated weight and bias of the l layer, Wl and θl indicate the original weight and bias of the l layer, ΔWl and Δθl represent the correction part of weight and bias of the l layer, and η is a fixed value that indicates the learning rate.The fundamental part of backpropagation is to minimize the loss through the update of weights and biases, using the gradient descent method (actually using the chain partial derivative). The specific derivation process is given below.The following assumptions are considered:The activation function from the input layer to the hidden layer isfx=1/1+e−x (sigmoid).The activation function from the input layer to the hidden layer isfx=x.Loss function isLoss=1/2×∑k=1nyak∗−yak2=1/2×∑k=1neak2, eak=yak∗−yak.The weights update process of the hidden-output layer:(4)∂Loss∂wjk=∂Loss∂eak×∂eak∂yak×∂yak∂Yak×∂Yak∂wjk=eak×−1×1×baj=−bajeak,wjk∗=wjk−η∂Loss∂wjk=wjk+ηbajeak.The weight update process of the input-hidden layer:(5)∂Loss∂wij=∂Loss∂eak×∂eak∂yak×∂yak∂Yak×∂Yak∂baj×∂baj∂Bij×∂Bij∂wij=eak×−1×wjk×baj1−baj×xai=−baj1−bajwjkxaieak,wij∗=wij−η∂Loss∂wij=wij+ηbaj1−bajeakwjkxai.In the same way, the updated value of the bias can be obtained:Hidden-output layer:γk∗=γk+ηeakInput-hidden layer:θj∗=θj+ηbaj1−bajwjkeak


## 3.2.2. Key Parameters

According to the model principle previously mentioned, there are 5 key parameters that determine the learning effect of the BP: activation function, loss function, learning rate, the number of hidden layer nodes, and gradient descent algorithm.(1)
Activation function:For the BP both, the hidden layer and the output layer need to use an activation function.For the hidden layer, the activation function is generally a nonlinear function. The reason for this is that, if the activation function is a linear function, the output is a linear combination of the input, which is equivalent to the effect of the no hidden layer (the hidden layer is invalid). The introduction of a nonlinear function as the hidden layer activation function makes the network more powerful, increases its ability to learn complex data, and reflects the nonlinear relationship between input and output. Therefore, we introduced four nonlinear activation functions that are widely used (Tables3 and 4).Although ReLU has two problems, it is currently the most commonly used activation function for BP. In addition, it is the default activation function used by most feedforward neural networks.For the output layer, the choice of its activation function depends on whether the problem is a regression problem or a classification problem. In the event, it is a classification problem, the sigmoid activation function represents a good choice; for regression problems, a linear activation function is more appropriate.(2)
Loss function:With regards to the problem and the output layer activation function to match different loss functions,(1)
Cross-entropy function: it is suitable for binary classification problems, and the output layer activation function is sigmoid(2)
Log-likelihood cost: it is suitable for multiclassification problems, and the output layer activation function is softmax(3)
Mean square error (MSE): it is suitable for regression problems, and the output layer activation function is a linear function(3)
Learning rate:The learning rate value is an important part of the BP, which represents the speed of information accumulation in the neural network over time, and its value is between [0, 1]. Under ideal circumstances, we would start with a large learning rate and gradually reduce the speed until the loss value no longer diverges (if the learning rate is set too low, the training progress will be very slow because only very few adjustments to the weight of the network are made. However, if the learning rate is set too high, it may bring undesirable consequences on the loss function (Figure1)).(4)
The number of hidden layer nodes:The number of hidden layer nodes has a great influence on the prediction accuracy of BP; if the number of nodes is too small, the network cannot perform a proper learning process, and it will need more times to train. In addition, the training accuracy is also affected. When the number of nodes is too large, the training time increases and the network will result in overfitting. However, there is no conventional formula for determining the number of nodes. Some empirical formulas are given below for reference:l<n−1l<m+n+a, l=log2n,l≥k×n/n+m, where n indicates the number of input layer nodes, l indicates the number of hidden layer nodes, m is the number of output layer nodes, a represents a constant between 0–10, and k corresponds to the number of samples.In fact, the number of hidden layer nodes can be roughly calculated according to the reference formula. Later, trial and error is used to find the optimal number of nodes. Generally speaking, the BP error shows a trend where it first decreases and later increases with the increase of hidden layer nodes.(5)
Gradient descent algorithm:We introduce six well-known gradient descent algorithms (Tables5 and 6).According to this analysis, there are no perfect key parameters that can suit all conditions. The appropriate selection of key parameters depends on the specific problem of study. In the present research, we studied a regression problem. Thus, we chose ReLU andfx=x as the activation function and MSE as the loss function. Also, Adam may be appropriate as the gradient descent algorithm; however, the learning rate and the number of hidden layer nodes cannot be determined in advance. In summary, the final determination of all key parameters needs BP training. For this reason, the determination of key parameters is provided in Section 4.Table 3
Activation functions and figures.
NameFunctionDerivative functionSigmoidSx=1/1+e−xS′x=Sx1−SxSoftmaxSoftmaxZi=eZi/∑c=1CeZc=pi∂softmaxZi/∂Zj=pi1−pi,j=i,−pj⋅pi,j≠i.tanhtanhx=ex−e−x/ex+e−xT′x=1−Tx2ReLUReLUx=max0,xReLU′x=0,x≤0,1,x>0.Table 4
Advantages and disadvantages of activation functions.
NameAdvantagesDisadvantagesSigmoid(1) It can smoothly map the real-number field to [0, 1](1) Gradient vanishing, that is, in the process of backpropagation, the derivative will gradually become 0; thus, the parameters cannot be updated and the neural network cannot be optimized(2) Monotonically increasing, continuous derivable, and its derivative form is very simple(2) Nonzero-centered: the output value of the function is always greater than 0, which will slow down the convergence speed of the model training(3) Suitable for handling binary classification problems(3) Exponentiation is relatively time-consumingSoftmax(1) It maps the output value to (0, 1), and the sum of the mapped output value is 1(1) The operation of Softmax involves the calculation of exponential function; in consequence, an “overflow problem” for computers occurs(2) It divides the entire hyperspace according to the number of classifications(3) Suitable for multiclassification problems(2) Not suitable for face recognition taskstanh(1) It can smoothly map the real-number field to [−1, 1](1) Gradient vanishing, that is, in the process of backpropagation, the derivative will gradually become 0; thus, the parameters cannot be updated and the neural network cannot be optimized(2) Solve nonzero-centered problem(3) Suitable for handling binary classification problemsReLU(1) Solve the gradient vanishing in the positive interval(1) Nonzero-centered(2) The calculation is simple, no exponential calculation is required, and the activation value can be obtained with only one value(2) Dead ReLU problem: it is “vulnerable” during training; whenx < 0, the gradient is 0; the gradient of these nodes and subsequent nodes are always 0, and no longer responds to any data, causing the corresponding parameters to never be updated(3) The convergence speed is much faster than sigmoid and tanhFigure 1
The loss function results of different learning rates.
(a)(b)(c)Table 5
Gradient descent functions and parameter description.
NameGradient descent functionParameter descriptionBGDθ∗=θ−η⋅g⌢=θ−η⋅∇θJθg⌢ denotes gradient estimateθ denotes weightSGDθ∗=θ−η⋅g⌢=θ−η⋅∇θJθ;xi;yig⌢ denotes gradient estimateθ denotes weightMomentumvt=γvt−1+η⋅∇θJθvt denotes gradient estimateθ∗=θ−vtθ denotes weightAdaGradg⌢t,i=∇θJθig⌢t,i denotes gradient estimateθt+1,i∗=θt,i−η/g⌢t,i⋅g⌢t,i+ε⋅g⌢t,iθt,i denotes weightRMSpropgt=∇θJθt−1γ=0.9Eg2t denotes Exponential weighted averageEg2t=γEg2t−1+1−γgt2gt denotes gradient estimateθt+1∗=θt−η/Eg2t+εgtθt denotes weightAdamgt=∇θJθt−1β1=0.9β2=0.999ε=10e−8mt=β1mt−1+1−β1gtgt denotes gradient estimate θt denotes weightvt=β2vt−1+1−β2gt2m⌢t=mt/1−β1t, v⌢t=vt/1−β2tθt+1∗=θt−η/v⌢t+εm⌢tTable 6
Advantages and disadvantages of gradient descent algorithms.
NameAdvantagesDisadvantagesBGDThe principle of gradient descent is simple(1) Calculation is very slow(2) Difficult to handle a large dataset(3) Cannot add new data to update the modelSGD(1) Compared with BGD, SGD training speed is faster(1) Frequent updates may cause severe oscillations in the loss function(2) New data can be added to update the modelMomentum(1) Consider the speed of the previous step and the new gradient(2) Can speed up the convergence and suppress the shockAdaGrad(1) Compared with SGD, it adds a denominator(1) If gradient update is frequent, it may cause the subsequent gradient updates be slow or disappear(2) Handle the case where the number of gradient updates is smallRMSprop(1) Similar to momentum, it can reduce fluctuations(2) Overcome the problem of the sharp decrease or disappearance of gradient in AdaGrad(3) It performs better than SGD, momentum, and AdaGrad, based on the nonstationary objective functionAdam(1) Combine momentum and RMSProp(2) Integrate the contents of gradient descent, momentum, Adagrad, and RMSprop with certain improvement(3) Easy to use, insensitive to gradient scaling, can be used for large data, processing sparse data, easy to adjust hyperparameters, etc.


## 4. Empirical Results and Analysis

We focused on the allocation of regional carbon quotas, which is a regression problem, and our goal was to minimize the loss rate while ensuring that the test set results meet the realistic expected range in China. Therefore, based on experience and BP training, we chose MSE as the loss function and Adam [37] as the gradient descent algorithm in the backpropagation process, the final learning rate was 0.009, and number of iterations were 5000. Other parameters are shown in Table 7.Table 7
Parameters of the benchmark model.
LayerFunctionActivationImport-hidden(8, 6)ReLUHidden-export(6, 1)fx=xThe total loss rate of the benchmark model was 0.02419. And the comparison between the results of the training set and the historical quotas are shown in Table8 (benchmark model in Table 8). The test results are shown in Table 9 (benchmark model in Table 9). After adding up the historical quotas in pilot regions and the estimated quotas in the nonpilot regions, national quotas were obtained. The resulting national quotas are displayed in Table 10 (benchmark model in Table 10).Table 8
Comparison between model results and historical quotas.
YearRegionHistorical quotasBenchmark modelNew modelTraining resultsLossTraining resultsLoss2014Beijing0.50.51140.01140.5020.002Tianjing1.61.73970.13971.584−0.016Shanghai1.51.77950.27951.4785−0.0215Hubei3.243.0159−0.22413.24850.0085Guangdong3.883.6819−0.19813.8754−0.0046Shenzhen0.330.36230.03230.3169−0.0131Chongqing1.31.32830.02831.2949−0.00512015Beijing0.50.54410.04410.51440.0144Tianjing1.61.73120.13121.61690.0169Shanghai1.61.67280.07281.62140.0214Hubei3.243.0162−0.22383.1589−0.0811Guangdong4.083.8977−0.18234.0595−0.0205Shenzhen0.330.36230.03230.3182−0.0118Chongqing1.251.4190.1691.30370.05372016Beijing0.50.4536−0.04640.4882−0.0118Tianjing1.61.5162−0.08381.5506−0.0494Shanghai1.51.4923−0.00771.50990.0099Hubei2.82.93840.13842.87150.0715Guangdong44.05520.05524.01410.0141Shenzhen0.30.36360.06360.32090.0209Chongqing1.31.41760.11761.2564−0.04362017Beijing0.50.3623−0.13770.4946−0.0054Tianjing1.61.363−0.2371.63990.0399Shanghai1.61.2582−0.34181.5741−0.0259Hubei2.52.58910.08912.4883−0.0117Guangdong4.24.50110.30114.2080.008Shenzhen0.30.36230.06230.30430.0043Chongqing1.31.2166−0.08341.2883−0.0117Unit of quotas: 100 million tons.Table 9
Test quotas of nonpilot regions.
RegionBenchmark modelNew model20142015201620172014201520162017Hebei4.47424.72134.66393.79584.24894.25123.93893.4056Shanxi8.78749.55149.91797.93216.75757.26397.48576.0955Inner Mongolia5.45415.31785.83676.05775.41355.26435.31965.8676Liaoning5.70355.21966.7685.49175.95955.50835.81795.1618Jilin3.32953.19443.22862.80114.57324.6514.51394.3528Heilongjiang4.56584.76665.31264.67115.11615.30885.56125.4386Jiangsu4.90334.93794.92764.50713.48513.69383.70963.6284Zhejiang4.38964.32164.05973.93794.21724.32144.01434.036Anhui3.62313.6733.58413.17433.644.08333.82083.1997Fujian2.91492.94272.79862.59973.5583.56713.34283.4086Jiangxi2.15242.31862.4082.14463.53343.55353.40323.2111Shandong7.087.07777.29476.7665.64715.63225.41545.0304Henan3.60613.78023.83083.22013.15823.26543.36322.937Hunan2.16932.30322.36782.18152.69932.80952.8942.9341Guangxi2.07991.9912.0832.03833.26283.15193.04733.0321Hainan7.0927.77617.76545.80296.17146.56586.61435.6134Sichuan3.03562.95892.84582.5231.96392.10272.03261.9373Guizhou3.39243.28823.59772.48393.47293.40673.44042.8913Yunnan1.18690.99861.19441.10382.6982.65552.72582.5445Shaanxi4.47164.49874.56983.71512.75072.91262.6892.6752Gansu3.99624.29834.3573.23072.30422.61872.69272.0599Qinghai0.36230.36230.87580.36232.07441.99192.14662.0002Ningxia8.99519.00658.8929.03854.96345.10765.04975.1049Xinjiang8.2868.20769.16256.55783.68023.94034.2363.8333Unit of quotas: 100 million tons.Table 10
National quota.
National quotaBenchmark model2014201520162017118.4012120.1122124.3424108.137New model2014201520162017107.6989110.2274109.2749102.3993Unit of quotas: 100 million tons.The results of the benchmark model indicated that, during the initial stage of China’s unified carbon emission market (3-4 years), the national quota will increase and, after an adaptation period, China’s total quota and regional quotas will begin to decrease. This will stimulate enterprises to accelerate emission reduction and prove China’s determination to achieve carbon peaks before 2030 and carbon neutrality by 2060.


## 5. Further Analysis: Build a Comprehensive Feature System

It is unreasonable to allocate regional quotas only considering the factor of corporate CO2 emissions. In addition to corporate production factors, regional CO2 emissions should also consider regional human activities and the role of forests in reducing those emissions. Thus, we believe that, in addition to the CO2 emissions reported by enterprises, quotas in China’s pilot regions should also take into account other features such as forest carbon sinks, population, and natural endowments. Based on the comprehensive factors of these three aspects, the regional quotas were determined. Therefore, we added other factors related to people and regions (see Section 3.1.2, for details) into the feature system. Subsequently, we chose MSE as the loss function, Adam as the gradient descent algorithm, and trained BP to obtain the final learning rate (0.003), number of iterations (5000), and other parameters (Table 11).Table 11
Parameters of the new model.
LayerFunctionActivationImport-hidden(21, 19)ReLUHidden-export(19, 1)fx=xThe total loss rate of the new model was 0.00089, and the results of the comparison between the training set and historical quotas are shown in Table8 (new model in Table 8). The test results are displayed in Table 9 (new model in Table 9). After adding up the historical quotas in pilot regions and the estimated quotas in the nonpilot regions, the national quota was obtained (new model in Table 10).While using the comprehensive feature system, the new model displayed a lower loss rate than that obtained with the benchmark model, and the calculated national quota was closer to the CO2 emissions reported by China. These results indicated that the feature system has a certain degree of rationality and accuracy. Similarly, the calculation results of the model presented a trend, where the amount of national carbon quotas initially increased and later began to decrease.


## 6. Conclusions and Future Work

The whole world is expecting China to lead the economic recovery and green development after the global epidemic. It also expects China’s 14th Five-Year Plan to become the Guide for green recovery. In the same year, China established a unified national carbon emissions market. This represents not only China’s further exploration of the carbon emissions trading system to achieve green development but also one of the important tools for China to achieve two low-carbon goals. In addition, quota allocation is an important factor that determines the functionality of Chinese carbon market. In order to calculate other regional quotas, we trained a BP benchmark model. For this purpose, we considered historical quota data of China’s 7-carbon emissions trading market pilot regions from 2014 to 2017 and selected suitable features that fit China’s “bottom-up” total control method. Later, we built a feature system that included human, corporate, and regional factors, retrained the model, and recalculated quotas for other regions. The results are presented herein. First, both, the benchmark model results and the results obtained using the comprehensive feature system showed that within the sample interval, the amount of China’s national carbon quotas displayed an initial increase to later decrease. Second, the model trained with the characteristic data of the feature system built in the present research displayed a lower loss rate as compared with the benchmark model. These results demonstrated that the feature system proposed in this paper fits not only the actual situation of China’s CO2 emissions and quotas but also that the framework of the system is reasonable and accurate. Third, the feature system and training model proposed in the present article combined with the original “bottom-up” total control and post adjustment method can be used by Chinese CO2 emission decision makers to obtain advanced predictions. We have provided the content of China’s carbon emissions trading quota system, which can promote the operation of China’s carbon emissions market, encourage participants in market transactions to reduce emissions, and accelerate China’s low-carbon development.Of course, we also admit that, in the future, the feature system and model proposed in this article can be further improved and perfected as follows. First, the indicators related to enterprises in the currently constructed feature system are substitute indicators because the specific transaction data and enterprise-related data of China’s carbon emission market have not been unified and officially announced. Therefore, once the data is available, this part of the indicators will increase or decrease. Second, at present, China’s national unified carbon emissions trading market has just started, and the main participants are enterprises, with less individual participation. At the same time, the central government has not issued a policy about people’s low-carbon life. Therefore, when China’s emission reduction program enters the critical stage in the future, the features related to people will increase or decrease. Third, the current BP neural network model has only three layers. In the future, with the improvement of feature data, a certain number of hidden layers may be further increased to train a model with low loss rate and stronger generalization ability.



---

*Source: 6659302-2021-01-06.xml*