漫谈格兰杰因果关系（Granger Causality） |
您所在的位置：网站首页 › 格兰杰的左轮有必要捡吗 › 漫谈格兰杰因果关系（Granger Causality） |
漫谈格兰杰因果关系（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 |
【本文地址】
今日新闻 |
推荐新闻 |