AnInconvenientStatisticCodeExample.r~

> require(forecast)
> require(vars)
> var.data = read.csv(file.choose())
> head(var.data)

> #put data into a time series
> carbon.ts = ts(CO2, frequency=1, start=c(1850), end=c(2010))
> temp.ts = ts(Temp, frequency=1, start=c(1850), end=c(2010))
  #subset the data from 1900 until 2010
> surfacetemp = window(temp.ts, start=c(1900), end=c(2010))
> co2 = window(carbon.ts, start=c(1900), end=c(2010))
> climate.ts = cbind(co2, surfacetemp)
> plot(climate.ts)
  
> #determine stationarity and number of lags to achieve stationarity
> ndiffs(co2, alpha = 0.05, test = c("adf"))
> ndiffs(surfacetemp, alpha = 0.05, test = c("adf"))

> #difference to achieve stationarity
> d.co2 = diff(co2)
> d.temp = diff(surfacetemp)
  
> #again, we need a mts class dataframe
> climate2.ts = cbind(d.co2, d.temp)
> plot(climate2.ts)

> #determine the optimal number of lags for vector autoregression
> VARselect(climate2.ts, lag.max=10) $selection

> #vector autoregression with lag1
> var = VAR(climate2.ts, p=1)

> serial.test(var, lags.pt=10, type="PT.asymptotic")

#The null hypothesis is no serial correlation, so we can reject it with extreme prejudice.on to var3
> var3 = VAR(climate2.ts, p=3)
> serial.test(var3, lags.pt=10, type="PT.asymptotic")

> summary(var3, equation="d.temp")

> #does co2 granger cause temperature
> grangertest(d.temp ~ d.co2, order=3)

> #Clearly the model is not significant, so we can say that carbon emissions do not granger-cause surface temperatures.
  
> #does temperature granger cause co2
> grangertest(d.co2 ~ d.temp, order =3)

> #try again using lag 7
> grangertest(d.temp ~ d.co2, order=7)

> predict(var3, n.ahead=6, ci=0.95)

> fcst = forecast(var3)
> plot(fcst)