[This article was first published on Fear and Loathing in Data Science, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.
When I first learned about Granger-causality this past February, I was bemused and quite skeptical of the whole procedure. I felt it belonged on the scrapheap of impractical academic endeavors, preferring to possibly use an ARIMA transfer function model for the same task. However, several contemporaries threw the red challenge flag and upon further review, my initial impressions have been overturned. Not only am I fascinated by the technique, in my attempt to discover its value I have became a raving R fan. As such, my first blog entry is to provide some simple code to allow anyone to utilize this obscure econometric technique, but first some background.< o:p>
Given two sets of time series data, x and y, granger-causality is a method which attempts to determine whether one series is likely to influence change in the other. This is accomplished by taking different lags of one series and using that to model the change in the second series. We create two models which predict y, one with only past values of y (Ω), and the other with past values of y and x (π). The models are given below where k is the number of lags in the time series:< o:p>
Let Ω = yt = β0 + β1yt-1 +…+ βkyt-k + e< o:p>
And π = yt = β0 + β1yt-1 +…+ βkyt-k + α1xt-1 +…+ αkxt-k + e< o:p>
The residual sums of squared errors are then compared and a test is used to determine whether the nested model (Ω) is adequate to explain the future values of y or if the full model (π) is better. The F-test, t-test or Wald test (used in R) are calculated to test the following null and alternate hypotheses:< o:p>
H0: αi = 0 for each i of the element [1,k]< o:p>
H1: αi ≠ 0 for at least 1 i of the element [1,k]< o:p>
Essentially, we are trying to determine whether we can say that statistically x provides more information about future values of y than past values of y alone. Under this definition it is clear that we are not trying to prove actual causation, only that the two values are related by some phenomenon. Along those lines, we must also run this model in reverse to verify that that y does not provide information about future values of x. If we find that this is the case, it is likely that there is some exogenous variable, z, which needs to be controlled or could be a better candidate for granger causation. < o:p>
For a detailed explanation, one can read the original paper on the subject:< o:p>
Granger, Clive W., “Investigating Causal Relations by Econometric Models and Cross-Spectral Methods”, Econometrica, 37(1969): 424-38< o:p>
The R package “lmtest” incorporates the granger causality procedure, including a data set to answer the age old question of what came first, “the chicken or the egg”. The data was presented by Walter Thurman and Mark Fisher in the American Journal of Agricultural Economics, May 1988, titled “Chickens, Eggs, and Causality, or Which Came First?” It consists of two time series from 1930 to 1983, one of U.S. egg production and the other the estimated U.S. chicken population.< o:p>
Let’s get some code going, first loading the data from a saved .csv file.< o:p>
> chickegg <- read.csv(file.choose())< o:p>
> head(chickegg)< o:p>
Year chicken egg< o:p>
1 1930 468491 3581< o:p>
2 1931 449743 3532< o:p>
3 1932 436815 3327< o:p>
4 1933 444523 3255< o:p>
5 1934 433937 3156< o:p>
6 1935 389958 3081< o:p>
> attach(chickegg)< o:p>
> # plot the time series< o:p>
> par(mfrow=c(2,1))< o:p>
> plot.ts(chicken)< o:p>
> plot.ts(egg)< o:p>
The plots provide little information other than the data is likely not stationary. I’ve just started using the forecast package, so let’s load it and test for what will achieve stationarity.< o:p>
> library(forecast)< o:p>
> # test for unit root and number of differences required, you can also test for seasonality with nsdiffs< o:p>
> ndiffs(chicken, alpha=0.05, test=c(“kpss”)) < o:p>
[1] 1< o:p>
> ndiffs(egg, alpha=0.05, test=c(“kpss”))< o:p>
[1] 1< o:p>
> # differenced time series< o:p>
> dchick <- diff(chicken)< o:p>
> degg <- diff(egg)< o:p>
> plot.ts(dchick)< o:p>
> plot.ts(degg)< o:p>
Much better!< o:p>
That’s pretty standard stuff, but this is where the magic happens! There are several ways to find the optimal lag, which I will skip in the interest of time, but let’s say four is the magic number.< o:p>
> # do eggs granger cause chickens? < o:p>
> grangertest(dchick ~ degg, order=4)< o:p>
Granger causality test< o:p>
Model 1: dchick ~ Lags(dchick, 1:4) + Lags(degg, 1:4)< o:p>
Model 2: dchick ~ Lags(dchick, 1:4)< o:p>
Res.Df Df F Pr(>F) < o:p>
1 40 < o:p>
2 44 -4 4.1762 0.006414**< o:p>
—< o:p>
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1< o:p>
Highly significant p-value, but what about the other direction?< o:p>
> # do chickens granger cause eggs, at lag 4?< o:p>
> grangertest(degg ~ dchick, order=4)< o:p>
Granger causality test< o:p>
Model 1: degg ~ Lags(degg, 1:4) + Lags(dchick, 1:4)< o:p>
Model 2: degg ~ Lags(degg, 1:4)< o:p>
Res.Df Df F Pr(>F)< o:p>
1 40 < o:p>
2 44 -4 0.2817 0.8881< o:p>
It is not significant, so we can say the eggs Granger-Cause chickens!< o:p>
This is just the tip of the iceberg, but should be enough to strike up your curiosity and to make you dangerous. I’m working on commodity prices, bond prices and the U.S. stock market, but that is better left for another day.< o:p>
Reference:< o:p>
Thurman W.N. & Fisher M.E. (1988), Chickens, Eggs, and Causality, or Which Came First?, American Journal of Agricultural Economics, 237-238.< o:p>
To leave a comment for the author, please follow the link and comment on their blog: Fear and Loathing in Data Science.
R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.