We've seen in the preceding section that temperature is stationary and CO2 requires a first order difference. Another simple way to show this is with the forecast package's ndiffs() function. It provides an output that spells out the minimum number of differences needed to make the data stationary. In the function, you can specify which test out of the three available ones you would like to use: Kwiatkowski, Philips, Schmidt & Shin (KPSS), Augmented Dickey Fuller (ADF), or Philips-Peron (PP). I will use ADF in the following code, which has a null hypothesis that the data is not stationary:
> ndiffs(climate[, 1], test = "adf")
[1] 1
> ndiffs(climate[, 2], test = "adf")
[1] 1
We see that both require a first order difference to become stationary. To get started, we will create the difference. Then, we will complete the traditional approach, where both series are stationary. Let's also load our packages for this exercise:
> library(vars)
> library(aod)
> climateDiff <- diff(climate)
> climateDiff <- window(climateDiff, start = 1946)
> head(climateDiff)
CO2 Temp
[1,] 78 -0.099
[2,] 154 0.034
[3,] 77 0.001
[4,] -50 -0.035
[5,] 211 -0.100
[6,] 137 0.121
It is now a matter of determining the optimal lag structure based on the information criteria using vector autoregression. This is done with the VARselect function in the vars package. You only need to specify the data and number of lags in the model using lag.max = x in the function. Let's use a maximum of 12 lags:
> lag.select <- VARselect(climateDiff, lag.max = 12)
> lag.select$selection
AIC(n) HQ(n) SC(n) FPE(n)
5 1 1 5
We called the information criteria using lag$selection. Four different criteria are provided including AIC, Hannan-Quinn Criterion (HQ), Schwarz-Bayes Criterion (SC), and FPE. Note that AIC and SC are covered in Chapter 2, Linear Regression - The Blocking and Tackling of Machine Learning, so I will not go over the criterion formulas or differences here. If you want to see the actual results for each lag, you can use lag$criteria. We can see that AIC and FPE have selected lag 5 and HQ and SC lag 1 as the optimal structure to a VAR model. It seems to make sense that the 5-year lag is the one to use. We will create that model using the var() function. I'll let you try it with lag 1:
> fit1 <- VAR(climateDiff, p = 5)
The summary results are quite lengthy as it builds two separate models and would take up probably two whole pages. What I provide is the abbreviated output showing the results with temperature as the prediction:
> summary(fit1)
Residual standard error: 0.1006 on 52 degrees of freedom
Multiple R-Squared: 0.4509, Adjusted R-squared: 0.3453
F-statistic: 4.27 on 10 and 52 DF, p-value: 0.0002326
The model was significant with a resulting adjusted R-square of 0.35.
As we did in the previous section, we should check for serial correlation. Here, the VAR package provides the serial.test() function for multivariate autocorrelation. It offers several different tests, but let's focus on Portmanteau Test, and also please note that the DW test is for univariate series only. The null hypothesis is that autocorrelations are zero and the alternate is that they are not zero:
> serial.test(fit1, type = "PT.asymptotic")
Portmanteau Test (asymptotic)
data: Residuals of VAR object fit1
Chi-squared = 35.912, df = 44, p-value = 0.8021
With p-value at 0.3481, we do not have evidence to reject the null and can say that the residuals are not autocorrelated. What does the test say with 1 lag?
To do the Granger causality tests in R, you can use either the lmtest package and the Grangertest() function or the causality() function in the vars package. I'll demonstrate the technique using causality(). It is very easy as you just need to create two objects, one for x causing y and one for y causing x, utilizing the fit1 object previously created:
> x2y <- causality(fit1, cause = "CO2")
> y2x <- causality(fit1, cause = "Temp")
It is now just a simple matter to call the Granger test results:
> x2y$Granger
Granger causality H0: CO2_diff do not Granger-cause
climate2.temp
data: VAR object fit1
F-Test = 2.2069, df1 = 5, df2 = 104, p-value = 0.05908
> y2x$Granger
Granger causality H0: climate2.temp do not Granger-cause
CO2_diff
data: VAR object fit1
F-Test = 0.66783, df1 = 5, df2 = 104, p-value = 0.6487
The p-value value for CO2 differences of Granger causing temperature is 0.05908 and not significant in the other direction. So what does all this mean? The first thing we can say is that Y does not cause X. As for X causing Y, we cannot reject the null at the 0.05 significance level and therefore conclude that X does not Granger cause Y. However, is that the relevant conclusion here? Remember, the p-value evaluates how likely the effect is if the null hypothesis is true. Also, remember that the test was never designed to be some binary Yay or Nay. If this was a controlled experiment, it would be unlikely for us to hesitate to say we had insufficient evidence to reject the null, like the Food and Drug Administration (FDS) would do for a phase 3 clinical trial. Since this study is based on observational data, I believe we can say that it is highly probable that "CO2 emissions Granger cause Surface Temperature Anomalies". But, there is a lot of room for criticism on that conclusion. I mentioned upfront the controversy around the quality of the data. The thing that still concerns me is what year to start the analysis. I chose 1945 because it looked about right; you could say I applied proc eyeball in SAS terminology. What year is chosen has a dramatic impact on the analysis, changing the lag structure and also leading to insignificant p-values.
However, we still need to model the original CO2 levels using the alternative Granger causality technique. The process to find the correct number of lags is the same as before, except we do not need to make the data stationary:
> climateLevels <- window(climate, start = 1946)
> level.select <- VARselect(climateLevels, lag.max = 12)
> level.select$selection
AIC(n) HQ(n) SC(n) FPE(n)
10 1 1 6
Let's try the lag 6 structure and see whether we can achieve significance, remembering to add one extra lag to account for the integrated series. A discussion on the technique and why it needs to be done is available at http://davegiles.blogspot.de/2011/04/testing-for-granger-causality.html:
fit2 <- VAR(climateLevels, p = 7)
> serial.test(fit2, type = "PT.asymptotic")
Portmanteau Test (asymptotic)
data: Residuals of VAR object fit2
Chi-squared = 35.161, df = 36, p-value = 0.5083
Now, to determine Granger causality for X causing Y, you conduct a Wald test, where the coefficients of X and only X are 0 in the equation to predict Y, remembering to not include the extra coefficients that account for integration in the test.
The Wald test in R is available in the aod package we've already loaded. We need to specify the coefficients of the full model, its variance-covariance matrix, and the coefficients of the causative variable.
First, let's see how CO2 does Granger causing temperature:
> CO2terms <- seq(1, 11, 2)
> Tempterms <- seq(2, 12, 2)
We are now ready to run the wald test, described in the following code:
> wald.test(b = coef(fit2$varresult$Temp),
Sigma = vcov(fit2$varresult$Temp),
Terms = c(CO2terms))
Wald test:
----------
Chi-squared test:
X2 = 11.5, df = 6, P(> X2) = 0.074
How about that? We are close to the magical 0.05 p-value. Let's test the other direction causality with the following code:
> wald.test(b = coef(fit2$varresult$CO2),
Sigma = vcov(fit2$varresult$CO2),
Terms = c(Tempterms))
Wald test:
----------
Chi-squared test:
X2 = 3.9, df = 6, P(> X2) = 0.69
The last thing to show here is how to use a vector autoregression to produce a forecast. A predict function is available, so let's autoplot() it for a 25-year period and see what happens:
> autoplot(predict(fit2, n.ahead = 25, ci = 0.95))
It seems dark days lie ahead, perhaps, to coin a phrase, "Winter is Coming" from the popular TV series Game of Thrones. Fine by me as my investment and savings plan for a long time has consisted of canned goods and ammunition. What else am I supposed to do, ride a horse to work? I'll do that the day Al Gore does. In the meantime, I am going to work on my suntan.
If nothing else, I hope it has stimulated your thinking on how to apply the technique to your own real-world problems or maybe even to examine the climate change data in more detail. There should be a high bar when it comes to demonstrating causality, and Granger causality is a great tool for assisting in that endeavor.