基于小波變換多維時間序列最新變化點檢測

陳雪文

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

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

中圖分類號：O213???????文獻(xiàn)標(biāo)識碼：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)用中，會產(chǎn)生大量多維時間序列數(shù)據(jù)。多維時間序列數(shù)據(jù)中變化點的檢測至關(guān)重要，現(xiàn)有的變化點檢測方法大致可以分為兩類：回顧性（后驗）方法和在線檢測方法。后驗方法使用所有的觀測值，檢測已經(jīng)發(fā)生的變化點[1];在線檢測方法對新到達(dá)的數(shù)據(jù)進(jìn)行自適應(yīng)分析[2]。

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

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

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

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

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

1?單變量時間序列

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

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

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

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

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

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

β的值越高，意味著檢測到的變化點越少。根據(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]中的算法有效計算：

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

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

2?多維時間序列

2.1?小波變換

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

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

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

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

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

為檢測面板數(shù)據(jù)中的最近變化點集，令Gi（η）是第i個序列最近變化點在η的最小成本函數(shù)。假設(shè)有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）的優(yōu)化問題，選擇使得下式達(dá)到最小的K值：

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

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

2.3?最近變化點集

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

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

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

xi，η=1，如果序列i最近變化點在時間η0，其他

zη=1，如果存在序列在時間η有最近變化點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）

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

3?模擬研究

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

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

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

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

4?實例研究

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

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

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

5?結(jié)?論

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

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

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

參考文獻(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