陳雪文

摘?要：現(xiàn)代數(shù)據(jù)科學(xué)中存在大量的多維時(shí)間序列數(shù)據(jù)，檢測(cè)多維時(shí)間序列中的最新變化點(diǎn)對(duì)于短期預(yù)測(cè)很重要。一種改進(jìn)的方法被提出，以檢測(cè)此類多維時(shí)間序列數(shù)據(jù)中最新變化點(diǎn)。通過使用小波變換，將多維時(shí)間序列中的變化點(diǎn)檢測(cè)問題轉(zhuǎn)化為相對(duì)較容易的多維面板數(shù)據(jù)中的變化點(diǎn)檢測(cè)問題。該方法旨在跨時(shí)間序列合并信息，以便優(yōu)先推斷多個(gè)序列中同一時(shí)間點(diǎn)的最新變化。通過對(duì)每個(gè)時(shí)間序列的輸出進(jìn)行后處理，獲取最新變化點(diǎn)的時(shí)間集及該時(shí)間最近變化點(diǎn)的序列集。最后使用R軟件以模擬數(shù)據(jù)和S&P500數(shù)據(jù)實(shí)驗(yàn)證明了該方法的有效性。

關(guān)鍵詞：多維時(shí)間序列;小波變換;懲罰成本;最新變化點(diǎn)

中圖分類號(hào)：O213???????文獻(xiàn)標(biāo)識(shí)碼：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

在許多現(xiàn)代應(yīng)用中，會(huì)產(chǎn)生大量多維時(shí)間序列數(shù)據(jù)。多維時(shí)間序列數(shù)據(jù)中變化點(diǎn)的檢測(cè)至關(guān)重要，現(xiàn)有的變化點(diǎn)檢測(cè)方法大致可以分為兩類：回顧性（后驗(yàn)）方法和在線檢測(cè)方法。后驗(yàn)方法使用所有的觀測(cè)值，檢測(cè)已經(jīng)發(fā)生的變化點(diǎn)[1];在線檢測(cè)方法對(duì)新到達(dá)的數(shù)據(jù)進(jìn)行自適應(yīng)分析[2]。

在多個(gè)時(shí)間序列中檢測(cè)變化點(diǎn)會(huì)出現(xiàn)單個(gè)時(shí)間序列不存在變化點(diǎn)的問題，因此可以采用單變量變化點(diǎn)方法：一種是分別分析每個(gè)時(shí)間序列，另一種是匯總時(shí)間序列并分析所得的單變量序列。每種方法都有其缺點(diǎn)：前者在檢測(cè)變化點(diǎn)時(shí)忽略了不同時(shí)間序列可能在相似時(shí)間具有變化點(diǎn)的信息;如果來自影響少量序列的變化點(diǎn)的信號(hào)在匯總時(shí)被其余序列中的噪聲淹沒，則后者性能可能較差。

分析多維時(shí)間序列數(shù)據(jù)的另一種方法是將數(shù)據(jù)視為具有多變量觀測(cè)值的單個(gè)時(shí)間序列。對(duì)分段中的多元數(shù)據(jù)建模，并允許該模型以適當(dāng)?shù)姆绞皆诓煌糠种g有變化。如Lavielle[3]將數(shù)據(jù)建模為均值可能因分段而異的多元高斯模型，Matteson[4]提出一種非參數(shù)方法以檢測(cè)多元數(shù)據(jù)中的多個(gè)變化。但是，就像匯總數(shù)據(jù)一樣，如果變化僅影響時(shí)間序列的一小部分，則這些方法可能表現(xiàn)較差。

針對(duì)僅影響序列子集的變化點(diǎn)，Cho[5]基于CUSUM統(tǒng)計(jì)量[6]，提出一種檢測(cè)多個(gè)變化點(diǎn)存在和位置信息的方法。保留所有序列的CUSUM值能夠證明在給定時(shí)間點(diǎn)是否發(fā)生變化，但降低了其他時(shí)間序列值的權(quán)重。同樣，Xie[7]引入廣義似然比檢測(cè)僅影響序列子集的單個(gè)公共變化點(diǎn)。該檢驗(yàn)需要估計(jì)受變化點(diǎn)影響的序列比例，此估算會(huì)影響每個(gè)序列發(fā)生變化可能性權(quán)重，如給予變化點(diǎn)可能性較大的序列較大的權(quán)重，而給予其他序列較低的權(quán)重。

為估計(jì)每個(gè)時(shí)間序列的最新變化點(diǎn)，采用了不同的方法。該方法旨在將多維時(shí)間序列通過小波變換，轉(zhuǎn)化為多維面板數(shù)據(jù)，對(duì)面板數(shù)據(jù)序列進(jìn)行懲罰成本分析，最后輸出每個(gè)序列的最新變化點(diǎn)。

某些因素會(huì)影響一個(gè)時(shí)間序列的結(jié)構(gòu)變化，也可能會(huì)影響其他時(shí)間序列的結(jié)構(gòu)變化。但是，并非所有時(shí)間序列都可能在同一時(shí)間存在一個(gè)變化點(diǎn)。因此，對(duì)這些分析的輸出進(jìn)行后處理，將時(shí)間序列分組，同一組共享公共變化點(diǎn)。該方法能夠準(zhǔn)確估計(jì)每個(gè)時(shí)間序列的最新變化點(diǎn)的位置，這有助于提高序列短期預(yù)測(cè)準(zhǔn)確性。

1?單變量時(shí)間序列

假設(shè)時(shí)間序列為y1：T，變化點(diǎn)的數(shù)量為m，位置為τ=（τ1，…，τm），且τ0=0，τm+1=T。將時(shí)間序列在變化點(diǎn)處分段，對(duì)段內(nèi)數(shù)據(jù)建模，每一段的成本要與極大似然值成比例[8-9]，以將該模型擬合到數(shù)據(jù)段[10]，從而得出成本。假設(shè)段中數(shù)據(jù)的模型具有一定密度f（yθ）且獨(dú)立同分布，其中θ是段特定參數(shù)，那么段ys：t的成本[11]定義為：

為檢測(cè)均值的變化，一個(gè)簡(jiǎn)單的模型是分段內(nèi)數(shù)據(jù)獨(dú)立同高斯分布，方差均為σ2，分段特定均值為θ，那么成本函數(shù)如下：

若分段內(nèi)的均值是關(guān)于時(shí)間的線性函數(shù)，且這個(gè)線性模型在不同的段中不同。表示為θ=θ1，θ2，其中θ1是截距，θ2是斜率。假設(shè)此模型的噪聲是獨(dú)立同高斯分布，那么段成本函數(shù)為：

某些時(shí)間序列存在明顯的異常值，為對(duì)異常值具有魯棒性，采用Fearnhead（2017）[12]中方法，使用如下分段成本函數(shù)：

若異常值的殘差與分段均值的距離大于2倍標(biāo)準(zhǔn)偏差，則此成本會(huì)限制異常值的影響。對(duì)于剩余方差σ2，基于差分時(shí)間序列的中值絕對(duì)偏差，使用一個(gè)簡(jiǎn)單而魯棒的σ2估計(jì)量[13]。

為分段數(shù)據(jù)并找到變化點(diǎn)，希望將所有分段的成本總和降至最低。為避免過度擬合，為每個(gè)分段添加一個(gè)懲罰項(xiàng)β（β>0）。因此，對(duì)數(shù)據(jù)進(jìn)行分段，即解決以下優(yōu)化問題：

β的值越高，意味著檢測(cè)到的變化點(diǎn)越少。根據(jù)BIC準(zhǔn)則，若段特定參數(shù)的維數(shù)為p，則β=（p+1）log?T。定義F（t），t=1，2，…，T如下：

其中0=τ0<τ1<…<τm<τm+1=t，因此，F(xiàn)（t）是分段y1：t的最小成本。函數(shù)F（t）能夠使用PELT[14]中的算法有效計(jì)算：

假設(shè)最新的變化點(diǎn)在時(shí)間r時(shí)，最低成本為G（r）。G（r）與F（r）有關(guān)，因?yàn)镕（r）是分段y1：r的最小成本。

可知G（0）=C（y1：T），因?yàn)橐呀?jīng)對(duì)最近變化點(diǎn)之前的變化點(diǎn)的數(shù)量和位置進(jìn)行了優(yōu)化，因此可認(rèn)為G（r）的大小與邊側(cè)似然函數(shù)有關(guān)。最近變化點(diǎn)的估計(jì)由arg?min?r∈0，…，T-1G（r）得出。若最近變化點(diǎn)在r=0，則該時(shí)間序列沒有變化點(diǎn)。

2?多維時(shí)間序列

2.1?小波變換

Nason[15]最初提出使用小波作為非平穩(wěn)時(shí)間序列的基礎(chǔ)，類似于平穩(wěn)過程經(jīng)典Cramer表示中的傅立葉指數(shù)。使用最簡(jiǎn)單的小波系統(tǒng)，Haar小波：

其中j∈{-1，-2，…}是小波尺度因子，l∈Z是位置參數(shù)。

假設(shè)時(shí)間序列的維度為N，長(zhǎng)度為T，yi，t表示第i條時(shí)間序列，截至?xí)r間t的最近變化點(diǎn)為η（t），那么yηti，t小波系數(shù)為：

其中Lj是尺度為j的離散小波向量ψj=（ψj，0，…，ψj，Lj-1）T的支持長(zhǎng)度，Lj=M2-j（M是定值，取決于小波族）。序列yηti，t的小波周期圖序列和交叉周期圖序列分別為：D20F8BCA-7CEA-40EA-A250-51F7BF3BCEA4

2.2?多維面板數(shù)據(jù)檢測(cè)最近變化點(diǎn)

為檢測(cè)面板數(shù)據(jù)中的最近變化點(diǎn)集，令Gi（η）是第i個(gè)序列最近變化點(diǎn)在η的最小成本函數(shù)。假設(shè)有K個(gè)公共最近變化點(diǎn)，η1：K=（η1，…，ηK）。對(duì)于第k個(gè)最近變化點(diǎn)ηk，存在一個(gè)集合Ιk1，2，…，N，使得對(duì)于所有i∈Ιk的序列，最近變化點(diǎn)都是ηk。由此，集合I1：K可將多維序列分組。

為得到η1：K的估值，需最小化每個(gè)序列的成本總和：

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

對(duì)于K值的選取，使用式（13）的懲罰形式，針對(duì)一個(gè)范圍中不同的K值，解決式（13）的優(yōu)化問題，選擇使得下式達(dá)到最小的K值：

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

使用BIC選擇方法的懲罰項(xiàng)不利于添加最新的變化點(diǎn)，因此，在計(jì)算Gi（η）時(shí)，添加一個(gè)比BIC選擇略低的懲罰項(xiàng)β=（p+1/2）log?T。在沒有變化點(diǎn)的模擬數(shù)據(jù)中，隨著η的增大，較小的β值會(huì)產(chǎn)生更小的G（η），表明選擇較小的β會(huì)添加錯(cuò)誤的最近變化點(diǎn)。與之相比，使用β=（p+1/2）log?T產(chǎn)生的Gi（η）函數(shù)的平均值會(huì)隨著η趨于常數(shù)。

2.3?最近變化點(diǎn)集

令G是條件損失矩陣，Gi，η=Gi（η），i=1，2，…，N，η=0，1，…，T-1，Gi（η）是指第i個(gè)序列最近變化點(diǎn)在η的最優(yōu)成本。目的是尋找G的K列，使得這些列的每一行成本達(dá)到最小，且匯總這些列N行總和達(dá)到最小。即將受同一個(gè)變化點(diǎn)影響的多個(gè)序列劃分為一組，N個(gè)序列劃分為K個(gè)不相交的組。優(yōu)化問題如下：

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

其中S0，1，…，T-1，S=K。這個(gè)優(yōu)化問題在數(shù)學(xué)上等價(jià)于K中值問題，可以用具有二值變量xi，η和zη的整數(shù)規(guī)劃來解決，其中：

xi，η=1，如果序列i最近變化點(diǎn)在時(shí)間η0，其他

zη=1，如果存在序列在時(shí)間η有最近變化點(diǎn)0，其他

原始問題就被轉(zhuǎn)化為求如下的優(yōu)化問題：

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

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

xi，η≤zη，i，η

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

其中，第一個(gè)約束保證每個(gè)序列只有一個(gè)最近變化點(diǎn)，第二和第三個(gè)約束確?？偣灿蠯個(gè)不同的最近變化點(diǎn)。Reese[17]討論了解決K中值問題的相關(guān)方法，使用R包tbart中的方法解決該問題。

3?模擬研究

模擬數(shù)據(jù)集一維數(shù)據(jù)長(zhǎng)度為500，多維數(shù)據(jù)維數(shù)為100，長(zhǎng)度為500。對(duì)于給定K值，首先從集合300，320，…，480中模擬產(chǎn)生K個(gè)不同的最近變化點(diǎn)，以確保每個(gè)最近變化點(diǎn)和其他最近變化點(diǎn)至少間隔20個(gè)點(diǎn)的位置。將100個(gè)序列按照最近變化點(diǎn)分為K組，針對(duì)每條序列，在最近變化點(diǎn)之前的時(shí)間點(diǎn)以0.02的概率獨(dú)立模擬產(chǎn)生更早時(shí)間的變化點(diǎn)。然后，從均勻分布中模擬產(chǎn)生一個(gè)概率，模擬一個(gè)變化點(diǎn)以該概率獨(dú)立出現(xiàn)在每個(gè)時(shí)間序列中。每個(gè)分段中的觀測(cè)值都是獨(dú)立同高斯分布，方差固定為σ2=1，變化服從N（0，22）分布，每個(gè)變點(diǎn)躍度為1。結(jié)果如圖2和圖3所示，其中實(shí)線表示真實(shí)的最近變化點(diǎn)，虛線表示檢測(cè)的最近變化點(diǎn)。

3.1?一維數(shù)據(jù)（見圖2）

3.2?多維數(shù)據(jù)（見圖3）

由圖2、圖3可知，對(duì)于一維數(shù)據(jù)和最近變化點(diǎn)較少的多維數(shù)據(jù)，改進(jìn)方法對(duì)最近變化點(diǎn)的檢測(cè)較準(zhǔn)確。在多維數(shù)據(jù)最近變化點(diǎn)較多時(shí)，改進(jìn)方法會(huì)檢測(cè)出少量多余的變化點(diǎn)，但是可以準(zhǔn)確檢測(cè)出真實(shí)變化點(diǎn)的位置。

4?實(shí)例研究

使用美國(guó)標(biāo)準(zhǔn)普爾500指數(shù)2007年1月1日至2011年12月31日期間股票每日收盤價(jià)，選擇371家具有完整數(shù)據(jù)的機(jī)構(gòu)，因此時(shí)間序列的維數(shù)是371，長(zhǎng)度是1260。其累積對(duì)數(shù)收益如圖4（a）所示，由圖4（b）累積對(duì)數(shù)平方收益圖可知，總體上的最近變化點(diǎn)在2011年下半年。分別使用原始方法和改進(jìn)方法檢測(cè)最近變化點(diǎn)，小波周期圖和交叉周期圖由尺度j=-1的Haar小波計(jì)算，如式（17）和式（18）所示，將其作為2.2部分的輸入。

使用改進(jìn)方法檢測(cè)到的最近變化點(diǎn)位置及對(duì)應(yīng)公司數(shù)如表1所示，使用原始方法檢測(cè)的結(jié)果如表2所示。對(duì)比表1、表2可知，原始方法檢測(cè)出的最近變化點(diǎn)均靠近數(shù)據(jù)末端且密集，檢測(cè)結(jié)果受噪聲影響大，不符合實(shí)際情況。而小波具有去噪、自適應(yīng)等優(yōu)良性質(zhì)，特別適合非平穩(wěn)、非線性的時(shí)間序列。小波變換的變化點(diǎn)檢測(cè)結(jié)果表明，改進(jìn)方法是有效和準(zhǔn)確的。D20F8BCA-7CEA-40EA-A250-51F7BF3BCEA4

由表1，將所有數(shù)據(jù)按照最近變化點(diǎn)分為四組，分別作出分段線性回歸圖，如圖5所示，其中豎直線表示檢測(cè)到的最近變化點(diǎn)，可看出改進(jìn)方法檢測(cè)最近變化點(diǎn)效果顯著。將具有相同最近變化點(diǎn)的序列歸為同一組，有助于尋找變化產(chǎn)生的原因且有利于后期的短期預(yù)測(cè)。

5?結(jié)?論

結(jié)合小波變換與懲罰成本的變化點(diǎn)檢測(cè)方法，將多維時(shí)間序列中的變化點(diǎn)檢測(cè)問題轉(zhuǎn)化為相對(duì)較容易的多維面板數(shù)據(jù)中的變化點(diǎn)檢測(cè)問題。實(shí)驗(yàn)結(jié)果表明，改進(jìn)方法對(duì)最近變化點(diǎn)的檢測(cè)準(zhǔn)確有效，檢測(cè)的特定變化點(diǎn)是序列中不同且不相交子集的最新變化，且能識(shí)別出受不同變化影響的序列。

對(duì)于多維時(shí)間序列數(shù)據(jù)中的變化，若單獨(dú)分析每個(gè)序列，效率低且信息冗余;若匯總分析，則個(gè)別序列的變化點(diǎn)信號(hào)易被其余序列中的噪聲淹沒。而本方法可根據(jù)不同的最近變化點(diǎn)，對(duì)多維序列進(jìn)行分組分析，有助于更好地了解變化產(chǎn)生的原因。此外，根據(jù)最近變化點(diǎn)將數(shù)據(jù)分段做擬合和預(yù)測(cè)，有利于降低模型擬合誤差，提高預(yù)測(cè)精度。

但是改進(jìn)后的方法也存在一些問題，如增加了計(jì)算量，精確度還有待提高。此外，還有許多問題需要更深入的研究，如使用不同的小波基函數(shù)會(huì)有什么效果等。

參考文獻(xiàn)

[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，F(xiàn)EARNHEAD?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