漫谈格兰杰因果关系（Granger Causality）——第一章 野火烧不尽，春风吹又生
2017年7月9日上午6点10分，先师胡三清同志——新因果关系的提出者、植入式脑部电极癫痫治疗法的提出者、IEEE高级会员，因肺癌医治无效于杭州肿瘤医院去世，享年50岁。余蒙先师厚恩数载，一朝忽闻先师驾鹤西归，悲痛不已。瘁心之余，遂决意传先师之道，以慰先师在天之灵。如此，先师盖以瞑目矣！
格兰杰因果关系作为一种可以衡量时间序列之间相互影响关系的方法，最近十几年备受青睐。无论是经济学[1]，气象科学[2]，神经科学[3]都有广泛的应用，尽管后两者（气象和神经科学）连格兰杰自己都反对（格兰杰反对将格兰杰因果关系用在除经济学以外的其他领域，这就是本文题目所谓的“野火”）[4]。鉴于笔者从未在气象学有过半分建树，所以不敢妄谈。不过庆幸的是，经过神经科学家数十载的辛苦“洗地”，他们纷纷找到了自己‘合法’使用格兰杰因果关系的理由[5]（Anil Seth是英国皇家科学院院士，笔者首推的‘地表最强洗地王’，也即本文题目所谓的“春风”）。除了由克里夫·格兰杰本人提出的格兰杰因果关系之外，还有数种围绕格兰杰因果关系方法产生的变体，本文也将对这些变体分门别类，做出一些简介。当然，也正是因为胡先生认为应当遵从格兰杰爵士的论文，所以才创立了‘新因果关系’专门解决神经科学中的因果关系使用问题，不过这是后话了。本文的讲述将围绕以下几个方面进行：
本文的框架：
第一章 野火烧不尽，春风吹又生
一：格兰杰因果关系究竟是什么？
二：格兰杰因果关系原理是什么？
三：求解格兰杰因果关系的标准化流程是什么？
第二章 星星因果，可以燎原
四：格兰杰因果关系变种、原理及其应用领域？
五：正宗的格兰杰因果关系怎样用代码实现？（MATLAB版）
第三章 胡氏因果关系（新因果关系）（Innovation by Sanqing Hu）
一：格兰杰因果关系究竟是什么？
一般口头上所称的格兰杰因果关系全称是“格兰杰因果关系检验”（Granger causality test），看到这里，诸位看官大可以感叹一句“啊哈！原来格兰杰因果关系是统计假设检验的一种！”。那么“因果”二字又要从何谈起呢？这就要从时间序列间的关系开始说起。格兰杰是这样定义“因果”的。现在我们假设有一个时间序列X，它是由不同时间的采样点{x1,x2,x3,……,xn}共同组成的。同X一样，我们假设有时间序列Y，它形如X，由{y1,y2,y3,……,yn}共同组成。现在我们利用X的过去预测X的未来，比如用x1~xn-j（这就是X的过去）预测xn-j+1~xn（这就是X的未来），预测的过程中产生一个误差δ1（请先忽略预测的具体方法和误差产生的具体方法），然后把这个误差视为我们得到的第一个结果。再然后利用X和Y共同的过去预测X的未来，比如用{x1~xn-j（这就是X的过去）|y1~yn-j（这就是Y的过去）}预测xn-j+1~xn（这就是X的未来），预测的过程中产生一个误差δ2（请先忽略“X和Y是怎样共同的 ?”这个问题） ，然后把这个误差视为我们得到的第二个结果。如果δ1小于δ2，也就是说X和Y的联合预测误差小于X自身的预测误差，那么，必然是因为Y对X的预测起到了帮助，所以才减小了预测误差。在这种情况下，我们称Y对X有格兰杰因果关系。
二：格兰杰因果关系原理是什么？
这一节里，我们要重点解决上一节中出现的几个问题。即1）.预测究竟是怎样实现的？；2）.误差究竟是怎样产生的？；3）X和Y是怎样共同联合预测的 ? 在讲解这些问题之前。有一个非常重要的概念——回归问题，必须预先铺垫，否则，这三个问题无从下手。为了方便理解，我们在二维空间中表述这个问题，当然这个问题也可以被拓展到高维空间。当然，如果某位看官确保自己已经弄清楚了回归问题，那么可以直接跳到回归问题讲完之后的部分。
自回归问题讲解开始：
现假设存在一系列二维空间点的集合S={x1,x2,x3,……,xn }。如图1所示。现在我们想找到一条线，这条线最好穿过点集S中的每一个点。这些点包含一个隐变量T。这样，当找到这条线的时候，我们实际上找到了代表这条二维空间中的线的一个表示函数x = f(t)。
图1
那么怎样找到这条线呢？
STEP1，你需要假定这条线来源于一个阶数为m的函数。如果你的假设是阶数为3，那么这个函数就是一个形如x = f(t) = a0+a1t1+a2t2+a3t3 的函数。集合S中的每一个点都必须要满足（或尽力满足）x = f(t)这个公式。就拿第一个点[t1,x1]举个栗子：x1 = f(t) = a0+a1t11+a2t12+a3t13。当然，我们知道这条线y=f(x)完美的经过每一个点近乎是不可能的，因此我们允许误差的存在。只要y = f(x)+ε 即可，其中ε代表误差项，并且ε服从正态分布（有时也称高斯白噪声，高斯分布）。
STEP2，尝试求出系数a0~am满足STEP1的要求，使得对于点集S中所有的点满足x = f(t)+ε。一般情况下，求解系数a0~am 最常见的方法就是最小二乘法。最小二乘法的具体步骤这里不再描述。一般可以通过matlab函数polyfit进行求解，代码如下。求解结果如图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就是拟合参数，它负责根据时间变化1:9以及自变量X将参数a求解出来，注意，a在此处是个数组，由a0~a3组成。polyval函数根据参数数组a和点集T_new求解X_new。而绘制[T_new, X_new]就是图2中的蓝色曲线。
图2
可能细心的同学能发现，图1和图2中，时间轴T的值域不一样。在图1中，T轴的值域是1~9；图2中，x轴的值域1~10。而这多余出来的X10其实可以被视为我们根据已有数据集[X1-9]预测出来的未来的数据[X10]。
总之，自回归问题的一般形式就是给出一系列点集[T,X]，然后想办法找到一个函数X = f(T) 。我们要想法找到一系列参数a，使得函数代表的这条线尽可能经过点集。
自回归问题讲解结束。
讲完了自回归问题，以下两个问题就可以得到完美的解答：
同本文的开始相同，现在我们假定拥有两个时间序列X：{x1,x2,x3,……,xn }下标代表时间T从1到n，n是实数。
1）.预测究竟是怎样实现的？
读到这里时，我假设各位看官已经弄清楚了我在上面所述的回归问题。但是格兰杰因果关系中的回归问题同原问题的原理相同但过程有所不同。且听我慢慢道来他们的不同之处：
不同之处一：在回归的原问题中，我们以时间轴T为横坐标尽力求解自变量X。试图找到时间轴T同自变量X之间的关系。即X = f(T)。比如我们在上式中举出的例子里X = f(T) = a0+a1T1+a2T2+a3T3（请注意：这个问题实际上是三阶问题，'阶' 是指除了a0之外，还有几个参数a。例子中的回归问题其实是一个非线性回归（次数为一才是线性的））。在回归问题的末尾，我们说过，回归问题的本质是要找到一系列的参数a，因此，我们在例子中举到的3阶形式，存在的具体形式如下：
X = f(T) = a0+a1t11+a2t12+a3t13（式一）······普通回归问题（三阶一次问题（线性回归））
在格兰杰因果关系中，自回归模型就是类似于这种的一阶形式。但是它不再是T映射到X上的函数，而是过去的X映射到现在的X上的函数，即
X_now = f(X_past)，在考虑过去三个点的情况下，这个问题就变成了
xt =a1xt-1+a2xt-2 +a3xt-3（式二）······格兰杰回归问题（注意到没有常数项）
既然在式一中，知晓T和X，进而求解一些列参数a不成问题，那么在这里，知晓时间序列上的X，进而求解一系列参数a也不是问题。
不同之处二：
既然是回归问题就要涉及到回归阶数order的选择问题。通过不同之处一，我们知晓，在格兰杰因果关系检验中，阶数的作用更像是在指明采用过去的多少个点对当前点的点求解回归问题。因此阶数order有了一个崭新的名字：lag（迟延，迟滞）。在阐述回归问题的时候，我们假设迟滞lag=3，然后使用最小二乘法拟合了点集S。实际上，为了简化最小二乘拟合的问题，笔者武断的设置了迟滞lag=3而并没有采取适当地方法对lag应有的值做出选择。在格兰杰因果关系的计算过程中，针对阶数的选择方法大致可以分为三类：I. Akaike information criterion(赤池信息准则，简称AIC)；II. Bayes information criterion(贝叶斯信息准则，简称BIC)；III Rule of thumb(凭经验选择lag)。虽然最后一种给人感觉似乎不正规，但实际上，最后一种方法也是获得专家认可的（详见Anil Seth的论文对格兰杰因果关系的描述 ）。为了保证讲解的条理清晰性，这三种lag选择方法的具体计算流程这里不再详表，而是放到第五章代码实现中直接给出。这里请把阶数的获得方法当做一个黑盒过程。我们通过某种阶数选择方法得到迟延lag = m（在这里注意一个隐含条件，迟延m一定小于数据段长度n，且m是实数），那么针对这个迟延lag，预测流程即为：
利用时间范围T:1~lag上[X1~lag]的点预测[lag+1]上的点[Xlag+1]，这里Xplag+1是通过预测得出来的Xlag+1 ，Xlag+1-Xplag+1产生误差ε1。
利用时间范围T:2~lag+1上[X2~lag+1]的点预测[lag+2]上的点[Xlag+2]，这里Xplag+2是通过预测得出来的Xlag+2 ，Xlag+2-Xplag+2产生误差ε2。
利用时间范围T:3~lag+2上[X3~lag+2]的点预测[lag+3]上的点[Xlag+3]，这里Xplag+3是通过预测得出来的Xlag+3 ，Xlag+3-Xplag+3产生误差ε3。
…………
利用时间范围T:n-lag~n-1上[Xn-lag~n-1]的点预测[n]上的点[Xn ]，这里Xpn是通过预测得出来的Xn ，Xn-Xpn产生误差εn-lag。
以上的过程可以看作是一个预测窗函数不断滑动实现的，如图三：
图3
到这里，我们便知道，所谓的预测，就是通过最小二乘法注意预测紧挨着滑动窗口之后的那个点的值。与此同时产生预测值的预测误差ε
2）.误差究竟是怎样产生的？
在上一节里，我们了解了预测的来源。本节里，我们解释总体误差的产生方式。从第一个问题我们可以得知，对一个长度为n的数据段，我们需要进行n-lag次预测，每一次预测都会产生一个误差项。而在文章的一开始，我们提到过对一个长度为n的数据段，我们只需要得出一个误差δ1而不是n-lag个误差。那么总体误差δ1和这n-lag个误差究竟有什么关系呢。答案是，δ1 是自回归误差ε1 ~εn-lag的无偏估计。他们之间的关系表述：
3）.X和Y是怎样共同联合预测的 ?
联合回归问题讲解开始：
有了自回归问题的基础，相信理解联合回归问题不会是一件难事。为了不失一般性，我们仍然从高阶高次的联合回归讲起，然后回到高阶一次的形式。与自回归问题不同，现在我们假定拥有三个时间序列X：{x1,x2,x3,……,xn }，Y：{y1,y2,y3,……,yn }，Z：{z1,z2,z3,……,zn } 下标代表时间T从1到n，n是实数。问题从自回归问题的X = f(T)变成了Z = f(X,Y)。我们试图找到自变量X，Y同另一个自变量Z的关系。在限定问题为2次的情况下,针对点集中的第一个点，我们的问题要求解的具体形式为：z1 = a1x12+a2y12+a3x1+a4y1+a5x1y1+a6+ε，拟合函数的目的就是要求出参数
a1~a6。具体的求解过程不再细说。求解可以通过matlab函数regress实现。该问题在一次状态下的方程为
z1 = a1x12+a2y12+a3x1+a4y1+a5x1y1+a6+ε （式三） ······联合回归问题
在格兰杰因果关系中，我们将状态方程转化为X_new = f(X_past,Y_past)，假设lag = 3 ,则具体的形式是
xt = a1xt-1+a2xt-2+a3xt-3+a4yt-1+a5yt-2+ε（式四）······格兰杰联合回归问题（注意到没有常数项）
注意到，式四中y的lag不一定等于x的lag。（备注：在我们漫谈格兰杰因果关系的最后一章所推荐的GCCN工具箱(Anil Seith著，matlab版)里，y的lag等于x的lag。）
联合回归问题讲解结束。
这里，我们给出格兰杰因果关系下的一个一般性的回归公式：
式五即为式二的一般化形式，而式六就是式四的一般化形式。
最后我们采用union-regression替代autoregression方法重复‘不同之处二’一节中叙述过的过程：
利用时间范围T:1~lag上[X1~lag ,Y1~lag]的点预测[lag+1]上的点[Xlag+1]，这里Xplag+1是通过预测得出来的Xlag+1 ，Xlag+1-Xplag+1产生误差ε1。
利用时间范围T:2~lag+1上[X2~lag+1 ,Y2~lag+1]的点预测[lag+2]上的点[Xlag+2]，这里Xplag+2是通过预测得出来的Xlag+2 ，Xlag+2-Xplag+2产生误差ε2。
利用时间范围T:3~lag+2上[X2~lag+2 ,Y2~lag+2]的点预测[lag+3]上的点[Xlag+3]，这里Xplag+3是通过预测得出来的Xlag+3 ，Xlag+3-Xplag+3产生误差ε3。
…………
利用时间范围T:n-lag~n-1上[Xn-lag~n-1 ,Yn-lag~n-1]的点预测[n]上的点[Xn ]，这里Xpn是通过预测得出来的Xn ，Xn-Xpn产生误差εn-lag。
在得到一系列的误差后，再一次利用前文中提到的无偏估计方法求解联合回归产生的无偏误差估计δ2。
最后，一切仿佛又回到了开始：
假如δ2<δ1则认为，变量Y对变量X有格兰杰因果关系。否则则Y对变量X没有格兰杰因果关系。看到这一步，我们已经弄懂了格兰杰因果关系的绝大部分基本原理。大部分看官可以长吁一口气了。然而，事情真的完了吗？答案是：没有完成，还差最后一步。在得到δ2和δ1的值并比较大小之后，我们必须做自回归和联合回归误差的F检验，否则无法判断δ2和δ1值的大小是否有意义。因此，千万不要忘记最后一步F_test(δ2，δ1 )。至于为何要做F检验和F检验的具体方法，请参见笔者之前写的另一篇博文。我确信，如果你是一名研究人员，那篇博文值得一读。
三：求解格兰杰因果关系的标准化流程是什么 ？
在讲过了格兰杰因果关系的核心思想后，我们就要开始了解一下格兰杰因果关系的标准化流程[6]。可能有人会发出疑问：“我们不是已经在上面介绍了格兰杰因果关系处理的基本流程了吗？这里的标准化过程又是指什么？”这里的标准化流程是指为了求解时间序列间的因果关系而对时间序列数据所做一些基本预处理过程。而这个过程主要是为了保证格兰杰因果关系检验结果的有效性。下面给出格兰杰因果关系的标准化流程的框架图，每一步的原因以及相关的知识点放在之后解释。
图4
首先，普及一个概念。什么叫做平稳时间序列？
粗略的讲，平稳时间序列需要满足以下两个条件
1）时间序列没有变化趋势（含周期性变化趋势），换言之，它不可以是那种看起来有明显增加或减少的线条。像以下两种就绝非平稳：
2）在时间序列中，任意取出若干段数据并对该若干段数据求方差。则方差间不应该有显著性差异。
STEP1：数据去趋势化
由于格兰杰因果关系要求所处理的时间序列必须是平稳时间序列，所以必须要去趋势化。一个去除了趋势化的数据应该如下图所示。可以注意到，该时间序列围绕y=3轴上下波动。
去趋势化的方法有很多，这里不再详细介绍，第五节提供的matlab工具箱会有相关工具。
STEP2：数据去均值
格兰杰因果关系要求时间序列数据必须围绕y=0轴波动而不是像上图中围绕y=3波动，一个去过均值的时间序列应该如下图所示
在讲解STEP3之前，率先普及一个概念
平稳时间序列里不含有单位根，非平稳时间序列里一定含有单位根。
STEP3：ADF检验和KPSS检验
ADF检验全称 Augmented-Dickey-Fuller test（增广迪肯富乐检验），KPSS检验全称 Kwiatkowski-Phillips-Schmidt-Shin test。这两种检验均称为单位根检验。如果对时间序列数据做单位根检验得出的结论是该时间序列有单位根，则可以肯定该序列一定不是平稳的。
这里产生了四个问题：
1）我们不是经过STEP1和STEP2了吗，为何序列仍是不平稳的？
答：STEP1和STEP2并不保证序列的平稳性，只能使序列看似平稳。因此在完成STEP1和STEP2后，必须进行STEP3以便找到时间序列平稳数学上的证据。
2）如果时间序列经过STEP1和STEP2后仍然检验不平稳该怎么办？
答：问得好！解决这个问题正是STEP3的任务。
3）ADF检验和KPSS检验有什么区别？
KPSS 是右侧单边检验，原假设是序列是平稳的（不存在单位根）。
ADF是双边检验，原假设是序列是非平稳的（存在单位根）。
4）ADF检验和KPSS检验结果不一致怎么办？
不得不说，这事儿常发生。笔者在网上也好好搜索了一番（2017年7月25日搜索）。网友给出的结论基本分为两派：
A."看情况自己选"
B."有一个通过就算通过"
STEP4：一阶差分
正如同STEP3中所提到的一样，只有当经过STEP1和STEP2后仍然检验不平稳时才会对时间序列数据做差分。那么怎么做差分呢？
给定时间序列X{x1,x2,x3,……,xn }。则定义差分过的时间序列X_diff为{x2-x1, x3-x2, x4-x3, ......, xn-xn-1}。说白了就是后项减前项，这就是差分。差分带来的唯一问题是经过差分后的时间序列长度会缩减1个单位。另一个问题是，如果时间序列仍然不是平稳的该怎么办？答案是再差分一次知道平稳为止。根据笔者的经验，一般差分一次就足够了，很少有出现差分两次的。
STEP5：通过AIC准则求阶数lag
AIC信息准则的公式是：AIC=2*lag+nln(RSS/n)。其中，RSS是拟合的残差和（就是之前的误差无偏估计量）, n是数据段的长度。我们的目标是搜索在1~n的范围内搜索lag，使AIC的值取得最小。而那个取得AIC最小值的lag就是我们想要的阶数lag
STEP6：杜尔滨怀特讯检验
英文名称为Durbin-Watson test。相必大家还记得我们采用最小二乘法求回归问题时曾假设回归后理论值与实际值的误差为ε，且ε服从正态分布。实际上，误差服从正态分布是使用最小二乘法求解回归问题的先决条件。而杜尔滨怀特讯检验的目的就是检测回归完成后的残差是否服从正态分布。如果不服从则该段数据不满足使用最小二乘法的先决条件，也就不满足求解格兰杰因果关系的基础。
STEP7：一致性检验
当对时间序列的数据值点完成回归以后。并不能确定回归得出的理论值和实际值是否来自于同一分布。此时，应采用一致性检验。如果一致性检验的结论表明理论值和实际值差距较大，则回归结果失败，需要重新确定回归。
到这里，格兰杰因果关系的基本原理部分相信大家已经明白了。如果还有不明白的问题以及其他需求，请在站内给笔者留言。笔者将尽快给予相应的回复。后续部分笔者也将抽时间尽快完成。除此之外，本文中如有其它不妥之处，希望诸位看官不吝赐教，小生愿洗耳恭听。
[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] 维基百科，“格兰杰因果关系 ”词条中格兰杰原话：“ 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