河南科技大学数学与统计学院，河南 洛阳

收稿日期：2021年11月27日；录用日期：2021年12月11日；发布日期：2021年12月27日

摘要

聚类是时间序列数据挖掘的重要任务之一。本文基于有限混合MS-GARCH模型，提出一种时间序列聚类方法。利用贝叶斯马尔科夫链蒙特卡洛模拟方法，克服路径依赖的困难，给出了模型参数的估计。最后，选取23家中国上市公司股票数据进行实证分析，验证了所提方法的有效性。

关键词

时间序列聚类，有限混合模型，MCMC算法，MS-GARCH模型

Time Series Clustering with MS-GARCH Mixtures

Lin Wang, Xiaoquan Ding^*

School of Mathematics and Statistics, Henan University of Science and Technology, Luoyang Henan

Received: Nov. 27^th, 2021; accepted: Dec. 11^th, 2021; published: Dec. 27^th, 2021

ABSTRACT

Clustering is one of the important tasks of time series data mining. In this paper, we propose a novel time series clustering method based on the finite mixture MS-GARCH model. By using Bayesian Markov chain Monte Carlo simulation methods to overcome the difficulty of full path dependence, we estimate the model parameters. Finally, the empirical analysis of stock data of 23 Chinese listed companies verifies the effectiveness of our proposed method.

Keywords:Time Series Clustering, Finite Mixture Model, MCMC Algorithm, MS-GARCH Model

Copyright © 2021 by author(s) and Hans Publishers Inc.

This work is licensed under the Creative Commons Attribution International License (CC BY 4.0).

http://creativecommons.org/licenses/by/4.0/

1. 引言

时间序列数据是一类非常重要的复杂数据，广泛存在于经济、金融、管理、医疗卫生、气象、水文和工程等领域。由于其高维、高冗余和随时间动态变化的特点，相比直接分析原始数据，合理利用前人建立的时间序列模型，在研究中更能容易刻画序列的变化规律。

广义自回归条件异方差(GARCH)模型是由Bollerslev [1] 在Engle [2] 研究的基础上对自回归条件异方差(ARCH)模型推广得到的。该模型以可变条件方差描述序列波动特征，条件异方差依赖于给定信息和滞后条件方差，具有更灵活的滞后结构，能更好地体现序列的记忆特征。1990年，Lamoureux和Lastrages [3] 指出条件方差可能存在结构转换，如果不考虑这种结构变化，可能造成波动持续性的错误估计。受其启发，Cai [4] 和Hamilton，Sumsel [5] 将Markov状态切换和ARCH效应结合，提出MS-ARCH模型；Gray [6] 则进一步把马尔可夫状态切换引入GARCH模型，建立MS-GARCH模型。随后，许多学者对MS-GARCH模型做了不同的改进，并开展了深入研究 [7] - [13]。

聚类是把一个数据对象的集合划分成若干簇(子集)，使簇内对象彼此相似、簇间对象不相似的过程。面对庞杂而日益增长的时间序列数据，聚类是时间序列数据挖掘的重要任务之一。近年来，涌现出许多时间序列聚类方法，大体上可以分为基于原始数据的聚类、基于数据特征的聚类和基于数据模型的聚类等三种 [14] [15] [16]。

基于混合模型的聚类是一种有效、可解释、适应性强的统计分析方法 [17] [18] [19]，它假定要聚类的数据来自一个有限混合的概率分布，通过计算每个数据在此概率分布下的后验概率，取后验概率最大的成分作为聚类结果。近年来，该方法也被广泛用于时间序列数据的聚类 [20] - [27]。

在基于混合模型的聚类研究中，参数估计是关键工作，通常有期望最大化(EM)和Markov Chain Monte Carlo (MCMC)两种算法。EM算法基于极大似然估计，利用迭代优化估计模型的参数；MCMC算法则基于Bayes估计，利用MCMC模拟，解决了时间序列数据的高维、高冗余、路径依赖等不适用于EM算法的问题。研究表明，MCMC算法能很好地解决高维数据计算复杂的难题，拟合效果优于EM算法 [12] [13] [28] [29] [30]。

本文基于有限混合MS-GARCH模型，利用Bayes理论的数据增补方法和MCMC模拟，探讨时间序列数据的聚类问题。文 [21] 利用有限混合GARCH模型对美国上市公司股票收益数据进行聚类，最终确定按波动率的高、中、低分为三类，但并没有考虑机制切换问题。就我们所知，未见文献把有限混合MS-GARCH模型用于时间序列数据的聚类分析。

2. 基于MS-GARCH模型的时间序列聚类

2.1. 基于有限混合模型的聚类

设数据集 $\left\{y_1,y_2,\cdots ,y_I\right\}$ 是来自由K个成分构成的有限混合概率分布

$f\left(y|\rho ,\Theta \right)={\displaystyle \underset{k=1}{\overset{K}{\sum }}{\rho }_{k}f\left(y|{\Theta }_k\right)}$ (2.1)

的样本，其中 $\rho =\left({\rho }_1,{\rho }_2,\cdots ,{\rho }_K\right)$ 是混合权重， $\Theta =\left({\Theta }_1,{\Theta }_2,\cdots ,{\Theta }_K\right)$ 是成分参数， $f\left(y|{\Theta }_k\right)$ 是第k个成分的概率分布，权重 $\rho$ 满足：

$0\le {\rho }_k\le 1\left(k=1,2,\cdots ,K\right),{\displaystyle \underset{k=1}{\overset{K}{\sum }}{\rho }_k}=1.$

基于模型的聚类方法把混合模型的成分与聚类的簇相对应，成分数就是簇数。设 $z_i$ 代表数据 $y_i$ 的簇标记， $z_i=k$ 表示 $y_i$ 属于第k个簇。若模型的参数 $\rho ,\Theta$ 已知，则由Bayes定理

$\mathbb{P}\left\{z_i=k|y_i\right\}=\frac{{\rho }_{k}f\left(y_i|{\Theta }_k\right)}{{\displaystyle \underset{j=1}{\overset{K}{\sum }}{\rho }_{j}f\left(y_i|{\Theta }_j\right)}}.$

若

${\lambda }_i=\underset{k\in \left\{1,2,\cdots ,K\right\}}{\arg \min }\mathbb{P}\left\{z_i=k|y_i\right\},$

则 $y_i$ 归属于第 ${\lambda }_i$ 个簇。

2.2. MS-GARCH模型

称

$\{\begin{array}{l}{y}_t=\sqrt{h_t}\cdot {\epsilon }_t,\\ h_t={\omega }_{s_t}+{\displaystyle \underset{i=1}{\overset{p}{\sum }}{\alpha }_{i,s_t}{y}_{t-i}^2}+{\displaystyle \underset{j=1}{\overset{q}{\sum }}{\beta }_{j,s_t}{h}_{t-j}},\end{array}$

为具有Markov状态切换的GARCH(p, q)模型，简记为MS-GARCH，其中 ${\left\{s_t\right\}}_{t\in \mathbb{N}}$ 是状态空间为 $\left\{1,2,\cdots ,\mathcal{S}\right\}$ 、转移概率矩阵为 $\eta =\left({\eta }_{ij}\right)$ 的遍历Markov链， $\left\{{\epsilon }_t\right\}$ 独立同分布，对任意 $s=1,2,\cdots ,\mathcal{S}$，$i=1,2,\cdots ,p$，$j=1,2,\cdots ,q$，有 ${\omega }_s>0$，${\alpha }_{i,s}\ge 0$，${\beta }_{j,s}\ge 0$。

2.3. 基于MS-GARCH模型的时间序列聚类

设 $\left\{y_{i,t},i=1,2,\cdots ,I;t=1,2,\cdots ,T\right\}$ 是I个长度为T的时间序列组成的观察数据，记

$y_i=\left(y_{i,1},y_{i,2},\cdots ,y_{i,T}\right),i=1,\cdots ,I.$

假设 $\left\{y_1,y_2,\cdots ,y_I\right\}$ 来自有限混合概率分布(2.1)， $f\left(y|{\Theta }_k\right)$ 是第k个MS-GARCH(p, q)模型

$\{\begin{array}{l}{y}_t=\sqrt{h_t}\cdot {\epsilon }_t,\\ h_t={\omega }_{s_t^k}^k+{\displaystyle \underset{i=1}{\overset{p}{\sum }}{\alpha }_{i,s_t^k}^k{y}_{t-i}^2}+{\displaystyle \underset{j=1}{\overset{q}{\sum }}{\beta }_{j,s_t^k}^k{h}_{t-j}},\end{array}$ (2.2)

的概率密度，其中 ${\left\{s_t^k\right\}}_{t\in \mathbb{N}}$ 是状态空间为 $\left\{1,2,\cdots ,\mathcal{S}\right\}$ 、转移概率矩阵为 ${\eta }^k=\left({\eta }_{ij}^k\right)$ 的遍历Markov链。

为叙述方便，下面只对两状态的MS-GARCH(1, 1)模型讨论，其它情况类似。此时，模型(2.2)的参数为 ${\Theta }_k=\left({\theta }_k,{\eta }^k\right)$，其中

${\theta }_k=\left({\omega }_1^k,{\omega }_2^k,{\alpha }_1^k,{\alpha }_2^k,{\beta }_1^k,{\beta }_2^k\right),{\eta }^k=\left({\eta }_{11}^k,{\eta }_{12}^k,{\eta }_{21}^k,{\eta }_{22}^k\right),$

对 $j=1,2$，有

${\omega }_j^k>0,{\alpha }_j^k\ge 0,{\beta }_j^k\ge 0,{\eta }_{j1}^k\ge 0,{\eta }_{j2}^k\ge 0,{\alpha }_j^k+{\beta }_j^k<1,{\eta }_{j1}^k+{\eta }_{j2}^k=1.$

当 ${\epsilon }_t$ 服从标准正态分布时，联合概率密度

$f\left(y|{\Theta }_k\right)={\displaystyle \underset{t=1}{\overset{T}{\prod }}\frac{1}{\sqrt{2\pi h_t}}\exp }\left(-\frac{y_t^2}{2h_t}\right),$ (2.3)

其中 $h_t$ 由(2.2)的第二式递推给出。

定义随机变量 $z_i\in \left\{1,2,\cdots ,K\right\}$，$i=1,2,\cdots ,I$，其概率分布

$\mathbb{P}\left\{z_i=k|\rho \right\}={\rho }_k.$

记

$z=\left(z_1,z_2,\cdots ,z_I\right),S_T^k=\left(s_1^k,s_2^k,\cdots ,s_T^k\right),S_{-t}^k=\left(s_1^k,s_2^k,\cdots ,s_{t-1}^k,s_{t+1}^k,\cdots ,s_T^k\right).$

由于GARCH模型中条件方差存在路径依赖问题，参数估计的极大似然方法不再适用，我们利用MCMC方法对隐变量z和参数 $\rho ,\Theta$ 进行Gibbs抽样。为此，需要计算 $\left(z,\rho ,\Theta \right)$ 的后验分布。取先验分布 $\pi \left(z,\rho ,\Theta \right)$，使得

$\pi \left(z,\rho ,\Theta \right)=\pi \left(z|\rho \right)\pi \left(\rho \right)\pi \left(\Theta \right),$ (2.4)

其中

$\pi \left(z|\rho \right)={\displaystyle \underset{i=1}{\overset{I}{\prod }}{\rho }_{z_i}}={\displaystyle \underset{k=1}{\overset{K}{\prod }}{\rho }_k^{x_k}},x_k=\#\left\{z_i=k,i=1,2,\cdots ,I\right\},$ (2.5)

$\pi \left(\Theta \right)={\displaystyle \underset{k=1}{\overset{K}{\prod }}\pi }\left({\Theta }_k\right).$ (2.6)

利用Bayes定理，由(2.4)~(2.6)可得 $\left(z,\rho ,\Theta \right)$ 的后验分布

$p\left(z,\rho ,\Theta |y\right)\propto f\left(y|z,\rho ,\Theta \right)\pi \left(z,\rho ,\Theta \right)=\pi \left(\rho \right){\displaystyle \underset{k=1}{\overset{K}{\prod }}{\rho }_k^{x_k}\pi \left({\Theta }_k\right)}{\displaystyle \underset{i=1}{\overset{I}{\prod }}f}\left(y_i|{\Theta }_{z_i}\right).$ (2.7)

下面依次讨论z、 $\rho$ 和 $\Theta$ 的抽样。

2.3.1. 簇标记z的抽样

由(2.7)可得z的后验分布

$p\left(z_1,z_2,\cdots ,z_I|\rho ,\Theta ,y\right)\propto {\displaystyle \underset{i=1}{\overset{I}{\prod }}f\left(y_i|{\Theta }_{z_i}\right){\rho }_{z_i}}={\displaystyle \underset{i=1}{\overset{I}{\prod }}p}\left(z_i{|\rho ,\Theta ,y}_i\right).$ (2.8)

于是， $z_i|\rho ,\Theta ,y_i$ 服从多项分布，满足

$P\left(z_i=k|\rho ,\Theta ,y_i\right)\propto f\left(y_i|{\Theta }_k\right){\rho }_k,k=1,2,\cdots ,K.$ (2.9)

2.3.2. 混合概率ρ的抽样

由(2.7)可得 $\rho$ 的后验分布

$p\left(\rho |z,\Theta ,y\right)=p\left(\rho |z\right)\propto \pi \left(\rho \right){\displaystyle \underset{k=1}{\overset{K}{\prod }}{\rho }_k^{x_k}}.$ (2.10)

取 $\pi \left(\rho \right)$ 为Dirichlet分布 $\mathcal{D}\left(a_1,a_2,\cdots ,a_K\right)$ ：

$\pi \left(\rho |a\right)=\frac{{\displaystyle \underset{k=1}{\overset{K}{\prod }}\Gamma \left(a_k\right)}}{\Gamma \left({\displaystyle \underset{k=1}{\overset{K}{\sum }}{a}_k}\right)}{\displaystyle \underset{k=1}{\overset{K}{\prod }}{\rho }_k^{a_k-1}}\text{,}{a}_k>0,{\rho }_i\ge 0,{\displaystyle \underset{k=1}{\overset{K}{\sum }}{\rho }_k=1},$

则由(2.10)可知，后验分布 $p\left(z|\rho \right)$ 为Dirichlet分布 $\mathcal{D}\left(a_1+x_1,a_2+x_2,\cdots ,a_K+x_K\right)$。

2.3.3. 簇参数Θ的抽样

由(2.6)和(2.7)可知， $\Theta$ 的后验分布

$p\left(\Theta |z,\rho ,y\right)=p\left(\Theta |z,y\right)={\displaystyle \underset{k=1}{\overset{K}{\prod }}p\left({\Theta }_k|{\tilde{y}}_k\right)}\propto {\displaystyle \underset{k=1}{\overset{K}{\prod }}\pi \left({\Theta }_k\right)}{\displaystyle \underset{i\in I_k}{\prod }f}\left(y_i|{\Theta }_k\right),$ (2.11)

其中 $I_k=\left\{i|z_i=k\right\}$，${\tilde{y}}_k=\left\{y_i|i\in I_k\right\}$。于是，对 $\Theta$ 的抽样可通过对每个 ${\Theta }_k$ 分别独立抽样来实现。取先验分布 $\pi \left({\Theta }_k\right)$，使得

$\pi \left({\Theta }_k\right)=\pi \left({\theta }_k\right)\pi \left({\eta }^k\right).$

下面依次讨论 $S^k$ 、 ${\eta }^k$ 和 ${\theta }_k$ 的抽样。

1) 状态标记 $S_T^k$ 的抽样

状态变量 $s_t^k$ 是一个隐Markov过程，根据Markov链的性质， $s_t^k$ 的概率需要考虑 $s_{t-1}^k$ 和 $s_{t+1}^k$ 的影响。由

$\begin{array}{c}p\left(s_t^k|S_{-t}^k,{\theta }_k,{\eta }^k,{\tilde{y}}_k\right)=\frac{p\left(s_t^k,S_{-t}^k,{\theta }_k,{\eta }^k,{\tilde{y}}_k\right)}{p\left(S_{-t}^k,{\theta }_k,{\eta }^k,{\tilde{y}}_k\right)}=\frac{p\left({\tilde{y}}_k|S_T^k,{\theta }_k,{\eta }^k\right)p\left(S_T^k,{\theta }_k,{\eta }^k\right)}{p\left({\tilde{y}}_k|S_{-t}^k,{\theta }_k,{\eta }^k\right)p\left(S_{-t}^k,{\theta }_k,{\eta }^k\right)}\\ =\frac{f\left({\tilde{y}}_k|S_T^k,{\theta }_k,{\eta }^k\right)p\left(S_T^k|S_{-t}^k,{\theta }_k,{\eta }^k\right)}{f\left({\tilde{y}}_k|S_{-t}^k,{\theta }_k,{\eta }^k\right)},\end{array}$

可知

$\mathbb{P}\left(s_t^k=s|S_{-t}^k,{\theta }_k,{\eta }^k,{\tilde{y}}_k\right)\propto {\eta }_{s_{t-1},s}^k{\eta }_{s,s_{t+1}}^k{\displaystyle \underset{i\in I\left(k\right)}{\prod }{\displaystyle \underset{\tau =t}{\overset{T}{\prod }}{h}_{i\tau }^{-0.5}\exp }}\left(-\frac{y_{i\tau }^2}{2h_{i\tau }}\right).$ (2.12)

由于 $s_t^k$ 只有两个状态，由(2.12)可知对 $s_t^k$ 的抽样相当于从一个Bernoulli分布中抽样。

2) 转移概率 ${\eta }^k$ 的抽样

给定 ${\eta }^k$ 的先验分布 $\pi \left({\eta }^k\right)$，当 ${\eta }^k$ 独立于观察变量 ${\tilde{y}}_k$ 和参数 ${\theta }^k$ 时， ${\eta }^k$ 的后验分布

$p\left({\eta }^k|S_T^k,{\theta }^k,{\tilde{y}}_k\right)=p\left({\eta }^k|S_T^k\right)\propto \pi \left({\eta }^k\right){\displaystyle \underset{t=1}{\overset{T}{\prod }}{\eta }_{s_{t-1},s_t}^k}.$ (2.13)

取 $\pi \left({\eta }^k\right)=\pi \left({\eta }_{11}^k,{\eta }_{12}^k\right)\pi \left({\eta }_{21}^k,{\eta }_{22}^k\right)$，当 $\pi \left({\eta }_{11}^k,{\eta }_{12}^k\right)$ 和 $\pi \left({\eta }_{21}^k,{\eta }_{22}^k\right)$ 分别为Dirichlet分布 $\mathcal{D}\left(b_1,b_2\right)$ 和 $\mathcal{D}\left(c_1,c_2\right)$ 时，由(2.13)可知后验分布 $p\left({\eta }_{11}^k,{\eta }_{12}^k|S_T^k\right)$ 和 $p\left({\eta }_{21}^k,{\eta }_{22}^k|S_T^k\right)$ 分别为 $\mathcal{D}\left(b_1+n_{11},b_2+n_{12}\right)$ 和 $\mathcal{D}\left(c_1+n_{21},c_2+n_{22}\right)$，其中 $n_{gl}$ 为从状态g转移到状态l的次数。

3) GARCH参数 ${\theta }^k$ 的抽样

由于GARCH模型中当期条件方差受到前期条件方差影响的递推结构，各参数之间不独立，参数的后验分布没有合适的共轭先验，因此不能直接使用常规的Gibbs抽样。若取 ${\theta }^k$ 的先验分布为常数，则 ${\theta }^k$ 的后验分布

$p\left({\theta }^k|S_T^k,{\eta }^k,{\tilde{y}}_k\right)\propto {\displaystyle \underset{i\in I\left(k\right)}{\prod }{\displaystyle \underset{t=1}{\overset{T}{\prod }}{h}_{it}^{-0.5}\exp }}\left(-\frac{y_{it}^2}{2h_{it}}\right).$ (2.14)

利用griddy-Gibbs抽样方法 [12] [28]，根据(2.14)可对 ${\theta }^k$ 进行抽样。

3. 实证分析

本文选取美的集团、长安汽车、格力电器、海信家电、三全食品、比亚迪、东风汽车、宇通客车、上汽集团、江淮汽车、光明乳业、海尔智家、伊利股份、长城汽车、飞科电器、中信证券等23家上市公司2017年12月至2020年7月的所有日收益数据进行模拟实验(数据来源：网易财经数据库)，以验证本算法的有效性。

3.1. 数据预处理

1) 由于上市公司的股权变更、违规调查等外部因素造成的停牌现象，使数据中存在少数NULL值。对此执行删除操作，不影响原序列的特征信息。

2) 实验中选取上市公司股票的对数收益率数据进行计算，以消除时间不连续等负面影响，

3) 本文的MS-GARCH模型只考虑条件方差的波动特征，取对数收益率后，将序列减去各自样本均值以消除模型中均值参数的影响，有助于简化迭代计算过程。

3.2. 数据概览

我们从中选取两家公司的收益率序列数据，画出它们的时序图，如图1所示。可以看出各个序列都有明显的波动聚集特征，即某些时间段有较稳定的波动，某些时间段有较剧烈的震荡。

图1. 对数收益率图

3.3. 模型选择

选用两状态的MS-GARCH(1,1)混合模型进行聚类分析。在聚类数的选择方面，参照文 [21] 的结论，我们把股票收益率数据按照波动率的高、中、低划分为三类，分别对应投资风险高、投资风险中等和投资风险低的三类股票。

3.4. 聚类结果

首先，我们计算得出了23条股票数据(包含均值、标准差等)的描述性统计，如表1所示。这些数据确实存在一定差异性。例如，全部数据的平均峰度为6.3084，但范围从3.6448到11.9634。

表1. 描述性统计

算法程序采用Matlab软件实现。抽样算法设定循环3000次，其中前1000次作为预烧阶段，保留后2000次作为抽样结果。

对于参数初始化，我们根据每条数据方差的差异性，对混合概率和簇类归属作了简单的初始化设置。其中，11家公司(美的集团、格力电器、上汽集团、三元股份、海尔智家、伊利股份、长江电力、中国平安、工商银行、中国石油、建设银行)归属第一类，8家公司(海信家电、三全食品、比亚迪、宇通客车、光明乳业、长城汽车、飞科电器、中信证券)归属第二类，4家公司(长安汽车、贝因美、东风汽车、江淮汽车)归属第三类。初始混合概率的设置如表2所示。

表2. 初始混合概率

23条时间序列各自包含两个状态的参数，即共有46组不相同的参数。我们在程序中对这些模型参数分别设置了不同的初始值，同时它们又必须满足GARCH模型的平稳性条件。由于Gibbs抽样的收敛性质，只需要粗略估计其大致范围即可。除此之外，griddy-Gibbs计算中设置格点数为50 (格点的数量过大或过小都会对估计结果和运行时间产生负面影响)。最后，要注意计算过程中数值不能超过Matlab所允许的范围。对于GARCH参数初始化，我们以美的集团为例给出了详细情况，如表3所示。

表3. GARCH参数初始化信息

通过计算，我们得到了该组数据两个状态的参数后验概率密度，抽样情况良好。根据Gibbs抽样原理，经过迭代后的抽样值会在真实值附近波动。美的集团的GARCH参数的后验概率密度图，如图2所示。

图2. GARCH参数后验概率密度图

由实验结果中状态转移概率后验均值与状态转移概率初始值的对比，可以看到股票处于低波动状态的持续性长于高波动状态，表示该股票大多数时间处于低幅度变动，不会频繁出现剧烈震荡。这类相对稳定的股票对于投资者来说风险偏低，若近期走势良好，值得投资。相反，若实验所得的后验概率中，高波动状态所占比重较高，那么说明股票长期处于数值较大的震荡，不利于进行投资，投资者需结合实际情况慎重分析。状态转移概率初始值与后验均值，如表4所示。

表4. 状态转移概率

在最终实验结果中，根据混合概率后验均值的情况可以看出在我们提供数据的基础上，三个混合概率非常接近。其中，第二类混合后验概率与初始簇类混合概率相比变化极小，也就意味着我们对第二类的初始分类状况较好。而第一类和第三类经过迭代抽样后，都产生了不同程度的变化，如表5所示。

表5. 混合概率

根据GARCH参数可对三个簇进行解释。三类中持久性最高的为第一组(1状态 ${\alpha }_1^1+{\beta }_1^1$ 估计为0.4217，2状态 ${\alpha }_2^1+{\beta }_2^1$ 估计为0.4621)，第一组与第三组的持久性更加接近(第三组1状态 ${\alpha }_1^3+{\beta }_1^3$ 估计为0.4062，2状态 ${\alpha }_2^3+{\beta }_2^3$ 估计为0.4320)，即第一组代表持久性最高的股票类型，第二组代表持久性最低的股票类型，第三组代表持久性中等的股票类型。其中三组的滞后因子彼此相近(1状态为0.0273左右，2状态为0.2160左右)，而差异最明显的地方则是条件方差的自回归参数(第一组两状态分别为0.3942和0.2462，第二组分别估计为0.2974和0.2375，第三组估计分别为0.3790和0.2148)。其次，三类的无条件方差具有一定差异，第一类与第二类相对接近(第一组两状态估计分别为3.05e−05和1.61e−03，第二组估计分别为1.75e−04和1.67e−03)，第三类最小(第三类估计分别为4.61e−05和9e−04)。但由于数值较小，还不足以对各自模型的持久性产生决定性的影响。GARCH参数的后验均值结果，如表6所示。

表6. GARCH参数的后验均值

在23家上市公司股票数据聚类结果中，有5家公司属于第一组，9家公司属于第二组，9家公司属于第三组。与初始化值相比较，有9家公司在初始化中就被正确区分，其余14家公司经过迭代抽样产生了类别变化，说明我们基于Gibbs抽样的聚类算法有明显的纠正效果。除此之外，我们还展示了各股票数据各自模型参数 $\hat{\omega }$ 、 $\hat{\alpha }$ 和 $\hat{\beta }$ 的后验均值，在取值上均满足GARCH模型的平稳性条件，如表7所示。

对于Gibbs抽样，先验分布的设定可能影响算法对于参数的估计、状态识别和聚类结果。在实验过程中，我们尝试修改各参数初始值，观察结果变化，发现初始值的改变并不会对最终结果产生太大影响，这与理论分析是一致的。

表7. 实验结果

4. 结论

本文基于有限混合MS-GARCH模型，提出一种时间序列聚类方法。在新息服从标准正态分布的情况下，利用MCMC方法，给出了模型参数的估计，再通过计算后验概率得到聚类结果。最后，选取23家中国上市公司股票的收益率数据进行实证分析，验证了所提方法的有效性。本文的方法也可用于新息服从t-分布、广义误差分布等其它分布的情况。本文的聚类数目是根据经验和理论分析事先指定的，聚类数目未知的聚类问题将在未来的工作中讨论

文章引用

王 琳,丁孝全. 基于MS-GARCH模型的时间序列聚类
Time Series Clustering with MS-GARCH Mixtures[J]. 统计学与应用, 2021, 10(06): 1071-1082. https://doi.org/10.12677/SA.2021.106114

参考文献