Granger causality ordering variable

Question

I am trying to compute Granger causality from a VAR using vars package. I know that order of variable is important in a VAR to compute IRF, but here I have different result for Granger causality. Here is a reproducible example:

test <-
structure(c(324L, 444L, 456L, 451L, 484L, 535L, 625L, 622L, 532L, 
623L, 619L, 642L, 550L, 643L, 629L, 709L, 781L, 792L, 858L, 979L, 
974L, 1039L, 1062L, 1167L, 1345L, 1600L, 1789L, 1912L, 1965L, 
2086L, 2350L, 2507L, 2186L, 2224L, 2235L, 2435L, 1923L, 2035L, 
2253L, 2365L, 2214L, 2554L, 2714L, 2773L, 2777L, 2897L, 3060L, 
3053L, 2426L, 2714L, 2953L, 2925L, 3233L, 3494L, 3724L, 4137L, 
4147L, 4275L, 4726L, 4910L, 114L, 131L, 137L, 136L, 157L, 174L, 
207L, 165L, 168L, 189L, 176L, 187L, 138L, 153L, 176L, 183L, 196L, 
229L, 263L, 281L, 263L, 310L, 310L, 313L, 283L, 291L, 324L, 367L, 
410L, 423L, 470L, 452L, 465L, 529L, 559L, 602L, 657L, 729L, 721L, 
747L, 781L, 815L, 861L, 863L, 817L, 839L, 875L, 855L, 825L, 947L, 
1038L, 1086L, 1190L, 1235L, 1347L, 1435L, 1495L, 1482L, 1645L, 
1769L, 61L, 84L, 106L, 100L, 123L, 112L, 131L, 151L, 105L, 121L, 
127L, 167L, 166L, 164L, 153L, 212L, 228L, 236L, 264L, 316L, 223L, 
245L, 267L, 301L, 288L, 276L, 284L, 351L, 378L, 395L, 460L, 530L, 
479L, 483L, 472L, 528L, 399L, 456L, 477L, 480L, 500L, 526L, 550L, 
636L, 513L, 567L, 607L, 595L, 565L, 593L, 632L, 614L, 719L, 747L, 
874L, 1011L, 905L, 857L, 947L, 956L), .Dim = c(60L, 3L), .Dimnames = list(
    NULL, c("H", "M", "S")), index = structure(c(1996, 1996.08333333333, 
1996.16666666667, 1996.25, 1996.33333333333, 1996.41666666667, 
1996.5, 1996.58333333333, 1996.66666666667, 1996.75, 1996.83333333333, 
1996.91666666667, 1997, 1997.08333333333, 1997.16666666667, 1997.25, 
1997.33333333333, 1997.41666666667, 1997.5, 1997.58333333333, 
1997.66666666667, 1997.75, 1997.83333333333, 1997.91666666667, 
1998, 1998.08333333333, 1998.16666666667, 1998.25, 1998.33333333333, 
1998.41666666667, 1998.5, 1998.58333333333, 1998.66666666667, 
1998.75, 1998.83333333333, 1998.91666666667, 1999, 1999.08333333333, 
1999.16666666667, 1999.25, 1999.33333333333, 1999.41666666667, 
1999.5, 1999.58333333333, 1999.66666666667, 1999.75, 1999.83333333333, 
1999.91666666667, 2000, 2000.08333333333, 2000.16666666667, 2000.25, 
2000.33333333333, 2000.41666666667, 2000.5, 2000.58333333333, 
2000.66666666667, 2000.75, 2000.83333333333, 2000.91666666667
), class = "yearmon"), class = "zoo")

I follow Toda-Yamamoto procedure with the following code:

VARselect(test, lag.max=14)
var<-VAR(test, p=2, season=12)
serial.test(var)
var<-VAR(test, p=3, season=12)

aTester<-c(1)
#test causality of Industry on all others variables
wald.test(b=coef(var$varresult[[1]]), Sigma=vcov(var$varresult[[1]]), Terms=aTester)
#test causality of Finance on all others variables
wald.test(b=coef(var$varresult[[1]]), Sigma=vcov(var$varresult[[1]]), Terms=aTester+1)
#test causality of Consulting on all others variables
wald.test(b=coef(var$varresult[[1]]), Sigma=vcov(var$varresult[[1]]), Terms=aTester+2)

My result are different if a change the order of variable in VAR and i really don't understand why.

EDIT : Apparently I made a huge mistake and i had to change something in my code :

 wald.test(b=coef(var$varresult[[1]]), Sigma=vcov(var$varresult[[1]]), Terms=aTester)
    #test causality of Finance on all others variables
    wald.test(b=coef(var$varresult[[**2**]]), Sigma=vcov(var$varresult[[**2**]]), Terms=aTester+1)
    #test causality of Consulting on all others variables
    wald.test(b=coef(var$varresult[[**3**]]), Sigma=vcov(var$varresult[[**3**]]), Terms=aTester+2)

Answer 1

It looks to me like you specify your Terms= argument wrong. Here follows a brief explanation with examples by aod::wald.test() and the vars::causality():

If you want to test that "M,S do not Granger-cause H" then

wald.test(b=coef(var$varresult[[1]]), Sigma=vcov(var$varresult[[1]]), Terms=c(2,3,5,6,8,9))

which you could also obtain directly from the vars-package by:

causality(var,c("M","S"))

If you want to test that "S do not Granger-cause H,M" then you could do something along:

wald.test(
    b=Reduce(rbind,coef(var))[c(1:9,22:30),1], 
    Sigma=vcov(var)[c(2:10,23:31),c(2:10,23:31)], 
    Terms=c(3,6,9,12,15,18)
)

where it requieres a little more work because of a difference in variable ordering. In vars you could directly specify:

causality(var,"S")

At last if you want bivariate Granger causality tests, then you could use the function in MSBVAR:

library(MSBVAR)
granger.test(test, p=3) 

Hope this helps.