基于两阶段多变点估计对动力煤价格的实证分析

Empirical Analysis of Thermal Coal Price Based on Two-Stage Multiple Change-Point Estimation

摘要: 将变量选择技术应用到变点估计中可以同时得到多变点数量和位置的估计结果，大大降低了计算复杂度。文章首先介绍了两阶段多变点检测与估计方法的基本概念，以2021年12月22日至2024年11月6日的秦皇岛动力煤价格作为研究对象，并对其中的训练集建立线性回归模型。其次，利用此变点估计方法估计线性回归模型中系数的多变点，根据变点个数和位置的估计结果建立更精准的分段线性回归模型，并运用最后一个子段的线性回归模型对测试集进行短期预测。实证结果表明，分段线性回归模型能更精准地刻画和预测秦皇岛动力煤价格的走势。

Abstract: Applying variable selection techniques to multiple change-point estimation can yield estimates of the number and location of multiple change-point at the same time, significantly reducing computational complexity. The article first introduces the basic concepts of two-stage multiple change-point detection and estimation methods, takes Qinhuangdao thermal coal prices from December 22, 2021 to November 6, 2024 as the research object, and builds a linear regression model for the training set therein. Secondly, this change-point estimation method is utilized to estimate the multiple change-point of the coefficients in the linear regression model, establish more accurate segmented linear regression modeling based on the estimation of the number and location of change-point, and apply the linear regression model of the last subsection to make short-term predictions about the test set. The empirical results show that the segmented linear regression model can more accurately portray and predict the trend of Qinhuangdao thermal coal prices.

1. 引言

变点是一种广泛存在于自然界和人类社会的现象，经济政策的改变、突发事件、自然灾害等原因都可能导致模型的结构发生变化，现实中数据符合的模型结构通常不具有稳定性，在进行统计建模之前进行变点的检测与估计可以帮助我们建立更精准的模型。变点问题自1954年由统计学家Page (1954) [1]在质量控制领域提出后，其理论研究和应用不断发展，还被广泛应用于气象学[2]、金融学[3]和医学[4]等多个领域。利用变量选择方法检测和估计变点是目前十分流行且有效的方法。Jin等(2011) [5]通过分割数据并引入一个特殊的设计矩阵将多变点检测与高维线性回归模型的变量选择联系起来，利用极小极大凹惩罚(Minimax Concave Penalty, MCP)和带平滑削边绝对偏离法(Smoothly Clipped Absolute Deviation, SCAD)检测分段平稳自回归过程(Piecewise Stationary Autoregressive Process, PSAR)中的多变点。Jin等(2016) [6]提出了一种快速且精准的两阶段方法，应用MCP、SCAD、最小绝对值收敛(Least Absolute Shrinkage and Selection Operator, LASSO)和自适应Lasso (Adaptive Lasso)四种变量选择方法检测线性回归模型中的多变点。Sun和Wu (2020) [7]针对广义线性回归模型提出了一种两阶段多变点估计方法，并证明了估计量的相合性。吕丽和金百锁(2021) [8]基于两阶段多变点估计方法，应用随机加权自助法结合高斯混合模型得到线性回归模型中多变点的置信区间。

煤炭是中国能源的重要组成部分，动力煤是煤炭的一个重要类型，在电力、钢铁、建材和化工方面都起着重要作用。近些年，国内学者对于动力煤价格的分析与预测做了大量研究，唐静和王艳洁等(2024) [9]利用长短期记忆网络(Long Short Term Memory, LSTM)模型结合混合核密度估计方法给出了煤炭价格指数的波动区间以及概率评估，并利用环渤海动力煤价格验证了所提方法的有效性。赵修茗等(2020) [10]通过建立向量自回归(Vector Autoregression, VAR)模型研究动力煤期货价格与其影响因素的动态关系，得到了动力煤期货价格与煤炭行业失业率呈负相关的结论。云小鹏(2018) [11]从统计学和煤炭产业两个角度解释了动力煤价格的波动。刘威等(2020) [12]建立了时间序列模型、多元线性回归模型和神经网络模型，通过多种模型结合实现了秦皇岛动力煤价格的预测。本文选取2021年12月22日至2024年11月6日的秦皇岛动力煤价格作为研究对象，对其中的训练集建立线性回归模型，并利用两阶段多变点检测与估计方法估计回归系数的多变点，建立了更精准的分段线性回归模型来刻画和预测秦皇岛动力煤价格的走势。

2. 两阶段多变点检测与估计方法的介绍

Jin等(2016) [6]针对线性回归模型提出了两阶段多变点检测方法，第一阶段是分割阶段(cutting stage)，将数据序列分割为$p_n+1$ 段，该分割方法降低了变点估计的计算复杂度，并通过引入一个特殊的设计矩阵将线性回归模型的变点检测问题转化为高维线性回归模型的变量选择问题，利用MCP、SCAD、Lasso和自适应Lasso方法得到高维线性回归模型中回归系数的估计量。第二阶段是精炼阶段(refining stage)，基于拟似然比检验确定分割数据中准确的变点位置，并证明了变点个数估计和位置估计的相合性，具体如下：

Jin等考虑一个含有s个变点的线性回归模型，其中$s\ge 2$ ，$1<a_1<\cdots <a_s<n$ ，模型如下[6]：

$\begin{array}{c}{y}_i={\displaystyle \sum }_{j=1}^q{x}_{i,j}{\beta }_{j,0}+{\displaystyle \sum }_{l=1}^s{\displaystyle \sum }_{j=1}^q{x}_{i,j}{\delta }_{j,0}I\left(a_l<i\le n\right)+{\epsilon }_i\\ =x_i^{\text{T}}\left[{\beta }_0+{\displaystyle \sum }_{l=1}^s{\delta }_{l}I\left(a_l<i\le n\right)\right]+{\epsilon }_i,i=1,\cdots ,n\\ =\{\begin{array}{ll}{x}_i^{\text{T}}{\beta }_0+{\epsilon }_i\hfill & 1\le i\le a_1,\hfill \\ x_i^{\text{T}}\left({\beta }_0+{\delta }_1\right)+{\epsilon }_i\hfill & a_1<i\le a_2,\hfill \\ ⋮\hfill & ⋮\hfill \\ x_i^{\text{T}}\left({\beta }_0+{\displaystyle \sum }_{l=1}^s{\delta }_l\right)+{\epsilon }_i\hfill & a_s<i\le n,\hfill \end{array}\end{array}$ (1)

其中，$\left\{y={\left(y_1,\cdots ,y_n\right)}^{\text{T}}\right\}$ 是n维观测值，$\left\{x_i={\left(x_{i,1},\cdots ,x_{i,q}\right)}^{\text{T}}\right\}$ 是q维解释变量序列，${\beta }_0={\left({\beta }_{1,0},\cdots ,{\beta }_{q,0}\right)}^{\text{T}}\ne 0$ 是q维回归系数，s是变点数量，$a_1,\cdots ,a_s$ 是变点位置，${\delta }_l={\left({\delta }_{1,l},\cdots ,{\delta }_{q,l}\right)}^{\text{T}}$ ，$\left(l=1,\cdots ,s\right)$ 是变点处回归系数的增量，且${\beta }_0$ 、s、$a_j$ 、${\delta }_l$ 都是未知的，${\epsilon }_1,\cdots ,{\epsilon }_n$ 是随机误差。

2.1. 第一阶段(分割阶段)

2.1.1. 将线性回归模型转化为高维线性回归模型

先将观测数据分割为$p_n+1$ 段，其中$p_n+1>s$ ，第一段的长度为$n-p_{n}m$ ，其他所有段的长度都为m (段长m的选取方法参见文献[6])，$m=n/\left(p_n+1\right)$ (这里c满足$0\le c-c<1$ )，记第一段观测点的指标集为${ℐ}_1=\left\{1,\cdots ,n-p_{n}m\right\}$ ，第j段观测点的指标集为${ℐ}_j=\left\{n-\left(p_n-j+2\right)m+1,\cdots ,n-\left(p_n-j+1\right)m\right\}$ ($j=2,\cdots ,p_n+1$ )，令$X$ 是一个矩阵，

$\begin{array}{c}X={\left(X^{\left(1\right)},X^{\left(2\right)},\cdots ,X^{\left(p_n+1\right)}\right)}_{n\times q\left(p_n+1\right)}\\ =\left(\begin{array}{ccccc}{X}_{\left(1\right)}& 0& 0& \cdots & 0\\ X_{\left(2\right)}& X_{\left(2\right)}& 0& \cdots & 0\\ X_{\left(3\right)}& X_{\left(3\right)}& X_{\left(3\right)}& \cdots & 0\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ X_{\left(p_n+1\right)}& X_{\left(p_n+1\right)}& X_{\left(p_n+1\right)}& \cdots & X_{\left(p_n+1\right)}\end{array}\right)\end{array}$ (2)

上式中，$X^{\left(1\right)}={\left(X_{\left(1\right)}^{\text{T}},X_{\left(2\right)}^{\text{T}},\cdots ,X_{\left(p_n+1\right)}^{\text{T}}\right)}^{\text{T}}$ ，$X^{\left(j\right)}={\left(0_{q\times \left(n-\left(p_n-j+2\right)m\right)},X_{\left(j\right)}^{\text{T}},X_{\left(j+1\right)}^{\text{T}},\cdots ,X_{\left(p_n+1\right)}^{\text{T}}\right)}^{\text{T}}$ ，$X_{\left(1\right)}={\left(x_1,\cdots ,x_{n-p_{n}m}\right)}^{\text{T}}$ 为序列$\left\{x_i,i=1,2,\cdots ,n\right\}$ 第${ℐ}_1$ 个子段的解释变量，$X_{\left(j\right)}={\left(x_{n-\left(p_n-j+2\right)m+1},\cdots ,x_{n-\left(p_n-j+1\right)m}\right)}^{\text{T}}$ 为第${ℐ}_j$ 个子段的解释变量$\left(j=2,3,\cdots ,p_n+1\right)$ ，这样原来$n\times q$ 维的自变量矩阵就被拓展为$n\times q\left(p_n+1\right)$ 维。假设每个子段中至多有一个变点(At Most One Change-point, AMOC)，若变点$a_j$ 位于第$k_j$ 个子段中，即$a_j\in {ℐ}_{k_j}=\left\{n-\left(p_n-k_j+2\right)m+1,\cdots ,n-\left(p_n-k_j+1\right)m\right\}$ ，$\left(j=1,2,\cdots ,s\right)$ 。记${\mathcal{A}}_n=\left\{1,k_1,k_1+1,\cdots ,k_s,k_s+1\right\}$ ，即${\mathcal{A}}_n$ 表示由第一个子段和变点$a_j$ 所在的子段以及其后一子段的指标集构成的集合。

基于以上数据序列的分割及下三角矩阵(2)的构造，模型(1)的矩阵形式为：

$y=X{\theta }^{*}+X_{\omega }\overrightarrow{\omega }+\epsilon$ (3)

其中，$y={\left(y_1,y_2,\cdots ,y_n\right)}^{\text{T}}$ ，$X$ 是由(2)式定义的下三角矩阵，${\theta }^{*}={\left({\theta }_1^{\text{T}},{\theta }_2^{\text{T}},\cdots ,{\theta }_{p_n+1}^{\text{T}}\right)}^{\text{T}}$ ，关于$X_{\omega }$ 、$\overrightarrow{\omega }$ 、${\theta }^{*}$ 的定义以及模型(3)的推导详见文献[6]。这样，模型(1)中多变点的估计问题就转化为高维线性回归模型(3)的变量选择问题。

2.1.2. 应用SCAD方法对高维线性回归模型进行变量选择

由文献[13]和[14]知高维线性回归模型(3)中回归系数${\theta }^{*}$ 的SCAD估计为：

${\widehat{\theta }}^{*}=\underset{{\theta }^{*}}{\arg \min }\left\{\frac{1}{n}{\Vert y-X{\theta }^{*}\Vert }^2+{\displaystyle \sum }_{j=1}^{p_n+1}{\displaystyle \sum }_{i=1}^q{p}_{\lambda ,\gamma }\left(\left|{\theta }_{ji}\right|\right)\right\},$ (4)

其中，对于$u\in \left[0,\infty \right)$ ，SCAD的惩罚项为：

$p_{\lambda ,\gamma }\left(u\right)=\lambda u{I}_{\left[0,\lambda \right]}\left(u\right)+\frac{\gamma \lambda u-0.5\left(u^2+{\lambda }^2\right)}{\gamma -1}{I}_{\left(\lambda ,\gamma \lambda \right]}\left(u\right)+\frac{{\lambda }^2\left({\gamma }^2-1\right)}{2\left(\gamma -1\right)}{I}_{\left(\gamma \lambda ,\infty \right)}\left(u\right),$

(4)式中，$\gamma >0$ ，正则化参数$\lambda >0$ 且可通过贝叶斯信息准则(BIC)得到。并且SCAD的惩罚项满足以下假设：如果$u>\gamma \lambda$ ，则$p_{\lambda ,\gamma }\left(0\right)=0$ ，${p^{\prime }}_{\lambda ,\gamma }\left(u\right)=0$ ，${p^{\prime }}_{\lambda ,\gamma }\left(0\right)=\lambda$ 。

2.2. 第二阶段(精炼阶段)

令${\widehat{\theta }}^{*}={\left({\widehat{\theta }}_1^{\text{T}},{\widehat{\theta }}_2^{\text{T}},\cdots ,{\widehat{\theta }}_{p_n+1}^{\text{T}}\right)}^{\text{T}}$，其中${\widehat{\theta }}_j={\left({\widehat{\theta }}_{j1},\cdots ,{\widehat{\theta }}_{jq}\right)}^{\text{T}}$ ($j=1,\cdots ,p_n+1$ )。记${\widehat{\mathcal{A}}}_n=\left\{j:{\widehat{\theta }}_j\ne 0,j=1,\cdots ,p_n+1\right\}$ 为系数不为零的子段的指标集构成的集合，则${\widehat{\mathcal{A}}}_n^{\text{*}}=\left\{j:j\in {\widehat{\mathcal{A}}}_n,j-1\notin {\widehat{\mathcal{A}}}_n,j=2,\cdots ,p_n+1\right\}=\left\{{\widehat{k}}_1,\cdots ,{\widehat{k}}_{\widehat{s}}\right\}$为变点可能存在的子段的指标集构成的集合(${\widehat{k}}_1<\cdots <{\widehat{k}}_{\widehat{s}}$ )。为了得到变点$\left\{a_j\right\}$ 的估计值$\left\{{\widehat{a}}_j\right\}$ ，应用文献[15]中的定理3.1.1，利用拟似然比检验去检测${\widehat{ℐ}}_{\left(j\right)}={ℐ}_{{\widehat{k}}_j-1}\cup {ℐ}_{{\widehat{k}}_j}=\left\{n-\left(p_n-{\widehat{k}}_j+3\right)m+1,\cdots ,n-\left(p_n-{\widehat{k}}_j+1\right)m\right\}$ $\left(j=1,\cdots ,\widehat{s}\right)$ 内是否存在变点，具体如下：

考虑以下模型，其可能的变点位置为${\zeta }_j$ $\left(j=1,\cdots ,\widehat{s}\right)$

$y_i=x_i^{\text{T}}{\beta }_{j}I\left(n_j^{\left(l\right)}\le i\le {\zeta }_j\right)+x_i^{\text{T}}{\beta }_{j+1}I\left({\zeta }_j<i\le n_j^{\left(r\right)}\right)+{\epsilon }_i,$ (5)

其中，$n_j^{\left(l\right)}=n-\left(p_n-{\widehat{k}}_j+3\right)m+1$ ，$n_j^{\left(r\right)}=n-\left(p_n-{\widehat{k}}_j+1\right)m$ ，${\beta }_j$ 和${\beta }_{j+1}$ 是未知的q维回归系数，$\left(i=n_j+1,\cdots ,n_{j+1}\right)$ 。区间${\widehat{ℐ}}_{\left(j\right)}={ℐ}_{{\widehat{k}}_j-1}\cup {ℐ}_{{\widehat{k}}_j},\left(j=1,\cdots ,\widehat{s}\right)$ 中是否存在变点等价于假设检验问题$H_0:{\beta }_j={\beta }_{j+1}\leftrightarrow H_1:{\beta }_j\ne {\beta }_{j+1}$ 。

记$N_j=n_j^{\left(r\right)}-n_j^{\left(l\right)}$ 表示区间${\widehat{ℐ}}_{\left(j\right)}={ℐ}_{{\widehat{k}}_j-1}\cup {ℐ}_{{\widehat{k}}_j}$ 的长度，$b_j={\left(2\log \log N_j+q\left(\log \log \log N_j\right)/2-\log \Gamma \left(q/2\right)\right)}^2/2\log \log N_j$ ，$c_j={\left(b_j/\left(2\log \log N_j\right)\right)}^{1/2}$ ，${\widehat{\sigma }}_j^2=\underset{\beta }{\min }{\displaystyle {\sum }_{i=n_j^{\left(l\right)}}^{n_j^{\left(r\right)}}{\left(y_i-x_i^{\text{T}}\beta \right)}^2}$为回归模型中残差平方和的估计值，其中$\Gamma \left(x\right)={\displaystyle {\int }_0^{+\infty }{t}^{x-1}{\text{e}}^{-t}\text{d}t}\left(x>0\right)$ 是Gamma函数。由文献[15]中的引理3.1.9可知，模型(5)中的回归系数是否存在变点的拟似然比检验统计量为：

$T_j=N_j\left({\widehat{\sigma }}_j^2-\underset{\beta }{\text{min}}{\displaystyle \sum }_{i=n_j^{\left(l\right)}}^{{\widehat{\zeta }}_j}{\left(y_i-x_i^{\text{T}}\beta \right)}^2-\underset{\beta }{\text{min}}{\displaystyle \sum }_{i={\widehat{\zeta }}_j+1}^{n_j^{\left(r\right)}}{\left(y_i-x_i^{\text{T}}\beta \right)}^2\right)/{\widehat{\sigma }}_j^2 \text{,}$ (6)

且当$T_j>b_j+2c_j\log \left(-2/\log \left(1-\alpha \right)\right)$ 拒绝原假设$H_0$ ，说明模型(5)中的回归系数存在一个变点${\zeta }_j\in \left\{n_j^{\left(l\right)},\cdots ,n_j^{\left(r\right)}\right\}$ 使得${\beta }_j\ne {\beta }_{j+1}$ ，其中$1-\alpha$ 是置信水平。变点${\zeta }_j$ 的相合估计量为：

${\widehat{\zeta }}_j=\arg \underset{n_j^{\left(l\right)}+q<k<n_j^{\left(r\right)}-q}{\min }\left[\underset{\beta }{\text{min}}{\displaystyle \sum }_{i=n_j^{\left(l\right)}}^k{\left(y_i-x_i^{\text{T}}\beta \right)}^2+\underset{\beta }{\text{min}}{\displaystyle \sum }_{i=k+1}^{n_j^{\left(r\right)}}{\left(y_i-x_i^{\text{T}}\beta \right)}^2\right].$

若接受原假设$H_0$ ，则删除第j段。记所有能够检测出的变点为${\tilde{\zeta }}_1<\cdots <{\tilde{\zeta }}_{\tilde{s}}$ 。由于每个${\tilde{\zeta }}_j,\left(j=1,\cdots ,\tilde{s}\right)$ 都仅用了2m个观测值得到，为了避免“过度估计(overestimation)”，提高变点估计的精准度，令${\tilde{n}}_0=0$ ，${\tilde{n}}_j={\tilde{\zeta }}_j$ ，$\left(j=1,\cdots ,\tilde{s}\right)$ ，${\tilde{n}}_{\tilde{s}+1}=n$ ，再次使用拟似然比检验(6)式去检测区间$\left\{{\tilde{n}}_{j-1}+1,\cdots ,{\tilde{n}}_{j+1}\right\}$ $\left(j=1,\cdots ,\tilde{s}\right)$ 中是否存在变点。若存在变点，变点估计量为：

${\widehat{a}}_j=\arg \underset{{\tilde{n}}_{j-1}+q<k<{\tilde{n}}_{j+1}-q}{\min }\left[\underset{\beta }{\text{min}}{\displaystyle \sum }_{i={\tilde{n}}_{j-1}+1}^k{\left(y_i-x_i^{\text{T}}\beta \right)}^2+\underset{\beta }{\text{min}}{\displaystyle \sum }_{i=k+1}^{{\tilde{n}}_{j+1}}{\left(y_i-x_i^{\text{T}}\beta \right)}^2\right]$ (7)

($j=1,\cdots ,\tilde{s}$ )

3. 实证分析

3.1. 指标选取

动力煤在电力、钢铁和建材等行业中发挥着重要作用，本文选取2021年12月22日至2024年11月6日秦皇岛动力煤(Q5500K)平仓价格(元/吨)的日数据作为研究对象(后续简述为秦皇岛动力煤价格)，共690个数据，其中前678个数据(2021年12月22日~2024年10月21日)作为训练集，后12个数据(2024年10月22日~2024年11月6日)作为测试集。基于以往学者的研究成果[16] [17]，综合考虑数据的公开性和可获取性，从替换能源、需求、运输成本和库存量4个维度选取5个指标的日数据作为影响动力煤价格的因素，数据均来源于前瞻数据库(https://d.qianzhan.com/)，具体如表1所示。

Table 1. Table of comprehensive influencing factors of thermal coal prices in Qinhuangdao

表1. 秦皇岛动力煤价格综合影响因素表

目标层	一级指标	二级指标	单位	符号表示
秦皇岛动力煤价格(y)	替换能源	WTI原油期货收盘价	美元/桶	$x_1$
		液化天然气市场价	元/吨	$x_2$
	需求	螺纹钢期货收盘价	元/吨	$x_3$
	运输成本	中国沿海煤炭运价综合指数	2011年9月1日=1000点	$x_4$
	库存量	秦皇岛港煤炭库存量	万吨	$x_5$

3.2. 变点检测与模型建立

由于线性回归模型的解释性强、应用性广，下面建立秦皇岛动力煤价格(y)与表1中5个指标的全局线性回归模型：

$y=4.6946x_1+0.0649x_2-0.0850x_3+0.6620x_4-0.3120x_5+363.2168$ (8)

图1给出了表1中各变量之间的相关系数，并将结果通过相关性热图可视化。

Figure 1. Correlation heat map between variables

图1. 各变量间的相关性热图

表2给出了全局线性回归模型(8)中各变量$x_i$ $\left(i=1,\cdots ,5\right)$ 的方差膨胀系数。

Table 2. Table of variance inflation factor for each variable in the global linear regression model (8)

表2. 全局线性回归模型(8)中各变量方差膨胀系数表

变量	$x_1$	$x_2$	$x_3$	$x_4$	$x_5$
方差膨胀系数VIF值	1.9848	2.2737	1.9797	1.8368	1.1211

由图1可以看出各变量之间不存在显著相关性，且由表2可以看出方差膨胀系数$VI{F}_i$ $\left(i=1,\cdots ,5\right)$ 均小于10，即全局线性回归模型(8)的各变量间不存在多重共线性，则利用这5个自变量能够较好地解释因变量的变化。接下来应用两阶段多变点检测与估计方法(SCAD)对模型(8)中的回归系数进行多变点检测与估计，根据文献[6]中段长的选取方法，这里选取段长$m=91$ ，由(7)式得到回归系数的变点位置为164和451，对应的时间分别为2022年8月25日和2023年11月7日，如图2中红色虚线所示。

Figure 2. Qinhuangdao thermal coal price time series diagram

图2. 秦皇岛动力煤价格时间序列图

图2给出了2021年12月22日至2024年10月21日秦皇岛动力煤价格的时间序列图。第一个变点(2022年8月25日)出现的原因是2022年7月28日中共中央政治局会议上提出“保交楼、稳民生”政策，使得房地产复工建设量增加，导致螺纹钢的价格上升，由于动力煤是生产螺纹钢的重要能源，进而影响动力煤价格的升高；第二个变点(2023年11月7日)出现的原因是2023年10月WTI原油价格下跌，使得替换能源动力煤的价格下降。

上述变点将整个观测数据分割为3个子段(第一子段为1~163，第二子段为164~450，第三子段为451~678)。应用R软件对每个子段分别建立线性回归模型，得到以下分段线性回归模型：

$y=\{\begin{array}{l}8.5740x_1+0.0096x_2+0.0917x_3+0.7592x_4+0.2729x_5-905.20001\le t<164,\\ 5.3133x_1+0.0815x_2-0.1572x_3+0.9313x_4-0.4152x_5+485.7751164\le t<451,\\ -4.8286x_1-0.0080x_2+0.1133x_3-0.0500x_4+0.2218x_5+781.0473451\le t\le 678,\end{array}$ (9)

表3给出了全局线性回归模型(8)与分段线性回归模型(9)的调整后拟合优度(Adjusted R²)的结果。

Table 3. Comparison table of adjusted R² between model (8) and model (9)

表3. 模型(8)与模型(9)调整后的拟合优度(Adjusted R²)对比表

模型	模型(8)全局线性回归模型	模型(9)第一段的线性回归模型	模型(9)第二段的线性回归模型	模型(9)第三段的线性回归模型
调整后的拟合优度(Adjusted R^2)	0.6481	0.8459	0.8573	0.7863

图3是由全局线性回归模型(8)和分段线性回归模型(9)得到的2021年12月22日至2024年10月21日秦皇岛动力煤价格拟合值与真实值的对比图，其中蓝色实线(-)表示真实值曲线，绿色虚线(--)表示由全局线性回归模型(8)得到的拟合值曲线，红色点状线(·*·)表示由分段线性回归模型(9)得到的拟合值曲线。

Figure 3. Comparison chart of the fitted value and the real value of Qinhuangdao thermal coal price obtained from model (8) and model (9)

图3. 由模型(8)和模型(9)得到的秦皇岛动力煤价格拟合值与真实值的对比图

由表3可以看出全局线性回归模型(8)的拟合优度$R^2<0.65$ ，而分段线性回归模型(9)三个子段的拟合优度R²均大于0.78。并且结合图3可以得到，分段线性回归模型(9)能利用$x_i$ $\left(i=1,\cdots ,5\right)$ 更精准地刻画秦皇岛动力煤价格(y)的走势。

3.3. 短期预测

分别利用全局线性回归模型(8)和分段线性回归模型(9)的最后一个子段对2024年10月22日~2024年11月6日的测试集进行短期预测。

Table 4. Statistical index table of the prediction results of model (8) and the last subsegment of model (9)

表4. 模型(8)与模型(9)最后一个子段预测结果的统计指标表

统计指标	由全局线性回归模型(8)得到的预测值	由分段线性回归模型(9)得到的预测值
均方根误差(RMSE)	188.6133	34.4454
平均绝对误差(MAE)	188.3771	32.1677
平均绝对百分比误差(MAPE)	22.1268%	3.7803%
对称平均绝对百分比误差(SMAPE)	19.9183%	3.7002%

图4是由全局线性回归模型(8)和分段线性回归模型(9)最后一个子段得到的预测值与真实值的对比图，其中蓝色实心圆表示真实值散点，绿色实心正方形表示由模型(8)得到的预测值散点，红色实心三角表示由模型(9)最后一个子段得到的预测值散点。

Figure 4. Comparison chart of the predicted value and the true value obtained from the model (8) and the last subsegment of the model (9)

图4. 由模型(8)和模型(9)最后一个子段得到的预测值与真实值的对比图

从表4和图4可以看出在对秦皇岛动力煤价格的短期预测方面，分段线性回归模型(9)最后一个子段的预测效果优于全局线性回归模型(8)。说明应用线性回归模型刻画秦皇岛动力煤价格的走势前进行变点检测与估计是至关重要的，应用分段线性回归模型能够更精准地刻画秦皇岛动力煤价格的走势。

4. 结论

本文将两阶段多变点检测与估计方法应用到秦皇岛动力煤价格(y)与WTI原油期货收盘价($x_1$ )、液化天然气市场价($x_2$ )、螺纹钢期货收盘价($x_3$ )、中国沿海煤炭运价综合指数($x_4$ )、秦皇岛港煤炭库存量($x_5$ )的全局线性回归模型中，检测并估计了回归系数的多变点，通过变点数量和位置的估计结果建立了更精准的分段线性回归模型来刻画秦皇岛动力煤价格的走势，说明在进行数据分析和统计建模之前进行变点检测与估计是至关重要的，可以帮助我们建立更精准的模型。

基金项目

新疆维吾尔自治区自然科学基金项目(2023D01A37)。

NOTES

^*通讯作者。