基于區(qū)域結(jié)構(gòu)探測與邊緣辨識發(fā)現(xiàn)時(shí)間序列因果關(guān)系轉(zhuǎn)換

中圖分類號：TP391 文獻(xiàn)標(biāo)志碼：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，F(xiàn)uzhou 35O117，F(xiàn)ujian，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

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

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

本文中研究單個(gè)因果關(guān)系區(qū)域問題，即只考慮某時(shí)間段 T 內(nèi)有且僅有一個(gè)因果關(guān)系區(qū)域的情況，目標(biāo)是發(fā)現(xiàn)時(shí)間段 T 內(nèi)的因果關(guān)系轉(zhuǎn)換點(diǎn)；提出左、右區(qū)結(jié)構(gòu)探測與邊緣辨識（SDEI）方法識別時(shí)間序列中具有因果關(guān)系的區(qū)域。

1相關(guān)方法

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

給定時(shí)間點(diǎn)數(shù)量 T 的時(shí)間序列 ， x₂，… 和 ，擬合時(shí)間序列 X 和 Y 的2種回歸模型如下：

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

采用 F 檢驗(yàn)法判別完整模型是否比約減模型的預(yù)測更準(zhǔn)確（更優(yōu)），為此構(gòu)造 F 統(tǒng)計(jì)量

其中，預(yù)測誤差平方和為

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

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

CCM方法[4]以流形理論為基礎(chǔ)，其理論依據(jù)是：若時(shí)間序列 X 和時(shí)間序列 Y 共享一個(gè)流形，則認(rèn)為 X 和 Y 存在因果關(guān)系。如果 X 是 Y 的原因，則X 的歷史信息能夠很好地估計(jì) Y

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

首先定位 M_x 上的時(shí)滯坐標(biāo)向量 x（t） ，然后找到 x（t） 的 E+1 個(gè)最近鄰，并按從最近到最遠(yuǎn)將時(shí)間點(diǎn) t₁?t₂?…，t_E+1 排序。這些時(shí)間點(diǎn)被用于識別 M_y 中假定的鄰點(diǎn) Y（t_i） ，通過式（3）加權(quán)平均估計(jì) Y（t） 。

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

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

最后，計(jì)算預(yù)測值 與觀測值 Y（t） 的Pearson相關(guān)系數(shù)，將其定義為CCM 映射相關(guān)性，記為 ρ_XY 號

轉(zhuǎn)換點(diǎn)檢測算法[]的基本原理是將相似的觀測序列數(shù)據(jù)分割在相同區(qū)域，轉(zhuǎn)換點(diǎn)檢測問題可以看作一個(gè)優(yōu)化問題，優(yōu)化目標(biāo)如下：

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

交叉相關(guān)函數(shù)[12]是用來衡量2個(gè)時(shí)間序列之間關(guān)聯(lián)程度的統(tǒng)計(jì)量，可以用于估計(jì)2個(gè)時(shí)間序列之間的相對時(shí)滯。交叉相關(guān)函數(shù) C_XY（p） 的取值表示序列 X 與序列 Y 相對移動 p 個(gè)時(shí)間點(diǎn)后序列 X 和序列 Y 的相關(guān)性，計(jì)算公式為

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

2左、右區(qū)結(jié)構(gòu)探測與邊緣辨識法

現(xiàn)有的時(shí)間序列因果關(guān)系轉(zhuǎn)換點(diǎn)識別方法通常不考慮因果關(guān)系的區(qū)域結(jié)構(gòu)特點(diǎn)，但是因果關(guān)系的區(qū)域結(jié)構(gòu)特點(diǎn)提供的信息能夠提升轉(zhuǎn)換點(diǎn)識別的準(zhǔn)確率。以下利用因果關(guān)系的區(qū)域結(jié)構(gòu)特點(diǎn)，針對不同區(qū)域結(jié)構(gòu)特點(diǎn)設(shè)計(jì)不同的因果關(guān)系轉(zhuǎn)換點(diǎn)識別方法。

2.1 問題分析

給定時(shí)間序列 和時(shí)間序列 ，若在 t₁ 時(shí)刻到 t₂ 時(shí)刻范圍內(nèi)序列 X 的變化受序列 Y 的變化影響或序列 Y 的變化受序列 X 的變化影響，稱在 t₁ 時(shí)刻到 t₂ 時(shí)刻范圍內(nèi)變量 X 和 Y 之間存在因果關(guān)系， 稱為因果關(guān)系區(qū)域。設(shè)時(shí)間序列 X，Y 之間的關(guān)系狀態(tài)隨時(shí)間發(fā)展而發(fā)生改變，一種關(guān)系狀態(tài)將持續(xù)一段時(shí)間，關(guān)系狀態(tài)在有因果關(guān)系和無因果關(guān)系之間變化，時(shí)間序列中僅包含一個(gè)因果關(guān)系區(qū)域，目標(biāo)是發(fā)現(xiàn)給定時(shí)間序列 X，Y 之間因果關(guān)系狀態(tài)發(fā)生變化的時(shí)刻（簡稱因果關(guān)系狀態(tài)轉(zhuǎn)換點(diǎn)）。

對于僅包含一個(gè)因果關(guān)系區(qū)域的時(shí)間序列，按因果關(guān)系區(qū)域位置不同設(shè)計(jì)將其分類為3種結(jié)構(gòu)序列：左區(qū)結(jié)構(gòu)序列、中區(qū)結(jié)構(gòu)序列、右區(qū)結(jié)構(gòu)序列。圖1為3種結(jié)構(gòu)序列示意圖，其中陰影區(qū)域表示變量間存在因果關(guān)系。左區(qū)結(jié)構(gòu)序列是因果關(guān)系區(qū)域起始于序列左端點(diǎn)，終止于時(shí)刻 χ_t ；右區(qū)結(jié)構(gòu)序列是因果關(guān)系區(qū)域起始于時(shí)刻 χ_t ，終止于序列右端點(diǎn)；中區(qū)結(jié)構(gòu)序列是因果關(guān)系區(qū)域起始于時(shí)刻 t₁ ，終止與時(shí)刻 t₂，t₁ 和 t₂ 不在序列的左、右端點(diǎn)。

2.2 求解方法

針對上述3種結(jié)構(gòu)序列的特點(diǎn)，分別設(shè)計(jì)識別3種結(jié)構(gòu)序列因果關(guān)系狀態(tài)轉(zhuǎn)換點(diǎn)的方法以及甄別策略，選擇合適的識別方法。

2.2.1 左、右區(qū)結(jié)構(gòu)識別法

左區(qū)結(jié)構(gòu)序列的特點(diǎn)是因果關(guān)系區(qū)域起始于序列左端點(diǎn)，因果關(guān)系狀態(tài)轉(zhuǎn)換點(diǎn)左側(cè)有因果關(guān)系，轉(zhuǎn)換點(diǎn)右側(cè)無因果關(guān)系。轉(zhuǎn)換點(diǎn)出現(xiàn)的時(shí)間點(diǎn) χ_t 和序

列起始時(shí)間點(diǎn)0的間隔等于因果關(guān)系持續(xù)時(shí)長 V，t- 0=V ，即 t=V ，于是可以通過估計(jì) U 來識別轉(zhuǎn)換點(diǎn) χ_t 。

定義窗口為給定時(shí)間范圍內(nèi)的一個(gè)子區(qū)域，記作[f₁，f₂] ，其中 f₁，f₂ 分別為窗口左、右端點(diǎn)。記窗口

[f₁，f₂] 內(nèi)數(shù)據(jù)因果分?jǐn)?shù)為 S（f₁，f₂） ，因果分?jǐn)?shù)反映序列在窗口內(nèi)的因果關(guān)系強(qiáng)度，設(shè)計(jì)因果分?jǐn)?shù)滿足如下2個(gè)性質(zhì)：

性質(zhì)1若窗口內(nèi)全部時(shí)間序列之間存在因果關(guān)系，則窗口寬度越大，因果分?jǐn)?shù)越大。

性質(zhì)2若窗口內(nèi)因果關(guān)系持續(xù)時(shí)間不變，則窗口寬度越大，因果分?jǐn)?shù)越小。

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

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

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

式中 W₁ 1 W₂ 為假設(shè)的2個(gè)窗口的寬度。

當(dāng)左區(qū)結(jié)構(gòu)序列左端窗口寬度等于因果關(guān)系持續(xù)時(shí)間時(shí)，窗口得到的因果分?jǐn)?shù)最大，表達(dá)式如下：

依據(jù)上述思路和原理可以設(shè)計(jì)左區(qū)結(jié)構(gòu)識別法。將寬度為 W 的窗口的左端點(diǎn)固定于時(shí)間序列的左端，通過改變窗口寬度 W 獲得對應(yīng)不同 W 的窗口內(nèi)數(shù)據(jù)的因果分?jǐn)?shù)，估計(jì)左區(qū)結(jié)構(gòu)序列的因果關(guān)系持續(xù)時(shí)長 V ，計(jì)算公式如下：

其中 為因果分?jǐn)?shù)最小有效長度， 太小可能出現(xiàn)誤判， W₀ 太大可能忽略小因果關(guān)系區(qū)域。得到估計(jì)值 V 后可得 t=V ，算法如下。

算法1 左區(qū)結(jié)構(gòu)識別法（leftRegionDetect）

輸入：屬于左區(qū)結(jié)構(gòu)的時(shí)間序列 X 和 Y ，因果分?jǐn)?shù)函數(shù) S（f₁，f₂） ，因果分?jǐn)?shù)最小有效長度 ：

輸出：因果關(guān)系轉(zhuǎn)換點(diǎn) Ψ_t ：

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

3：end for

4：

5：return t （204號

右區(qū)結(jié)構(gòu)序列與左區(qū)結(jié)構(gòu)序列是對稱關(guān)系，右區(qū)結(jié)構(gòu)識別法（rightRegionDetect）是類似的，不同之處是轉(zhuǎn)換點(diǎn) t=T-V 0

2.2.2 中區(qū)結(jié)構(gòu)識別法

中區(qū)結(jié)構(gòu)序列的特點(diǎn)是在序列的左、右端點(diǎn)無因果關(guān)系轉(zhuǎn)換點(diǎn)。從右區(qū)結(jié)構(gòu)和左區(qū)結(jié)構(gòu)的視角看，中區(qū)結(jié)構(gòu)的左轉(zhuǎn)換點(diǎn)在右區(qū)結(jié)構(gòu)子序列中，右轉(zhuǎn)換點(diǎn)在左區(qū)結(jié)構(gòu)子序列中，于是須要從中區(qū)結(jié)構(gòu)序列中提取左區(qū)結(jié)構(gòu)子序列和右區(qū)結(jié)構(gòu)子序列

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

算法2 中區(qū)結(jié)構(gòu)識別法（midRegionDetect）

輸入：屬于中區(qū)結(jié)構(gòu)的時(shí)間序列 X 和 Y ，因果分?jǐn)?shù)函數(shù) S（f₁，f₂） ，因果分?jǐn)?shù)最小有效長度 ：

輸出：左轉(zhuǎn)換點(diǎn) t₁ 和右轉(zhuǎn)換點(diǎn) 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 因果分?jǐn)?shù)的設(shè)計(jì)

在合適的數(shù)據(jù)下滿足因果分?jǐn)?shù)性質(zhì)的統(tǒng)計(jì)量有Granger因果檢驗(yàn)的 F 值，即式（2），以及CCM方法的CCM映射相關(guān)性，即式（5）。

在使用Granger因果分?jǐn)?shù)前須要檢驗(yàn)序列間是否存在相對時(shí)滯，具體做法是利用交叉相關(guān)函數(shù)找到序列間相關(guān)性最高的相對時(shí)滯 p ：若 p=0 ，則認(rèn)為序列間不存在相對時(shí)滯；若 p≠0 ，則須要調(diào)整數(shù)據(jù)消除相對時(shí)滯；若 pgt;0 ，則將序列 Y 左移 |p| 個(gè)時(shí)間步;若plt;0 ，則將序列 X 右移 |p| 個(gè)時(shí)間步

2.2.4區(qū)域結(jié)構(gòu)探測方法的設(shè)計(jì)

根據(jù)3種結(jié)構(gòu)序列的不同特點(diǎn)，甄別策略的思路是利用選擇的因果分?jǐn)?shù)檢驗(yàn)給定序列左端是否有因果關(guān)系：若左端有因果關(guān)系則為左區(qū)結(jié)構(gòu)序列，使用左區(qū)結(jié)構(gòu)識別法識別因果關(guān)系轉(zhuǎn)換點(diǎn)；若左端無因果關(guān)系，則檢驗(yàn)序列右端是否有因果關(guān)系，若右端有因果關(guān)系則為右區(qū)結(jié)構(gòu)序列，使用右區(qū)結(jié)構(gòu)識別法識別因果關(guān)系轉(zhuǎn)換點(diǎn)；若不存在上述2種情況，則使用中區(qū)結(jié)構(gòu)識別法識別因果關(guān)系轉(zhuǎn)換點(diǎn)。

根據(jù)選擇的因果分?jǐn)?shù)對窗口內(nèi)時(shí)間序列進(jìn)行因果關(guān)系檢驗(yàn)。當(dāng)選擇Granger因果檢驗(yàn) F 值作為因果分?jǐn)?shù)時(shí)，檢驗(yàn)方法為Granger因果檢驗(yàn)；當(dāng)選擇CCM映射相關(guān)性作為因果分?jǐn)?shù)時(shí)，則采用CCM法檢驗(yàn)因果關(guān)系。為了減小因果關(guān)系檢驗(yàn)誤判的概率，使用3個(gè)不同尺度的窗口，分別計(jì)算窗口內(nèi)的因果分?jǐn)?shù)并獲得3個(gè)窗口的識別結(jié)果，以多數(shù)票勝出的方式判斷窗口內(nèi)數(shù)據(jù)是否存在因果關(guān)系區(qū)域[9]

3 實(shí)驗(yàn)結(jié)果

3.1 實(shí)驗(yàn)數(shù)據(jù)

第1個(gè)模擬數(shù)據(jù)是可分離系統(tǒng)。參考 F 界檢測法[和差異區(qū)域平衡法產(chǎn)生時(shí)間長度 T 的平穩(wěn)時(shí)間序列 X，Y ，記因果關(guān)系區(qū)域?yàn)? ，序列 X 的產(chǎn)生公式為 ，其中 ε_Xt 為均值 μ_x 與方差 σ_X² 的高斯噪聲。序列 Y 的產(chǎn)生公式如下：

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

第2個(gè)模擬數(shù)據(jù)是弱耦合系統(tǒng)。參考 CCM^[4] 和

KCP-CCM[1o]產(chǎn)生時(shí)間長度 T 的耦合非線性差分系統(tǒng)。時(shí)間序列 X，Y 的因果關(guān)系區(qū)域?yàn)? ，序列X，Y 的產(chǎn)生公式如下：

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

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

3.2 性能評價(jià)指標(biāo)

實(shí)驗(yàn)使用調(diào)整蘭德指數(shù)（ARI）作為檢測準(zhǔn)確率評價(jià)指標(biāo)。設(shè) n 個(gè)真實(shí)的因果關(guān)系轉(zhuǎn)換點(diǎn)將區(qū)間[1， T] 為 n+1 個(gè)子區(qū)間 ，測試方法判斷得到 ?_m 個(gè)因果關(guān)系轉(zhuǎn)換點(diǎn)將區(qū)間 [1，T] 分為m+1 個(gè)子區(qū)間 。從時(shí)間范圍[1，T] 中任取2個(gè)時(shí)間點(diǎn)構(gòu)成時(shí)間對 ，則有 C_T² 個(gè)不同的時(shí)間對，調(diào)整蘭德指數(shù) I_AR 計(jì)算公式如下：

式中： C_T² 為從 T 個(gè)時(shí)間點(diǎn)中選擇2個(gè)時(shí)間點(diǎn)的組合數(shù)， C_T²=T（T-1）/2 n_a 為時(shí)間對 Ξ（t_p1，t_p2） 在子區(qū)間 中的相同子區(qū)間且在子區(qū)間 P 中的相同子區(qū)間的數(shù)量; n_b 為時(shí)間對 （t_p1，t_p2） 在子區(qū)間 中的相同子區(qū)間且在子區(qū)間 P 中的不同子區(qū)間的數(shù)量； n_c 為時(shí)間對 Ξ（t_p1，t_p2） 在子區(qū)間 中的不同子區(qū)間且在子區(qū)間 P 中的相同子區(qū)間的數(shù)量； n_d 為時(shí)間對 （2號在子區(qū)間 中的不同子區(qū)間且在子區(qū)間 P 中的不同子區(qū)間的數(shù)量。

3.3 實(shí)驗(yàn)方法

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

KCP-CCM和KCP-GC分別用CCM映射相關(guān)性作為因果分?jǐn)?shù)和用Granger因果檢驗(yàn)的 F 值作為因果分?jǐn)?shù)的轉(zhuǎn)換點(diǎn)識別方法。SDEI-CCM和SDEI-GC分別用CCM映射相關(guān)性作為因果分?jǐn)?shù)和用Granger因果檢驗(yàn)的 F 值作為因果分?jǐn)?shù)的區(qū)域邊緣探測方法。對于常規(guī)滑動窗口和DrB法，設(shè)窗口寬度為30，滑動步長為 10^[9] ;對于KCP-CCM 和 KCP-GC，設(shè)窗口寬度為 50^[10] ;對于基于自舉的方法，設(shè)最小窗口為90[11]

3.4 實(shí)驗(yàn)結(jié)果

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

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

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

注： ①KCP-GC 為基于Granger因果檢驗(yàn)的核變化點(diǎn)檢測法。 ② SDEI-GC為基于Granger因果檢驗(yàn)的左、右區(qū)結(jié)構(gòu)探測與邊緣辨識法。 ③ B-rolling為自舉滑動窗口法。 ④ B-forward為自舉前向擴(kuò)展窗口法。 ⑤B. -recursive為自舉遞歸窗口法。

從上述實(shí)驗(yàn)結(jié)果可以看出，SDEI法能夠在多數(shù)情況下較準(zhǔn)確地識別因果關(guān)系的轉(zhuǎn)換點(diǎn)，原因是其通過辨識因果關(guān)系區(qū)域結(jié)構(gòu)得到了因果關(guān)系區(qū)域的結(jié)構(gòu)信息，這些結(jié)構(gòu)信息有利于提高轉(zhuǎn)換點(diǎn)的檢測準(zhǔn)確性。

4結(jié)語

本文中將包含一個(gè)因果關(guān)系區(qū)域的二元時(shí)間序列根據(jù)因果關(guān)系區(qū)域位置設(shè)計(jì)為3種不同的區(qū)域結(jié)構(gòu)序列，針對每種區(qū)域結(jié)構(gòu)序列的特點(diǎn)設(shè)計(jì)對應(yīng)的辨識方法，設(shè)計(jì)漸增的探測窗口及其因果關(guān)系強(qiáng)度指標(biāo)。對于包含無因果關(guān)系區(qū)域和一個(gè)因果關(guān)系區(qū)域的二元可分離時(shí)間序列，SDEI-CCM能準(zhǔn)確識別出因果關(guān)系區(qū)域及其轉(zhuǎn)換點(diǎn)。對于包含無因果關(guān)系區(qū)域和一個(gè)因果關(guān)系區(qū)域的二元可分離時(shí)間序列，SDEI-GC能準(zhǔn)確識別出因果關(guān)系區(qū)域及其轉(zhuǎn)換點(diǎn)。如何將邊緣辨識方法應(yīng)用于包含多個(gè)因果關(guān)系區(qū)域的二元時(shí)間序列值得今后進(jìn)一步研究。

參考文獻(xiàn)：

[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］王開軍，曾元鵬，繆忠劍．差異區(qū)域平衡法探索時(shí)間序列變 化的因果關(guān)系[J]．電子與信息學(xué)報(bào)，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.

（責(zé)任編輯：劉飚）