差異區域平衡法探索時間序列變化的因果關系

王開軍 曾元鵬 繆忠劍

(福建師范大學數學與信息學院 福州 350117)

(福建師范大學數字福建環境監測物聯網實驗室 福州 350117)

1 引言

發現和識別事物之間的因果關系是許多研究領域的重要主題[1,2]，而探索時間序列之間的因果關系對于研究動態發展的系統具有重要的理論和實際意義[3]。已有一些挖掘時間序列間因果關系的方法，例如狀態空間方法[4]通過建立狀態空間模型描述具有時間序列的動態系統，然后根據測量方程來判別因果關系，文獻[5]在回歸模型中并入時滯序列數據，然后用條件概率和期望傳播方法求解模型參數來獲得因果關系。針對2元時間序列的因果關系分析問題，經典的Granger因果分析方法[6]的解決思路是：當某個序列變量加入Granger因果模型有助于獲得更好的預測結果時，可認為模型中的變量之間具有統計意義上的因果關系。基于時間序列的Granger因果分析具有好的可解釋性，在許多領域得到了應用，例如Granger因果分析方法用于分析阻塞性睡眠呼吸暫停患者的多頻道睡眠時間序列來檢測腦電波和心率變異性之間的關聯性[7]，用于探究我國經濟發展和社會發展的因果關系[8]。針對多元時間序列，文獻[9]結合蒙特卡羅方法提出了基于多變量的Granger因果分析方法，文獻[10]提出了基于Hilbert-Schmidt獨立準則的Lasso Granger因果分析方法來分析多元時間序列之間的非線性因果關系，更多關于多變量的Granger因果分析方法參見文獻[3]。

對于隨時間而發展變化的系統，其時間序列(序列變量)之間的因果關系可能隨時間而發生變化。探尋時間序列中發展變化的因果關系是尚待解決的困難問題，有關的研究文獻較少。不失一般性，這里討論兩條時間序列之間的線性因果關系隨著時間改變的情況，側重于因果關系隨時間的變化問題，例如它們之間的因果關系在某個時間區間存在，而在下一個時間區間則變化為不存在。常規的滑動窗口方法[11,12]采用滑動時間窗口搜索時間序列，并用Granger因果方法檢測當前時間窗口內時間序列之間的因果關系；但是，窗寬和步長的不同會導致結果和性能的差異?；瑒哟翱诜椒ㄈ舾F舉所有可能的時間區間(遍歷方式)進行因果關系探尋時，對于長度為T、時滯為L的兩條時間序列，其時間復雜度高達O(T3L2)[13]。F界檢測法[13]采用遍歷方式對包含時刻t的所有不同長度的時間區間進行因果關系檢測，其中對已檢測過的區間的相鄰區域使用F檢驗的上下界判別因果關系；對每個t計算包含t的因果區間數量占全部區間數量的比值作為因果分，依據因果分可判別不同時間區間上序列間因果關系是否存在。文獻[14]采用轉折點檢測算法找出時間序列上的轉折點(作為分割點)，然后在分割區域上進行線性Granger因果檢測，獲得隨時間演化的關系網絡(簡稱轉折點檢測法)；該方法認為由轉折點分割出的不同時間區域在產生時間序列的機制上對應著因果關系結構的改變。但是，復雜多變的時間序列會干擾這種轉折點的正確識別，進而導致該方法性能的下降?？梢钥闯?，常規滑動窗方法的性能對不同窗寬和步長比較敏感而帶來性能不穩定；F界檢測法遍歷不同長度的時間區間的方式偏重全局信息，局部信息傾向于被弱化，則在信號弱的局部區域上容易出現誤判導致性能降低；轉折點檢測法的性能很容易被時間序列數據的波動干擾而降低；因而上述方法在探尋時間序列間發展變化的因果關系上的性能有待進一步改進。

針對滑動窗方法的不足，本文提出一種差異區域平衡方法，設計了前向探索策略、不同長度區域的平衡檢測方案和多數票贏的結果綜合等措施，探索時間序列之間的因果關系隨時間發展而發生變化的情況，期待在合理的窗寬步長時新方法具有更好的性能且穩定。

2 時間序列因果發現的相關方法

常規滑動窗方法探索時間序列隨時間變化的關系的基本框架為：寬wid e的滑動窗Wt按步長step沿時間t方向移動，同時在當前窗口Wt內對時間序列間關系進行探查，并采用Granger因果檢測方法探查窗口內時間序列之間是否存在因果關系，記錄探查結果。

Granger因果分析的基本思想是：如果時間序列X引起了時間序列Y 的變化，相較于僅使用Y 的歷史數據預測Y 的未來值(常用自回歸模型描述，稱約減模型Mr)，使用X和Y 的歷史數據對Y 的未來值進行預測(常用回歸模型描述，稱完整模型Mf)，將提高其預測準確性。下文簡介實現上述思想的Granger因果檢測方法。

給定長度為T 的時間序列X={x1,x2,···,x T}和Y={y1,y2,···,y T}，擬合時間序列X和Y 的兩種回歸模型[13]為(時滯參數L、回歸參數a和b)

采用F檢驗法判別Mf是否比Mr的預測更準確(更優)，為此構造F統計量

F界檢測法[13]主要設計有關于時間的外循環迭代每個時刻t(t=1,2,···,T )，內循環迭代時間區間[i,i+j](其中i=t, j=1,2,···,T–t)：首先執行起步檢測，即對起步區間[i,k](k–i+1的初值可設置為由時滯L 限定的最小長度)內的時間序列采用上述Granger因果檢測方法進行因果關系的檢測；然后，計算擴展區Ut=[i,k+j](其中j=1,2,···,T–k)內F統計量的上下界，使用該上下界判別Ut內的因果關系，若滿足起步條件則k=k+j并執行起步檢測；繼續下一次迭代。最后對每個t，計算包含t的因果區間數量占全部區間數量的比值作為因果分，用因果分判別是否存在因果關系。F界檢測法的主要貢獻在于使用了已檢測過區域的擬合模型和誤差項，替代計算擴展區Ut內的擬合模型和誤差項，來構建F統計量的上下界，節省了計算量。

3 時間序列因果發現的差異區域平衡方法

常規滑動窗方法的主要弱點是其性能對窗口寬度敏感，不同寬度的窗口有可能帶來不同的性能，不合適的窗寬很可能導致低性能。性能改進的新思路是：將這種單一滑動窗口(單一區域)的因果檢測方式改進為3個有差異區域(即當前滑動窗口、當前窗口擴展出的前向和后向區域)的因果檢測方式；于是，當一個區域的數據體現的序列間關系有偏差，而其余2個區域的數據體現的序列間關系為正確時，以多數票贏方式綜合3個區域的結果仍為正確。這樣，與單一區域相比，3區域方式可降低性能對窗口寬度的敏感性。次要弱點是性能對滑動窗的移動步長也敏感。考慮到移動步長的主要作用是控制窗口移動的快慢和計算量，改進性能的思路是：將步長限制在較小的范圍內，從而限制步長敏感性在一個較低的范圍。

具體地設計與當前時間窗口Wt緊相鄰的長wide/2的前向區域作為前向相鄰窗口Ut；與窗口Wt緊相鄰的長wide/4的后向區域作為后向相鄰窗口Vt。在探索時間序列間關系的常規滑動窗方法的框架基礎上，差異區域平衡方法的設計包括如下措施：

(1)前向探索策略：計算窗口Wt內序列的波動程度Sw，并作為局部波動界；計算窗口Ut內序列的波動程度Su；若Su≤Sw，實施差異區域的平衡檢測；若Su>Sw，前向區域出現序列間關系轉換的可能性增大，則實施對稱的差異平衡檢測。

(2)差異區域的平衡檢測：對窗口Wt(長wide)內序列間關系檢測一次，對窗口Wt與Vt的合并區域(長1.25 wide)進行一次關系檢測，對窗口Wt與Ut的合并區域(長1.5 wide)內序列間關系檢測一次。將這3次不同長度區域的檢測結果進行綜合作為當前窗口上的序列間關系，有利于降低單一長度區域時檢測結果可能的偏差影響，提高因果檢測的準確性和穩定性。

(3)對稱區域的平衡檢測：采用上述平衡檢測方案，但其中窗口Vt的長度改為wide/2，此時對窗口Wt而言，窗口Ut與Vt是前后對稱的。當前窗口出現序列間關系的轉換點時，這種對稱方式體現了轉換點前后檢測區域的平衡性，有利于降低轉換點區域的檢測偏差。

(4)多次檢測結果的綜合：多次檢測可能出現不一致的結果，則由多數票贏的投票法來決策有無因果關系。

差異區域平衡法需要時間窗口內數據波動程度的衡量方法，為此設計如下方法。

定義1(平均單間隔差分)設時間軸上在時刻t的滑動時間窗口Wt包含時間序列Y 的m+1個元素{y t?m,···,y t?1,y t}，衡量窗口Wt中元素波動幅度的平均單間隔差分定義為

差異區域平衡法的實現過程由3部分組成：主框架部分是滑動窗的迭代移動過程，即寬wide的滑動窗口Wt以步長step沿時間t方向移動；第2部分是在移動到當前窗口Wt時，依據前向探索策略檢測時間序列X和Y 之間的關系；第3部分是綜合多次檢測結果作為輸出，即每個時間點所記錄的多次檢測結果(有或無關系)采用多數票贏的投票法。

設兩個長度為T 的時間序列X和Y，僅評估目標變量Y的波動界，降低波動界為0.95Sw以加強分析。差異區域平衡法的實現步驟設計如下：

步驟1確定滑動時間窗口W=[i,j]的寬度wide、滑動步長step_slide，序列的自回歸時滯參數L；固定參數有最大探索區域長度uwide=0.5wide，區域長度增量du=0.5uwide；初始參數i=1,j=w id e,j_t ail=j，時長增量d t=0，后向增量dback=0，檢測結果Q={q k=0},k=1:T；

步驟2窗口W=[i,j]上的波動界計算：按前文定義計算Y的平均單、雙間隔差分sf[i,j]與df[i,j]作為波動界，dback=du；

步驟3前向探索波動程度：d t=d t+d u,j_tail=j+dt，前向區域U=[j,j_tail]，計算U上Y的sf[j,j_tail]與df[j,j_tail]；若sf[j,j_tail]>0.95×sf[i,j]或df[j,j_tail]>0.95×df[i,j]，則令dback=uwide，轉下一步驟；若dt≤du，重復本步驟；

步驟4平衡檢測之后向檢測：i_head=max(0, i–dback)，對后向合并區域[i_head, j]內時間序列X和Y 采用Granger因果檢測方法進行檢測；若有因果關系X→Y 或X←Y，記錄qk=q k+1,k=i_head:j；若無因果關系，記錄qk=q k?1,k=i_head:j;dt=0；

步驟5平衡檢測之前向檢測(包括窗口W及前向合并區域W+U)：j_tail=j+dt，對區域[i,j_tail]內時間序列X和Y 采用Granger因果檢測方法進行檢測；若有因果關系X→Y或X←Y，記錄q k=q k+1,k=i:j_tail；若無因果關系，q k=q k?1,k=i:j_tail；若dt

步驟6窗口W移動：若i+step_slide+wide+uwide≤T，令i=i+step_slide,j=i+wide,dt=0，轉步驟2；否則，若j

步驟7輸出結果：Q中q k>0的區域為有因果關系，qk≤0的區域為無因果關系；停止。

現在分析本文方法比現有的常規滑動窗方法、F界檢測法和轉折點方法具有更好的性能穩定性。考慮圖1所示的一般情況，設很少數據點的時間區域上的序列間關系是無意義的；設兩條時間序列X和Y 在時間區域[t1,t2]具有近似線性關系，而在區域[t2,t7]沒有線性關系，但在區域[t4,t5]出現干擾性波動(有線性關系的假象)。

對于常規滑動窗方法，由于區域[t4,t5]中體現出近似線性關系的數據占比大約為60%，該方法的窗寬小于該區域寬度時的Granger因果檢測結果是存在因果關系的概率大，即出現誤判的概率大。而本文方法雖然在區域[t4,t5]出現誤判的概率大，但在區域[t4,t6]和[t3,t5]體現出近似線性關系的數據占比明顯小于50%，使得Granger因果檢測結果是存在因果關系的概率小，即出現誤判的概率??；于是，3個差異區域[t4,t5],[t4,t6]和[t3,t5]的多數票贏帶來正確結果的概率大，即本文方法比常規滑動窗方法具有更穩定的性能。

F界檢測法的遍歷方式偏重全局信息，局部信息傾向于被抑制或弱化。該方法出現弱點的情形是，圖1中當區域[t1,t2]相對于[t2,t7]比較短時(情形1)，短區域[t1,t2]在全區域上占比小，則其(以區間數量占比為基礎的)因果分會比較小，因大多數遍歷區域屬于區域[t2,t7]；相反地，當區域[t1,t2]相對于[t2,t7]比較長時(情形2)，區域[t2,t7]上的因果分比較大，因大多數遍歷區域屬于區域[t1,t2]。這兩種情形會帶來因果分判別準則設計的矛盾性，即取大的因果分閾值會使情形1的區域[t1,t2]出現誤判的概率大，取小的因果分閾值會使情形2的區域[t2,t7]出現誤判的概率大；于是，上述復雜情形容易導致F界檢測法性能下降。而本文方法以窗口局部信息為主，不受全局信息干擾，故比F界檢測法具有更穩定的性能。

轉折點方法出現弱點的情形是，當序列數據波動多時找出的轉折點也多，但時間序列波動的轉折點也可能不是序列間關系的轉折點，例如圖1中時刻t5位置就不是序列間關系的轉折點，這種誤識別會使該方法性能下降。而本文方法不涉及序列波動的轉折點，上述情形不影響本文方法的性能，故本文方法比轉折點方法具有更穩定的性能。

圖1 兩條序列的有近似線性關系區域[t1,t2]和無線性關系區域[t2,t7]

4 實驗結果

本節對差異區域平衡法(簡記為平衡)與常規滑動窗方法(簡記為常規)、F界檢測法(簡記為F界)、轉折點檢測法(簡記為轉折)進行實驗對比，各方法的程序均采用如下相同的參數和設置：運行環境為pycharm和python 3.8.3；主框架為滑動窗移動過程，相鄰區域如Ut與Vt的長度不超過窗口Wt長度的1/2；時滯參數L給出約束條件，即滑動窗口或檢測區域的長度均不小于2L+2；F檢驗的顯著性水平為0.05。對于轉折法，奇異譜分析的嵌入維度為6，對前3條奇異譜進行距離計算，距離計算所用的序列長度為窗口Wt長度的1/2，轉折點的判別閾值為2倍的歷史平均距離(任一條奇異譜序列滿足即可)，相鄰的2個轉折點確定一個檢測區域。對于F界檢測法，以步長1移動的同時遍歷包含時刻t的所有不同長度的時間區域；判別閾值取最高因果分的90%，因果分大于該閾值的區域為有因果關系。

性能評價指標采用正確率，即統計各個時間點上算法的輸出結果與已知因果關系相同的數量(m)，計算該數量占總時間點數量(T )的百分比(100%×m/T )。

模擬數據集參照文獻[13]產生長1000的平穩時間序列X和Y，時滯參數L=2，有關系區間[s*,t*]=[450,650],X的產生公式為x t=1+sin(0.08t)+εt，其中εt為零均值與方差σ2的高斯噪聲。每次模擬實驗都產生一個隨機模擬數據集，包括噪聲εt隨機產生，序列Y 的初值y0和y1從0～0.1的均勻分布中隨機采樣，關系區間的邊界s*和t*分別添加區間[–20,20]中的隨機整數，使得s*和t*每次實驗不相同。使用不同的窗口寬度wide和步長step_slide組合(但不適用于F界檢測法)，對各種方法重復進行50次性能實驗，獲得的平均結果列于表1中。

表1 不同方法在模擬數據集上發掘因果關系的正確率(%)

真實數據集Dropoff-tweet來自文獻[13]的圖4，是2012年10月在紐約會議中心下客的出租車數量與觀眾評論帖數量隨時間變化的時間序列，樣本數量為745(t=1～745)。在紐約會議中心有喜劇演出的11～14日，下客的出租車數量(Dropoff序列X)是引起該演出的評論帖數量(Tweet序列Y )變化的原因。對序列X和Y 的樣本計算平方根作為數據預處理；時滯參數L=5。這樣，時滯對齊后的序列X和Y 在演出區間[248,334]具有近似線性關系(因果關系)，其中每天的夜眠時間(凌晨1～7點)不存在因果關系。實驗任務是檢測31天中在何時段時間序列之間存在因果關系?？紤]一天24小時內序列X和Y之間的影響關系，采用不同的窗口寬度wide和步長step_slide的組合(但不適用于F界檢測法)，其中F界檢測法遍歷的最大區間長度為72(與文獻[13]一致)，實驗結果列于表2中。

表2 在數據集Dropoff-tweet上發掘因果關系的正確率(%)

真實數據集Tweet-pickup來自文獻[13]的圖4，是2012年10月在紐約會議中心上客的出租車數量與觀眾評論帖數量隨時間變化的時間序列。演出結束前后時間，關于演出的評論帖數量(Tweet序列X)是引起上客的出租車數量(Pickup序列Y )變化的原因。其余信息與數據集Dropoff-tweet相同，實驗結果列于表3中。

表3 在數據集Tweet-pickup上發掘因果關系的正確率(%)

真實數據集Fish-school[15]是未訓練的魚(游動方向序列Y)跟隨被訓練的魚(游動方向序列X)游向喂食點的時間序列集中第2個子集，樣本數量為596(t=1～596)。實驗任務是檢測該時間序列之間在何時段存在因果關系。數據做如下預處理：對序列X計算每20個(類似地對Y是每10個)連續的序列樣本的均值(滾動均值)，獲得在t=20～596時的平滑樣本，而t=1～20時的樣本取t=20時的值；然后對平滑樣本計算平方根。時滯參數L=65。這樣，在魚游向喂食點期間t=300～479，時滯對齊后的序列X和Y之間具有近似線性關系(因果關系)。考慮最小窗寬需大于132，實驗采用不同的窗口寬度wide和步長step_slide組合(但不適用于F界檢測法)，實驗結果列于表4中。

表4 在數據集Fish-school上發掘因果關系的正確率(%)

真實數據集Baboon-troop[15]是一隊狒狒(移動方向序列Y)跟隨狒狒ID3(移動方向序列X)移動的時間序列集中第1個子集，樣本數量為599(t=1～599)。實驗任務是檢測該時間序列之間在何時段存在因果關系。對序列X和Y的數據均做如下預處理：計算每25個連續的序列樣本的均值(滾動均值)，獲得在t=2 5 ～5 9 6 時的平滑樣本，而t=1～24時的樣本取t=25時的值。時滯參數L=12。這樣，在一隊狒狒跟隨狒狒ID3移動期間t=92～445，時滯對齊后的序列X和Y之間具有近似線性關系(因果關系)?？紤]是近似線性關系，在較短時間區間上序列間關系不穩定，故取較大窗寬進行實驗；實驗采用不同的窗口寬度wide和步長step_slide組合(但不適用于F界檢測法)，實驗結果列于表5中。

表5 在數據集Baboon-troop上發掘因果關系的正確率(%)

從上面各表中可以看出，差異區域平衡法成功地揭示出時間序列之間隨時間變化的因果關系。在數據集Dropoff-tweet,Tweet-pickup和模擬數據集上，平衡方法的性能與對比方法相當，其原因是：出現因果關系的區域剛好位于全部時間區域的中部，且序列X和Y 出現顯著增大/變化的值(即體現因果關系的信號強)，而在其他時間區域序列X的值為0或接近0，這種區分性極其顯著的數據集使得各方法容易獲得好的性能。在數據集Fish-school和Baboon-troop上，絕大部分情況下平衡方法的正確率都較大幅度高于對比方法；對比方法性能較低的原因是，數據集出現上節提及的因果區間長于非因果區間、序列數據波動較多這樣導致對比方法性能下降的情形。雖然在少量情況下平衡方法的正確率略低于對比方法，但在高正確率且性能穩定的綜合性能上，差異區域平衡法優于對比方法。

實驗中還獲得關于參數的觀察體會：對于轉折點方法，從表1可以看出，其性能隨著噪聲的增大下降很多，即容易受數據波動的影響；該方法更適合時間序列波動的波峰波谷差別顯著且與序列間關系轉換對應的情況；當序列波動較小時應當設置較小的判別閾值，但小閾值很可能找出過多的波峰波谷(轉折點)，而設置大閾值則可能丟失轉折點；因此，閾值參數本身就是個矛盾體。對于F界檢測法，實驗中使用的閾值是依據多次實驗經驗選取的最好閾值，如何更好地設置因果分閾值仍需進一步的研究。對于常規滑動窗方法，文獻[12]建議窗寬不超過20，而上文的實驗結果表1—表3顯示，窗寬超過20反而出現更好的性能。因此，如何調整參數以提升性能仍需要更深入的研究。

上述實驗結果顯現出窗口寬度和滑動步長對正確率有一定的影響。故而，需要根據實際問題相對合理地選取窗口寬度和滑動步長。在相對合理的窗寬和步長下，正確率高和性能穩定都是我們選用方法的考量指標。因此，從正確率和穩定性這兩方面的綜合角度考量，差異區域平衡法勝過對比方法。

5 結束語

在常規滑動窗方法和Granger因果檢測的基礎上，本文方法新設計了3個差異區域替代單一滑動窗口，新設計了3個差異區域的平衡檢測方式，設計了多次檢測結果不一致時按多數票贏的投票法確定最終是否存在因果關系的方法。這種差異區域平衡策略有利于降低單一區域檢測結果可能的偏差影響，有利于降低方法性能對窗寬的敏感性，提高因果關系檢測的準確性和穩定性?；瑒哟胺椒ǖ拇翱趯挾群突瑒硬介L之間是否存在某種影響性能的關聯性值得進一步研究。