When it comes to causality tests, the typical Granger-causality test can be problematic. Testing for Granger-causality using F-statistics when one or both time series are non-stationary can lead to spurious causality (He & Maekawa, 1999).
Professor Giles gives an excellent example of how the Toda-Yamamoto (TY) method can be implemented. More formal explanations can be found in the original TY (1995) paper.
In this post, I will show how Professor Giles' example can be implemented in R. The procedure is based on the following steps:
- Test for integration (structural breaks need to be taken into account). Determine the maximum order of integration (m). If none of the series is integrated, the usual Granger-causality test can be done.
- Set up a VAR model in the levels (do not difference the data).
- Determine the lag length p. The VAR model is thus VAR(p).
- Carry out tests for misspecification, especially for residual serial correlation.
- Add the maximum order of integration to the number of lags. This is the augmented VAR model, VAR(p+m).
- Carry out a Wald test for the first p variables only, with p degrees of freedom.
You may want to run a test of cointegration. If the series are cointegrated, there must be causality. However, Toda and Yamamoto (1995) noted that one advantage of the TY method is that you don't have to test for cointegration — a pretest bias can therefore be avoided.
The example is about causalities between the prices of Robusta and Arabica coffee. You can download the data (CSV) to follow along. The script below is annotated; let me know if I can clarify anything or if there is room for improvement.
library(fUnitRoots)
library(urca)
library(vars)
library(aod)
library(zoo)
library(tseries)
#Load data
cof <- read.csv("https://christophpfeiffer.org/data/coffee_data.csv", header=T,sep=";")
names(cof)
#Adjust Date format
cof["Date"]<-paste(sub("M","-",cof$Date),"-01",sep="")
#Visualize
plot(as.Date(cof$Date),cof$Arabica,type="l",col="black",lwd=2)
lines(as.Date(cof$Date),cof$Robusta,col="blue",lty=2,lwd=1)
legend("topleft",c("Arabica","Robusta"),col=c("black","blue"),lty=c(1,2),lwd=c(2,1),bty="n")
#Restrict the sample if there is an early structural break (adjust to your data)
cof1<-cof[193:615,]
#Visualize
plot(as.Date(cof1$Date),cof1$Arabica,type="l",col="black",lwd=2,ylim=range(cof1$Robusta))
lines(as.Date(cof1$Date),cof1$Robusta,col="blue",lty=2,lwd=1)
legend("topright",c("Arabica","Robusta"),col=c("black","blue"),lty=c(1,2),lwd=c(2,1),bty="n")
#Test for unit roots
adf.test(cof$Arabica)
adf.test(cof$Robusta)
kpss.test(cof$Arabica)
kpss.test(cof$Robusta)
adf.test(diff(cof$Arabica,1))
adf.test(diff(cof$Robusta,1))
kpss.test(diff(cof$Arabica,1))
kpss.test(diff(cof$Robusta,1))
# Since first order differencing eliminates the unit root, the maximum order of integration
# is concluded to be I(1).
#Set up VAR-Model
#select lag order // either 2 or 6
VARselect(cof1[,2:3],lag=20,type="both")
#VAR Model, lag=2
V.2<-VAR(cof1[,2:3],p=2,type="both")
serial.test(V.2)
#VAR-Model, lag=6
V.6<-VAR(cof1[,2:3],p=6,type="both")
serial.test(V.6)
#Stability analysis (thanks to Erdogan Cevher)
1/roots(V.6)[[1]] # ">1"
1/roots(V.6)[[2]] # ">1"
# Model with p=6 is less likely to be serially correlated. Thus model with p=6 is selected.
# Wald-test for the first 6 lags
# VAR model is seperately set up as a linear model; makes the wald test easier
#lag variables
arab<-zoo(cof1["Arabica"])
robu<-zoo(cof1["Robusta"])
arab.l<-lag(arab,-(0:7),na.pad=T)
robu.l<-lag(robu,-(0:7),na.pad=T)
lm1<-lm(arab~arab.l[,2:8]+robu.l[,2:8]+index(arab))
lm2<-lm(robu~arab.l[,2:8]+robu.l[,2:8]+index(arab))
#Wald-test (H0: Robusta does not Granger-cause Arabica)
vcov(lm1)
wald.test(b=coef(lm1), Sigma=vcov(lm1), Terms= c(9:14),df=6)
# Could not be rejected (X2=8.6; p=0.2)
#Wald.test (H0: Arabica does not Granger-cause Robusta)
vcov(lm2)
wald.test(b=coef(lm2), Sigma=vcov(lm2), Terms= c(2:7),df=6)
# Could be rejected at 10% (X2=12.3; p=0.056)
# It seems that Arabica Granger-causes Robusta prices, but not the other way around.
Let me know if you have any suggestions.
— C
References
He, Z.; Maekawa, K. (1999). On spurious Granger causality. Economics Letters, 73(3), 307–313.
Toda, H. Y.; Yamamoto, T. (1995). Statistical inference in vector autoregressions with possibly integrated processes. Journal of Econometrics, 66, 225–250.
Discussion
A few questions and answers from the comments:
Dave Giles — whose example this builds on:
This is just great! Thanks for sharing this.
Cindy M. — Why is the index(arab) term in the regression?
index(arab)captures the trend.index(robu)would work just as well and yields the same result.
"london" — I ran this on a small data set (22 observations, 9 variables) and get NAs — am I overfitting?
With so few observations, statistical analysis is difficult. Any chance to obtain more data? And have you considered bootstrapping?
Mike — wald.test also reports an F-statistic, but you interpret the χ² result. How should the F output be read here?
Toda and Yamamoto showed that if you add an additional lag to a correctly specified VAR model for which at least one series is integrated, the parameters asymptotically follow a chi-squared distribution. Under those conditions, looking at the F-statistic for the augmented model would simply yield a meaningless answer.
Rick — Why Terms = c(9:14) and c(2:7), rather than c(8:13) / c(1:6), given the last lag is exogenous just for the asymptotics?
For the first test we assess whether the Robusta coefficients are significant for the price of Arabica; for the second, whether the Arabica coefficients are significant for Robusta. The Wald test targets the relevant block of lag coefficients in each model accordingly.