1. 引言

全国碳排放权交易市场的建设是实现“碳达峰、碳中和”目标的核心政策工具之一。碳交易把原本游离在资产负债表外的气候变化因素纳入进来，联通了金融资本和实体经济，改变了企业的收支结构。同时，碳定价权的争夺也将成为国际碳市场中的竞争核心。我国于2011年末开展“两省五市”碳排放权交易试点工作，并于2022年7月启动上线全国碳排放权交易市场交易。经过多年试点和推广，我国碳市场初步形成，但各试点地区碳交易价格存在着价格悬殊、涨跌波动幅度大等特点，这不利于碳市场的长期稳定运行，而合理的碳价对于我国实现“双碳目标”，以及激励碳减排企业能够起到积极作用。

目前针对碳交易价格预测及波动特征的研究方法主要有以下两类：一是基于时间序列分析方法。学者通过构建ARIMA模型、GARCH族模型 [1] [2] [3] 、VAR模型及相关组合模型(夏睿瞳，2018 [4] 和白强等，2022 [5] )，对国内碳排放权价格的波动和风险特征进行了诸多分析，如：辛姜和赵春艳(2018) [6] 将马尔可夫转移模型与VAR模型结合，朱若瑾(2022) [7] 对新冠疫情前后两个情景使用GARCH模型、刘红琴和胡淑慧(2022) [8] 对所设置的4种交易情境使用TGARCH模型，实证发现我国碳市场的以下一些特征：碳价格受到的影响多维时变、疫情冲击能够一定程度修正碳价格的未来波动、不同情境下碳排放交易市场风险存在差异等。二是机器学习和深度学习方法。使用单一机器学习和深度学习方法，如：神经网络 [9] 、经验模态分解 [10] 、Transformers [11] 等。使用组合模型方法，包括：组合不同机器学习和深度学习模型、与ARMA模型组合等。崔焕影和窦祥胜(2018) [12] 使用五种算法模型(EMD、GA、BP、PSO和LSSVM)做短期预测和长期预测、闫梦和王聪(2020) [13] 使用四种算法模型(EMD、BP、DNN和SVM)预测，均发现组合模型的预测误差比单一模型小。吕靖烨等(2020) [14] 考虑到碳交易价格时间序列同时具有线性和非线性特征，构建ARIMA-SVM组合模型。谢旭升和严思屏(2021) [15] 分别从最优权重和残差优化角度，建立ARMA-BP组合模型。公维凤等(2022) [16] 考虑到碳交易价格收益率的波动聚集特征、风险溢价及非对称性，构建GARCH-BP组合模型。

尽管机器学习和深度学习方法能够有效提高预测精确度，但由于其黑盒特性解释性较差。因此，本文将聚焦于我国碳排放权交易价格的波动特征，使用单变量和多元时间序列分析方法建模和预测，为完善我国碳排放定价机制和交易制度提供理论意义。

2. 模型与方法

2.1. ARMA模型

自回归滑动平均模型(ARMA模型)，由自回归模型(AR模型)与移动平均模型(MA模型)为基础混合构成。ARMA(p,q)模型表达式如下，

$r_t={\varphi }_0+\underset{i=1}{\overset{p}{{\displaystyle \sum }}}{\varphi }_i{r}_{t-i}+a_t-\underset{i=1}{\overset{q}{{\displaystyle \sum }}}{\theta }_i{a}_{t-i}$(1)

其中， $\left\{a_t\right\}$是白噪声序列，AR模型和MA模型都是ARMA(p,q)的特殊情形。

2.2. GARCH族模型

2.2.1. 简单GARCH模型

GARCH模型是自回归条件异方差模型(ARCH模型)的拓展，能准确地模拟收益率时间序列变量波动性的变化，使人们能更加准确地把握风险。对于一个对数收益率序列 $r_t$，令 $a_t=r_t-{\mu }_t$为t时刻的新息，称 $a_t$服从GARCH(m,s)模型，若 $a_t$满足下式，

$a_t={\sigma }_t{\epsilon }_t,{\sigma }_t^2=\omega +\underset{i=1}{\overset{m}{{\displaystyle \sum }}}{\alpha }_i{a}_{t-i}^2+\underset{j=1}{\overset{s}{{\displaystyle \sum }}}{\beta }_j{\sigma }_{t-j}^2$(2)

其中， $\omega >0$， ${\alpha }_i\ge 0$， ${\beta }_j\ge 0$， $\left\{{\epsilon }_t\right\}$是均值为0、方差为1的独立同分布随机变量序列，通常假定 ${\epsilon }_t$是标准正态的或标准化的学生t分布。在许多情形下，只用到低阶的GARCH模型，通常GARCH(1,1)模型就足够了。

2.2.2. EGRCH模型

在简单GARCH模型中，对正的和负的扰动给予相同的反应，为了克服在处理金融时序时的这一弱点，Nelson (1991) [17] 提出了指数GARCH (EGARCH)模型，通过考虑加权的新息来体现正的和负的资产收益率的非对称效应：

$g\left({\epsilon }_t\right)=\alpha {\epsilon }_t+\gamma \left[\left|{\epsilon }_t\right|-E\left(\left|{\epsilon }_t\right|\right)\right]=\{\begin{array}{l}\left(\alpha +\gamma \right){\epsilon }_t-\gamma E\left(\left|{\epsilon }_t\right|\right){\epsilon }_t\ge 0\\ \left(\alpha -\gamma \right){\epsilon }_t-\gamma E\left(\left|{\epsilon }_t\right|\right){\epsilon }_t<0\end{array}$(3)

其中， $\alpha$和 $\gamma$是实常数。使用滞后算子形式，EGARCH(m,s)模型可以表示为

$a_t={\sigma }_t{\epsilon }_t,\ln \left({\sigma }_t^2\right)=\omega +\frac{1+{\beta }_{1}B+\cdots +{\beta }_{s-1}{B}^{s-1}}{1-{\alpha }_{1}B-\cdots -{\alpha }_m{B}^m}g\left({\epsilon }_{t-1}\right)$(4)

其中， $\omega$是常数，B是延迟算子， $Bg\left({\epsilon }_t\right)=g\left({\epsilon }_{t-1}\right)$。更一般地，模型表达式为

$a_t={\sigma }_t{\epsilon }_t,\ln \left({\sigma }_t^2\right)=\omega +{\displaystyle {\sum }_{i=1}^m\left[{\alpha }_i\left(\left|{\epsilon }_{t-j}\right|-E\left|{\epsilon }_{t-j}\right|\right)+{\gamma }_i{\epsilon }_{t-j}\right]}+{\displaystyle {\sum }_{j=1}^s{\beta }_j\ln \left({\sigma }_{t-j}^2\right)}$(5)

其中， ${\gamma }_i$能够反映对数收益率的正负扰动对波动率的不同影响。

本文将对均值方程和波动率方程进行联合估计，对资产收益率序列建立ARMA-GARCH和ARMA-EGARCH模型。

2.3. VAR模型

向量自回归模型(VAR模型)是AR模型的推广，可以估计联合内生变量的动态关系。时间序列 $r_t$服从VAR(p)模型，如果它满足

$r_t={\varphi }_0+{\Phi }_1{r}_{t-1}+\cdots +{\Phi }_p{r}_{t-p}+a_t$(6)

其中， ${\varphi }_0$是一个k维向量， ${\Phi }_j$是一个k × k矩阵， $\left\{a_t\right\}$是一个均值为0、协方差矩阵为 $\Sigma$的序列不相关的随机向量序列，通常假定其为多元正态。

将VAR(p)模型用过去的新息表示，可以改写为

$r_t=\mu +a_t+{\psi }_1{a}_{t-1}+{\psi }_2{a}_{t-2}+\cdots$(7)

其中， $\mu ={\left[\Phi \left(1\right)\right]}^{-1}{\varphi }_0$， ${\psi }_i$为过去新息 $a_{t-i}$对 $r_t$的影响，称为 $r_t$的脉冲响应函数。脉冲响应分析，可以描述某一内生变量对于残差冲击的反应，即在随机误差项上施加一个标准差大小的某一内生变量冲击后，对感兴趣的内生变量的当期值和未来值所产生的影响。

VAR模型的应用，还可以采用方差分解方法研究模型的动态特征，分析预测残差的标准差由不同新息的冲击影响的比例，进一步评价各内生变量对预测方差的贡献度。

3. 实证分析

3.1. 数据来源和选取

考虑到地区分布、碳交易规模和活跃性等因素，本文选取了广州、湖北、北京和上海四个我国典型的碳交易市场为研究对象，将这四个市场的碳排放权交易价格的日收盘价作为样本数据，计算其对数收

益率(百分数) $r_t=\ln \left(\frac{P_t}{P_{t-1}}\right)\ast 100$，以此来研究我国碳排放权交易价格的波动特征。选取样本区间为2017

年4月至2022年12月，数据来源于各碳排放权交易所，分别得到1368、1344、786和860条原始数据。

3.2. 碳排放权交易价格及对数收益率分布特征

四个省市的碳排放权交易价格和收益率走势如图1所示，通过时序图可以看到，四个碳市场的收益率序列均围绕0上下波动，且呈现“波动聚集”特征。其中，广州碳排放权交易价格大致可以分为两个阶段，在2019年以前对数收益率波动较大，2019年以后波动减小但价格上涨趋势明显，均价为34.25元/吨；湖北碳排放权交易价格对数收益率波动比较平稳，在2019年年中和2022年年初出现过价格明显抬升，均价为30.29元/吨；北京碳排放权交易价格则相对较高，均价为72.82元/吨，对数收益率波动也较大，在2021年存在两次大幅下跌和2022年下半年连续上涨趋势；上海碳排放权交易价格与湖北类似，对数收益率波动较小，均价为43.49元/吨。

碳交易价格的对数收益率序列分布特征如表1所示，由描述性统计表可以看到：四个碳市场的收益率序列均平稳且不服从正态分布，其中只有上海呈现为轻微右偏且超额峰度较大，其他三个碳市场均为轻微左偏和尖峰态，由此可见碳交易价的对数收益率普遍存在“尖峰厚尾”的特征；另外，北京的收益率标准差最大，表明该序列的价格波动最大，风险较高。

3.3. 基于单变量时间序列分析的碳排放权交易价格波动特征研究

3.3.1. ARMA模型构建与模型检验

首先，由单位根ADF检验结果表明，广州、湖北、北京和上海四个碳排放权交易市场的日收益率序列均不存在单位根，满足平稳序列条件。然后，分别计算各市场日收益率的自相关系数(ACF)、偏相关系数(PACF)及EACF，并结合R语言系统定阶结果，对滞后阶数进行调整，在保证各回归系数统计显著和通过模型拟合充分检验的基础上，选择对数似然值最大、AIC值最小的模型作为最优拟合模型，得到各碳市场的最优ARMA模型及相关残差检验结果如表2所示。

由ARMA模型的残差检验结果可知，四个模型的Ljung-Box统计量均不能拒绝对数收益率没有显著的前后相关性这一原假设，即模型残差为白噪声，模型拟合充分。ARCH LM卡方检验p值均显著小于0.05，说明对数收益率的ARMA模型均存在ARCH效应。因此，后续考虑在上述模型基础上加入GARCH族模型来描述条件异方差性，更好地拟合收益率波动特征。

3.3.2. ARMA-GARCH族模型构建与模型检验

分别在四个ARMA模型基础上加入GARCH(1,1)和EGARCH(1,1)进行联合估计，除北京碳市场外，ARMA模型的滞后阶数与上文一致，因为在对北京碳市场建模时，考虑到系数的统计显著性进行了修正，改用ARMA(1,1)模型与GARCH族模型结合。根据AIC准则，对广州、湖北和上海碳市场选取学生t分布，对北京碳市场用有偏学生t分布来刻画条件均值方程残差分布的有偏、肥尾等特征，模型结果如表3所示。由各模型的标准化残差和标准化残差平方的序列相关Ljung-Box检验可知，八个模型的拟合正确且充分。

分析各市场的拟合模型发现：1) 从接近于1的ARCH项与GARCH项系数 ${\alpha }_1$、 ${\beta }_1$，可以发现当期收益率的波动程度会受到前期波动的显著影响，且波动的持续时间普遍较长。2) 由ARMA-EGARCH模型中的非对称项系数 ${\gamma }_1$，得到广州和上海碳市场的 ${\gamma }_1$大于0，说明正冲击的影响大于负冲击，而在湖北和北京碳市场中则负冲击的影响大于正冲击，且只有上海的 ${\gamma }_1$系数统计显著。3) 由信息冲击曲线图(图2)，也可以得到相同的结论，广州、湖北和北京三个碳市场在利空和利好消息时曲线均走势陡峭，上海碳市场则只对利好消息敏感。4) 从似然值和AIC准则角度来看，对广州、北京和上海市场来说，ARMA-GARCH模型更优；对湖北市场来说，ARMA-EGARCH模型更优。

注：1.^***、^**、^*分别表示1%、5%、10%水平上的显著，括号中为系数显著性检验统计量的p值。2. Q(20)和Q²(20)分别表示标准化残差和标准化残差平方的序列相关Ljung-Box检验的p值。

3.3.3. 模型预测

基于所构建的ARMA模型、ARMA-GARCH模型和ARMA-EGARCH模型对四个市场分别进行向前20步预测，如图3所示，其中黑色圆点表示对数收益率的实际值，红色×表示预测值，蓝色虚线表示置信水平为95%的预测区间。可以看到：在ARMA模型下，预测区间能较好的包含四个市场的实际值，其中在北京市场下的效果稍差些；在ARMA-GARCH(1,1)模型下，预测区间能包含全部实际值，但是其在广州和湖北碳市场中，当预测步长较大时的预测区间跨度较大，ARMA-EGARCH(1,1)模型则能够克服以上两种模型的缺点，在合适的预测区间大小内基本包含了所有实际值；而在北京和上海碳市场中，从预测效果来看，ARMA-GARCH(1,1)模型则更优。

3.4. 基于多元时间序列分析的碳排放权交易价格波动特征研究

3.4.1. VAR模型构建与模型检验

考虑到四个碳市场间可能存在相互影响，本节将通过对四个碳市场的对数收益率序列构建VAR模型进行考察。由于四个序列长度不等，保留四个碳市场都存在数据的交易日数据，得到502 × 4条观测。通过R语言系统定阶，最终选择VAR(6)进行拟合。

此时，模型残差的Portmanteau检验p值为0.2305，不拒绝残差为白噪声的原假设，说明模型拟合充分。另外，通过模型稳定性检验得到该模型稳定，可以进行脉冲响应分析，脉冲响应结果如图4所示(y1、y2、y3、y4分别代表广州、湖北、北京和上海)。分别给四个碳市场一个标准差的冲击，可以看到：广州、湖北和上海碳市场对北京碳市场均有一定影响，广州和湖北碳市场对上海也有一定影响，而北京碳市场对其他三个市场影响很小。进一步采用方差分解方法研究各市场对预测方差的贡献度，如图5所示，可以看到三个碳市场均受自身市场影响大，而受其他市场影响小，说明我国碳交易市场流动性较低。

3.4.2. 模型预测

基于VAR模型对四个碳市场做20步预测(如图6所示)，可以看到：当预测步长较小时对四个碳市场的收益率走势预测与实际基本一致，同时，预测区间能够较好地包含所有实际值，且预测区间范围更为合理。因此，使用多元时间序列方法能够更好的捕捉各碳市场交易价格的对数收益率波动特征。

4. 总结

本文基于我国四个典型的碳交易市场(广州、湖北、北京和上海)，以其碳排放权交易价格的对数收益率(%)为研究对象，通过构建ARMA-GARCH模型、ARMA-EGARCH模型和VAR模型，对我国碳交易市场的波动特征进行了实证研究。主要结论如下：

首先，上述四个碳交易市场的交易活跃度、交易量和价格存在一定的区域差异。其中，广州和湖北的碳交易较为活跃，拥有着较大的交易量；北京的碳交易的平均价格最高，波动幅度也最大；上海碳交易价格较为平稳，价格走势与湖北存在一定相似性。

其次，四个碳交易市场的日对数收益率均存在显著的自相关性和条件异方差性，且存在一定的相互影响。当期收益率会同时受前期收益率大小和波动的显著影响，且波动的持续时间普遍较长，使用GARCH模型联合估计能够帮助拟合波动特征。基于多元时间序列分析的VAR模型则能够更好地捕捉市场间的影响，在预测步长较小时，对四个碳市场的收益率走势预测与实际基本一致，且预测区间范围比单变量分析更为合理。

另外，四个碳市场的抗风险能力和流动性仍存在提高空间。广州、湖北和北京三个碳市场对利空和利好消息均敏感，上海碳市场则只对利好消息敏感，其中，在湖北和北京碳市场中负冲击的影响大于正冲击，而在广州和上海碳市场中则相反。同时，四个碳交易市场均受自身市场影响大，而受其他市场影响小。其中，广州、湖北和上海碳市场对北京碳市场均有一定影响，广州和湖北碳市场对上海也有一定影响，而北京对其他市场基本没有影响。

由此可见，我国碳交易市场的价格对数收益率有很强的记忆性，该收益率特征与有效市场的随机游走相悖，说明我国碳交易市场尚未成熟，市场效率偏低，价格缺乏市场机制决定的特点，并且对消息较为敏感，市场流动性也较低。

由此，本文建议“坚持市场导向 + 政府部门监管和引导”：1) 坚持市场导向，提高市场流动性，建立合理的价格体系，积极发挥碳价格信号的激励作用。一方面，适当降低企业、机构和个人投资者的准入门槛，增加碳产品种类，进一步提高碳交易的规模和活跃度。另一方面，建立完善全国统一的碳交易市场，解决市场分割问题，帮助实现地区间要素、技术和产品按照价格信号进行自由流动。2) 政府部门加强监管和引导，帮助碳排放权交易价格趋于稳定。一方面，完善立法、监管等政策制度，形成全国统一的管理标准，完善信息披露制度，监控利空和利好消息面的影响，建立风险预警机制。另一方面，积极借鉴和学习国外先进市场的管理经验，结合试点经验完善我国碳市场制度设计，建立符合我国国情的全国性碳排放权交易市场。

参考文献