基于小波变换多维时间序列最新变化点检测

陈雪文

摘 要：現代数据科学中存在大量的多维时间序列数据，检测多维时间序列中的最新变化点对于短期预测很重要。一种改进的方法被提出，以检测此类多维时间序列数据中最新变化点。通过使用小波变换，将多维时间序列中的变化点检测问题转化为相对较容易的多维面板数据中的变化点检测问题。该方法旨在跨时间序列合并信息，以便优先推断多个序列中同一时间点的最新变化。通过对每个时间序列的输出进行后处理，获取最新变化点的时间集及该时间最近变化点的序列集。最后使用R软件以模拟数据和S&P500数据实验证明了该方法的有效性。

关键词：多维时间序列;小波变换;惩罚成本;最新变化点

中图分类号：O213       文献标识码：A

Latest Change Point Detection of Multidimensional

Time Series Based on Wavelet Transform

CHEN Xuewen

（College of Science， Hohai University， Nanjing， Jiangsu 211100，China）

Abstract：There are a lot of multidimensional time series data in modern data science， and detecting the latest change points in multidimensional time series is very important for shortterm forecasting. An improved method is proposed to detect the latest change points in such multidimensional time series data. By using wavelet transform， the problem of detecting change points in the secondorder structure of multidimensional time series is transformed into a relatively easy problem of detecting change points in multidimensional panel data. This method aims to merge information across time series in order to prioritize the latest changes at the same time point in multiple series. By postprocessing the output of each time series， the time set of the latest change point and the sequence set of the latest change point in time are obtained. Finally， the use of R software to simulate data and S&P500 data experiments proved the effectiveness of the method.

Key words：multidimensional time series; wavelet transform; penalty cost; latest change pointD20F8BCA-7CEA-40EA-A250-51F7BF3BCEA4

在许多现代应用中，会产生大量多维时间序列数据。多维时间序列数据中变化点的检测至关重要，现有的变化点检测方法大致可以分为两类：回顾性（后验）方法和在线检测方法。后验方法使用所有的观测值，检测已经发生的变化点[1];在线检测方法对新到达的数据进行自适应分析[2]。

在多个时间序列中检测变化点会出现单个时间序列不存在变化点的问题，因此可以采用单变量变化点方法：一种是分别分析每个时间序列，另一种是汇总时间序列并分析所得的单变量序列。每种方法都有其缺点：前者在检测变化点时忽略了不同时间序列可能在相似时间具有变化点的信息;如果来自影响少量序列的变化点的信号在汇总时被其余序列中的噪声淹没，则后者性能可能较差。

分析多维时间序列数据的另一种方法是将数据视为具有多变量观测值的单个时间序列。对分段中的多元数据建模，并允许该模型以适当的方式在不同部分之间有变化。如Lavielle[3]将数据建模为均值可能因分段而异的多元高斯模型，Matteson[4]提出一种非参数方法以检测多元数据中的多个变化。但是，就像汇总数据一样，如果变化仅影响时间序列的一小部分，则这些方法可能表现较差。

针对仅影响序列子集的变化点，Cho[5]基于CUSUM统计量[6]，提出一种检测多个变化点存在和位置信息的方法。保留所有序列的CUSUM值能够证明在给定时间点是否发生变化，但降低了其他时间序列值的权重。同样，Xie[7]引入广义似然比检测仅影响序列子集的单个公共变化点。该检验需要估计受变化点影响的序列比例，此估算会影响每个序列发生变化可能性权重，如给予变化点可能性较大的序列较大的权重，而给予其他序列较低的权重。

为估计每个时间序列的最新变化点，采用了不同的方法。该方法旨在将多维时间序列通过小波变换，转化为多维面板数据，对面板数据序列进行惩罚成本分析，最后输出每个序列的最新变化点。

某些因素会影响一个时间序列的结构变化，也可能会影响其他时间序列的结构变化。但是，并非所有时间序列都可能在同一时间存在一个变化点。因此，对这些分析的输出进行后处理，将时间序列分组，同一组共享公共变化点。该方法能够准确估计每个时间序列的最新变化点的位置，这有助于提高序列短期预测准确性。

1 单变量时间序列

假设时间序列为y1：T，变化点的数量为m，位置为τ=（τ1，…，τm），且τ0=0，τm+1=T。将时间序列在变化点处分段，对段内数据建模，每一段的成本要与极大似然值成比例[8-9]，以将该模型拟合到数据段[10]，从而得出成本。假设段中数据的模型具有一定密度f（yθ）且独立同分布，其中θ是段特定参数，那么段ys：t的成本[11]定义为：

为检测均值的变化，一个简单的模型是分段内数据独立同高斯分布，方差均为σ2，分段特定均值为θ，那么成本函数如下：

若分段内的均值是关于时间的线性函数，且这个线性模型在不同的段中不同。表示为θ=θ1，θ2，其中θ1是截距，θ2是斜率。假设此模型的噪声是独立同高斯分布，那么段成本函數为：

某些时间序列存在明显的异常值，为对异常值具有鲁棒性，采用Fearnhead（2017）[12]中方法，使用如下分段成本函数：

若异常值的残差与分段均值的距离大于2倍标准偏差，则此成本会限制异常值的影响。对于剩余方差σ2，基于差分时间序列的中值绝对偏差，使用一个简单而鲁棒的σ2估计量[13]。

为分段数据并找到变化点，希望将所有分段的成本总和降至最低。为避免过度拟合，为每个分段添加一个惩罚项β（β>0）。因此，对数据进行分段，即解决以下优化问题：

β的值越高，意味着检测到的变化点越少。根据BIC准则，若段特定参数的维数为p，则β=（p+1）log T。定义F（t），t=1，2，…，T如下：

其中0=τ0<τ1<…<τm<τm+1=t，因此，F（t）是分段y1：t的最小成本。函数F（t）能够使用PELT[14]中的算法有效计算：

假设最新的变化点在时间r时，最低成本为G（r）。G（r）与F（r）有关，因为F（r）是分段y1：r的最小成本。

可知G（0）=C（y1：T），因为已经对最近变化点之前的变化点的数量和位置进行了优化，因此可认为G（r）的大小与边侧似然函数有关。最近变化点的估计由arg min r∈0，…，T-1G（r）得出。若最近变化点在r=0，则该时间序列没有变化点。

2 多维时间序列

2.1 小波变换

Nason[15]最初提出使用小波作为非平稳时间序列的基础，类似于平稳过程经典Cramer表示中的傅立叶指数。使用最简单的小波系统，Haar小波：

其中j∈{-1，-2，…}是小波尺度因子，l∈Z是位置参数。

假设时间序列的维度为N，长度为T，yi，t表示第i条时间序列，截至时间t的最近变化点为η（t），那么yηti，t小波系数为：

其中Lj是尺度为j的离散小波向量ψj=（ψj，0，…，ψj，Lj-1）T的支持长度，Lj=M2-j（M是定值，取决于小波族）。序列yηti，t的小波周期图序列和交叉周期图序列分别为：D20F8BCA-7CEA-40EA-A250-51F7BF3BCEA4

2.2 多维面板数据检测最近变化点

为检测面板数据中的最近变化点集，令Gi（η）是第i个序列最近变化点在η的最小成本函数。假设有K个公共最近变化点，η1：K=（η1，…，ηK）。对于第k个最近变化点ηk，存在一个集合Ιk1，2，…，N，使得对于所有i∈Ιk的序列，最近变化点都是ηk。由此，集合I1：K可将多维序列分组。

为得到η1：K的估值，需最小化每个序列的成本总和：

CK=min I1，…，IK，η1，…，ηK∑Kk=1∑i∈IkGi（ηk）     （13）

对于K值的选取，使用式（13）的惩罚形式，针对一个范围中不同的K值，解决式（13）的优化问题，选择使得下式达到最小的K值：

CK+Nlog 2K+Klog 2T      （14）

使用BIC选择方法的惩罚项不利于添加最新的变化点，因此，在計算Gi（η）时，添加一个比BIC选择略低的惩罚项β=（p+1/2）log T。在没有变化点的模拟数据中，随着η的增大，较小的β值会产生更小的G（η），表明选择较小的β会添加错误的最近变化点。与之相比，使用β=（p+1/2）log T产生的Gi（η）函数的平均值会随着η趋于常数。

2.3 最近变化点集

令G是条件损失矩阵，Gi，η=Gi（η），i=1，2，…，N，η=0，1，…，T-1，Gi（η）是指第i个序列最近变化点在η的最优成本。目的是寻找G的K列，使得这些列的每一行成本达到最小，且汇总这些列N行总和达到最小。即将受同一个变化点影响的多个序列划分为一组，N个序列划分为K个不相交的组。优化问题如下：

minS∑Ni=1min η∈SGi，η          （15）

其中S0，1，…，T-1，S=K。这个优化问题在数学上等价于K中值问题，可以用具有二值变量xi，η和zη的整数规划来解决，其中：

xi，η=1，如果序列i最近变化点在时间η0，其他

zη=1，如果存在序列在时间η有最近变化点0，其他

原始问题就被转化为求如下的优化问题：

min ∑Ni=1∑T-1η=0Gi，ηxi，η

s.t. ∑T-1η=0xi，η=1，i

xi，η≤zη，i，η

∑T-1η=0zη=K        （16）

其中，第一个约束保证每个序列只有一个最近变化点，第二和第三个约束确保总共有K个不同的最近变化点。Reese[17]讨论了解决K中值问题的相关方法，使用R包tbart中的方法解决该问题。

3 模拟研究

模拟数据集一维数据长度为500，多维数据维数为100，长度为500。对于给定K值，首先从集合300，320，…，480中模拟产生K个不同的最近变化点，以确保每个最近变化点和其他最近变化点至少间隔20个点的位置。将100个序列按照最近变化点分为K组，针对每条序列，在最近变化点之前的时间点以0.02的概率独立模拟产生更早时间的变化点。然后，从均匀分布中模拟产生一个概率，模拟一个变化点以该概率独立出现在每个时间序列中。每个分段中的观测值都是独立同高斯分布，方差固定为σ2=1，变化服从N（0，22）分布，每个变点跃度为1。结果如图2和图3所示，其中实线表示真实的最近变化点，虚线表示检测的最近变化点。

3.1 一维数据（见图2）

3.2 多维数据（见图3）

由图2、图3可知，对于一维数据和最近变化点较少的多维数据，改进方法对最近变化点的检测较准确。在多维数据最近变化点较多时，改进方法会检测出少量多余的变化点，但是可以准确检测出真实变化点的位置。

4 实例研究

使用美国标准普尔500指数2007年1月1日至2011年12月31日期间股票每日收盘价，选择371家具有完整数据的机构，因此时间序列的维数是371，长度是1260。其累积对数收益如图4（a）所示，由图4（b）累积对数平方收益图可知，总体上的最近变化点在2011年下半年。分别使用原始方法和改进方法检测最近变化点，小波周期图和交叉周期图由尺度j=-1的Haar小波计算，如式（17）和式（18）所示，将其作为2.2部分的输入。

使用改进方法检测到的最近变化点位置及对应公司数如表1所示，使用原始方法检测的结果如表2所示。对比表1、表2可知，原始方法检测出的最近变化点均靠近数据末端且密集，检测结果受噪声影响大，不符合实际情况。而小波具有去噪、自适应等优良性质，特别适合非平稳、非线性的时间序列。小波变换的变化点检测结果表明，改进方法是有效和准确的。D20F8BCA-7CEA-40EA-A250-51F7BF3BCEA4

由表1，將所有数据按照最近变化点分为四组，分别作出分段线性回归图，如图5所示，其中竖直线表示检测到的最近变化点，可看出改进方法检测最近变化点效果显著。将具有相同最近变化点的序列归为同一组，有助于寻找变化产生的原因且有利于后期的短期预测。

5 结 论

结合小波变换与惩罚成本的变化点检测方法，将多维时间序列中的变化点检测问题转化为相对较容易的多维面板数据中的变化点检测问题。实验结果表明，改进方法对最近变化点的检测准确有效，检测的特定变化点是序列中不同且不相交子集的最新变化，且能识别出受不同变化影响的序列。

对于多维时间序列数据中的变化，若单独分析每个序列，效率低且信息冗余;若汇总分析，则个别序列的变化点信号易被其余序列中的噪声淹没。而本方法可根据不同的最近变化点，对多维序列进行分组分析，有助于更好地了解变化产生的原因。此外，根据最近变化点将数据分段做拟合和预测，有利于降低模型拟合误差，提高预测精度。

但是改进后的方法也存在一些问题，如增加了计算量，精确度还有待提高。此外，还有许多问题需要更深入的研究，如使用不同的小波基函数会有什么效果等。

参考文献

[1] BAI P，SAFIKHANI A，MICHAILIDIS G.Multiple change point detection in low rank and sparse high dimensional vector autoregressive models[J].IEEE Transactions on Signal Processing，2020，68：3074-3089.

[2] ZAMENI M，SADRI A，GHAFOORI Z，et al.Unsupervised online change point detection in highdimensional time series[J].Knowledge and Information Systems，2019，62（2）： 719-750.

[3] LAVIELLE M，TEYSSIRE G.Detection of multiple change points in multivariate time series[J].Lithuanian Mathematical Journal，2006，46：287-306.

[4] MATTESON D S，JAMES N A.A nonparametric approach for multiple change point analysis of multivariate data [J].Journal of the American Statistical Association，2014， 109：334-345.

[5] CHO H.Changepoint detection in panel data via double  CUSUM statistic[J].Electronic Journal of Statistics，2016， 10：2000-2038.

[6] LEE S，KIM C K，LEE S.Hybrid CUSUM change point test for time series with timevarying volatilities based on suport vector regression[J].Entropy，2020，22（5）：578.

[7] XIE Y，SIEGMUND D.Sequential multisensor changepoint detection[J].Annals of Statistics，2013，41：670-692.

[8] MILLER D J，GHALYAN N F，MONDAL S，et al.HMM conditionallikelihood based change detection with strict delay tolerance[J].Mechanical Systems and Signal Processing，2020，147（15）.D20F8BCA-7CEA-40EA-A250-51F7BF3BCEA4

[9] CHEN Z，XU Q，LI H.Inference for multiple change points in heavytailed time series via rank likelihood ratio scan statistics[J].Economics Letters，2019，179（1）：53-56.

[10]MA L，GRANT A J，SOFRONOV G.Multiple change point detection and validation in autoregressive time series data[J]. Statal Papers，2020，61（4）：1507-1528.

[11]LAWRENCE B，PAUL F，IDRIS A，et al.Most recent change point detection in panel data[J].Technometrics，2019，61：88-98.

[12]FEARNHEAD P，RIGAILL G.Changepoint detection in the presence of outliers[J].Journal of the American Statistical Association，2019，114（525）：169-183.

[13]FRYZLEWICZ P. Wild binary segmentation for multiple change point detection[J].The Annals of Statistics，2014， 42：2243-2281.

[14]KILLICK R，FEARNHEAD P，ECKLEY I A.Optimal detection of changepoints with a linear computational cost[J].Journal of the American Statistical Association，2012，107：1590-1598.

[15]NASON G P，VON SACHS R，KROISANDT G.Wavelet processes and adaptive estimation of the evolutionary wavelet spectrum[J].Journal of the Royal Statistical Society， Series B.2000，62：271-292.

[16]MATTEO B， HAERAN C， PIOTR F.Simultaneous multiple changepoint and factor analysis for highdimensional time series[J].Journal of Econometrics，2018，206：187-225.

[17]REESE J.Solution methods for the pmedian problem：an  annotated bibliography[J].Networks，2006，48：125-142.D20F8BCA-7CEA-40EA-A250-51F7BF3BCEA4