中国煤炭价格指数的多元时间序列分析

Multivariate Time Series Analysis of China Coal Price Index

摘要: 煤炭是我国目前的主要能源，煤炭价格指数是煤炭价格现状的重要呈现。对煤炭价格指数的预测可以帮助企业计划生产、国家调整政策。本文通过ARIMAX模型进行预测发现，WTI对我国煤炭价格指数有显著影响，当没有突发事件发生时，2020年初煤炭价格指数可能会出现小幅上涨。与ARIMA模型拟合对比显示，ARIMA模型和ARIMAX模型均能对我国煤炭价格指数进行良好的预测，但ARIMAX模型优于ARIMA模型。

Abstract: Coal is currently the main energy source in my country, and the coal price index is an important representation of the current status of coal prices. Forecasting the coal price index can help companies plan production and adjust national policies. This article predicts through the ARIMAX model that WTI has a significant impact on my country’s coal price index. When there are no emergencies, the coal price index may rise slightly in early 2020. Comparison with the ARIMA model shows that both the ARIMA model and the ARIMAX model can make a good prediction of my country’s coal price index, but the ARIMAX model is superior to the ARIMA model.

1. 引言

我国是煤炭大国，目前煤炭消费是我国主要的能源消费。据统计，2019年我国能源消费中，煤炭消费量约39.4亿吨，占比57.7%，超过半数 [1]。作为最经济可靠的能源，在可预见的未来，煤炭在我国能源结构中的地位不会发生较大改变。

煤炭价格指数反映了煤炭市场的供需关系，展示了煤炭价格波动的真实情况，其变化发展不仅影响着煤炭的开采与挖掘，对全国各行各业的生产经营都起着十分重要的作用。通过对煤炭价格指数的预测，可以降低生产者的风险，也可以帮助国家在宏观调控上制定政策、采取措施，使我国经济保持稳定持续的发展。

中国煤炭价格指数(China Coal Price Index，简称CCPI)是一种可以全面、客观、及时地描述全国以及各区域、各品种煤炭市场价格变化的数据，它创建于2006年，以2006年1月1日为基期(100点)，在“中国煤炭市场网”按周发布。

申万(2014) [2] 等通过对澳大利亚煤炭价格和迪拜石油价格进行回归分析，证实石油价格与煤炭价格存在正相关关系。王延伟、景晓真(2016) [3] 采用灰色综合评价法对我国煤炭价格进行研究，发现其受经济发展影响大。王帮俊、朱超群(2017) [4] 对环渤海动力煤价格进行短期分析，建立ARIMA模型，得到了良好的预测结果。张同功、赵德志(2018) [5] 对秦皇岛煤炭价格运用逐步回归法发现，对煤炭价格影响较大的是煤炭产量、煤炭进出口量和国际石油价格等。赵修茗、张甜、邹绍辉(2020) [6] 对我国煤炭价格指数等交易数据构建VAR模型，发现我国煤炭经济变动与动力煤期货价格呈均衡关系。

煤炭价格指数受多种因素影响，国际原油价格、原煤产量、原煤进口量、煤炭销量等都会影响其波动，但由于ARIMAX模型要求输入变量为平稳时间序列，通过检验，本文最终选取美国西德克萨斯轻质原油(West Texas Intermediate，简称WTI)价格为ARIMAX模型的输入变量。

使用求和自回归移动平均模型(ARIMA)和传递函数模型(ARIMAX)对2013~2019年我国煤炭价格指数和WTI价格周数据进行分析，研究两种模型在分析我国煤炭价格指数方面的优劣性，并使用两种方法分别对CCPI进行短期预测。

2. 模型理论

2.1. ARIMA模型简介

求和自回归移动平均(ARIMA)模型是在进行一元时间序列分析时的经典模型。其结构为：

$\{\begin{array}{l}\Phi \left(B\right){\nabla }^d{x}_t=\Theta \left(B\right){\epsilon }_t\\ E\left({\epsilon }_t\right)=0,Var\left({\epsilon }_t\right)={\sigma }_{\epsilon }^2,E\left({\epsilon }_t{\epsilon }_s\right)=0,s\ne t\\ E\left({\epsilon }_t{x}_s\right)=0,\forall s<t\end{array}$

式中， ${\nabla }^d={\left(1-B\right)}^d$ ； $\Phi \left(B\right)=1-{\varphi }_{1}B-{\varphi }_2{B}^2-\cdots -{\varphi }_p{B}^p$，为平稳可逆ARMA(p,q)模型的自回归系数多项式；，为平稳可逆ARMA(p,q)模型的移动平滑系数多项式。

通过差分将非平稳时间序列转化为平稳，可以解决绝大多数单变量时间序列的分析问题；但因为很多时间序列的变化都会受到其他序列的影响，若仅考虑自身过去对现在的影响，得出的模型可能并不精确，这就涉及多元时间序列分析。

2.2. ARIMAX模型简介

带有输入序列的ARIMA模型称为传递函数(ARIMAX)模型，可以对多元时间序列进行分析。传递函数模型基于假定：输入序列x_t、输出序列y_t均为零均值平稳时间序列。若序列不平稳，可以通过差分运算或方差齐性变换使其平稳，非零均值序列可通过序列中心变换使序列均值变为零。

实际应用的传递函数模型为：

$Y_t=\frac{\Omega \left(B\right)B^b}{E\left(B\right)}{X}_t+\frac{\Theta \left(B\right)}{\Phi \left(B\right)}{a}_t$

其中， $\Omega \left(B\right)={\omega }_0-{\omega }_{1}B-{\omega }_2{B}^2-\cdots -{\omega }_s{B}^s,{\omega }_0\ne 1$， $E\left(B\right)=1-{\epsilon }_{1}B-{\epsilon }_2{B}^2-\cdots -{\epsilon }_r{B}^r$。

2.3. ARIMAX建模步骤

① 对x_t、y_t进行平稳性检验和纯随机性检验：若不平稳，可通过取对数、求差分等方法使其平稳。

② 预白噪化输入序列x_t：对x_t进行ARMA建模，得到 ${\Phi }_x\left(B\right)x_t={\Theta }_x\left(B\right){\alpha }_t$，即 ${\alpha }_t=\frac{{\Phi }_x\left(B\right)}{{\Theta }_x\left(B\right)}{x}_t$，此时 ${\alpha }_t$ 为白噪声序列。

③ 预白噪化输出序列y_t：令 ${\beta }_t=\frac{{\Phi }_x\left(B\right)}{{\Theta }_x\left(B\right)}{y}_t$。

④ 估计脉冲响应函数 ${\hat{v}}_k$ ：计算 ${\alpha }_t$ 和 ${\beta }_t$ 的互相关系数，脉冲响应函数 ${\hat{v}}_k=\frac{{\sigma }_{\beta }}{{\sigma }_{\alpha }}{\hat{\rho }}_{\alpha \beta }\left(k\right)$。

⑤ 初步判定传递函数模式：识别延迟参数b，衰减参数r和参数s，基于 $v\left(B\right)=\frac{\Omega \left(B\right)}{E\left(B\right)}{B}^b$ 估计参数 ${\hat{\epsilon }}_j$， ${\hat{\omega }}_j$，从而得到传递函数的估计： $\hat{v}\left(B\right)=\frac{\hat{w}\left(B\right)}{\hat{\delta }\left(B\right)}{B}^b$。

⑥ 估计噪声序列： ${\hat{\eta }}_t=y_t-\hat{v}\left(B\right)x_t=y_t-\frac{\Omega \left(B\right)}{E\left(B\right)}{B}^b{x}_t$。

⑦ 噪声序列ARMA建模：对噪声序列构建ARMA模型： ${\Phi }_{\hat{\eta }}\left(B\right){\hat{\eta }}_t={\Theta }_{\hat{\eta }}\left(B\right)a_t$，即，进而得到传递函数模型： $y_t=\frac{\hat{w}\left(B\right)}{\hat{\delta }\left(B\right)}{B}^b{x}_t+\frac{{\Theta }_{\hat{\eta }}\left(B\right)}{{\Phi }_{\hat{\eta }}\left(B\right)}{a}_t$。

⑧ 参数估计：基于以上过程所得相应延迟阶数，对传递函数模型： $y_t=\frac{w\left(B\right)}{\delta \left(B\right)}{B}^b{x}_t+\frac{{\Theta }_{\hat{\eta }}\left(B\right)}{{\Phi }_{\hat{\eta }}\left(B\right)}{a}_t$ 所有未知参数重新进行最大似然估计。

⑨ 模型检验：1) 残差序列 $a_t$ 纯随机性 $Q_{LB}$ 检验；2) 残差序列 $a_t$ 与输入序列 $x_t$ 纯随机性 $Q_{LB}$ 检验。

3. 实例分析

虽然Eviews软件可以实现ARIMAX模型建模的部分步骤，但相较之下，SAS软件能程序化实现整个建模过程，所以本文选用了SAS软件进行实现。

3.1. 数据选择

本文选用2013年1月至2019年12月CCPI和WTI价格进行分析。所用WTI价格来源于美国能源信息网站，CCPI来源于wind数据库。数据时序图如图1所示，CCPI均值为150.2761，标准差为14.0821，最高值出现在2013年1月3日为171.1，最低值出现在2015年11月13日为124.7。

Figure 1. CCPI and WTI price timing chart

图1. CCPI与WTI价格时序图

3.2. 数据检验

由图1可知WTI价格和我国煤炭价格指数均为非平稳时间序列。分别对其进行一阶差分，得到的差分后序列时序图如图2、图3所示。

Figure 2. After difference WTI price sequence diagram

图2. 差分后WTI价格时序图

Figure 3. Differential CCPI sequence diagram

图3. 差分后CCPI时序图

由图2、图3可初步判别两数据差分后为平稳序列。进一步通过单位根检验证实了此结论。

二者的白噪声检验结果显示两序列均为非白噪声序列，可以进行建模。记差分后WTI价格序列为x，差分后CCPI序列为y。

3.3. 建立模型

首先，对 $y_t$ 进行预白造化处理，通过AIC准则对比最小信息量结果，得到最优模型AR(1)，去掉不显著的常数项后拟合结果为

$x_t=\left(1+0.2996B\right)a_t$ ； $a_t\stackrel{iid}{~}N\left(0,{\sigma }_{\epsilon }^2\right)$。

预白噪化CCPI差分序列 $y_t$ 得到序列 ${\tilde{y}}_t$ ：

${\tilde{y}}_t=\frac{1}{1+0.2996B}{y}_t$。

计算 $a_t$ 与 ${\tilde{y}}_t$ 的互相关系数 ${\rho }_{a_t{\tilde{y}}_t}\left(k\right)$，结果见图4，基于该互相关系数初步判定脉冲响应函数的相关参数：b = 1，r = 1，s = 0。

Figure 4. Intercorrelation function of input-output pre-white sequence

图4. 输入输出预白造化后序列的互相关函数

估计结果常数项不显著，去掉常数项后，各个参数显著，估计结果如表1所示。

Table 1. Parameter estimation

表1. 参数估计

基于参数估计结果，脉冲响应函数的估计结果为：

$\hat{v}\left(B\right)=\frac{0.0119B+\mathrm{0.0.0180}{B}^2}{1-0.7904B}$

对噪声序列 ${\eta }_t=y_t-\left(\frac{0.0119B+\mathrm{0.0.0180}{B}^2}{1-0.7904B}\right)x_t$ 尝试ARMA(2,0)模型进行参数估计，结果如表2所示。

Table 2. Parameter estimation

表2. 参数估计

参数通过了参数显著性检验，得到最终传递函数模型：

$y_t=\left(\frac{0.0242B+0.02986B^2}{1+0.4971B}\right)x_t+\left(\frac{1}{1-0.5521B+0.1578B^2}\right){\epsilon }_t$

该估计结果中的参数都通过了显著性检验，模型残差通过了白噪声检验，但模型残差与输入变量之间仍存留一定相关性，对该模型的预测精度有一定影响。

可以看出，输入变量WTI价格对滞后1期CCPI有显著影响，随着时间的发展，当期国际原油价格对滞后期煤炭价格指数的影响将变得越来越小。

3.4. 预测

通过以上模型拟合，对我国2019年12月至2020年2月煤炭价格指数进行预测，预测结果如表3。

Table 3. Forecast of China’s coal price index

表3. 我国煤炭价格指数预测

结果显示，我国煤炭价格指数在2019年末和2020年初有递增趋势，且浮动不大。这可能是因为进入冬季，天气寒冷，供暖消耗大量煤炭，所以导致价格上涨。

3.5. 对比结论

通过以上模型拟合，可以看出，国际原油价格对我国煤炭价格指数有十分重要的影响。输入变量为国际原油价格时，可以使用ARIMAX模型对CCPI进行预测。

使用ARIMA模型和ARIMAX模型的预测效果图如图5、图6所示，可见两方法的拟合效果均良好。

Figure 5. After difference WTI price sequence diagram

图5. 基于ARIMA模型的CCPI预测

Figure 6. Differential CCPI sequence diagram

图6. 基于ARIMAX模型的CCPI预测

将两模型的检验结果对比，对比结果如表4所示。相比之下，ARIMAX模型的标准误差较小，不同准则下的最小信息量也较小，整体优于ARIMA模型。

Table 4. ARIMA and ARIMAX model test results

表4. ARIMA与ARIMAX模型检验结果对比

4. 小结

4.1. 结论

分别使用ARIMA模型和ARIMAX模型对我国煤炭价格指数预测的过程中发现，两种预测方法均能得到较好的拟合效果，且ARIMAX模型的结果优于ARIMA模型。

WTI对我国煤炭价格指数有显著影响。在我国能源结构中，煤炭和石油是主要消费的一次能源，二者互为替代品，当其中一种价格上升时，对另一种的需求就会增大，从而也就带动了另一种价格的上涨。

但是由于ARIMAX模型对输入变量的平稳性有较高要求，所以当其无法满足时也可使用ARIMA模型。ARIMAX模型作为一种限制条件较多的模型，应用受到很大的限制，但其结果十分可靠。使用更多解释变量的ARIMAX模型可能会得到更好的预测效果。

4.2. 不足与展望

我国煤炭价格指数受诸多因素影响，本文主要研究了CCPI与WTI价格之间的关系，还可以尝试寻找其他平稳或差分平稳时间序列对CCPI进行ARIMAX模型建模。

现有模型普遍在短期预测时效果较好，后续需根据实际情况对模型进行检验与调整。

虽然通过ARIMA模型和ARIMAX模型均得到了较好的拟合效果，且二者预测结果较为一致，但对于突发事件的影响是无法预测的。如受今年新冠肺炎的影响，物流收阻，人员流动受限，CCPI的上涨趋势可能也会受到影响。

随着经济全球化趋势不断发展，国内外联系日趋紧密，我国煤炭市场与国际原油市场的关联将日益增加，相关企业与国家在进行决策时都应考虑国际原油价格的长远影响。

基金项目

本文受中国矿业大学(北京)“课程思政”示范课程建设项目(时间序列分析)的支持。

NOTES

^*通讯作者。