R version 2.9.0 (2009-04-17) Copyright (C) 2009 The R Foundation for Statistical Computing ISBN 3-900051-07-0 R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit R. > y <- c(98.2,98.7,113.3,104.6,99.3,111.8,97.3,97.7,115.6,111.9,107.0,107.1,100.6,99.2,108.4,103.0,99.8,115.0,90.8,95.9,114.4,108.2,112.6,109.1,105.0,105.0,118.5,103.7,112.5,116.6,96.6,101.9,116.5,119.3,115.4,108.5,111.5,108.8,121.8,109.6,112.2,119.6,104.1,105.3,115.0,124.1,116.8,107.5,115.6,116.2,116.3,119.0,111.9,118.6,106.9,103.2,118.6,118.7,102.8,100.6) > x <- c(95.1,97.0,112.7,102.9,97.4,111.4,87.4,96.8,114.1,110.3,103.9,101.6,94.6,95.9,104.7,102.8,98.1,113.9,80.9,95.7,113.2,105.9,108.8,102.3,99.0,100.7,115.5,100.7,109.9,114.6,85.4,100.5,114.8,116.5,112.9,102.0,106.0,105.3,118.8,106.1,109.3,117.2,92.5,104.2,112.5,122.4,113.3,100.0,110.7,112.8,109.8,117.3,109.1,115.9,96.0,99.8,116.8,115.7,99.4,94.3) > par8 = '3' > par7 = '0' > par6 = '0' > par5 = '2.0' > par4 = '12' > par3 = '0' > par2 = '0' > par1 = '2.0' > #'GNU S' R Code compiled by R2WASP v. 1.0.44 () > #Author: Prof. Dr. P. Wessa > #To cite this work: Wessa P., (2008), Bivariate Granger Causality (v1.0.0) in Free Statistics Software (v$_version), Office for Research Development and Education, URL http://www.wessa.net/rwasp_grangercausality.wasp#output/ > #Source of accompanying publication: Office for Research, Development, and Education > #Technical description: > library(lmtest) Loading required package: zoo Attaching package: 'zoo' The following object(s) are masked from package:base : as.Date.numeric > par1 <- as.numeric(par1) > par2 <- as.numeric(par2) > par3 <- as.numeric(par3) > par4 <- as.numeric(par4) > par5 <- as.numeric(par5) > par6 <- as.numeric(par6) > par7 <- as.numeric(par7) > par8 <- as.numeric(par8) > ox <- x > oy <- y > if (par1 == 0) { + x <- log(x) + } else { + x <- (x ^ par1 - 1) / par1 + } > if (par5 == 0) { + y <- log(y) + } else { + y <- (y ^ par5 - 1) / par5 + } > if (par2 > 0) x <- diff(x,lag=1,difference=par2) > if (par6 > 0) y <- diff(y,lag=1,difference=par6) > if (par3 > 0) x <- diff(x,lag=par4,difference=par3) > if (par7 > 0) y <- diff(y,lag=par4,difference=par7) > x [1] 4521.505 4704.000 6350.145 5293.705 4742.880 6204.480 3818.880 4684.620 [9] 6508.905 6082.545 5397.105 5160.780 4474.080 4597.905 5480.545 5283.420 [17] 4811.305 6486.105 3271.905 4578.745 6406.620 5606.905 5918.220 5232.145 [25] 4900.000 5069.745 6669.625 5069.745 6038.505 6566.080 3646.080 5049.625 [33] 6589.020 6785.625 6372.705 5201.500 5617.500 5543.545 7056.220 5628.105 [41] 5972.745 6867.420 4277.625 5428.320 6327.625 7490.380 6417.945 4999.500 [49] 6126.745 6361.420 6027.520 6879.145 5950.905 6715.905 4607.500 4979.520 [57] 6820.620 6692.745 4939.680 4445.745 > y [1] 4821.120 4870.345 6417.945 5470.080 4929.745 6249.120 4733.145 4772.145 [9] 6681.180 6260.305 5724.000 5734.705 5059.680 4919.820 5874.780 5304.000 [17] 4979.520 6612.000 4121.820 4597.905 6543.180 5853.120 6338.880 5950.905 [25] 5512.000 5512.000 7020.625 5376.345 6327.625 6797.280 4665.280 5191.305 [33] 6785.625 7115.745 6658.080 5885.625 6215.625 5918.220 7417.120 6005.580 [41] 6293.920 7151.580 5417.905 5543.545 6612.000 7699.905 6820.620 5777.625 [49] 6681.180 6750.720 6762.345 7080.000 6260.305 7032.480 5713.305 5324.620 [57] 7032.480 7044.345 5283.420 5059.680 > (gyx <- grangertest(y ~ x, order=par8)) Granger causality test Model 1: ~ Lags(, 1:3) + Lags(, 1:3) Model 2: ~ Lags(, 1:3) Res.Df Df F Pr(>F) 1 50 2 53 -3 9.3078 5.393e-05 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 > (gxy <- grangertest(x ~ y, order=par8)) Granger causality test Model 1: ~ Lags(, 1:3) + Lags(, 1:3) Model 2: ~ Lags(, 1:3) Res.Df Df F Pr(>F) 1 50 2 53 -3 9.014 7.1e-05 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 > postscript(file="/var/www/html/rcomp/tmp/1cj6a1260794230.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556) > op <- par(mfrow=c(2,1)) > (r <- ccf(ox,oy,main='Cross Correlation Function (raw data)',ylab='CCF',xlab='Lag (k)')) Autocorrelations of series 'X', by lag -14 -13 -12 -11 -10 -9 -8 -7 -6 -5 -4 -0.216 0.096 0.604 0.018 -0.171 0.065 -0.021 0.199 0.333 0.193 0.082 -3 -2 -1 0 1 2 3 4 5 6 7 0.138 -0.164 0.239 0.957 0.034 -0.136 0.208 0.033 0.244 0.349 0.147 8 9 10 11 12 13 14 -0.005 0.011 -0.171 0.203 0.617 -0.063 -0.180 > (r <- ccf(x,y,main='Cross Correlation Function (transformed and differenced)',ylab='CCF',xlab='Lag (k)')) Autocorrelations of series 'X', by lag -14 -13 -12 -11 -10 -9 -8 -7 -6 -5 -4 -0.224 0.086 0.602 0.028 -0.171 0.064 -0.007 0.200 0.329 0.199 0.098 -3 -2 -1 0 1 2 3 4 5 6 7 0.144 -0.165 0.227 0.962 0.041 -0.147 0.212 0.049 0.238 0.343 0.158 8 9 10 11 12 13 14 0.006 0.013 -0.162 0.196 0.612 -0.056 -0.198 > par(op) > dev.off() null device 1 > postscript(file="/var/www/html/rcomp/tmp/2vo1x1260794230.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556) > op <- par(mfrow=c(2,1)) > acf(ox,lag.max=round(length(x)/2),main='ACF of x (raw)') > acf(x,lag.max=round(length(x)/2),main='ACF of x (transformed and differenced)') > par(op) > dev.off() null device 1 > postscript(file="/var/www/html/rcomp/tmp/3xj8d1260794230.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556) > op <- par(mfrow=c(2,1)) > acf(oy,lag.max=round(length(y)/2),main='ACF of y (raw)') > acf(y,lag.max=round(length(y)/2),main='ACF of y (transformed and differenced)') > par(op) > dev.off() null device 1 > > #Note: the /var/www/html/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab > load(file="/var/www/html/rcomp/createtable") > > a<-table.start() > a<-table.row.start(a) > a<-table.element(a,'Granger Causality Test: Y = f(X)',5,TRUE) > a<-table.row.end(a) > a<-table.row.start(a) > a<-table.element(a,'Model',header=TRUE) > a<-table.element(a,'Res.DF',header=TRUE) > a<-table.element(a,'Diff. DF',header=TRUE) > a<-table.element(a,'F',header=TRUE) > a<-table.element(a,'p-value',header=TRUE) > a<-table.row.end(a) > a<-table.row.start(a) > a<-table.element(a,'Complete model',header=TRUE) > a<-table.element(a,gyx$Res.Df[1]) > a<-table.element(a,'') > a<-table.element(a,'') > a<-table.element(a,'') > a<-table.row.end(a) > a<-table.row.start(a) > a<-table.element(a,'Reduced model',header=TRUE) > a<-table.element(a,gyx$Res.Df[2]) > a<-table.element(a,gyx$Df[2]) > a<-table.element(a,gyx$F[2]) > a<-table.element(a,gyx$Pr[2]) > a<-table.row.end(a) > a<-table.end(a) > table.save(a,file="/var/www/html/rcomp/tmp/46w3h1260794230.tab") > a<-table.start() > a<-table.row.start(a) > a<-table.element(a,'Granger Causality Test: X = f(Y)',5,TRUE) > a<-table.row.end(a) > a<-table.row.start(a) > a<-table.element(a,'Model',header=TRUE) > a<-table.element(a,'Res.DF',header=TRUE) > a<-table.element(a,'Diff. DF',header=TRUE) > a<-table.element(a,'F',header=TRUE) > a<-table.element(a,'p-value',header=TRUE) > a<-table.row.end(a) > a<-table.row.start(a) > a<-table.element(a,'Complete model',header=TRUE) > a<-table.element(a,gxy$Res.Df[1]) > a<-table.element(a,'') > a<-table.element(a,'') > a<-table.element(a,'') > a<-table.row.end(a) > a<-table.row.start(a) > a<-table.element(a,'Reduced model',header=TRUE) > a<-table.element(a,gxy$Res.Df[2]) > a<-table.element(a,gxy$Df[2]) > a<-table.element(a,gxy$F[2]) > a<-table.element(a,gxy$Pr[2]) > a<-table.row.end(a) > a<-table.end(a) > table.save(a,file="/var/www/html/rcomp/tmp/53uho1260794231.tab") > > try(system("convert tmp/1cj6a1260794230.ps tmp/1cj6a1260794230.png",intern=TRUE)) character(0) > try(system("convert tmp/2vo1x1260794230.ps tmp/2vo1x1260794230.png",intern=TRUE)) character(0) > try(system("convert tmp/3xj8d1260794230.ps tmp/3xj8d1260794230.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 0.912 0.472 1.667