0
$\begingroup$

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.

$\endgroup$

migrated from stackoverflow.com Jul 1 '15 at 8:46

This question came from our site for professional and enthusiast programmers.

  • $\begingroup$ 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. $\endgroup$ – MrFlick Jul 1 '15 at 0:05
  • $\begingroup$ It can be easily extended as you can see my answer. $\endgroup$ – user227710 Jul 1 '15 at 1:51
2
$\begingroup$

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);
}
$\endgroup$
1
$\begingroup$
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))
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.