CO2 and Global Temperature

CO2 and Global Temperature:

Trying to follow this RPubs article by Ivo Welch.

Data Set:

From the quoted paper:

One of the most prominent graphs in climate activism — the correlation between CO2 and global temperature over the last 400,000 years, based on the Vostok Ice Cores — is misleading when it is shown in order to bolster the case that CO2 increases have caused global warming. The graph intermingles two effects. Temperature can predict CO2. CO2 can predict temperature. Both effects are well known and one does not preclude the other. When presented in a CO2-temperature context, this intermixing renders the graph deceptive. In truth, the prime of the analysis is the high auto-correlation of changes in CO2. It makes it difficult to pin down the role of CO2 shocks on temperature innovations. Put into econometric-analysis terms, the x variable is poorly identified. It is itself endogenous most of the time. Perhaps most importantly, this joint 400,000 year history is not relevant today. Although the driving forces of CO2 were not well-known for most of this long period, scientists know that humanity has injected about 130 ppm of CO2 into the atmosphere, unrelated to other factors, during the last 200 years.

Endogenous: Price and quantity are related by some linear function \(P = a + b \times Q\). One depends on the other: an increase in quantity results in a decrease in price. Therefore, price is endogenous.

Exogenous: Income is exogenous. An increase in disposable income will shift up the price regardless of the quantity: the line will be translated upwards.

The dataset can be found in here and contains 4,117 observations at 100 year intervals.

options(scipen=999)
library(data.table) # Provides the function shift
dat <- read.csv("Vostok_temp_co2.csv")
dat <- dat[complete.cases(dat),]
names(dat) <- c('year', 'temp', 'co2', 'solar')
summary(dat)

##       year              temp             co2            solar      
##  Min.   :-414000   Min.   :-9.222   Min.   :182.2   Min.   :388.0  
##  1st Qu.:-311100   1st Qu.:-7.001   1st Qu.:203.5   1st Qu.:424.5  
##  Median :-208200   Median :-5.530   Median :224.9   Median :440.7  
##  Mean   :-208200   Mean   :-4.926   Mean   :227.2   Mean   :440.8  
##  3rd Qu.:-105300   3rd Qu.:-3.466   3rd Qu.:248.0   3rd Qu.:455.9  
##  Max.   :  -2400   Max.   : 3.187   Max.   :298.3   Max.   :496.7

dat[1:2,]

##       year  temp     co2   solar
## 82 -414000 0.837 285.260 443.192
## 83 -413900 0.827 284.977 443.480

The canonical graph that is being challenged is

plot(dat$year/1000, scale(dat$temp), col="blue", xlab="Relative Year, in Thousands", ylab="Z-Score",
       type="l", lwd=3, main="Misleading Association Graph of 400,000 Years")

plot(dat$year/1000, scale(dat$temp), col="blue", xlab="Relative Year, in Thousands", ylab="Z-Score",
     type="l", lwd=3, main="Misleading Association Graph of 400,000 Years")
lines(dat$year/1000, scale(dat$co2), col="brown", lwd=3)
lines(dat$year/1000, scale(dat$solar), col="orange")
legend("topleft", col=c("brown","blue"), lty=1, box.col="white", lwd=3, legend=c("CO2","Temp"))

Regression Analysis:

Again from the source:

Clearly, CO2 and temperature move together. The associated verbal suggestion — sometimes explicit, sometimes left to the audience’s imagination — is often that CO2 determines temperature. The visually implied X-Y association can also be graphed:

plot(dat$co2, dat$temp, xlab="CO2", ylab="Temperature", 
     main="Misleading Association Graph of 400,000 Years, XY", cex=0.4)
lines(dat$co2,dat$temp, col="gray")
abline(lm(dat$temp ~ dat$co2), col="blue", lwd=2)

The blue line is the best-fit OLS line, \(\text{temp}_t=−24+0.084\times\text{CO2}_t\). The suggestion in both plots is that an increase in \(100\) ppm of CO2 (about \(33\%\) of the mean CO2 in the data) predicts warming of \(8.4\)°C. This is of course absurd.

The current level of CO2 may be around \(410\) ppm. An increase of \(100\) ppm would result, according to the regression coefficients, in a temperature increase of \(9.6\) °C!

Once again from the source:

This naive interpretation, that this graph even suggests that CO2 drove global warming, is misleading and therefore wrong.

Anyone who understands basic data analysis should understand this plot is misleading as far as establishing a (causal) link from CO2 to temperature is concerned. Indeed, the only point of my analysis is to plead not to use this graph any longer. This plotted relationship is a misleading and classic example of a spurious relation. A classic example is the association between ice cream sales and murders. Both are higher in summer, and the two plots between ice-cream sales and murder would look just like two plots of CO2 and temperature above.

There are better ways to analyze the CO2, temperature, and solar data, shown below. These better ways address the facts that the graph misleads with respect to two s:

Confounding: Could a third variable — such as trends, volcanos, solar radiation, or anything else — have caused (co-)variation in both CO2 and temperature?

Confounding: A third variable, not the dependent (outcome) or main independent (exposure) variable of interest, that distorts the observed relationship between the exposure and outcome.

Causality: Is CO2 causing warming or is warming causing CO2, or are both causing one another?

From Wikipedia:

Ordinarily, regressions reflect “mere” correlations, but Clive Granger argued that causality in economics could be tested for by measuring the ability to predict the future values of a time series using prior values of another time series. Since the question of “true causality” is deeply philosophical, and because of the post hoc ergo propter hoc fallacy of assuming that one thing preceding another can be used as a proof of causation, econometricians assert that the Granger test finds only “predictive causality”.

The remedy to the first is to work in changes of variables, not in levels of variables. The remedy to the second is to work with lead-lag associations. I am not the first to have noticed that temperature changes can also anticipate CO2 changes. However, some climate-change critics have jumped to the equally incorrect conclusion that such feedback effects then reject the hypothesis that CO2 drives temperature. Feedbacks are not mutually exclusive with respect to the hypothesis of interest, which is whether CO2 changes anticipate temperature changes.

Models:

Confounders \((1)\):

If the model is correct:

\[y_t = a + b\times x_t\] it follows that

\[\Delta y_t=0+b \times \Delta x_t\] > where Δ is the change (also called the first difference). A good test uncovering many spurious correlations is to estimate regressions both in levels and differences. If the coefficient \(b\) is not the same (or at least similar) in both regression models, it suggests that the level correlation is spurious. In such a case, we would have learned that \(y_t = a + b\times x_t\) was not the correct model to begin with.

The original data set is altered using the defined functions chg and lag that calculate along columns the differences between one measurement and the preceding (chg), and slide the values one step along the time line:

# chg is a function with three lines separated by semicolons:
# shift(x) simply moves the x[1] to x[2], x[2] to x[3], etc.
# the resultant vector is subtracted from the original x
# the result is stored in a local variable "o" inside the function
# the second line assigns names to all the vectors in x (x is a df)
# the names are the original names of x with a "d" in front.
# the third line says that "o" is returned.
# lag simply slides all the entries in the df by 1 step and relabels.
chg <- function(x,...) { o <- x - shift(x,...); names(o) <- paste0("d",names(x)); o }
lag <- function(x,...) { o <- shift(x,...); names(o) <- paste0("l",names(x)); o }
   ## the above use data.table's shift function
ds <- cbind(dat, chg(dat)); rownames(ds) <- NULL  ## combine levels and changes into one data set
# define some variables used later as controls
# first some relabeling of variables:
ds <- within(ds, {lagsolar <- lag(solar); lagco2 <- lag(co2); lagtemp <- lag(temp)})
ds <- ds[complete.cases(ds),]
head(ds)[1:2,]

##      year  temp     co2   solar dyear  dtemp   dco2 dsolar lagtemp  lagco2
## 2 -413900 0.827 284.977 443.480   100 -0.010 -0.283  0.288   0.837 285.260
## 3 -413800 0.818 284.694 443.778   100 -0.009 -0.283  0.298   0.827 284.977
##   lagsolar
## 2  443.192
## 3  443.480

ds2 <- subset(ds, select=c(dco2,dtemp)) 
# This dataframe can be converted into a time series:
tsds <- ts(ds2, start = -414000, deltat = 100)
head(tsds)[1:2,]

##        dco2  dtemp
## [1,] -0.283 -0.010
## [2,] -0.283 -0.009

Running the original regression:

# The function specified below, using the update() function for the model is explained in here:
# https://www.r-bloggers.com/2010/05/using-the-update-function-during-variable-selection/
showcoef <- function(formula, controls = (~.), data=ds){
  round(coef(summary(lm(update(formula, controls),  data=data)))[,c(1,3,4)],4)
} 

showcoef(temp ~ co2)

##             Estimate   t value Pr(>|t|)
## (Intercept) -23.9062 -138.7603        0
## co2           0.0835  110.9774        0

Comparing to the model specified with first differences:

\[\Delta \text{temp}_t=a+b\times (\Delta \text{CO}_t)\]

showcoef(dtemp ~ dco2)

##             Estimate t value Pr(>|t|)
## (Intercept)  -0.0004 -0.1018    0.919
## dco2          0.0331  6.2510    0.000

The coefficient estimate of 0.03 is much smaller than 0.08, suggesting spurious level trend correlation in the prior regression.

The same coefficients should show with second differences:

showcoef(chg(dtemp) ~ chg(dco2))

##             Estimate t value Pr(>|t|)
## (Intercept)  -0.0003 -0.0545   0.9565
## chg(dco2)     0.0177  1.4155   0.1570

but again, the coefficients are different.

Causality \((2)\):

Note that the above regression inputs were still contemporaneous. They only solve the spurious correlation issue with respect to trends. They do not address the question of whether CO2 drives warming or vice-versa.

A better test is based on the idea that if CO2 really changes temperature, then (unexpected) changes in CO2 should anticipate (unexpected) changes in temperature. Econometricians call this Granger-Sims causality (GSC). GSC is a necessary but not a sufficient condition to read potential causality into data, even if there are no important omitted variables. If there is no GSC, then the data cannot show that CO2 causally influences temperature. GSC tests may be too weak in actual data to find an association that exists.

This is addressed first with simple OLS, and later with Sims’ Vector Auto-Regression (VAR)

Prelimary OLS for CO2 levels:

First, consider changes in CO2. The strongest association in the data, by far, is that they are highly autocorrelated:

par(mfrow = c(1,2))
acf(tsds[,1],
     main="ACF", 
     xlab="Lag",     
     ylab="ACF", 
     ylim=c(-0.1,1), # this sets the y scale to -0.2 to 0.2
     las=1,
     xaxt="n")
pacf(diff(tsds[,1]), main= 'PACD')

The geometric decay in the ACF suggests and auto-regressive process.

This can be seen also with a simple OLS between changes in values of CO2 at one time regressed over lagged values:

showcoef(dco2 ~ lag(dco2))

##             Estimate t value Pr(>|t|)
## (Intercept)   0.0002  0.0306   0.9756
## lag(dco2)     0.8384 98.5864   0.0000

This is stable and applies even to changes in changes in CO2:

showcoef(lag(dco2) ~ lag(lag(dco2)))

##                Estimate t value Pr(>|t|)
## (Intercept)      0.0002  0.0254   0.9797
## lag(lag(dco2))   0.8383 98.5548   0.0000

The data inform us that once CO2 changes got going into a particular direction, they tended to continue the same direction. CO2 changes had strong internal dynamics, almost random-walk like: whenever CO2 increased over 100 years, it strongly tended to increase over the next 100 years again, and by almost as much. Ergo, when a shock to CO2 occurs, it has long-term effects, far beyond a century. The half-life of shocks to changes in CO2 is about \(−\log(2)/\log(0.83)≈3.7\) centuries.

Empirically, the autocorrelation of CO2 is what will make it difficult to determine how CO2 changes influenced temperature changes — CO2 changes over the last 100 years tend to look very similar to CO2 changes from, say, 500 to 600 years ago.

Predicting Changes in Temperature:

The true variable of interest are changes in temperatures, not changes in CO2. Rather than start with the simplest regression, the following regression already includes a host of controls:

# In what follows, the regression equation dtemp ~ lag(solar) we are looking at 
# differences in temperature explained by solar radiation 
# Why? Because the lag solar means that the independent variable
# is data from the past: in x[2] there is now the previous entry from x[1].
# The last value of the original solar vector is thrown out.

# The model specified above, using the update() function for the model as
# https://www.r-bloggers.com/2010/05/using-the-update-function-during-variable-selection/
# ends up just printing up the full model with coefficients and t values:
# m = lm(dtemp ~ lag(dco2) + lag(dsolar) + lag(dtemp) +
# lag(temp) + lag(dtemp)*lag(temp) + lag(co2) + lag(solar) + dsolar, ds)
# coef(summary(m))[,c(1,3)]

showcoef(dtemp ~ lag(dco2) , ~ . + lag(dsolar) + lag(dtemp) +
            lag(temp) + lag(dtemp)*lag(temp) + lag(co2) + lag(solar) + dsolar)

##                      Estimate t value Pr(>|t|)
## (Intercept)           -0.7789 -5.2031   0.0000
## lag(dco2)              0.0305  5.7282   0.0000
## lag(dsolar)            0.0354  0.3695   0.7118
## lag(dtemp)            -0.0959 -3.3599   0.0008
## lag(temp)             -0.0215 -6.2229   0.0000
## lag(co2)               0.0015  4.6059   0.0000
## lag(solar)             0.0008  3.3917   0.0007
## dsolar                -0.0271 -0.2827   0.7774
## lag(dtemp):lag(temp)  -0.0471 -8.6498   0.0000

And here are the main points:

This regression suggests the following:

When solar radiation was high lag(solar), temperature tended to increase (\(0.00076\)). Changes in solar radiation dsolar were fairly unimportant (\(-0.02708\)), suggesting a very slow response of temperature to solar radiation.

When lagged CO2 lag(co2) was high (\(0.0015\)) and when lagged CO2 increased recently lag(dco2) (\(0.03046\)), temperature tended to increase. This is good evidence that changes in CO2 influenced future warming, both short-term and long-term. (The relation is robust to inclusion or exclusion of many control variables, such as the state variables included here, too.)

Earth has a strong thermostat: (a) when temperature lag(temp) was high, it tended to go down (\(-0.021\)); (b) when temperature recently has gone up lag(dtemp), it tended to go down again (\(-0.096\)); and (c) temperature really wanted to go down again if both temperature was high and temperature had recently gone up lag(dtemp):lag(temp) (\(-0.047\)).

There are two interesting earth-science questions beyond my expertise related to this regression output.

The first concerns earth’s thermostat. It seems to work, even controlling for solar forcing and CO2. What is its cause? Will it continue?

The second concerns the magnitude of the coefficient estimate on lag(dco2). It’s still way too big. It suggests that on the margin, an increase of \(250\) ppm (about doubling) in CO2 predicts global warming of nearly \(0.03046×250≈7.6 °C.\) Standard climate change models would suggest increases of less than half this much, about \(3\) °C. The \(400,000\) association between lagged CO2 changes and temperature changes was far too strong. (This is even more worrisome because Archer suggests that about \(3/4\) of the CO2 disappears in the carbon cycle rapidly before it has much opportunity to drive the greenhouse effect.)

Vector Autoregressions (VARs)

The improved statistical estimation method are Sims’ vector autoregressions. They estimate equations on both ΔCO2 dynamics and ΔTemp dynamics together to disentangle better how they influence one another (in innovations, too). The specification explicitly allows for the two variables to influence one another, too.

library(vars)
library(ggplot2)
plot.ts(tsds, col=2)

(author)We begin with a one-lag VAR analysis. The format of the coefficient-test output is first the dependent variable, then the independent variables, then an indicator of the lag.

To determine the optimal lag:

VARselect(tsds, lag.max = 10, type='const')$selection

## AIC(n)  HQ(n)  SC(n) FPE(n) 
##     10      9      2     10

So let’s try \(10\):

Arguments to the VAR() function:

y Data item containing the endogenous variables p Integer for the lag order (default is p=1). type Type of deterministic regressors to include. season Inlusion of centered seasonal dummy variables (integer value of frequency). exogen Inlusion of exogenous variables. lag.max Integer, determines the highest lag order for lag length selection according to the choosen ic. ic Character, selects the information criteria, if lag.max is not NULL. x Object with class attribute ‘varest’. digits the number of significant digits to use when printing

var1 <- VAR(tsds, type="const", lag.max=10, season=NULL, exogen=NULL)
print(var1)

## 
## VAR Estimation Results:
## ======================= 
## 
## Estimated coefficients for equation dco2: 
## ========================================= 
## Call:
## dco2 = dco2.l1 + dtemp.l1 + dco2.l2 + dtemp.l2 + dco2.l3 + dtemp.l3 + dco2.l4 + dtemp.l4 + dco2.l5 + dtemp.l5 + dco2.l6 + dtemp.l6 + dco2.l7 + dtemp.l7 + dco2.l8 + dtemp.l8 + dco2.l9 + dtemp.l9 + dco2.l10 + dtemp.l10 + const 
## 
##       dco2.l1      dtemp.l1       dco2.l2      dtemp.l2       dco2.l3 
##  0.8486433333  0.0583904532  0.0207034282  0.0130353230 -0.0384191281 
##      dtemp.l3       dco2.l4      dtemp.l4       dco2.l5      dtemp.l5 
##  0.0500279087  0.0009740159  0.0845627997 -0.0152519079  0.0003582399 
##       dco2.l6      dtemp.l6       dco2.l7      dtemp.l7       dco2.l8 
##  0.0003833884  0.0977294343 -0.1183632879  0.0661121849  0.0688685879 
##      dtemp.l8       dco2.l9      dtemp.l9      dco2.l10     dtemp.l10 
##  0.0690137940  0.0337005917  0.0726874924  0.0103238450  0.0090563637 
##         const 
##  0.0005004690 
## 
## 
## Estimated coefficients for equation dtemp: 
## ========================================== 
## Call:
## dtemp = dco2.l1 + dtemp.l1 + dco2.l2 + dtemp.l2 + dco2.l3 + dtemp.l3 + dco2.l4 + dtemp.l4 + dco2.l5 + dtemp.l5 + dco2.l6 + dtemp.l6 + dco2.l7 + dtemp.l7 + dco2.l8 + dtemp.l8 + dco2.l9 + dtemp.l9 + dco2.l10 + dtemp.l10 + const 
## 
##       dco2.l1      dtemp.l1       dco2.l2      dtemp.l2       dco2.l3 
##  0.0218955830  0.1116157807 -0.0059990359 -0.2423288442  0.0157072166 
##      dtemp.l3       dco2.l4      dtemp.l4       dco2.l5      dtemp.l5 
## -0.0256826278  0.0015560968 -0.0367748120 -0.0125179205  0.0139528635 
##       dco2.l6      dtemp.l6       dco2.l7      dtemp.l7       dco2.l8 
##  0.0229467028  0.0133577315 -0.0162688651 -0.0052952902  0.0370506907 
##      dtemp.l8       dco2.l9      dtemp.l9      dco2.l10     dtemp.l10 
##  0.0277687735 -0.0115820798  0.0391833528  0.0081247837 -0.0460107469 
##         const 
## -0.0005536633

# No exogenous variables.
# No seasonality
# No trend: type = c("const", "trend", "both", "none")
summary(var1)

## 
## VAR Estimation Results:
## ========================= 
## Endogenous variables: dco2, dtemp 
## Deterministic variables: const 
## Sample size: 4106 
## Log Likelihood: -2876.973 
## Roots of the characteristic polynomial:
## 0.9037 0.824 0.824 0.7911 0.7911 0.7848 0.7848 0.7539 0.7539 0.7512 0.7512 0.711 0.711 0.7082 0.6718 0.6718 0.6556 0.6556 0.3226 0.3226
## Call:
## VAR(y = tsds, type = "const", exogen = NULL, lag.max = 10)
## 
## 
## Estimation results for equation dco2: 
## ===================================== 
## dco2 = dco2.l1 + dtemp.l1 + dco2.l2 + dtemp.l2 + dco2.l3 + dtemp.l3 + dco2.l4 + dtemp.l4 + dco2.l5 + dtemp.l5 + dco2.l6 + dtemp.l6 + dco2.l7 + dtemp.l7 + dco2.l8 + dtemp.l8 + dco2.l9 + dtemp.l9 + dco2.l10 + dtemp.l10 + const 
## 
##             Estimate Std. Error t value             Pr(>|t|)    
## dco2.l1    0.8486433  0.0156698  54.158 < 0.0000000000000002 ***
## dtemp.l1   0.0583905  0.0258210   2.261             0.023789 *  
## dco2.l2    0.0207034  0.0205232   1.009             0.313138    
## dtemp.l2   0.0130353  0.0260615   0.500             0.616979    
## dco2.l3   -0.0384191  0.0204984  -1.874             0.060967 .  
## dtemp.l3   0.0500279  0.0267767   1.868             0.061787 .  
## dco2.l4    0.0009740  0.0204296   0.048             0.961976    
## dtemp.l4   0.0845628  0.0267952   3.156             0.001612 ** 
## dco2.l5   -0.0152519  0.0204259  -0.747             0.455292    
## dtemp.l5   0.0003582  0.0268052   0.013             0.989338    
## dco2.l6    0.0003834  0.0204144   0.019             0.985017    
## dtemp.l6   0.0977294  0.0268062   3.646             0.000270 ***
## dco2.l7   -0.1183633  0.0204116  -5.799        0.00000000718 ***
## dtemp.l7   0.0661122  0.0268337   2.464             0.013789 *  
## dco2.l8    0.0688686  0.0204893   3.361             0.000783 ***
## dtemp.l8   0.0690138  0.0268357   2.572             0.010155 *  
## dco2.l9    0.0337006  0.0205348   1.641             0.100844    
## dtemp.l9   0.0726875  0.0261141   2.783             0.005403 ** 
## dco2.l10   0.0103238  0.0156443   0.660             0.509347    
## dtemp.l10  0.0090564  0.0260073   0.348             0.727690    
## const      0.0005005  0.0068940   0.073             0.942132    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Residual standard error: 0.4418 on 4085 degrees of freedom
## Multiple R-Squared: 0.7126,  Adjusted R-squared: 0.7112 
## F-statistic: 506.4 on 20 and 4085 DF,  p-value: < 0.00000000000000022 
## 
## 
## Estimation results for equation dtemp: 
## ====================================== 
## dtemp = dco2.l1 + dtemp.l1 + dco2.l2 + dtemp.l2 + dco2.l3 + dtemp.l3 + dco2.l4 + dtemp.l4 + dco2.l5 + dtemp.l5 + dco2.l6 + dtemp.l6 + dco2.l7 + dtemp.l7 + dco2.l8 + dtemp.l8 + dco2.l9 + dtemp.l9 + dco2.l10 + dtemp.l10 + const 
## 
##             Estimate Std. Error t value             Pr(>|t|)    
## dco2.l1    0.0218956  0.0095372   2.296              0.02174 *  
## dtemp.l1   0.1116158  0.0157156   7.102     0.00000000000144 ***
## dco2.l2   -0.0059990  0.0124911  -0.480              0.63107    
## dtemp.l2  -0.2423288  0.0158620 -15.277 < 0.0000000000000002 ***
## dco2.l3    0.0157072  0.0124761   1.259              0.20811    
## dtemp.l3  -0.0256826  0.0162973  -1.576              0.11513    
## dco2.l4    0.0015561  0.0124342   0.125              0.90041    
## dtemp.l4  -0.0367748  0.0163085  -2.255              0.02419 *  
## dco2.l5   -0.0125179  0.0124320  -1.007              0.31404    
## dtemp.l5   0.0139529  0.0163146   0.855              0.39247    
## dco2.l6    0.0229467  0.0124249   1.847              0.06484 .  
## dtemp.l6   0.0133577  0.0163152   0.819              0.41299    
## dco2.l7   -0.0162689  0.0124233  -1.310              0.19042    
## dtemp.l7  -0.0052953  0.0163320  -0.324              0.74578    
## dco2.l8    0.0370507  0.0124706   2.971              0.00299 ** 
## dtemp.l8   0.0277688  0.0163332   1.700              0.08918 .  
## dco2.l9   -0.0115821  0.0124982  -0.927              0.35414    
## dtemp.l9   0.0391834  0.0158940   2.465              0.01373 *  
## dco2.l10   0.0081248  0.0095217   0.853              0.39355    
## dtemp.l10 -0.0460107  0.0158290  -2.907              0.00367 ** 
## const     -0.0005537  0.0041960  -0.132              0.89503    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Residual standard error: 0.2689 on 4085 degrees of freedom
## Multiple R-Squared: 0.08596, Adjusted R-squared: 0.08148 
## F-statistic: 19.21 on 20 and 4085 DF,  p-value: < 0.00000000000000022 
## 
## 
## 
## Covariance matrix of residuals:
##           dco2    dtemp
## dco2  0.195143 0.006539
## dtemp 0.006539 0.072288
## 
## Correlation matrix of residuals:
##          dco2   dtemp
## dco2  1.00000 0.05506
## dtemp 0.05506 1.00000

Notice that all the roots are inside the unit circle: errors withing contour lines of Gaussian distribution in 3D.

Test for serial correlation:

serial.test(var1, lags.pt=-4, type="PT.asymptotic")

## 
##  Portmanteau Test (asymptotic)
## 
## data:  Residuals of VAR object var1
## Chi-squared = 0.90972, df = -24, p-value = NA

The low p value is bad news: there is serial correlation.

Heteroskadisticity:

arch.test(var1, multivariate.only=T)

## 
##  ARCH (multivariate)
## 
## data:  Residuals of VAR object var1
## Chi-squared = 1937.2, df = 45, p-value < 0.00000000000000022

Bad news again for heteroskadisticity.

Normality test residuals:

normality.test(var1, multivariate.only=T)

## $JB
## 
##  JB-Test (multivariate)
## 
## data:  Residuals of VAR object var1
## Chi-squared = 14724668, df = 4, p-value < 0.00000000000000022
## 
## 
## $Skewness
## 
##  Skewness only (multivariate)
## 
## data:  Residuals of VAR object var1
## Chi-squared = 6584.7, df = 2, p-value < 0.00000000000000022
## 
## 
## $Kurtosis
## 
##  Kurtosis only (multivariate)
## 
## data:  Residuals of VAR object var1
## Chi-squared = 14718083, df = 2, p-value < 0.00000000000000022

The model fails the normality test.

Trying Granger causality (from here):

Granger <- causality(var1, cause='dco2')
Granger$Granger

## 
##  Granger causality H0: dco2 do not Granger-cause dtemp
## 
## data:  VAR object var1
## F-Test = 7.318, df1 = 10, df2 = 8170, p-value = 0.00000000001223

Impulse-response function:

tempIRF <- irf(var1, impulse='dco2', response='dtemp', n.ahead=10, boot=T)
plot(tempIRF, ylab='temp', xlab='CO2', main='Shock from sudden change in Co2')

In the other direction:

Granger2 <- causality(var1, cause='dtemp')
Granger2$Granger

## 
##  Granger causality H0: dtemp do not Granger-cause dco2
## 
## data:  VAR object var1
## F-Test = 4.5873, df1 = 10, df2 = 8170, p-value = 0.000001581

Granger test from lmtest (from here):

library(lmtest)
grangertest(dtemp ~ dco2, order=10, data=tsds)

## Granger causality test
## 
## Model 1: dtemp ~ Lags(dtemp, 1:10) + Lags(dco2, 1:10)
## Model 2: dtemp ~ Lags(dtemp, 1:10)
##   Res.Df  Df     F          Pr(>F)    
## 1   4085                              
## 2   4095 -10 7.318 0.0000000000139 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# H0: Time series x does not cause time series y to grange cause itself
# H1: Time series x does cause time series y to grange cause itself
# A low P VALUE: Knowing the values of CO2 is valuable for forecasting future temps.

# Compared to the test from `vars`:

causality(VAR(tsds, p = 3), cause = 'dco2')$Granger

## 
##  Granger causality H0: dco2 do not Granger-cause dtemp
## 
## data:  VAR object VAR(tsds, p = 3)
## F-Test = 16.018, df1 = 3, df2 = 8212, p-value = 0.0000000002219

In both cases we have to conclude that CO2 predicts future temperatures.

Impulse-response function:

tempIRF <- irf(var1, impulse='dtemp', response='dco2', n.ahead=10, boot=T)
plot(tempIRF, ylab='CO2', xlab='temp', main='Shock from sudden change in temp')

Variance decomposition:

var_dec <- fevd(var1, n.ahead=10)
plot(var_dec, col=c('orange','darkblue'), border=F)

The shock effects on each variable is mainly due to itself: shocks in CO2 do not change the future values of temperature.

Forecasting:

library(forecast)
library(tseries)
fc <- predict(var1, n.ahead=100, ci=0.95, col=2)
fanchart(fc, colors=heat.colors(9))

What about temperatures over CO2 levels (as opposed to changes in these values):

ds3 <- subset(ds, select=c(co2, temp)) 
# This dataframe can be converted into a time series:
tsds2 <- ts(ds3, start = -414000, deltat = 100)
head(tsds2)[1:2,]

##          co2  temp
## [1,] 284.977 0.827
## [2,] 284.694 0.818

vartempco2 <- VAR(tsds2, type='both', lag.max=10, season=NULL, exogen=NULL)
Granger2 <- causality(vartempco2, cause='co2')
Granger2$Granger

## 
##  Granger causality H0: co2 do not Granger-cause temp
## 
## data:  VAR object vartempco2
## F-Test = 8.2333, df1 = 10, df2 = 8168, p-value = 0.0000000000002069

Granger3 <- causality(vartempco2, cause='temp')
Granger3$Granger

## 
##  Granger causality H0: temp do not Granger-cause co2
## 
## data:  VAR object vartempco2
## F-Test = 7.2498, df1 = 10, df2 = 8168, p-value = 0.00000000001654

fc2 <- predict(vartempco2, n.ahead=100, ci=0.95, col=2)
fanchart(fc2, colors=heat.colors(9))

This suggests (as before) that

Carbon-dioxide changes, ΔCO2 are highly autocorrelated (\(0.837\)). When CO2 has increased, it wants to continue to increase. When CO2 has decreased, it wants to continue to decrease. Of all the association in this \(400,000\) year data, it is by far the strongest one. It is evidence of strong buffers and/or a strong CO2 feedback effect on earth.

When temperature has recently gone up, ΔCO2 wants to go up just a little more (\(0.046\)).

When temperature has recently gone up, then the temperature wants to go up just a little more (\(0.1\)). This disappears with better control for the level of temperature and recent temperature changes interacted. The reason can be inferred due to the acclerating shape of the \(xy\) graph above.

The association of most interest to us: When CO2 has recently gone up, then the temperature wants to go up just a little more. This is what we found before, and the magnitude of the coefficient remains troubling. The coefficient is still far too large, suggesting a warming effect of about \(2.4\times 250≈6\)°C for a doubling of CO2. And this is even more disconcerting, because it is not even for a recent 1-50 year increase but for a 100-year lagged increase in CO2.

Expected Decay of Autocoefficients With Lag:

The theory further predicts that the coefficient on lag CO2 should decrease with lag. A change in CO2 this century should have more ability to predict the temperature in the next 100 years than, say, in 100 years in five centuries. This is a quasi-placebo test. There should be very little association beyond the first one or two lags.

The following includes 10 lags of CO2 changes, i.e., the last ten centuries. The printouts are only for the coefficients on the ΔCO2 predictors in the Δtemp prediction, although the analysis itself remains based on the full VAR.

var10 <- VAR(tsds, type="none", lag.max=10)
## with controls: use VAR(ds2, lag.max=10, exog= subset(ds, T, select=c(lagsolar, lagco2, lagtemp)))
var10

## 
## VAR Estimation Results:
## ======================= 
## 
## Estimated coefficients for equation dco2: 
## ========================================= 
## Call:
## dco2 = dco2.l1 + dtemp.l1 + dco2.l2 + dtemp.l2 + dco2.l3 + dtemp.l3 + dco2.l4 + dtemp.l4 + dco2.l5 + dtemp.l5 + dco2.l6 + dtemp.l6 + dco2.l7 + dtemp.l7 + dco2.l8 + dtemp.l8 + dco2.l9 + dtemp.l9 + dco2.l10 + dtemp.l10 
## 
##       dco2.l1      dtemp.l1       dco2.l2      dtemp.l2       dco2.l3 
##  0.8486444688  0.0583891934  0.0207034312  0.0130321837 -0.0384193134 
##      dtemp.l3       dco2.l4      dtemp.l4       dco2.l5      dtemp.l5 
##  0.0500250232  0.0009739560  0.0845594964 -0.0152518932  0.0003554732 
##       dco2.l6      dtemp.l6       dco2.l7      dtemp.l7       dco2.l8 
##  0.0003834058  0.0977259933 -0.1183633612  0.0661098105  0.0688687199 
##      dtemp.l8       dco2.l9      dtemp.l9      dco2.l10     dtemp.l10 
##  0.0690125369  0.0337006275  0.0726852537  0.0103224585  0.0090555516 
## 
## 
## Estimated coefficients for equation dtemp: 
## ========================================== 
## Call:
## dtemp = dco2.l1 + dtemp.l1 + dco2.l2 + dtemp.l2 + dco2.l3 + dtemp.l3 + dco2.l4 + dtemp.l4 + dco2.l5 + dtemp.l5 + dco2.l6 + dtemp.l6 + dco2.l7 + dtemp.l7 + dco2.l8 + dtemp.l8 + dco2.l9 + dtemp.l9 + dco2.l10 + dtemp.l10 
## 
##      dco2.l1     dtemp.l1      dco2.l2     dtemp.l2      dco2.l3     dtemp.l3 
##  0.021894327  0.111617174 -0.005999039 -0.242325371  0.015707422 -0.025679436 
##      dco2.l4     dtemp.l4      dco2.l5     dtemp.l5      dco2.l6     dtemp.l6 
##  0.001556163 -0.036771158 -0.012517937  0.013955924  0.022946684  0.013361538 
##      dco2.l7     dtemp.l7      dco2.l8     dtemp.l8      dco2.l9     dtemp.l9 
## -0.016268784 -0.005292663  0.037050545  0.027770164 -0.011582119  0.039185829 
##     dco2.l10    dtemp.l10 
##  0.008126318 -0.046009848

plotvarcoef <- function(ctbl, vnm="CO2", yl=0.05) {
   plot( 1:nrow(ctbl), ctbl[,1], xlab=paste("100-Year Lag of",vnm,"Change"), ylab=paste("Coefficient on",vnm,"Change"), type="b", ylim=c(-yl,yl), main="Explaining 100-Year Ahead Temperature Changes")
   lines( 1:nrow(ctbl), ctbl[,1] + 1*ctbl[,2], col="gray", lty=2 )
   lines( 1:nrow(ctbl), ctbl[,1] - 1*ctbl[,2], col="gray", lty=2 )
   lines( 1:nrow(ctbl), ctbl[,1] + 2*ctbl[,2], col="gray", lty=3 )
   lines( 1:nrow(ctbl), ctbl[,1] - 2*ctbl[,2], col="gray", lty=3 )
   lines( c(0,20), c(0,0), lty=2, col="gray")
   points( 1, ctbl[1,1], col="blue", cex=2)
}

coef.dtemp <- coef(var10)$dtemp
coef.dtemp.dco2 <- coef.dtemp[grepl("dco2", rownames(coef.dtemp)),]
# print(coef.dtemp.dco2[,1])
plotvarcoef(coef.dtemp.dco2)

Archer explains that the theory says that the coefficients should be decaying. More specifically, theory predicts a long-run coefficient of about \(0.01\) (in this sample, a \(250\) ppm increase should induce a \(3\)°C increase). It should be split into about \(0.0075\) on lag0, \(0.002\) on lag1, and lower coefficients on further lags. The most recent CO2 change should be more powerful than more lagged CO2 changes.

If I were to claim that this data suggests that changes in CO2 have driven changes in temperature, then

why are CO2 changes from many centuries ago similarly powerful as the most recent CO2 changes? why are the coefficient estimates so large?

The statistical reason for the first part of this mess is the high correlation among CO2 changes. When CO2 increased in the last 100 years, it also likely increased before. As far as the regression is concerned, many recently past CO2 changes look somewhat alike in their power and could have been responsible for their influence on predicting temperature changes. 4,000 centuries should have been enough to uncover the relationship, but just weren’t. The reason for the second part of this mess, the terribly high coefficient estimates are a mystery to me.

The data absolutely do not reject the hypothesis that changes in CO2 drive changes in global warming. (I would go further and characterize this as “they hint at an association.”) The data just do not reject the hypothesis that the relationship is not strong enough to identify this relationship cleanly. Thus, the use of the visual graph at the outset is not only misleading (for ignoring the reverse association), it is badly misleading.

CO2 and Global Temperature

NOTES ON STATISTICS, PROBABILITY and MATHEMATICS