福建省2030年碳达峰前二氧化碳排放趋势研究——基于GM(1,1)、GM(2,1)与GM(1,1)邓聚龙灰色预测模型

Study of CO2 Emission Trends in Fujian Province before Carbon Peaking in 2030—Based on GM(1,1), GM(2,1) and GM(1,1) Deng Jurong Grey Forecast Models

摘要: 本文采用GM(1,1)、GM(2,1)和GM(1,1)邓聚龙灰色预测模型，选取了福建省1995年至2019年的二氧化碳排放量，用以拟合分析我国提出在2030年达到碳达峰前的福建省二氧化碳排放量数据增长趋势，实证分析三种灰色预测模型之间的差异性与一致性，得出GM(2,1)与GM(1,1)邓聚龙模型此次拟合效果具有一致性，均可以采用。并且二氧化碳排放预测结果表明：福建省从2019年开始至2030年，碳排放量呈现一种加速上升的趋势。本文也以此提出我国要在2060年达到碳中和目标的一些举措。

Abstract: This paper, the GM(1,1), GM(2,1) and GM(1,1) Deng Jurong gray prediction models are used to fit the CO2 emissions in Fujian Province from 1995 to 2019 to analyze the growth trend of CO2 emissions data in Fujian Province before China proposes to reach the carbon peak in 2030, and the differences and consistency between the three gray prediction models are empirically analyzed, and it is concluded that the GM(2,1) and GM(1,1) Deng Jurong models are consistent and can be used, and the prediction results of CO2 emissions show that the carbon emissions in Fujian Province will show an accelerated upward trend from 2019 to 2030. This paper also uses this to propose some initiatives for China to reach the carbon neutrality target in 2060.

1. 引言

伴随着社会经济的不断发展，世界各国逐渐将“环保”这一主题提上议程，这意味着各国不再追求单纯的经济发展，而是寻求环境保护与经济发展的平衡 [1]。气候变化是全人类的问题，中国作为GDP总量排名前几名的大国，更加具有责任和担当。在近年两会上，我国将碳达峰、碳中和写入了政府工作报告，并承诺在2030年碳达峰，这标志着一些高能耗、碳排放超标的产业可能需要寻求转型和技术升级，未来的经济驱动不能再仅仅依靠能源消耗来拉动 [2]。对于各级政府来说，工作的重心也要慢慢调整，向降低碳排放的目标努力。对于个人而言，行为要与国家发展的整体趋势相一致，这不仅是政府工作目标实现的保证，更与个人的生存和发展质量息息相关。要实现“碳达峰”和“碳中和”的目标，仅仅依靠经济发展的自然趋势是远远不够的，必须要结合相应的政策和规章制度 [3]。为了体现这一点，明确政策和规则制度调整的重要性，本文利用了GM(1,1)，GM(2,1)和GM(1,1)邓聚龙模型对福建省未来十年的碳排放量进行了预测 [4]，希望能够模拟在没有政策约束的情况下碳排放量增长的趋势，以证明为了实现控制碳排放量的目标，我国政府和人民还需要付出更多的努力。

2. 原理介绍

2.1. GM(1,1)模型理论基础

灰色系统模型(Grey System Model, GM)，实际上是利用数学公式中的指数函数构造的模型。GM(1,1)是此类模型中目前最简单也是应用范围最广的 [5]，其基础是解常微分方程，从而得到方程的解析解。

在这里我们假设有一个变量x(t)，其是关于时间变量t的函数，它满足一阶常微分方程条件如下：

$\frac{\text{d}x\left(t\right)}{\text{d}t}+ax\left(t\right)=b$ (1)

需要补充说明的是这里，参数a和参数b是两个常系数，并非是关于某个变量的函数，该微分方程也可以变为如下形式：

$\frac{\text{d}x\left(t\right)}{\text{d}t}+ax\left(t\right)-b=0$ (2)

如果我们假设参数C为任意常数，不限定其取值范围，那么，该微分方程的解析解(又称通解)为：

$x\left(0\right)=\frac{b+{\text{e}}^{aC}}{a}$ (3)

如果我们假设 $x\left(t\right)$ 具有一个初值，则易得知，在 $t=0$ 时，初值则为 $x\left(0\right)$，那么，当函数取值为 $x\left(0\right)$ 也需满足这个微分方程的通解：

$\frac{\text{d}x\left(t\right)}{\text{d}t}+ax\left(t\right)=b$ (4)

因此，常数参数C满足关系式：

${\text{e}}^{aC}=ax\left(0\right)-b$ (5)

代入微分方程的通解，就有：

$x\left(t\right)=\frac{a}{b}+\left(x\left(0\right)-\frac{b}{a}\right){\text{e}}^{-at}$ (6)

该函数具有指数函数的特征，由于其对于参数a和b的依赖性很强，而参数a和b又是未知的，且在正常情况下不易估计，因此，该模型也被称为灰色模型。由于只有一个变量 $x\left(t\right)$，且是一阶微分，所以记为GM(1,1)。

在实证分析中，我们发现对于最优解初值的选取，对我们最后得到最优解a，b有一定的影响。为了取得稳定的最优解，初值解取随机数较为合适。与GM(2,1)相比，由于GM(1,1)模型只有a，b两个系数，收敛速度较快，比较容易得到最优解a，b的值 [6]。

2.2. GM(2,1)模型理论基础

同GM(1,1)里所述，我们依然假设变量 $x\left(t\right)$ 是关于时间变量t的二阶微分方程的解，则在二阶情况下，需要满足如下条件：

$\frac{{\text{d}}^{2}x\left(t\right)}{\text{d}{t}^2}+a\frac{\text{d}x\left(t\right)}{\text{d}t}+bx\left(t\right)=0$ (7)

通过数学计算，此时容易解得，上述二阶微分方程的通解如下所示：

$x\left(t\right)=C_1{\text{e}}^{\frac{\left(-a-\sqrt{a^2-4b}\right)t}{2}}+C_2{\text{e}}^{\frac{\left(-a-\sqrt{a^2-4b}\right)t}{2}}$ (8)

为了容易在python程序里对其进行编程计算，我们可以对该通解关系式予以代数式的简化，使得：

$x\left(t\right)=C_1{\text{e}}^{mt}+C_2{\text{e}}^{kt}$ (9)

则此时归纳上述方程，可以知道，存在关系式：

$m=\frac{-a-\sqrt{a^2-4b}}{2}$ (10)

$k=\frac{-a+\sqrt{a^2-4b}}{2}$ (11)

$C_1=\frac{x\left(0\right){\text{e}}^k-x\left(1\right)}{{\text{e}}^k-{\text{e}}^m}$ (12)

$C_2=\frac{x\left(1\right)-x\left(0\right){\text{e}}^m}{{\text{e}}^k-{\text{e}}^m}$ (13)

可以注意到此时函数关系式里存在参数C₁，C₂，它的确定我们可以通过两个特解条件来求解，当 $t=0,1$ 时：

$X\left(0\right)=C_1+C_2$ (14)

$X\left(1\right)=C_1{\text{e}}^m+C_2{\text{e}}^k$ (15)

对于该二阶常微分方程通解的验证较为容易，因为：

$\frac{\text{d}x\left(t\right)}{\text{d}t}=C_{1}m{\text{e}}^{mt}+C_{2}k{\text{e}}^{kt}$ (16)

$\frac{{\text{d}}^{2}x\left(t\right)}{\text{d}{t}^2}=C_1{m}^2{\text{e}}^{mt}+C_2{k}^2{\text{e}}^{kt}$ (17)

$m^2+am+b={\left(m+\frac{a}{2}\right)}^2-\frac{a^2-4b}{4}=0$ (18)

$k^2+ak+b={\left(k+\frac{a}{2}\right)}^2-\frac{a^2-4b}{4}=0$ (19)

所以，很容易求证微分方程的关系式为零。由于只有一个变量 $x\left(t\right)$，但是微分方程为二阶，所以记为GM(2,1)。

2.3. GM(1.1)邓聚龙模型理论基础

灰色系统GM(1,1)参数可以通过一个近似矩阵求解的方式来估计。

假设变量 $x\left(t\right)$ 的原始离散序列为 $X^{\left(0\right)}$，是一个非负序列，具有形式：

$X^{\left(0\right)}=\left\{X_0^{\left(0\right)},\cdots ,X_{n-1}^{\left(0\right)}\right\}\left(n\ge 3;X_i^{\left(0\right)}\ge 0;i=0,1,2,\cdots ,n-1\right)$ (20)

假设离散序列 $X^{\left(\text{1}\right)}$ 为 $X^{\left(0\right)}$ 的一次累加序列，具有如下形式：

$X^{\left(1\right)}=\left\{X_0^{\left(1\right)},\cdots ,X_{n-1}^{\left(1\right)}\right\}\left(n\ge 3;X_i^{\left(1\right)}\ge X_{i-1}^{\left(1\right)};i=0,1,2,\cdots ,n-1\right)$ (21)

$X_i^{\left(1\right)}=X_0^{\left(0\right)}+X_1^{\left(0\right)}+X_2^{\left(0\right)}+\cdots +X_i^{\left(0\right)}={\displaystyle {\sum }_{l=0}^i{X}_l^{\left(0\right)}}$ (22)

$X_i^{\left(1\right)}=X_{i-1}^{\left(1\right)}+X_i^{\left(0\right)}$ (23)

$X_0^{\left(1\right)}=X_0^{\left(0\right)}$ (24)

假设离散序列 $Z^{\left(\text{1}\right)}$ 为一个均值序列：

$Z^{\left(1\right)}=\left\{Z_1^{\left(1\right)},\cdots ,Z_{n-1}^{\left(1\right)}\right\}\left(n\ge 3;i=0,1,2,\cdots ,n-1\right)$ (25)

$Z_i^{\left(1\right)}=\frac{1}{2}\left(X_{i-1}^{\left(1\right)}+X_i^{\left(1\right)}\right)$ (26)

那么，关于序列 $X^{\left(0\right)}\left(t\right)$ 的一阶常微分方程求解过程，就被转化成如下形式的一阶微分方程求解过程：

$\frac{\text{d}{x}^{\left(1\right)}\left(t\right)}{\text{d}t}+a{x}^{\left(0\right)}\left(t\right)=b$ (27)

其参数(a, b)估计问题，就等价于一个对于离散序列 $X_k^{\left(1\right)},\left(k=1,2,\cdots ,n-1\right)$ 的GM(1,1)灰色差分模型过程求解问题：

$X_k^{\left(0\right)}+a{Z}_k^{\left(1\right)}=b,\left(k=1,2,\cdots ,n-1\right)$ (28)

其参数(a, b)可以根据最小二乘法来估计：

$\hat{a}=\left(\begin{array}{c}a\\ b\end{array}\right)={\left(B^{\text{T}}B\right)}^{-1}{B}^{\text{T}}{y}_N$ (29)

$B=\left(\begin{array}{cc}-Z_1^{\left(1\right)}& 1\\ ⋮& ⋮\\ -Z_{n-1}^{\left(1\right)}& 1\end{array}\right),y_N=\left(\begin{array}{c}{X}_1^{\left(0\right)}\\ ⋮\\ X_{n-1}^{\left(0\right)}\end{array}\right)$ (30)

相应于等价微分方程，离散序列X⁽¹⁾的估计值解可以由下式定义：

${\hat{X}}_i^{\left(1\right)}=\left(X_1^{\left(0\right)}-\frac{b}{a}\right){\text{e}}^{-a\left(i-1\right)}+\frac{b}{a}$ (31)

${\hat{X}}_i^{\left(0\right)}={\hat{X}}_i^{\left(1\right)}-{\hat{X}}_{i-1}^{\left(1\right)}$ (32)

前面介绍的公式只有一个变量。假设有N个变量( $i=1,2,\cdots ,N$ )，那么，灰色模型就为G(1, N)，由于公式推导较为复杂，这里不再讨论。

3. 模型拟合过程及结果

3.1. 数据来源

根本文计算利用到的二氧化碳排放相关数据均选自历年由福建省统计局公布的《福建统计年鉴》(1995~2020)，历年数据如下表1所示。

Table 1. Carbon dioxide emissions data of Fujian Province over the years

表1. 福建省二氧化碳历年排放数据

3.2. 模型GM(1.1)、GM(2,1)与GM(1,1)邓聚龙拟合结果

我们将样本数据分别带入以Python为媒介建立GM(1,1)、GM(1,1)邓聚龙与GM(1,1)模型中进行回归分析与预测，在这里，我们在本篇文章中需要预测的是福建省2030年预计碳达峰时点时二氧化碳的排放量，因此在灰色预测模型中的预计时点为2030年，为了防止存在误差，我们将模型时点延至2031年，得到的的相关数据以及预测拟合程度结果表2与图1~3如下所示。

Table 2. Model data fitting effect

表2. 模型数据拟合效果

上表是对三种模型拟合效果的描述，由表可以看出相关三种模型的系数(R-Squared)除了GM(1,1)外均高于为0.83，说明GM(2,1)与GM(1,1)邓聚龙的相关性均非常好，拟合度非常高，其中GM(2,1)的相关系数最高，达到了0.990；与此同时GM(1,1)的相关系数只有0.071，小于0.1，说明呈现高度不相关，因此不适用于预测福建省二氧化碳的排放，从下图1中的曲线拟合效果来看也可以说明这点，因此提前将GM(1,1)排除考虑。剩余两种模型的均方根误差(RMSE)最高值为0.359，相对于以万吨为单位(图中python程序的单位为亿吨)的二氧化碳排放总量来说，实际值与预测值的误差是处于一个可接受的区间范围内，相比较而言GM(1,1)邓聚龙的拟合误差效果优于GM(2,1)。并且在GM(2,1)与GM(1,1)邓聚龙模型的拟合过程中，已经对参数进行适当的修正与调整，因此得到的模型与福建省的二氧化碳排放量的适用性非常好，比较适合这类数据的分析，拟合效果如图2，图3所示，从图中可以看出，曲线总体是上升的，但是曲线斜率是呈现一种逐渐增加的趋势，说明未来10年二氧化碳的排放量是呈现一种快速上升的势头，这比较符合接下来为达到碳达峰的排放值，各行业的排放量会逐步增加的预测趋势。

在这里笔者需要补充说明一点，RMSE表示均方根或者称为标准差，其值等于MSE的平方根，最主要的用途式用来评估我么所应用的模型存在的误差大小，通过对值进行开平方可以使得误差的单位级和我们所采用数据单位级保持一致 [7]。

GM(1,1)对应的公式(图1)：

$x\left(t\right)=\frac{81574.25}{123636.49}+\left(x\left(0\right)-\frac{81574.25}{123636.49}\right){\text{e}}^{-81574.25t}$ (33)

GM(2,1)对应的公式(这里取得值为R-Squared最优时的值) (图2)：

$x\left(t\right)=0.0005{\text{e}}^{-40.972t}+0.4535{\text{e}}^{0.0795t}$ (34)

GM(1,1)邓聚龙对应的公式(图3)：

(35)

Figure 1. Fitting graph of GM(1,1)

图1. GM(1,1)拟合图

Figure 2. Fitting graph of GM(2,1)

图2. GM(2,1)拟合图

Figure 3. Fitting graph of GM(1,1) Deng Jurong

图3. GM(1,1)邓聚龙拟合图

4. 结论

在对三种模型的对比分析中，由上图可以看出GM(2,1)与GM(1,1)邓聚龙的拟合趋势相当接近，图像走势中样本点与预测线交叉情况的也一致。总体上GM(2,1)的R-Squared比GM(1,1)邓聚龙小的同时，其RMSE值比GM(1,1)邓聚龙大，但是二者的值差距不大，说明GM(2,1)于GM(1,1)邓聚龙均符合我们用来预测福建省二氧化碳排放量的趋势预测。而GM(1,1)的误差对比前二者而言均更大，且拟合相关性很差，说明GM(1,1)灰色预测模型在此次的预测效果并没有GM(2,1)与GM(1,1)邓聚龙好，不予以采用。综合来说，GM(2,1)与GM(1,1)邓聚龙模型对这种值不太多的预测都有较好的效果，呈现的曲线走势也相对一致，均可以采用来预测福建省在2030年碳达峰前的二氧化碳排放量。

针对如何在2030年到2060年的30年间，如何在比西方国家更快的速度达到碳中和，我们也提出了几点建议以供参考：首先，应改变现有的能源发电方式，采用更多的节能方式，如风力发电，太阳能发电等方式来实现摆脱对煤炭等会大量排放二氧化碳的旧能源的依赖性；其次，针对快速发展的数字媒体技术，我们可以采用互联网来实现对能源的高效管理，例如针对某些高二氧化碳排放的企业，可以利用某些实时互联网技术监测设备来进行管控，实现能源技术与信息技术相结合。最后，碳交易市场的快速发展、碳关税可行性的研究以及碳外汇交易的实现，综合说明了我们可以采取一些政策工具与市场交易机制来控制二氧化碳的排放，以降低社会总成本为目标，实现我国对外界的承诺。

致谢

感谢厦门国家会计学院阎虎勤老师对本篇文章的悉心指导，以及厦门国家会计学院2020级同学的答疑解惑。

基金项目

本论文得到了厦门国家会计学院“云顶课题：Python财务数据分析”项目和大米(厦门)科技股份有限公司的支持