漫談格蘭杰因果關(guān)系（Granger Causality）——第一章 野火燒不盡，春風(fēng)吹又生
2017年7月9日上午6點(diǎn)10分，先師胡三清同志——新因果關(guān)系的提出者、植入式腦部電極癲癇治療法的提出者、IEEE高級(jí)會(huì)員，因肺癌醫(yī)治無(wú)效于杭州腫瘤醫(yī)院去世，享年50歲。余蒙先師厚恩數(shù)載，一朝忽聞先師駕鶴西歸，悲痛不已。瘁心之余，遂決意傳先師之道，以慰先師在天之靈。如此，先師蓋以瞑目矣！
格蘭杰因果關(guān)系作為一種可以衡量時(shí)間序列之間相互影響關(guān)系的方法，最近十幾年備受青睞。無(wú)論是經(jīng)濟(jì)學(xué)[1]，氣象科學(xué)[2]，神經(jīng)科學(xué)[3]都有廣泛的應(yīng)用，盡管后兩者（氣象和神經(jīng)科學(xué)）連格蘭杰自己都反對(duì)（格蘭杰反對(duì)將格蘭杰因果關(guān)系用在除經(jīng)濟(jì)學(xué)以外的其他領(lǐng)域，這就是本文題目所謂的“野火”）[4]。鑒于筆者從未在氣象學(xué)有過(guò)半分建樹(shù)，所以不敢妄談。不過(guò)慶幸的是，經(jīng)過(guò)神經(jīng)科學(xué)家數(shù)十載的辛苦“洗地”，他們紛紛找到了自己‘合法’使用格蘭杰因果關(guān)系的理由[5]（Anil Seth是英國(guó)皇家科學(xué)院院士，筆者首推的‘地表最強(qiáng)洗地王’，也即本文題目所謂的“春風(fēng)”）。除了由克里夫·格蘭杰本人提出的格蘭杰因果關(guān)系之外，還有數(shù)種圍繞格蘭杰因果關(guān)系方法產(chǎn)生的變體，本文也將對(duì)這些變體分門(mén)別類(lèi)，做出一些簡(jiǎn)介。當(dāng)然，也正是因?yàn)楹壬J(rèn)為應(yīng)當(dāng)遵從格蘭杰爵士的論文，所以才創(chuàng)立了‘新因果關(guān)系’專(zhuān)門(mén)解決神經(jīng)科學(xué)中的因果關(guān)系使用問(wèn)題，不過(guò)這是后話了。本文的講述將圍繞以下幾個(gè)方面進(jìn)行：
本文的框架：
第一章 野火燒不盡，春風(fēng)吹又生
一：格蘭杰因果關(guān)系究竟是什么？
二：格蘭杰因果關(guān)系原理是什么？
三：求解格蘭杰因果關(guān)系的標(biāo)準(zhǔn)化流程是什么？
第二章 星星因果，可以燎原
四：格蘭杰因果關(guān)系變種、原理及其應(yīng)用領(lǐng)域？
五：正宗的格蘭杰因果關(guān)系怎樣用代碼實(shí)現(xiàn)？（MATLAB版）
第三章 胡氏因果關(guān)系（新因果關(guān)系）（Innovation by Sanqing Hu）
一：格蘭杰因果關(guān)系究竟是什么？
一般口頭上所稱(chēng)的格蘭杰因果關(guān)系全稱(chēng)是“格蘭杰因果關(guān)系檢驗(yàn)”（Granger causality test），看到這里，諸位看官大可以感嘆一句“啊哈！原來(lái)格蘭杰因果關(guān)系是統(tǒng)計(jì)假設(shè)檢驗(yàn)的一種！”。那么“因果”二字又要從何談起呢？這就要從時(shí)間序列間的關(guān)系開(kāi)始說(shuō)起。格蘭杰是這樣定義“因果”的。現(xiàn)在我們假設(shè)有一個(gè)時(shí)間序列X，它是由不同時(shí)間的采樣點(diǎn){x1,x2,x3,……,xn}共同組成的。同X一樣，我們假設(shè)有時(shí)間序列Y，它形如X，由{y1,y2,y3,……,yn}共同組成。現(xiàn)在我們利用X的過(guò)去預(yù)測(cè)X的未來(lái)，比如用x1~xn-j（這就是X的過(guò)去）預(yù)測(cè)xn-j+1~xn（這就是X的未來(lái)），預(yù)測(cè)的過(guò)程中產(chǎn)生一個(gè)誤差δ1（請(qǐng)先忽略預(yù)測(cè)的具體方法和誤差產(chǎn)生的具體方法），然后把這個(gè)誤差視為我們得到的第一個(gè)結(jié)果。再然后利用X和Y共同的過(guò)去預(yù)測(cè)X的未來(lái)，比如用{x1~xn-j（這就是X的過(guò)去）|y1~yn-j（這就是Y的過(guò)去）}預(yù)測(cè)xn-j+1~xn（這就是X的未來(lái)），預(yù)測(cè)的過(guò)程中產(chǎn)生一個(gè)誤差δ2（請(qǐng)先忽略“X和Y是怎樣共同的 ?”這個(gè)問(wèn)題） ，然后把這個(gè)誤差視為我們得到的第二個(gè)結(jié)果。如果δ1小于δ2，也就是說(shuō)X和Y的聯(lián)合預(yù)測(cè)誤差小于X自身的預(yù)測(cè)誤差，那么，必然是因?yàn)閅對(duì)X的預(yù)測(cè)起到了幫助，所以才減小了預(yù)測(cè)誤差。在這種情況下，我們稱(chēng)Y對(duì)X有格蘭杰因果關(guān)系。
二：格蘭杰因果關(guān)系原理是什么？
這一節(jié)里，我們要重點(diǎn)解決上一節(jié)中出現(xiàn)的幾個(gè)問(wèn)題。即1）.預(yù)測(cè)究竟是怎樣實(shí)現(xiàn)的？；2）.誤差究竟是怎樣產(chǎn)生的？；3）X和Y是怎樣共同聯(lián)合預(yù)測(cè)的 ? 在講解這些問(wèn)題之前。有一個(gè)非常重要的概念——回歸問(wèn)題，必須預(yù)先鋪墊，否則，這三個(gè)問(wèn)題無(wú)從下手。為了方便理解，我們?cè)诙S空間中表述這個(gè)問(wèn)題，當(dāng)然這個(gè)問(wèn)題也可以被拓展到高維空間。當(dāng)然，如果某位看官確保自己已經(jīng)弄清楚了回歸問(wèn)題，那么可以直接跳到回歸問(wèn)題講完之后的部分。
自回歸問(wèn)題講解開(kāi)始：
現(xiàn)假設(shè)存在一系列二維空間點(diǎn)的集合S={x1,x2,x3,……,xn }。如圖1所示。現(xiàn)在我們想找到一條線，這條線最好穿過(guò)點(diǎn)集S中的每一個(gè)點(diǎn)。這些點(diǎn)包含一個(gè)隱變量T。這樣，當(dāng)找到這條線的時(shí)候，我們實(shí)際上找到了代表這條二維空間中的線的一個(gè)表示函數(shù)x = f(t)。
圖1
那么怎樣找到這條線呢？
STEP1，你需要假定這條線來(lái)源于一個(gè)階數(shù)為m的函數(shù)。如果你的假設(shè)是階數(shù)為3，那么這個(gè)函數(shù)就是一個(gè)形如x = f(t) = a0+a1t1+a2t2+a3t3 的函數(shù)。集合S中的每一個(gè)點(diǎn)都必須要滿足（或盡力滿足）x = f(t)這個(gè)公式。就拿第一個(gè)點(diǎn)[t1,x1]舉個(gè)栗子：x1 = f(t) = a0+a1t11+a2t12+a3t13。當(dāng)然，我們知道這條線y=f(x)完美的經(jīng)過(guò)每一個(gè)點(diǎn)近乎是不可能的，因此我們?cè)试S誤差的存在。只要y = f(x)+ε 即可，其中ε代表誤差項(xiàng)，并且ε服從正態(tài)分布（有時(shí)也稱(chēng)高斯白噪聲，高斯分布）。
STEP2，嘗試求出系數(shù)a0~am滿足STEP1的要求，使得對(duì)于點(diǎn)集S中所有的點(diǎn)滿足x = f(t)+ε。一般情況下，求解系數(shù)a0~am 最常見(jiàn)的方法就是最小二乘法。最小二乘法的具體步驟這里不再描述。一般可以通過(guò)matlab函數(shù)polyfit進(jìn)行求解，代碼如下。求解結(jié)果如圖2。
X=[4,3.5,3.6,2.1,4,6,5.7,5,4]; a = polyfit(1:9,X,3); T_new = 0:0.1:10; X_new = polyval(a,T_new); plot(1:9,X,'r*');hold on; plot(T_new,X_new,'b-');
在上述代碼中，polyfit就是擬合參數(shù)，它負(fù)責(zé)根據(jù)時(shí)間變化1:9以及自變量X將參數(shù)a求解出來(lái)，注意，a在此處是個(gè)數(shù)組，由a0~a3組成。polyval函數(shù)根據(jù)參數(shù)數(shù)組a和點(diǎn)集T_new求解X_new。而繪制[T_new, X_new]就是圖2中的藍(lán)色曲線。
圖2
可能細(xì)心的同學(xué)能發(fā)現(xiàn)，圖1和圖2中，時(shí)間軸T的值域不一樣。在圖1中，T軸的值域是1~9；圖2中，x軸的值域1~10。而這多余出來(lái)的X10其實(shí)可以被視為我們根據(jù)已有數(shù)據(jù)集[X1-9]預(yù)測(cè)出來(lái)的未來(lái)的數(shù)據(jù)[X10]。
總之，自回歸問(wèn)題的一般形式就是給出一系列點(diǎn)集[T,X]，然后想辦法找到一個(gè)函數(shù)X = f(T) 。我們要想法找到一系列參數(shù)a，使得函數(shù)代表的這條線盡可能經(jīng)過(guò)點(diǎn)集。
自回歸問(wèn)題講解結(jié)束。
講完了自回歸問(wèn)題，以下兩個(gè)問(wèn)題就可以得到完美的解答：
同本文的開(kāi)始相同，現(xiàn)在我們假定擁有兩個(gè)時(shí)間序列X：{x1,x2,x3,……,xn }下標(biāo)代表時(shí)間T從1到n，n是實(shí)數(shù)。
1）.預(yù)測(cè)究竟是怎樣實(shí)現(xiàn)的？
讀到這里時(shí)，我假設(shè)各位看官已經(jīng)弄清楚了我在上面所述的回歸問(wèn)題。但是格蘭杰因果關(guān)系中的回歸問(wèn)題同原問(wèn)題的原理相同但過(guò)程有所不同。且聽(tīng)我慢慢道來(lái)他們的不同之處：
不同之處一：在回歸的原問(wèn)題中，我們以時(shí)間軸T為橫坐標(biāo)盡力求解自變量X。試圖找到時(shí)間軸T同自變量X之間的關(guān)系。即X = f(T)。比如我們?cè)谏鲜街信e出的例子里X = f(T) = a0+a1T1+a2T2+a3T3（請(qǐng)注意：這個(gè)問(wèn)題實(shí)際上是三階問(wèn)題，'階' 是指除了a0之外，還有幾個(gè)參數(shù)a。例子中的回歸問(wèn)題其實(shí)是一個(gè)非線性回歸（次數(shù)為一才是線性的））。在回歸問(wèn)題的末尾，我們說(shuō)過(guò)，回歸問(wèn)題的本質(zhì)是要找到一系列的參數(shù)a，因此，我們?cè)诶又信e到的3階形式，存在的具體形式如下：
X = f(T) = a0+a1t11+a2t12+a3t13（式一）······普通回歸問(wèn)題（三階一次問(wèn)題（線性回歸））
在格蘭杰因果關(guān)系中，自回歸模型就是類(lèi)似于這種的一階形式。但是它不再是T映射到X上的函數(shù)，而是過(guò)去的X映射到現(xiàn)在的X上的函數(shù)，即
X_now = f(X_past)，在考慮過(guò)去三個(gè)點(diǎn)的情況下，這個(gè)問(wèn)題就變成了
xt =a1xt-1+a2xt-2 +a3xt-3（式二）······格蘭杰回歸問(wèn)題（注意到?jīng)]有常數(shù)項(xiàng)）
既然在式一中，知曉T和X，進(jìn)而求解一些列參數(shù)a不成問(wèn)題，那么在這里，知曉時(shí)間序列上的X，進(jìn)而求解一系列參數(shù)a也不是問(wèn)題。
不同之處二：
既然是回歸問(wèn)題就要涉及到回歸階數(shù)order的選擇問(wèn)題。通過(guò)不同之處一，我們知曉，在格蘭杰因果關(guān)系檢驗(yàn)中，階數(shù)的作用更像是在指明采用過(guò)去的多少個(gè)點(diǎn)對(duì)當(dāng)前點(diǎn)的點(diǎn)求解回歸問(wèn)題。因此階數(shù)order有了一個(gè)嶄新的名字：lag（遲延，遲滯）。在闡述回歸問(wèn)題的時(shí)候，我們假設(shè)遲滯lag=3，然后使用最小二乘法擬合了點(diǎn)集S。實(shí)際上，為了簡(jiǎn)化最小二乘擬合的問(wèn)題，筆者武斷的設(shè)置了遲滯lag=3而并沒(méi)有采取適當(dāng)?shù)胤椒▽?duì)lag應(yīng)有的值做出選擇。在格蘭杰因果關(guān)系的計(jì)算過(guò)程中，針對(duì)階數(shù)的選擇方法大致可以分為三類(lèi)：I. Akaike information criterion(赤池信息準(zhǔn)則，簡(jiǎn)稱(chēng)AIC)；II. Bayes information criterion(貝葉斯信息準(zhǔn)則，簡(jiǎn)稱(chēng)BIC)；III Rule of thumb(憑經(jīng)驗(yàn)選擇lag)。雖然最后一種給人感覺(jué)似乎不正規(guī)，但實(shí)際上，最后一種方法也是獲得專(zhuān)家認(rèn)可的（詳見(jiàn)Anil Seth的論文對(duì)格蘭杰因果關(guān)系的描述 ）。為了保證講解的條理清晰性，這三種lag選擇方法的具體計(jì)算流程這里不再詳表，而是放到第五章代碼實(shí)現(xiàn)中直接給出。這里請(qǐng)把階數(shù)的獲得方法當(dāng)做一個(gè)黑盒過(guò)程。我們通過(guò)某種階數(shù)選擇方法得到遲延lag = m（在這里注意一個(gè)隱含條件，遲延m一定小于數(shù)據(jù)段長(zhǎng)度n，且m是實(shí)數(shù)），那么針對(duì)這個(gè)遲延lag，預(yù)測(cè)流程即為：
利用時(shí)間范圍T:1~lag上[X1~lag]的點(diǎn)預(yù)測(cè)[lag+1]上的點(diǎn)[Xlag+1]，這里Xplag+1是通過(guò)預(yù)測(cè)得出來(lái)的Xlag+1 ，Xlag+1-Xplag+1產(chǎn)生誤差ε1。
利用時(shí)間范圍T:2~lag+1上[X2~lag+1]的點(diǎn)預(yù)測(cè)[lag+2]上的點(diǎn)[Xlag+2]，這里Xplag+2是通過(guò)預(yù)測(cè)得出來(lái)的Xlag+2 ，Xlag+2-Xplag+2產(chǎn)生誤差ε2。
利用時(shí)間范圍T:3~lag+2上[X3~lag+2]的點(diǎn)預(yù)測(cè)[lag+3]上的點(diǎn)[Xlag+3]，這里Xplag+3是通過(guò)預(yù)測(cè)得出來(lái)的Xlag+3 ，Xlag+3-Xplag+3產(chǎn)生誤差ε3。
…………
利用時(shí)間范圍T:n-lag~n-1上[Xn-lag~n-1]的點(diǎn)預(yù)測(cè)[n]上的點(diǎn)[Xn ]，這里Xpn是通過(guò)預(yù)測(cè)得出來(lái)的Xn ，Xn-Xpn產(chǎn)生誤差εn-lag。
以上的過(guò)程可以看作是一個(gè)預(yù)測(cè)窗函數(shù)不斷滑動(dòng)實(shí)現(xiàn)的，如圖三：
圖3
到這里，我們便知道，所謂的預(yù)測(cè)，就是通過(guò)最小二乘法注意預(yù)測(cè)緊挨著滑動(dòng)窗口之后的那個(gè)點(diǎn)的值。與此同時(shí)產(chǎn)生預(yù)測(cè)值的預(yù)測(cè)誤差ε
2）.誤差究竟是怎樣產(chǎn)生的？
在上一節(jié)里，我們了解了預(yù)測(cè)的來(lái)源。本節(jié)里，我們解釋總體誤差的產(chǎn)生方式。從第一個(gè)問(wèn)題我們可以得知，對(duì)一個(gè)長(zhǎng)度為n的數(shù)據(jù)段，我們需要進(jìn)行n-lag次預(yù)測(cè)，每一次預(yù)測(cè)都會(huì)產(chǎn)生一個(gè)誤差項(xiàng)。而在文章的一開(kāi)始，我們提到過(guò)對(duì)一個(gè)長(zhǎng)度為n的數(shù)據(jù)段，我們只需要得出一個(gè)誤差δ1而不是n-lag個(gè)誤差。那么總體誤差δ1和這n-lag個(gè)誤差究竟有什么關(guān)系呢。答案是，δ1 是自回歸誤差ε1 ~εn-lag的無(wú)偏估計(jì)。他們之間的關(guān)系表述：
3）.X和Y是怎樣共同聯(lián)合預(yù)測(cè)的 ?
聯(lián)合回歸問(wèn)題講解開(kāi)始：
有了自回歸問(wèn)題的基礎(chǔ)，相信理解聯(lián)合回歸問(wèn)題不會(huì)是一件難事。為了不失一般性，我們?nèi)匀粡母唠A高次的聯(lián)合回歸講起，然后回到高階一次的形式。與自回歸問(wèn)題不同，現(xiàn)在我們假定擁有三個(gè)時(shí)間序列X：{x1,x2,x3,……,xn }，Y：{y1,y2,y3,……,yn }，Z：{z1,z2,z3,……,zn } 下標(biāo)代表時(shí)間T從1到n，n是實(shí)數(shù)。問(wèn)題從自回歸問(wèn)題的X = f(T)變成了Z = f(X,Y)。我們?cè)噲D找到自變量X，Y同另一個(gè)自變量Z的關(guān)系。在限定問(wèn)題為2次的情況下,針對(duì)點(diǎn)集中的第一個(gè)點(diǎn)，我們的問(wèn)題要求解的具體形式為：z1 = a1x12+a2y12+a3x1+a4y1+a5x1y1+a6+ε，擬合函數(shù)的目的就是要求出參數(shù)
a1~a6。具體的求解過(guò)程不再細(xì)說(shuō)。求解可以通過(guò)matlab函數(shù)regress實(shí)現(xiàn)。該問(wèn)題在一次狀態(tài)下的方程為
z1 = a1x12+a2y12+a3x1+a4y1+a5x1y1+a6+ε （式三） ······聯(lián)合回歸問(wèn)題
在格蘭杰因果關(guān)系中，我們將狀態(tài)方程轉(zhuǎn)化為X_new = f(X_past,Y_past)，假設(shè)lag = 3 ,則具體的形式是
xt = a1xt-1+a2xt-2+a3xt-3+a4yt-1+a5yt-2+ε（式四）······格蘭杰聯(lián)合回歸問(wèn)題（注意到?jīng)]有常數(shù)項(xiàng)）
注意到，式四中y的lag不一定等于x的lag。（備注：在我們漫談格蘭杰因果關(guān)系的最后一章所推薦的GCCN工具箱(Anil Seith著，matlab版)里，y的lag等于x的lag。）
聯(lián)合回歸問(wèn)題講解結(jié)束。
這里，我們給出格蘭杰因果關(guān)系下的一個(gè)一般性的回歸公式：
式五即為式二的一般化形式，而式六就是式四的一般化形式。
最后我們采用union-regression替代autoregression方法重復(fù)‘不同之處二’一節(jié)中敘述過(guò)的過(guò)程：
利用時(shí)間范圍T:1~lag上[X1~lag ,Y1~lag]的點(diǎn)預(yù)測(cè)[lag+1]上的點(diǎn)[Xlag+1]，這里Xplag+1是通過(guò)預(yù)測(cè)得出來(lái)的Xlag+1 ，Xlag+1-Xplag+1產(chǎn)生誤差ε1。
利用時(shí)間范圍T:2~lag+1上[X2~lag+1 ,Y2~lag+1]的點(diǎn)預(yù)測(cè)[lag+2]上的點(diǎn)[Xlag+2]，這里Xplag+2是通過(guò)預(yù)測(cè)得出來(lái)的Xlag+2 ，Xlag+2-Xplag+2產(chǎn)生誤差ε2。
利用時(shí)間范圍T:3~lag+2上[X2~lag+2 ,Y2~lag+2]的點(diǎn)預(yù)測(cè)[lag+3]上的點(diǎn)[Xlag+3]，這里Xplag+3是通過(guò)預(yù)測(cè)得出來(lái)的Xlag+3 ，Xlag+3-Xplag+3產(chǎn)生誤差ε3。
…………
利用時(shí)間范圍T:n-lag~n-1上[Xn-lag~n-1 ,Yn-lag~n-1]的點(diǎn)預(yù)測(cè)[n]上的點(diǎn)[Xn ]，這里Xpn是通過(guò)預(yù)測(cè)得出來(lái)的Xn ，Xn-Xpn產(chǎn)生誤差εn-lag。
在得到一系列的誤差后，再一次利用前文中提到的無(wú)偏估計(jì)方法求解聯(lián)合回歸產(chǎn)生的無(wú)偏誤差估計(jì)δ2。
最后，一切仿佛又回到了開(kāi)始：
假如δ2<δ1則認(rèn)為，變量Y對(duì)變量X有格蘭杰因果關(guān)系。否則則Y對(duì)變量X沒(méi)有格蘭杰因果關(guān)系。看到這一步，我們已經(jīng)弄懂了格蘭杰因果關(guān)系的絕大部分基本原理。大部分看官可以長(zhǎng)吁一口氣了。然而，事情真的完了嗎？答案是：沒(méi)有完成，還差最后一步。在得到δ2和δ1的值并比較大小之后，我們必須做自回歸和聯(lián)合回歸誤差的F檢驗(yàn)，否則無(wú)法判斷δ2和δ1值的大小是否有意義。因此，千萬(wàn)不要忘記最后一步F_test(δ2，δ1 )。至于為何要做F檢驗(yàn)和F檢驗(yàn)的具體方法，請(qǐng)參見(jiàn)筆者之前寫(xiě)的另一篇博文。我確信，如果你是一名研究人員，那篇博文值得一讀。
三：求解格蘭杰因果關(guān)系的標(biāo)準(zhǔn)化流程是什么 ？
在講過(guò)了格蘭杰因果關(guān)系的核心思想后，我們就要開(kāi)始了解一下格蘭杰因果關(guān)系的標(biāo)準(zhǔn)化流程[6]。可能有人會(huì)發(fā)出疑問(wèn)：“我們不是已經(jīng)在上面介紹了格蘭杰因果關(guān)系處理的基本流程了嗎？這里的標(biāo)準(zhǔn)化過(guò)程又是指什么？”這里的標(biāo)準(zhǔn)化流程是指為了求解時(shí)間序列間的因果關(guān)系而對(duì)時(shí)間序列數(shù)據(jù)所做一些基本預(yù)處理過(guò)程。而這個(gè)過(guò)程主要是為了保證格蘭杰因果關(guān)系檢驗(yàn)結(jié)果的有效性。下面給出格蘭杰因果關(guān)系的標(biāo)準(zhǔn)化流程的框架圖，每一步的原因以及相關(guān)的知識(shí)點(diǎn)放在之后解釋。
圖4
首先，普及一個(gè)概念。什么叫做平穩(wěn)時(shí)間序列？
粗略的講，平穩(wěn)時(shí)間序列需要滿足以下兩個(gè)條件
1）時(shí)間序列沒(méi)有變化趨勢(shì)（含周期性變化趨勢(shì)），換言之，它不可以是那種看起來(lái)有明顯增加或減少的線條。像以下兩種就絕非平穩(wěn)：
2）在時(shí)間序列中，任意取出若干段數(shù)據(jù)并對(duì)該若干段數(shù)據(jù)求方差。則方差間不應(yīng)該有顯著性差異。
STEP1：數(shù)據(jù)去趨勢(shì)化
由于格蘭杰因果關(guān)系要求所處理的時(shí)間序列必須是平穩(wěn)時(shí)間序列，所以必須要去趨勢(shì)化。一個(gè)去除了趨勢(shì)化的數(shù)據(jù)應(yīng)該如下圖所示。可以注意到，該時(shí)間序列圍繞y=3軸上下波動(dòng)。
去趨勢(shì)化的方法有很多，這里不再詳細(xì)介紹，第五節(jié)提供的matlab工具箱會(huì)有相關(guān)工具。
STEP2：數(shù)據(jù)去均值
格蘭杰因果關(guān)系要求時(shí)間序列數(shù)據(jù)必須圍繞y=0軸波動(dòng)而不是像上圖中圍繞y=3波動(dòng)，一個(gè)去過(guò)均值的時(shí)間序列應(yīng)該如下圖所示
在講解STEP3之前，率先普及一個(gè)概念
平穩(wěn)時(shí)間序列里不含有單位根，非平穩(wěn)時(shí)間序列里一定含有單位根。
STEP3：ADF檢驗(yàn)和KPSS檢驗(yàn)
ADF檢驗(yàn)全稱(chēng) Augmented-Dickey-Fuller test（增廣迪肯富樂(lè)檢驗(yàn)），KPSS檢驗(yàn)全稱(chēng) Kwiatkowski-Phillips-Schmidt-Shin test。這兩種檢驗(yàn)均稱(chēng)為單位根檢驗(yàn)。如果對(duì)時(shí)間序列數(shù)據(jù)做單位根檢驗(yàn)得出的結(jié)論是該時(shí)間序列有單位根，則可以肯定該序列一定不是平穩(wěn)的。
這里產(chǎn)生了四個(gè)問(wèn)題：
1）我們不是經(jīng)過(guò)STEP1和STEP2了嗎，為何序列仍是不平穩(wěn)的？
答：STEP1和STEP2并不保證序列的平穩(wěn)性，只能使序列看似平穩(wěn)。因此在完成STEP1和STEP2后，必須進(jìn)行STEP3以便找到時(shí)間序列平穩(wěn)數(shù)學(xué)上的證據(jù)。
2）如果時(shí)間序列經(jīng)過(guò)STEP1和STEP2后仍然檢驗(yàn)不平穩(wěn)該怎么辦？
答：?jiǎn)柕煤茫〗鉀Q這個(gè)問(wèn)題正是STEP3的任務(wù)。
3）ADF檢驗(yàn)和KPSS檢驗(yàn)有什么區(qū)別？
KPSS 是右側(cè)單邊檢驗(yàn)，原假設(shè)是序列是平穩(wěn)的（不存在單位根）。
ADF是雙邊檢驗(yàn)，原假設(shè)是序列是非平穩(wěn)的（存在單位根）。
4）ADF檢驗(yàn)和KPSS檢驗(yàn)結(jié)果不一致怎么辦？
不得不說(shuō)，這事兒常發(fā)生。筆者在網(wǎng)上也好好搜索了一番（2017年7月25日搜索）。網(wǎng)友給出的結(jié)論基本分為兩派：
A."看情況自己選"
B."有一個(gè)通過(guò)就算通過(guò)"
STEP4：一階差分
正如同STEP3中所提到的一樣，只有當(dāng)經(jīng)過(guò)STEP1和STEP2后仍然檢驗(yàn)不平穩(wěn)時(shí)才會(huì)對(duì)時(shí)間序列數(shù)據(jù)做差分。那么怎么做差分呢？
給定時(shí)間序列X{x1,x2,x3,……,xn }。則定義差分過(guò)的時(shí)間序列X_diff為{x2-x1, x3-x2, x4-x3, ......, xn-xn-1}。說(shuō)白了就是后項(xiàng)減前項(xiàng)，這就是差分。差分帶來(lái)的唯一問(wèn)題是經(jīng)過(guò)差分后的時(shí)間序列長(zhǎng)度會(huì)縮減1個(gè)單位。另一個(gè)問(wèn)題是，如果時(shí)間序列仍然不是平穩(wěn)的該怎么辦？答案是再差分一次知道平穩(wěn)為止。根據(jù)筆者的經(jīng)驗(yàn)，一般差分一次就足夠了，很少有出現(xiàn)差分兩次的。
STEP5：通過(guò)AIC準(zhǔn)則求階數(shù)lag
AIC信息準(zhǔn)則的公式是：AIC=2*lag+nln(RSS/n)。其中，RSS是擬合的殘差和（就是之前的誤差無(wú)偏估計(jì)量）, n是數(shù)據(jù)段的長(zhǎng)度。我們的目標(biāo)是搜索在1~n的范圍內(nèi)搜索lag，使AIC的值取得最小。而那個(gè)取得AIC最小值的lag就是我們想要的階數(shù)lag
STEP6：杜爾濱懷特訊檢驗(yàn)
英文名稱(chēng)為Durbin-Watson test。相必大家還記得我們采用最小二乘法求回歸問(wèn)題時(shí)曾假設(shè)回歸后理論值與實(shí)際值的誤差為ε，且ε服從正態(tài)分布。實(shí)際上，誤差服從正態(tài)分布是使用最小二乘法求解回歸問(wèn)題的先決條件。而杜爾濱懷特訊檢驗(yàn)的目的就是檢測(cè)回歸完成后的殘差是否服從正態(tài)分布。如果不服從則該段數(shù)據(jù)不滿足使用最小二乘法的先決條件，也就不滿足求解格蘭杰因果關(guān)系的基礎(chǔ)。
STEP7：一致性檢驗(yàn)
當(dāng)對(duì)時(shí)間序列的數(shù)據(jù)值點(diǎn)完成回歸以后。并不能確定回歸得出的理論值和實(shí)際值是否來(lái)自于同一分布。此時(shí)，應(yīng)采用一致性檢驗(yàn)。如果一致性檢驗(yàn)的結(jié)論表明理論值和實(shí)際值差距較大，則回歸結(jié)果失敗，需要重新確定回歸。
到這里，格蘭杰因果關(guān)系的基本原理部分相信大家已經(jīng)明白了。如果還有不明白的問(wèn)題以及其他需求，請(qǐng)?jiān)谡緝?nèi)給筆者留言。筆者將盡快給予相應(yīng)的回復(fù)。后續(xù)部分筆者也將抽時(shí)間盡快完成。除此之外，本文中如有其它不妥之處，希望諸位看官不吝賜教，小生愿洗耳恭聽(tīng)。
[1] Testing for Linear and Nonlinear Granger Causality in the Stock Price-Volume Relation
[2] Spatial-temporal causal modeling for climate change attribution
[3] Mapping directed influence over the brain using Granger causality and fMRI
[4] 維基百科，“格蘭杰因果關(guān)系 ”詞條中格蘭杰原話：“ Of course, many ridiculous papers appeared”
[5] Causality connectivity of evolved neural networks during behavior
[6] RESTING-STATE BRAIN NETWORKS REVEALED BY GRANGER CAUSAL CONNECTIVITY IN FROGS
浙公網(wǎng)安備 33010602011771號(hào)