基于區域結構探測與邊緣辨識發現時間序列因果關系轉換

中圖分類號：TP391 文獻標志碼：A

Discovering Causal Transitions in Time Series Based on Region Structure Detection and Edge Identification

XIE Jiea， WANG Kaijuna，b， FANG Yinga，b， LUO Tianjiana， b （a.Collegeof Computer and Cyber Security，b.Digit Fujian Internet-Things Laboratory ofEnvironmental Monitoring， FujianNormal University，Fuzhou 35O117，Fujian，China）

Abstract：Toimprove theaccuracyof existing methods for mining causal relationships over time，forabinarytime series containing one causal region，amethod was proposed tominethecausal transition pointsof timeseries byidentifying differentregional structuresandtheedgesof structure.Themethod designedthe posiblepositions ofthe causal region in timeseries as left，right，and center structure，used the existing causalitydiscovery methods to detect therough causal region，anddistinguisheditasacertainregional structureaccording tothediferentcharacteristicsof theregionstructures.According to the characteristics of diffrentregional structures，coresponding edge identificationmeasures were designed，and gradualy increasing detection windows andcausality intensity indexes were setto identify regional structural edges as causal transition points，and improvethe identification accuracyofcausal transition points.Experimentson two simulated datasetsand tworeal datasetsverified theaccuracyof the proposed method inrecognizing causal transition points.The resultsshow thatthe average acuracyofcausal transition pointsobtained bythe proposed method using Grangercausality scores onseparable simulated datasets is higher than those of the comparison methods，the average auracyof causal transition points obtained by convergent cross mapping causality scores on weakly coupled simulated datasets is higher than those of thecomparison methodatcoupling degresof O.O1and O.5O，and the accuracyof causal transitionpoints obtained byusing Grangercausalityscoresontworeal datasetsis higher than thatof thecomparison method.

Keywords：timeseries；causalityelation；causalrelation transition；convergentcross maping；Grangercausalitydetection

挖掘時間序列的因果關系，理解變量之間的因果關系對于預測、制定決策和解決問題至關重要[①]分析時間序列因果關系的常用方法包括Granger因果檢驗[2-3]和收斂交叉映射因果檢測（CCM）[4]Granger因果檢驗和CCM方法在探索時間序列的因果關系方面的側重點不同：對于可分離系統和緊密度很高的耦合系統，Granger因果檢驗方法有效；對于緊密度較低的弱耦合系統，CCM方法有效。

在許多現實情況下，因果關系可能會隨著時間的推移而發生改變[5]。目前，發現時間序列之間動態變化的因果關系的主流方法有滑動時間窗口搜索法[5-7]、 F 界檢測法[8]、差異區域平衡法（DrB）[9]和基于CCM 的核變化點檢測法（KCP-CCM）[10]。 F 界檢測法計算檢測窗口的擴展區域的Granger檢驗值上、下界，利用上、下界判別擴展區域內的因果關系，可以減少計算量。差異區域平衡法綜合考慮前向擴展窗口、當前窗口和后向擴展窗口的Granger檢驗結果，以多數票取勝的方式輸出最終結果，降低了檢驗結果對窗口寬度的敏感性，確保最終結果的準確性和穩定性。KCP-CCM方法采用CCM映射相關性作為因果分數，在序列上用滑動窗口從左到右得到因果分數序列，發現因果分數序列的轉換點作為因果關系的轉換點。這些方法挖掘因果區域的性能仍有待提高。

本文中研究單個因果關系區域問題，即只考慮某時間段 T 內有且僅有一個因果關系區域的情況，目標是發現時間段 T 內的因果關系轉換點；提出左、右區結構探測與邊緣辨識（SDEI）方法識別時間序列中具有因果關系的區域。

1相關方法

Granger 因果檢驗的基本原理[8-9]是：如果時間序列 X 引起了時間序列 Y 的變化，與僅使用 Y 的歷史數據預測 Y 的未來值（常用自回歸模型描述，稱約減模型RM）、使用 X 和 Y 的歷史數據預測 Y 的未來值（常用回歸模型描述，稱完整模型FM）相比，將顯著提高預測的準確性。

給定時間點數量 T 的時間序列 ， x₂，… 和 ，擬合時間序列 X 和 Y 的2種回歸模型如下：

式中： L 為設定的最大時滯； a_l，b_l 為指定正整數 l 時的回歸參數； 為 χ_t 時刻約簡模型的預測值； 為 χ_t 時刻完整模型的預測值。

采用 F 檢驗法判別完整模型是否比約減模型的預測更準確（更優），為此構造 F 統計量

其中，預測誤差平方和為

式中 E_r，E_f 分別為約減模型和完整模型的預測誤差平方和。

在零假設（完整模型的預測不顯著優于約減模型）下，若由時間序列數據計算得出的 F 值大于顯著性水平 α 下的閾值 c_{L，T-2L-1}^* ，則拒絕上述零假設，判定 X 是引起 Y 的原因 XY ；否則，判定無因果關系。類似地，可以判別 Y 是否是引起 X 的原因 YX 。

CCM方法[4]以流形理論為基礎，其理論依據是：若時間序列 X 和時間序列 Y 共享一個流形，則認為 X 和 Y 存在因果關系。如果 X 是 Y 的原因，則X 的歷史信息能夠很好地估計 Y

將時間序列 X 和 Y 表示為 和 ，以 X（t） 為例， Φ_t 時刻的時滯向量為 式中： E 為嵌入維度； τ 為時間步長。把所有時滯向量 的集合看作一個重構流形，記作 M_x 。時間序列 Y 的重構流形 M_y 也用同樣的方法構造。

首先定位 M_x 上的時滯坐標向量 x（t） ，然后找到 x（t） 的 E+1 個最近鄰，并按從最近到最遠將時間點 t₁?t₂?…，t_E+1 排序。這些時間點被用于識別 M_y 中假定的鄰點 Y（t_i） ，通過式（3）加權平均估計 Y（t） 。

式中： 是 Y（t） 的估計值； i 為嵌入維度索引；w_i 為第 i 個嵌入維度下基于 與其在 M_x 中最近鄰的距離權重，由式（4）計算得到。

式中 d[?] 為2個向量的歐氏距離。從 Y 到 的交叉映射定義同理。

最后，計算預測值 與觀測值 Y（t） 的Pearson相關系數，將其定義為CCM 映射相關性，記為 ρ_XY 號

轉換點檢測算法[]的基本原理是將相似的觀測序列數據分割在相同區域，轉換點檢測問題可以看作一個優化問題，優化目標如下：

式中： k 為轉換點的索引， k=1，2，…，K;t_k 為第 k 次轉換的時刻，定義 t₀：=1 ， t_k+1：=T ， c（s_tk：tk+1） 是衡量子序列 不相似程度的代價函數。

交叉相關函數[12]是用來衡量2個時間序列之間關聯程度的統計量，可以用于估計2個時間序列之間的相對時滯。交叉相關函數 C_XY（p） 的取值表示序列 X 與序列 Y 相對移動 p 個時間點后序列 X 和序列 Y 的相關性，計算公式為

式中： 、 { Y } 分別為序列 X，Y 的均值； X_t 為時間序列 X 在 χ_t 時刻的值； Y_t+p 為時間序列 Y 在 t+p 時刻的值。該公式表示在相對時滯為 p 時序列 X 和 Y 之間的相關性。

2左、右區結構探測與邊緣辨識法

現有的時間序列因果關系轉換點識別方法通常不考慮因果關系的區域結構特點，但是因果關系的區域結構特點提供的信息能夠提升轉換點識別的準確率。以下利用因果關系的區域結構特點，針對不同區域結構特點設計不同的因果關系轉換點識別方法。

2.1 問題分析

給定時間序列 和時間序列 ，若在 t₁ 時刻到 t₂ 時刻范圍內序列 X 的變化受序列 Y 的變化影響或序列 Y 的變化受序列 X 的變化影響，稱在 t₁ 時刻到 t₂ 時刻范圍內變量 X 和 Y 之間存在因果關系， 稱為因果關系區域。設時間序列 X，Y 之間的關系狀態隨時間發展而發生改變，一種關系狀態將持續一段時間，關系狀態在有因果關系和無因果關系之間變化，時間序列中僅包含一個因果關系區域，目標是發現給定時間序列 X，Y 之間因果關系狀態發生變化的時刻（簡稱因果關系狀態轉換點）。

對于僅包含一個因果關系區域的時間序列，按因果關系區域位置不同設計將其分類為3種結構序列：左區結構序列、中區結構序列、右區結構序列。圖1為3種結構序列示意圖，其中陰影區域表示變量間存在因果關系。左區結構序列是因果關系區域起始于序列左端點，終止于時刻 χ_t ；右區結構序列是因果關系區域起始于時刻 χ_t ，終止于序列右端點；中區結構序列是因果關系區域起始于時刻 t₁ ，終止與時刻 t₂，t₁ 和 t₂ 不在序列的左、右端點。

2.2 求解方法

針對上述3種結構序列的特點，分別設計識別3種結構序列因果關系狀態轉換點的方法以及甄別策略，選擇合適的識別方法。

2.2.1 左、右區結構識別法

左區結構序列的特點是因果關系區域起始于序列左端點，因果關系狀態轉換點左側有因果關系，轉換點右側無因果關系。轉換點出現的時間點 χ_t 和序

列起始時間點0的間隔等于因果關系持續時長 V，t- 0=V ，即 t=V ，于是可以通過估計 U 來識別轉換點 χ_t 。

定義窗口為給定時間范圍內的一個子區域，記作[f₁，f₂] ，其中 f₁，f₂ 分別為窗口左、右端點。記窗口

[f₁，f₂] 內數據因果分數為 S（f₁，f₂） ，因果分數反映序列在窗口內的因果關系強度，設計因果分數滿足如下2個性質：

性質1若窗口內全部時間序列之間存在因果關系，則窗口寬度越大，因果分數越大。

性質2若窗口內因果關系持續時間不變，則窗口寬度越大，因果分數越小。

性質1、2分別由式（8）、（9）描述如下：

S（1，W₁）2），W₁2

S（1，W₁）gt;S（1，W₂），U12

式中 W₁ 1 W₂ 為假設的2個窗口的寬度。

當左區結構序列左端窗口寬度等于因果關系持續時間時，窗口得到的因果分數最大，表達式如下：

依據上述思路和原理可以設計左區結構識別法。將寬度為 W 的窗口的左端點固定于時間序列的左端，通過改變窗口寬度 W 獲得對應不同 W 的窗口內數據的因果分數，估計左區結構序列的因果關系持續時長 V ，計算公式如下：

其中 為因果分數最小有效長度， 太小可能出現誤判， W₀ 太大可能忽略小因果關系區域。得到估計值 V 后可得 t=V ，算法如下。

算法1 左區結構識別法（leftRegionDetect）

輸入：屬于左區結構的時間序列 X 和 Y ，因果分數函數 S（f₁，f₂） ，因果分數最小有效長度 ：

輸出：因果關系轉換點 Ψ_t ：

2： for W in range W₀ to T-W₀+1 doscores[W]=S（1，W）

3：end for

4：

5：return t （204號

右區結構序列與左區結構序列是對稱關系，右區結構識別法（rightRegionDetect）是類似的，不同之處是轉換點 t=T-V 0

2.2.2 中區結構識別法

中區結構序列的特點是在序列的左、右端點無因果關系轉換點。從右區結構和左區結構的視角看，中區結構的左轉換點在右區結構子序列中，右轉換點在左區結構子序列中，于是須要從中區結構序列中提取左區結構子序列和右區結構子序列

提取子序列的思路是，若能找到給定序列因果關系區域的一個子區域 ，則有包含左區結構的子序列 和包含右區結構的子序列 [1，t_b] 。設計利用核變化點檢測方法求給定序列的粗略因果關系區域，得到因果關系區域的左轉換點的近似值o₁ 和右轉換點的近似值 o₂ ，則 為近似的因果關系區域中點， σ_m 的鄰域 W₀/2] 為近似因果關系子區域，于是可得 [o_m-W₀/2 ，T] 為近似包含左區結構的子序列， 為近似包含右區結構的子序列，算法如下。

算法2 中區結構識別法（midRegionDetect）

輸入：屬于中區結構的時間序列 X 和 Y ，因果分數函數 S（f₁，f₂） ，因果分數最小有效長度 ：

輸出：左轉換點 t₁ 和右轉換點 t₂

1： ， O₂=KCP（X，Y）

2：

3： t₁= rightRegionDetect（1， O_m⁺+W₀/2）

4： t₂= leftRegionDetect（ o_m^-W₀/2 ，T）

5 ： return t₁ ， t₂

2.2.3 因果分數的設計

在合適的數據下滿足因果分數性質的統計量有Granger因果檢驗的 F 值，即式（2），以及CCM方法的CCM映射相關性，即式（5）。

在使用Granger因果分數前須要檢驗序列間是否存在相對時滯，具體做法是利用交叉相關函數找到序列間相關性最高的相對時滯 p ：若 p=0 ，則認為序列間不存在相對時滯；若 p≠0 ，則須要調整數據消除相對時滯；若 pgt;0 ，則將序列 Y 左移 |p| 個時間步;若plt;0 ，則將序列 X 右移 |p| 個時間步

2.2.4區域結構探測方法的設計

根據3種結構序列的不同特點，甄別策略的思路是利用選擇的因果分數檢驗給定序列左端是否有因果關系：若左端有因果關系則為左區結構序列，使用左區結構識別法識別因果關系轉換點；若左端無因果關系，則檢驗序列右端是否有因果關系，若右端有因果關系則為右區結構序列，使用右區結構識別法識別因果關系轉換點；若不存在上述2種情況，則使用中區結構識別法識別因果關系轉換點。

根據選擇的因果分數對窗口內時間序列進行因果關系檢驗。當選擇Granger因果檢驗 F 值作為因果分數時，檢驗方法為Granger因果檢驗；當選擇CCM映射相關性作為因果分數時，則采用CCM法檢驗因果關系。為了減小因果關系檢驗誤判的概率，使用3個不同尺度的窗口，分別計算窗口內的因果分數并獲得3個窗口的識別結果，以多數票勝出的方式判斷窗口內數據是否存在因果關系區域[9]

3 實驗結果

3.1 實驗數據

第1個模擬數據是可分離系統。參考 F 界檢測法[和差異區域平衡法產生時間長度 T 的平穩時間序列 X，Y ，記因果關系區域為 ，序列 X 的產生公式為 ，其中 ε_Xt 為均值 μ_x 與方差 σ_X² 的高斯噪聲。序列 Y 的產生公式如下：

式中： y_t 為時間序列 Y 在 Ψ_t 時刻的值； l 為 χ_t 時刻第 l 個滯后的時間點； a_l 控制自相關； b_l 控制Granger 因果關系； ε_Yt 為均值 μ_Y 與方差 σ_Y² 的高斯噪聲。

第2個模擬數據是弱耦合系統。參考 CCM^[4] 和

KCP-CCM[1o]產生時間長度 T 的耦合非線性差分系統。時間序列 X，Y 的因果關系區域為 ，序列X，Y 的產生公式如下：

式中 β_YX，β_XY 分別表示序列 Y 對序列 X 的影響、序列 X 對序列 Y 的影響。

真實數據集下客的出租車數量-推文量和推文量-上客的出租車數量8]對應的實際事件是關于2012年10月某幾日在美國紐約會議中心舉行喜劇演出（出現具有因果影響的事件），而其他時間沒有演出事件發生。該數據集是在該會議中心上、下客的出租車數量與觀眾推文數量的時間序列，樣本數量為745，即 T=745 ，因果關系區域為［247，337]，數據的分布特征見表1。

3.2 性能評價指標

實驗使用調整蘭德指數（ARI）作為檢測準確率評價指標。設 n 個真實的因果關系轉換點將區間[1， T] 為 n+1 個子區間 ，測試方法判斷得到 ?_m 個因果關系轉換點將區間 [1，T] 分為m+1 個子區間 。從時間范圍[1，T] 中任取2個時間點構成時間對 ，則有 C_T² 個不同的時間對，調整蘭德指數 I_AR 計算公式如下：

式中： C_T² 為從 T 個時間點中選擇2個時間點的組合數， C_T²=T（T-1）/2 n_a 為時間對 Ξ（t_p1，t_p2） 在子區間 中的相同子區間且在子區間 P 中的相同子區間的數量; n_b 為時間對 （t_p1，t_p2） 在子區間 中的相同子區間且在子區間 P 中的不同子區間的數量； n_c 為時間對 Ξ（t_p1，t_p2） 在子區間 中的不同子區間且在子區間 P 中的相同子區間的數量； n_d 為時間對 （2號在子區間 中的不同子區間且在子區間 P 中的不同子區間的數量。

3.3 實驗方法

為了驗證本文中提出的SDEI方法的性能，以下采用常規滑動窗口、DrB法、KCP-CCM、基于Granger因果檢驗（GC）的核變化點檢測法（KCP-GC）[10]、SDEIC-CM、SDEI-GC、自舉滑動窗口法（B-rolling）[11]、自舉前向擴展窗口法（B-forward）[]、自舉遞歸窗口法（B-recursive）[11]等方法進行對比。

KCP-CCM和KCP-GC分別用CCM映射相關性作為因果分數和用Granger因果檢驗的 F 值作為因果分數的轉換點識別方法。SDEI-CCM和SDEI-GC分別用CCM映射相關性作為因果分數和用Granger因果檢驗的 F 值作為因果分數的區域邊緣探測方法。對于常規滑動窗口和DrB法，設窗口寬度為30，滑動步長為 10^[9] ;對于KCP-CCM 和 KCP-GC，設窗口寬度為 50^[10] ;對于基于自舉的方法，設最小窗口為90[11]

3.4 實驗結果

對于第1個可分離模擬數據，令 μ_?X=μ_?Y=0，T= 1000，生成3組不同噪聲模擬時間序列 X 和 Y ，噪聲方差分別取 0.01，0.20，0.50， 。每一組實驗分別生成50個左區結構序列、50個右區結構序列、50個中區結構序列，計算每個方法在左區、右區、中區結構序列的平均ARI，然后計算每種方法在3種結構序列的平均ARI。對于左區結構序列，令 t₁= 1， 為 -100～100 的隨機數， v～V（-100 100）];對于右區結構序列，令 t₁=500+v ， t₂=T ；對于中區結構序列，令 t₁=300+v ， t₂=700+v 。由于Granger因果檢驗法適用于識別可分離系統中的因果關系，因此選擇Granger因果檢驗的 F 值作為因果分數，實驗結果見表2。由表中數據可以看出：當噪聲方差為0.01、0.20時，SDEI-GC、KCP-GC、差異區域平衡法均分別對于左區、右區、中區結構序列的檢測準確率最高，SDEI-GC的檢測準確率平均值最大，分別為0.907、0.904；當噪聲方差為0.50時，SDEI-GC、KCP-GC仍分別對于左區、右區結構序列的檢測準確率最高，SDEI-GC對于中區結構序列的檢測準確率最高，SDEI-GC的檢測準確率平均值最大，為0.899。綜上，SDEI-GC法的綜合性能優于對其他對比方法的。

對于第2個弱耦合模擬數據，令=0.2，=0.4，耦合程度 β=β_XX=β_XY ， T=10 00 ，生成3組不同耦合程度的模擬時間序列 X，Y ， β 分別取0.01、0.20， 0. 50 。由于CCM方法適用于識別弱耦合系統中的因果關系，因此選擇CCM映射相關性作為因果分數，實驗方法、因果關系區域類似第1個模擬數據，實驗結果見表3。結果顯示：當耦合程度為0.01時，SDEI-CCM對于左區、右區、中區結構序列的檢測準確率均最高，檢測準確率平均值為0.739；當耦合程度為0.20時，SDEI-CCM對于左區、右區結構序列的檢測準確率均最高，KCP-CCM對于中區結構序列的檢測準確率最高，且KCP-CCM的檢測準確率平均值最大，為0.741；當耦合程度為0.50時，SDEI-CCM對于左區、右區結構序列的檢測準確率均最高，KCP-CCM對于中區結構序列的檢測準確率最高，SDEI-CCM的檢測準確率平均值最大，為0.828。綜上，SDEI-CCM的綜合性能優于其他對比方法的。

對真實數據采用對數預處理，令新時間序列 ， ，并使用Granger因果檢驗的 F 值為因果分數。利用交叉相關函數[12]計算序列之間相對時滯，得到下客的出租車數量-推文量序列之間相對時滯為2，推文量-上客的出租車數量序列之間相對時滯為1。不同方法在真實數據集的檢測準確率如表4所示。從表中數據可以看出，SDEI-GC方法在2個真實數據的ARI均最大，分別為0.9758、0.9179，具有相對最優的轉換點檢測準確率。

注： ①KCP-GC 為基于Granger因果檢驗的核變化點檢測法。 ② SDEI-GC為基于Granger因果檢驗的左、右區結構探測與邊緣辨識法。 ③ B-rolling為自舉滑動窗口法。 ④ B-forward為自舉前向擴展窗口法。 ⑤B. -recursive為自舉遞歸窗口法。

從上述實驗結果可以看出，SDEI法能夠在多數情況下較準確地識別因果關系的轉換點，原因是其通過辨識因果關系區域結構得到了因果關系區域的結構信息，這些結構信息有利于提高轉換點的檢測準確性。

4結語

本文中將包含一個因果關系區域的二元時間序列根據因果關系區域位置設計為3種不同的區域結構序列，針對每種區域結構序列的特點設計對應的辨識方法，設計漸增的探測窗口及其因果關系強度指標。對于包含無因果關系區域和一個因果關系區域的二元可分離時間序列，SDEI-CCM能準確識別出因果關系區域及其轉換點。對于包含無因果關系區域和一個因果關系區域的二元可分離時間序列，SDEI-GC能準確識別出因果關系區域及其轉換點。如何將邊緣辨識方法應用于包含多個因果關系區域的二元時間序列值得今后進一步研究。

參考文獻：

[1]ASSAAD C K，EDVIJVER E，GAUSSIER E. Survey and evaluation of causal discovery methods for time series[J]． Journal of Artificial Intelligence Research，2022，73：767.

[2]DAS P，BABADI B.Non-asymptotic guarantees for reliable identification of Granger causality via the LASSO[J]. IEEE Transactions onInformation Theory，2023，69（11）：7439.

[3]GAO W， YANG H Z. Time-varying group LASSO Granger causality graph for high dimensional dynamic system[J]. Pattern Recognition，2022，130：108789.

[4]SUGIHARA G，MAY R，YE H，et al. Detecting causality in complex ecosystems[J]. Science，2012，338：496.

[5]FINKLEJD，WU JJ，BAGHERI N.Windowed Granger causal inference strategy improves discovery of gene regulatory networks [J].Proceedings of the National Academy of Sciences of the United States of America，2018，115（9）：2252.

[6]CHANG T， TSAI S L，HAGA K A. Uncovering the interrelationship between the U. S. stock and housing markets：a bootstrap rolling window Granger causality approach[J]. Applied Economics， 2017，49： 5841.

[7]MASNADI-SHIRAZI M， MAURYA M R， PAO G， et al. Time varying causal network reconstruction of a mouse cell cycle[J]. BMC Bioinformatics，2019，20：294.

[8] LI Z H，ZHENG G J，AGARWAL A，et al. Discovery of causal time intervals[C]//The 17th SIAM International Conference on Data Mining，April 27-29，2017，Houston，USA.Philadelphia： SIAM，2017：804.

［9］王開軍，曾元鵬，繆忠劍．差異區域平衡法探索時間序列變 化的因果關系[J]．電子與信息學報，2021，43（8）：2414.

[10] GEXL，LIN A J. Kernel change point detection based on convergent cross mapping[J]. Communications in Nonlinear Science and Numerical Simulation，2022，109：106318.

[11] SHI SP，HURN S，PHILLIPSPCB.Causal change detection inpossibly integrated systems：revisiting the money-income relationship[J].Journal of Financial Econometrics，2O2O，18（1）：158.

[12] TRUONG C，OUDRE L，VAYATIS N. Selective review of offline change point detection methods[J]. Signal Processing， 2020，167：107299.

（責任編輯：劉飚）