面向非平穩時間序列的因果關系發現算法

周嘉穎,周躍進

(安徽理工大學數學與大數據學院,安徽 淮南 232001)

0 引言

時間序列數據是指在相同間隔的時間段內,觀察某個研究對象的數據變化過程及將這些數據按照時間先后順序排列所形成的序列數據。在日常生活中,時間序列數據廣泛存在于農業、醫學、工業等多個領域[1-3]。隨著大數據時代的不斷發展,產生了海量、非平穩、非線性的時間序列數據,這為挖掘其有效信息增加了難度。因此,挖掘復雜時間序列數據的潛在信息,揭示未來發展規律,成為當前一個重要研究方向。

分析事物或現象之間的因果關系是現實中的常見問題,例如分析大腦不同區域間因果關系以構建大腦網絡[4]和分析商品價格與房產價格之間的雙向因果關系[5]等。由于時間序列的時序性,通過了解時序變量之間的因果關系,預測事物或現象的發展情況,因此因果關系分析方法得到廣大學者關注。GRANGER[6]首次提出了Granger因果關系分析法,是一種判別二元時間序列之間是否存在因果關系的方法,其關鍵假設是:(1)因果關系產生的過程可以用一組結構方程來表示;(2)任何時間點產生的因果效應均受到過去時間點的影響。由于傳統的Granger因果關系分析法只能用于判別二元線性時間序列,因此出現了大量改進模型。GEWEKE[7]提出了條件Granger因果模型,用于判別多元時間序列之間的因果關系。Granger因果關系通常是在線性系統的背景下研究的,隨著研究的深入,學者們開始探索在非線性系統中的因果關系。ANCONA等[8]在Granger因果關系和徑向基函數的基礎上提出了可用于判別非線性時間序列的RBF-Granger因果模型。

傳統的Granger因果關系分析法及其推廣分析法只能給出定性分析結果,分析高維時間序列時容易產生虛假的因果關系?；诖?SCHREIBER[9]首次提出了轉移熵的概念,一種基于信息理論的因果關系分析法,可以捕獲時間序列的非線性特征和定量分析因果關系的強弱。在Granger因果關系分析法和轉移熵因果關系分析法中,均假設時間序列的因果關系都受到原因的影響且時滯是固定的,受固定滯后時間序列影響仍然存在于Granger因果關系和轉移熵結果中。因果關系之間的時滯是固定的假設對于現實生活中自然和社會現象來說是過于絕對的,時間序列的時滯并不是固定的,不同情況下時滯會隨時間變化,一組時間序列可以由多組原因序列影響。因此,AMORNBUNCHORNVEJ[10]提出了一種時滯可變的轉移熵(Variable-Lag Transfer Entropy,VL-TE),可以分析具有任意時滯時間序列之間的因果關系,轉移熵值越大,其因果關系越強。

傳統的Granger因果關系分析法及其推廣分析法、轉移熵法及VL-TE法只能分析平穩時間序列間的因果關系,而不能直接分析非平穩時間序列之間的關系。為了解決此問題,本文提出一種基于分段聚合近似可變時滯轉移熵(Piecewise Aggregate Approximated Variable-Lag Transfer Entropy, PAAVL-TE)因果關系發現算法。該算法利用分段聚合近似法對時間序列進行轉換,提取時間序列的特征信息,使用動態時間彎曲距離尋找相似程度最高的時間序列進行轉移熵的計算,能夠實現在具有可變時滯的非平穩時間序列中的因果關系判定。最后在模擬數據集及真實數據集進行了實驗,并與Granger因果關系分析法、轉移熵法及VL-TE法對比,通過實驗驗證本文方法的有效性和應用性。

1 相關理論

1.1 信息理論基礎

SHANNON[11]提出了香農熵的概念來表示系統的混亂程度和隨機變量所含信息量的多少及其不確定性的程度。系統混亂程度越高,熵值越大,隨機變量的不確定性就越大。香農熵H(X)定義為:

(1)

為了衡量多個變量含有的共同信息量,提出了聯合熵H(X,Y)的概念,其定義為:

(2)

同時,為了衡量在一個條件下隨機變量的復雜程度,提出了條件熵H(Y|X)的概念,其定義為:

(3)

為了衡量兩個隨機變量的相關程度,提出了互信息(Mutual Information,MI)的概念,反映變量之間信息交互情況的度量?；バ畔⒅翟酱?變量之間的相關性也越高?；バ畔I定義為:

(4)

基于信息理論中香農熵、聯合熵、條件熵及互信息的概念,有學者提出了多種基于互信息理論的因果分析法,其中包括轉移熵(Transfer Entropy,TE)。

1.2 轉移熵與VL-TE

由于互信息無法得到信息的傳遞方向,Schreiber基于互信息理論提出了轉移熵的概念。轉移熵是一種基于兩個隨機過程的過去值及當前值來測量兩個過程的信息傳遞方向及信息傳遞量的因果關系分析方法。給定時間序列{Xt}和{Yt},從{Xt}到{Yt}的轉移熵定義為:

(5)

基于香農熵作信息量含量不確定性的度量時,香農轉移熵的定義為:

(6)

然而,這不符合實際中時間序列的時滯階數是可變的現實。AMORNBUNCHORNVEJ[10]提出了一種時滯階數可變的轉移熵(VL-TE),可以判定具有可變時滯時間序列之間的因果關系,VL-TE的定義為:

(7)

為了判定存在可變時滯時X是否是導致Y變化的原因,定義了可變時滯轉移熵比值:

(8)

當T(X,Y)VLr大于1時,表示在可變時滯轉移熵中X是導致Y變化的原因。比值越大,X導致Y變化的原因程度也越大。

由于實際中的時間序列不都是平穩的,而VL-TE只可發現平穩時間序列的因果關系。因此,本文的工作是探究非平穩的時滯階數可變時間序列的因果關系。

2 基于分段聚合近似可變時滯轉移熵的因果關系發現算法

2.1 時間序列預處理

對于缺失的數據,采用多重插補法進行填補,再對處理過后的時間序列進行Z-Score標準化,即將數據轉換成均值為0、標準差為1的數據:

(9)

其中,X′和X分別表示標準化后的數據和原始數據,μ和σ分別表示原始數據的均值和標準差。

對預處理過后的時間序列,采用分段聚合近似表示法(Piecewise Aggregate Approximation,PAA)進行處理。

2.2 分段聚合近似表示法

(10)

分段聚合近似表示法把長度為70的時間序列平均分成10段,如圖1所示。

圖1 基于分段聚合近似的時間序列降維示意圖

(11)

其中,c為最大壓縮比,θ為允許誤差,F為時間序列的變化頻率,k為常數,n為時間序列的數據總數。

2.3 動態時間彎曲距離

動態時間彎曲距離用于度量兩個時間序列之間的距離,是一種度量待分類的兩個時間序列間相似性的方法。通過對時間序列上不同時間點對應的元素進行動態彎曲調整,獲取一條最優路徑,以檢測不同時間序列的相似程度,從而最大限度地減少了時間移位和失真的影響。

給定兩個時間序列{Xt1,t1=1,2,…,m},{Yt2,t2=1,2,…,n}進行距離度量,構造一個矩陣D={d(i,j)}m×n,其中距離度量d(i,j)=(Xi-Yj)2。動態時間彎曲則是需要在矩陣D中所有彎曲路徑P={p1,p2,…,pK}找到一條連續的最優路徑P*,其中第k個元素為pk=(i,j)k且(max{m,n}≤K

(12)

(13)

圖2 動態時間彎曲路徑

動態時間彎曲距離的代碼如下所示:

輸入:時間序列{Xt1,t1=1,2,…,m}, {Yt2,t2=1,2,…,n}

建立元素距離矩陣D={d(i,j)}m×n,其中d(i,j)=(Xi-Yj)2

fort1=2∶m

end for

fort2=2∶n

end for

fort1=2∶m,t2=2∶n

end for

設定pK(i,j)=(m,n)

whilepK(i,j)≠(1,1)

end

記P的長度為K,長度為K的最優動態彎曲路徑P*=P(K-l+1),l∈{1,2,…,K}

輸入:時間序列{Xt}和{Yt}

由動態時間彎曲距離算法得到動態彎曲路徑P={p1,p2,…,pl,…,pK}

設定互相關時間序列{Xt}和{Yt}的結果為P0={p0,…,p0,…}

for alll

ifS(Xt-pl,Yt)

else if

end if

end for

本文提出的基于分段聚合近似可變時滯轉移熵因果關系發現算法(PAAVL-TE)的基本思想是:首先對各個時間序列進行分段聚合近似(PAA)計算,得到特征時間序列。然后對這些特征時間序列進行動態時間彎曲距離計算,找到相似程度最高的新的時間序列,最后計算轉移熵值分析因果關系。

基于分段聚合近似可變時滯轉移熵因果關系發現算法如下:

輸入:時間序列{Xt1,t1=1,2,…,m}, {Yt2,t2=1,2,…,n},可能的最大滯后階數δmax,分段數目nX和nY

輸出:轉移熵因果分析結果,TX→Y,TY→X

數據預處理:得到新時間序列{X′t1,t1=1,2,…,m}, {Y′t2,t2=1,2,…,n}

計算壓縮比cX′=m/nX,cY′=n/nY

fort1=1∶nX

fort2=1∶nY

分別計算轉移熵值TX→Y,TY→X,

ifT(X,Y)r>1

存在X→Y的因果關系

else if

不存在X→Y的因果關系

end if

3 仿真模擬

為了驗證本文提出的PAAVL-TE算法的有效性,使用具有可變時滯的模擬數據集對算法進行實驗評估,并與G、TE、VL-TE算法進行對比。

在模擬數據集中生成了固定時滯無法判別因果關系而可變時滯可以判別的時間序列,每次生成每組長度為500的時間序列,正常模型的時間序列值來自標準正態分布N(0,1),自回歸模型(Autoregressive Model,AR)的時間序列值來自:

Xt=0.5Xt-1+1.2Xt-2+εt,εt～N(0,1).

模擬生成5類時間序列數據集,具體來說:(1)15組存在明顯因果關系的時間序列{Xt}和{Yt}的數據集,其中{Yt}具有時滯階數p=5,Yt=Xt-p+0.1εt,εt～N(0,1)。為確保時滯的可變性,{Yt}在第200至250時間步長直接保持不變,模仿第180時間步長的{Xt},使{Yt}受到X的可變滯后的影響;(2)15組{Xt}和{Yt}獨立生成的不存在因果關系的時間序列,確保PAAVL-TE不會判別出錯誤的因果關系;(3)15組時間序列{Xt}和{Yt}存在因果關系的AR模型數據集;(4)15組時間序列{Xt}和{Yt}不存在因果關系的AR模型;(5)15組時間序列{Xt}來自正常模型數據集,{Yt}來自AR模型。

生成的模擬數據集的真實因果關系如圖3所示。圖3中的箭頭所指方向表示從原因時間序列(如{X1})到結果時間序列(如{Y1})的因果方向,{Yij}表示存在滯后的時間序列。{X1,X2,X3}由正常模型及AR模型生成,每次模擬生成15組時間序列,設置最大時滯階數δmax=12。此外,本文將F檢驗顯著性水平α設為0.05。當T(X,Y)r>1時,才判定存在因果關系。

圖3 模擬數據集的真實因果關系圖 圖4 不同方法的ROC曲線

對于算法評價,本文采用ROC曲線下的面積AUC(Area Under Curve)指標進行評價。ROC曲線是以假陽率(FPR)為橫坐標,真陽率(TPR)為縱坐標,由各閾值下的點坐標構成的曲線,其中假陽率和真陽率的計算公式如下:

其中,將預測因果關系的結果與實際因果關系的結果作對比,真陽性TP是實際有因果關系且預測也有因果關系的數量;真陰性TN是實際無因果關系且預測也無因果關系的數量;假陽性FP是實際無因果關系但預測有因果關系的數量;假陰性FN是實際有因果關系但預測無因果關系的數量,(FPR,TPR)為各閾值下的點坐標。AUC指數越高,算法預測效果越好,與實際因果關系越接近。

為了評估PAAVL-TE從非平穩時間序列中推斷出正確的因果圖邊,同時采用精確率(Precision)、召回率(Recall)和F1得分評價本文方法及其他對比方法的效果。

其中,P表示精確率,R表示召回率,F1表示F1得分,Tp表示事實與預測結果都有Xi→Yj因果圖邊的數量,Fp表示因果圖事實無邊但預測有邊的數量,FN表示事實有Xi→Yj因果圖邊但預測無邊的數量。

圖4給出了G、TE、VL-TE、PAAVL-TE四種方法的ROC曲線圖。從圖4可以看出,本文方法預測效果最好。相比已有方法,本文具有更好的性能。

表1給出了模擬數據集中G、TE、VL-TE、PAAVL-TE四種方法的實驗結果。從表1可以看出,PAAVL-TE在精確率及F1得分上明顯優于其他方法,這表明PAAVL-TE方法可以有效處理復雜時間序列的因果推理任務,分析出正確時滯可變時的因果關系,排除其他方法不能排除的無關因素,提高精確率。

表1 模擬數據集中不同方法的實驗結果

圖5給出了在改變最大時滯階數δmax時四種方法推斷因果方向的平均準確度。從圖5可以看出,無論最大時滯階數δmax為何種情況,在δmax變化的范圍內PAAVL-TE方法的平均準確度均高于其他方法,這表明在非平穩時間序列中提出的PAAVL-TE方法具有更好的表現?？傮w來說,PAAVL-TE算法在具有可變時滯的時間序列的因果關系發現問題上具有較好的實驗效果。

圖5 改變最大時滯階數δmax時不同方法推斷因果方向的平均準確度

4 應用實例

以2013年3月1日至2017年2月28日北京市昌平區的PM2.5濃度、污染物濃度及氣象時間序列為研究對象分析因果關系,數據來源于UCI數據庫的Beijing Multi-Site Air-Quality數據,一共有11維數據(PM2.5濃度、SO2濃度、NO2濃度、CO濃度、O3濃度、溫度、壓強、露點、降雨量、風速、風向),其中包含5維空氣污染物濃度數據和6維氣象數據。利用PAAVL-TE方法找出影響PM2.5濃度變化的主要因素,同時與其他方法進行對比。

首先,對數據進行缺失填補及Z-Score標準化,再進行PAA處理。然后,以PM2.5濃度為分析目標,利用PAAVL-TE尋找影響PM2.5濃度變化的因素,計算各因變量對應的可變時滯轉移熵比值,剔除無關和冗余變量。比值大于1則認為有因果關系,且比值越大,因果關系程度越強。最后,使用保留有因果關系的相關變量,進行預測建模。

四種因果分析方法顯示的影響PM2.5濃度的因素如表2所示。根據本文方法得到,對PM2.5濃度變化有因果關系的變量有SO2濃度、NO2濃度、CO濃度、溫度、降雨量、風速及風向。隨著北京市發展及汽車的普及程度增高,汽車尾氣排放、煤炭燃燒及工業排放等會造成大氣污染,汽車尾氣排放的污染物主要有SO2、NO2、CO等。由于在大氣中SO2、NO2、CO存在二次轉化過程,如NO2經過一系列復雜物理化學反應可轉化為硝酸鹽二次顆粒,使得PM2.5濃度上升。因此,PM2.5濃度變化會受到SO2、NO2、CO污染物濃度的影響。此外,北京市位于華北平原的西北緣,且地形呈簸箕狀。當風向為西北風時,容易生成更大風速的風,使得污染物更易擴散,空氣中PM2.5濃度也隨之下降。反之,東南風則會把污染物吹向簸箕地形,不利于污染物的擴散,從而PM2.5濃度升高。降雨量及溫度也影響污染物濃度,如當近地面大氣溫度較高時,大氣的對流作用加劇,使得PM2.5濃度降低,所以PM2.5濃度也受溫度及降雨量的影響。因此,運用PAAVL-TE法的因果關系分析所得到的結果,與北京市PM2.5濃度變化的影響關系一致,驗證了本文方法的有效性。

表2 PM2.5在不同方法下的預測結果

得出影響PM2.5濃度變化的因素后,本文建立長短時記憶網絡(Long Short Term Memory,LSTM)預測模型,通過這些因素進行建模預測并驗證分析的結果,針對每一種方法的預測結果,取不同參數下10折交叉驗證所得到的最好結果。采用平均絕對誤差(MAE)、對稱平均絕對百分比誤差(SMAPE)和均方根誤差(RMSE)三個指標來衡量預測的精度,定義如下:

表2給出了PM2.5在G、TE、VL-TE、PAAVL-TE四種方法下的預測結果。由表2可見,本文方法在三種誤差的預測評價指標均小于其他方法,表明其預測精度最高。同時用于預測的變量數目為7,小于數據集的原始變量數目,減小了預測過程中的計算難度及成本。相比之下,本文方法能夠有效應用于具有可變時滯的非平穩時間序列因果關系分析,為模型預測選擇合適的輸入變量,提高預測精度。這表明本文方法表現良好,具有現實應用價值。

5 結語

針對具有可變時滯的非平穩時間序列因果關系分析問題,本文提出了PAAVL-TE方法,利用分段近似聚合和動態時間彎曲距離的方法找到與原序列相似程度最高的時間序列,并計算可變時滯轉移熵,克服了傳統的Granger因果關系分析法及其推廣分析法、轉移熵法及VL-TE法難以運用于非平穩時間序列因果關系分析的不足。同時,有效分析出正確的因果關系和因果關系程度的高低,減少了計算時間,說明該算法具有良好的應用性。在未來的工作中將研究基于狀態空間的因果分析方法,在相空間重構時嘗試利用神經網絡以確定Granger因果模型的延遲階數和嵌入維數,降低模型的復雜度。