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

Question

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.

Answer 1

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))

Answer 2

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);
}