# How do I perform a Wald Test with multivariate Granger Causality Analysis

I am doing a Granger Causality Analysis for three economic variables (GDP, CO2 emissions and Total Energy Consumption) of Puerto Rico. I am using a Toda-Yamamoto Procedure implemented in R R. I am basing my code on example given by Christoph Pfieffer in this blog post. However, the example that is given in the blog post has two variables, which makes choosing the coefficient for the "Terms" argument fairly easy. I have to do a Wald Test to compare the causality potential of three variables, for example: "GDP Granger-cause Energy Consumption", "Energy Consumption Granger-cause CO2 emissions" and does "GDP Granger cause CO2 emissions"? My question is how do I perform a Wald with three variables and how I choose the coefficient for the "Terms" argument.

• This doesn't sound like a very specific programming question. If you need help with statistical methods, you should post your question to Cross Validated instead where such matters on on-topic. Otherwise create a reproducible example with sample code and data to make it clear what the programming question is. Jul 1 '15 at 0:05
• It can be easily extended as you can see my answer. Jul 1 '15 at 1:51

library(urca)
#library(fUnitRoots)
library(urca)
library(vars)
library(aod)
library(zoo)
library(tseries)
data(denmark)

lapply(2:5,function(i)adf.test(denmark[[i]])) # I haven't gone in details over these
lapply(2:5,function(i)adf.test(diff(denmark[[i]]))) # I haven't gone in details over these; you need to check but I assume that they are of I(1)
VARselect(denmark[,2:4],lag=20,type="both") #something like lag 11 is prefered
V.11<-VAR(denmark[,2:4],p=11,type="both")
serial.test(V.11) # I haven't gone in details over these
1/roots(V.11)[[1]] # I haven't gone in details over these
plot(stability(V.11)) # I haven't gone in details over these
V.12<-VAR(denmark[,2:4],p=12,type="both") # I haven't gone in details over these (add 1 lag as in the linked website)
V.12$varresult # I haven't gone in details over these summary(V.12) # I haven't gone in details over these #Wald-test (H0: LRY does not Granger-cause LRM) wald.test(b=coef(V.12$varresult[[1]]), Sigma=vcov(V.12$varresult[[1]]), Terms=seq(2,33,3)) #Wald-test (H0: LPY does not Granger-cause LRM) wald.test(b=coef(V.12$varresult[[1]]), Sigma=vcov(V.12$varresult[[1]]), Terms=seq(3,33,3)) #Wald-test (H0: LRM does not Granger-cause LRY) wald.test(b=coef(V.12$varresult[[2]]), Sigma=vcov(V.12$varresult[[2]]), Terms=seq(1,33,3)) #Wald-test (H0: LPY does not Granger-cause LRY) wald.test(b=coef(V.12$varresult[[2]]), Sigma=vcov(V.12$varresult[[2]]), Terms=seq(3,33,3)) #Wald-test (H0: LRM does not Granger-cause LPY) wald.test(b=coef(V.12$varresult[[3]]), Sigma=vcov(V.12$varresult[[3]]), Terms=seq(1,33,3)) #Wald-test (H0: LRY does not Granger-cause LPY) wald.test(b=coef(V.12$varresult[[3]]), Sigma=vcov(V.12$varresult[[3]]), Terms=seq(2,33,3))  I ran into the same problem and coded a specific function that takes a var object with two or possible more variables, adds the extra lag and then conducts the Wald test. It returns a data frame with the result. Note: This function only works for I(1) series, because it adds a single lag. For higher orders of integration, more lags are required. toda.yamamoto <- function(var) { # add the magic lag to the existing VAR model ty.df <- eval(var$call$y); ty.varnames <- colnames(ty.df); ty.lags <- var$p + 1;
ty.augmented_var <- VAR(ty.df, ty.lags, type=var$type); ty.results <- data.frame(predictor = character(0), causes = character(0), chisq = numeric(0), p = numeric(0)); for (current_variable in ty.varnames) { # construct the restriction matrix: to test if *current_variable* causes any of the others, # we test if the lagged values of current variable (ignoring the magic lag) are jointly insignificant ty.restrictions <- as.matrix(Bcoef(ty.augmented_var))*0+1; ty.coefres <- head(grep(current_variable, colnames(ty.restrictions), value=T), -1); ty.restrictions[which(rownames(ty.restrictions) != current_variable), ty.coefres] <- 0; # estimate restricted var ty.restricted_var <- restrict(ty.augmented_var, 'manual', resmat=ty.restrictions); for (k in 1:length(ty.varnames)) { if (ty.varnames[k] != current_variable) { my.wald <- waldtest(ty.augmented_var$varresult[[k]], ty.restricted_var$varresult[[k]], test='Chisq'); ty.results <- rbind(ty.results, data.frame( predictor = current_variable, causes = ty.varnames[k], chisq = as.numeric(my.wald$Chisq[2]),
p = my.wald\$Pr(>Chisq)[2])
);
}
}
}
return(ty.results);
}