グレンジャー因果検定

データ解析

グレンジャー因果検定：

  • 注意事項

＃本当の因果関係でないことに注意。
例えば，精度のよい天気予報と天気は因果関係ありと判定してしまう．
また，線形モデルを仮定するので，モデルが非線形の場合は正確ではない．

x→yの相関を考える．

y(t+1)を予測するときの予測誤差を考える．
このとき，予測モデルf1とf2に関して，
f1(y(1),y(2),…,y(t))の予測誤差σ1
とf2(y(1),y(2),…,y(t),x(1),x(2),…,x(t))の予測誤差σ2
を比較したとき，σ1>σ2なら相関（グレンジャー因果）ありとする．

  • やり方，

モデルはVARモデルとする．
モデル pラグ　T：データ数
y(t+1)=a1y(t)+a2y(t-1)+…apt(y-p+1)
y(t+1)=c1y(t)+c2y(t-1)+…cp(y-p+1)+b1y(t)+b2y(t-1)+…bpt(y-p+1)

の誤差合計から以下のF値を検定する．
F = *1

ここでSSEr=Σσ2
SSEu=Σσ1

  • 検定

帰無仮説：相関なし
仮説：相関あり

F値が0に近い→p値が大きい→相関なし！
F値が0から離れてる→p値が小さい→相関あり！

  • Rでの実装
grangertest(x,y)

とか

granger.test(cbind(x,y))

両方ともデフォルトでは次数１であることに注意．

import scala.math._
import org.apache.commons.math3.stat.regression.OLSMultipleLinearRegression;
import org.apache.commons.math3.distribution.FDistribution;
import scala.io.Source
object Granger_test{

        def granger_func(x:Array[Double], y:Array[Double],p:Int=1):Array[Double]={
                var n=x.length
                var x1=Array.ofDim[Double](n-p,p);
                var x2=Array.ofDim[Double](n-p,2*p);
                var y1=y.slice(p,n)

                (0 to (p-1)).foreach{i=>

                        ((0) to (n-p-1)).foreach{j=>


                                x1(j)(i)=y(j+i)

                                x2(j)(2*i)=y(j+i)

                                x2(j)(2*i+1)=x(j+i)
                        }
                }


                var reg1=new OLSMultipleLinearRegression();
                var reg2=new OLSMultipleLinearRegression();
                reg1.newSampleData(y1,x1)
                reg2.newSampleData(y1,x2)

                var resi1=reg1.estimateResiduals();
                var resi2=reg2.estimateResiduals();
                var ssqU=resi2.map{a=> a*a}.sum;
                var ssqR=resi1.map{a=> a*a}.sum;
              
                var ftest=((ssqR-ssqU)/p)/(ssqU/(y1.length-2*p-1))
                var f=new FDistribution(p,y1.length-2*p-1);
                var pvalue=1-f.cumulativeProbability(ftest);
                Array(ftest,pvalue)


        }

        def main(args: Array[String]){
                //var x=Array(1.0,2.0,3.0,4.0,5.0,4.0,5.0)
                //var y=Array(1.0,2.0,3.0,4.0,6.0,9.0,3.0)
                var j=Source.fromFile("test.txt","utf-8").getLines.toArray

                var x=j.map{s=> s.split(" ")(0).toDouble}
                var y=j.map{s=> s.split(" ")(1).toDouble}

                var j2=granger_func(x,y,1);
                println("f-statistic:"+j2(0).toString+" P-value:"+j2(1).toString);
        }





}

test.txtの例

1 1
2 2
3 3
4 4
5 6
4 9
5 3

参考Johansen rank 検定法　
http://web.iss.u-tokyo.ac.jp/~matsumur/UT5108A.ppt
http://blog.livedoor.jp/yukismd/archives/1096939.html

*1:SSEr - SSEu) / p )÷ (SSEu / (T-2p-1