Toda-Yamamoto implementation in 'R'

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 TY method can be implemented. 


More formal explanations can be found in the original TY (1995) paper or for example here


In this post, I will show how Professor Giles' example can be implemented in R. 


The procedure is based on the following steps: 


1. Test for integration (structural breaks need to be taken into account). Determine max order of integration (m). If none of the series in integrated, the usual Granger-causality test can be done. 


2. Set up a VAR-model in the levels (do not difference the data). 


3. Determine lag length. Let the lag length be p. The VAR model is thus VAR(p).


4. Carry out tests for misspecification, especially for residual serial correlation. 


5. Add the maximum order of integration to the number of lags. This is the augmented VAR-model, VAR(p+m). 


6. Carry out a Wald test for the first p variables only with p degrees of freedom. 


You may want to do a test of cointegration. If series are cointegrated, there must be a causality. However, Toda and Yamamoto (1995) noted that one advantage of the TY-method is that you don't have to test for cointegration and, therefore, a pretest bias can be avoided.




The example is about causalities between prices in Robusta and Arabica coffee. The excel-file can be downloaded here. But in order to be loaded into R, the data should be put in the csv. format. The csv. file is available here


Update: If you want examine the data interactively, have a look here.


The script below tests for causality between these two time series. The script is annotated, but let me know if I can clarify anything or if there is room for improvement. 




#Load data
cof <- read.csv("", header=T,sep=";")

#Adjust Date format


#Possible structural break in 1970s. Therefore only values from 1976:01 onwards are regarded


#Test for unit roots


# 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

#VAR Model, lag=2

#VAR-Model, lag=6

# Model with p=6 is less likely to be serially correlated. Thus model with p=6 is selected. 

# Model with additional lag is set up. 

# Wald-test for the first 6 lags 
# VAR model is seperately set up as a linear model; makes the wald test easier

#lag variables



#Wald-test (H0: Robusta does not Granger-cause Arabica)
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)
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. 

You can download the R-code as well as the csv. file in "Files".


Let me know if you have any suggestions. 


--- C





He, Z.; Maekawa, K. (1999). On spurious Granger causality. Economic 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. 



Write a comment

Comments: 33

  • #1

    Dave Giles (Wednesday, 07 November 2012 23:42)

    Christoph. This is just great! Thanks for sharing this.

  • #2

    Claudio D. Shikida (Thursday, 08 November 2012 16:10)

    Just change this for the tests.

    It's not adf.test, but ur.df. For example:


    Great job! Thanks!

  • #3

    Christoph (Thursday, 08 November 2012 17:38)

    Claudio, thanks for pointing out.

    adf.test(), kpss.test() work es well, but we need the "tseries" package loaded. Should be fine now.

  • #4

    Cindy M. (Thursday, 06 December 2012 17:53)

    Thanks for sharing this! Two questions for you.
    In lm1<-lm(arab~arab.l[,2:8]+robu.l[,2:8]+index(arab)) , can you please explain why you need to have the index(arab) term in the regression?

    Also in:
    >#Wald.test (H0: Arabica does not Granger-cause Robusta)
    >wald.test(b=coef(lm1), Sigma=vcov(lm1), Terms= c(2:7),df=6)

    Should this wald.test be using lm2 instead of lm1?

  • JimdoBusiness

    christophpfeiffer (Thursday, 06 December 2012 18:08)

    Hi Cindy,

    thanks for the catch, it is indeed lm2 and not lm1.

    index(arab) caputres the trend index(robu) would of course also work and yields the same result.

    -- Christoph

  • #6

    Cindy M. (Thursday, 06 December 2012 18:23)

    Thanks for the quick response!

  • #7

    london (Thursday, 21 March 2013 20:12)

    Dear Christoph,
    I used the code on a data set which has 22 observations and 9 variables all variables entering into the vAR model. I determined p=k+dmax=3. this gives me 27 coefficient estimates. I have NAs in the results. Just wondering if I am doing something wrong or overfitting?

  • JimdoBusiness

    christophpfeiffer (Thursday, 21 March 2013 21:52)


    it seems that you have relatively few observations which makes statistical analysis difficult. Any chance obtain more data? Have you considered bootstrapping?

    -- Christoph

  • #9

    london (Thursday, 21 March 2013 23:31)

    Dear Christoph,

    I am trying to run a bootstrap regression, however, the regression is producing NAs in place of coefficients for the last 2 coefficients after every bootstrap.

    would appreciate any suggestion!

  • JimdoBusiness

    christophpfeiffer (Friday, 22 March 2013 17:32)

    If you like, you can send me your R-code and the data and I'll have a look.

  • #11

    london (Friday, 22 March 2013 23:54)

    Dear Christoph,
    Many thanks, I will email you my code and data set.

  • #12

    Aviral Kumar Tiwari (Sunday, 19 May 2013 10:08)

    I tried to use your R-codes but in my computer it it shows "Error: could not find function "wald.test""
    when i used
    > wald.test(b=coef(lm1), Sigma=vcov(lm1), Terms= c(9:14),df=6)

    Please help me ..

  • JimdoBusiness

    christophpfeiffer (Sunday, 19 May 2013 12:14)

    Hi Aviral,
    make sure you have all packages installed. So you shouldn't get any error messages for the first 6 lines. You can install packages with the function


    where ... is the name of the package. Specifically, for the function wald.test() you need the package "aod" to be installed. For this you need to execute the following line:



    Hope this helps


  • #14

    Aviral Kumar Tiwari (Sunday, 19 May 2013 13:55)

    thankssssss it worked out ...

  • #15

    Mike (Friday, 07 June 2013 00:49)


    The output for wald.test also gives the F-statistic, but your comments show the test is determined by the X2 results. What is the interpretation of the F output results in this case? Thanks!

  • JimdoBusiness

    christophpfeiffer (Sunday, 16 June 2013 12:05)

    sorry for the late answer, I just got back from vacation.
    Toda and Yamamoto have shown that if you add an additional lag to a correctly specified VAR-model for which at least one time-series is integrated, the parameters asymptotically follow a chi-squared distribution. Assuming these two conditions are met, looking at the F-statistic for the augmented model would simply yield a meaningless answer.

  • #17

    Rick (Sunday, 14 July 2013 02:35)

    In your wald.test, why are you using Terms c(9:14) and c(2:7)? I would have thought it would be c(8:13) and (1:6) because the last lag is supposed to be exogenous and just there for the asymptotics?

  • #18

    Christoph (Monday, 15 July 2013 09:36)

    Hi Rick,

    for the first test we need to assess if the coefficients of Robusta are significant for the price of Arabica. For the second test we are testing if the coefficients of Arabica are significant for the price of Robusta. These are the models:

    Arabica = Intercept+A*X_1+B*X_2 + Trend # lm1
    No. of terms: (1) (7) (7) (1) = 16

    Robusta = Intercept+A*X_1+B*X_2 + Trend # lm2
    No. of terms: (1) (7) (7) (1) = 16

    With A being the coefficients for the values of Arabica, B coefficients for Robusta, X_1 and X_2 the respective lagged values. For each matrix of lagged values we have 7 terms but only want to test the first 6. Then, for the first Wald test we need to look at terms 9 to 14 and for the second model at terms 2 to 7. You can look at the coefficients of each model with coef(lm1) or coef(lm2).

  • #19

    Martin (Monday, 28 October 2013 13:21)

    I'm having trouble interpreting the serial.tests. What's the significance of the degrees of freedom and the p-values? (The p-value for serial.test(V.6) is .5.)

  • #20

    Bob (Wednesday, 11 December 2013 19:44)

    Any way to do this without the packages?

  • JimdoBusiness

    christophpfeiffer (Saturday, 14 December 2013 15:14)

    of course, but why would you do it without the packages? If you want to look deeper, you can just look at the functions provided by the packages.

  • JimdoBusiness

    christophpfeiffer (Saturday, 14 December 2013 15:20)

    sorry for the late reply. The serial test used here is a portmanteau test with the null hypothesis that residual errors are not serially correlated. A higher p-value indicates a favorable model which is why a model with 6 vs. 2 lags is chosen.

  • #23

    Yolande (Friday, 14 March 2014 20:18)

    I was wondering how the method and the R code would change if one of the series has unit root and the other does not.

  • #24

    Christoph (Friday, 14 March 2014 21:36)


    good question. The first step says:

    1. Test for integration (structural breaks need to be taken into account). Determine max order of integration (m). If none of the series in integrated, the usual Granger-causality test can be done.

    So, if one series has an order of integration of 1 and the other is not integrated (order of integration is 0), the maximum order of integration would be m=1.

    Jumping to step 5:

    5. Add the maximum order of integration to the number of lags. This is the augmented VAR-model, VAR(p+m).

    This means we would still have to use the augmented model even if one series is not integrated at all. It's the maximum order of integration that counts. The only case in which you can use the standard Granger test for non-causality is when neither of the series are integrated.

  • #25

    Cory (Monday, 17 March 2014 18:00)


    I was wondering how I could change the code in my blog to follow the procedure you outline above. I don't understand the idea of including an additional lag, but not testing for it.

  • #26

    Yolande (Monday, 17 March 2014 18:00)

    Thank you for the prompt reply.

  • #27

    Yolande (Tuesday, 18 March 2014 15:06)

    I do have an additional question: Where in the code and algorithm the model v.7 was used after it is set

  • #28

    Pablo (Thursday, 03 April 2014 22:38)

    I was wondering if you've ever used the function "causality" from the package "vars"? Are the results from this function compatible with the ones from your code?

    Thanks a lot

  • #29

    Christoph (Saturday, 05 April 2014 16:17)


    there are some issues with differencing the data first for integrated time series and then using the standard Granger causality test (e.g. Enders, W. Applied Econometric Times Series ).

    You can think of the idea to add an additional lag and then not using it for the test as a trick to get rid of the problems involved in testing integrated time series for causality.

  • #30

    Christoph (Saturday, 05 April 2014 16:18)


    you are very welcome. :)

    Good point, model v.7 is actually not necessary for the following tests.

  • #31

    Christoph (Saturday, 05 April 2014 16:20)

    Hi Pablo,

    as far as I know, the causality function from 'vars' only uses the standard Granger causality test and would not work in this context.

  • #32

    Yolande (Monday, 07 April 2014 21:33)

    Hi Christoph,
    This is a general question.
    Why would I got -Infinity in my varselect. The lag length is 3 based on all criteria.
    > VARselect(M,lag.max=20,type="both")
    AIC(n) HQ(n) SC(n) FPE(n)
    3 3 3 3

    1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
    AIC(n) 5.517709 4.919171 -Inf -Inf -Inf -Inf -Inf -Inf -Inf -Inf -Inf -Inf -Inf -Inf -Inf -Inf -Inf
    HQ(n) 5.139389 4.351692 -Inf -Inf -Inf -Inf -Inf -Inf -Inf -Inf -Inf -Inf -Inf -Inf -Inf -Inf -Inf
    SC(n) 5.693020 5.182138 -Inf -Inf -Inf -Inf -Inf -Inf -Inf -Inf -Inf -Inf -Inf -Inf -Inf -Inf -Inf
    FPE(n) 284.562815 237.788270 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
    18 19 20
    AIC(n) -Inf -Inf -Inf
    HQ(n) -Inf -Inf -Inf
    SC(n) -Inf -Inf -Inf
    FPE(n) 0 0 0

  • JimdoBusiness

    christophpfeiffer (Friday, 18 April 2014 19:31)

    Hi Yolande,
    this probably happens because the number of lags available in your dataset is too low to calculate the criteria for lags higher than two. If this does not explain it, you can send me the dataset and I will have a look.

  • loading
Profil von Christoph Pfeiffer auf LinkedIn anzeigen