On Spurious Causality, CO₂, and Global Temperature

Causality is fundamental to science. While causal statements can reasonably be made without hardship in controlled environments, things are far less straightforward when only observational data is available. Answering the question of what happens to $Y$ if one intervenes on $X$ is compromised by the simple fact that $X$ and $Y$ were not generated by exogenously modulating $X$, but by (often) endogenous interactions between the two variables.

Yet, for many scientific interrogations of capital importance – like the relationship between greenhouse gases and global temperature, only observational data is available. How does different forcings cause global mean surface temperature anomalies (GMTA)? In a popular article, Stips, Macias, Coughlan, Garcia-Gorriz, and Liang (2016) (henceforth, SMCGL) set up to answer the question using information flow (IF) from one variable to another for bivariate stochastic dynamical systems – a methodology developed in Liang2008; Liang2014; Liang2015; Liang2016. Authors make grand claims about the technique being able to extract "true" and "rigorous" causality, which clash with the usually somber tone of causal analysis.

We show that the IF methodology will not extract causality for the near-universe of bivariate stochastic dynamic systems estimated on real data – with those of SMCGL included. In a nutshell, this occurs because Liang2008 formulas being used throughout, assume that a certain statistical quantity is 0 when it is not. This can be tested, and it is rejected almost all the time. More precisely, Liang2014 formula assumes that when conditioning on the past state of the system ($X_{t-1}$ and $Y_{t-1}$), the remaining variation in $X_{t}$ and $Y_{t}$ (the innovations driving the system) are uncorrelated. In SMCGL, where the time unit $t$ is one year, this implies that forcings and GMTA are uncorrelated within the same year, conditional on last year’s values of both. That is unlikely to happen, and the data will testify to that. The assumption is presented as "linearity" (Liang2014), but it has little to do with it. Rather, it assumes that the system is identified, meaning that the joint dynamics of the time series data fall on the knife edge case where causality can indeed be claimed from the data without any external assumption (sims1980). The problem is that for discretely sampled time series, the assumption is almost always rejected by lower-frequency data.¹¹1In tawia2019time, IF are used on daily data, which can alleviate the problem if there are no intra-day relationships. This last condition is something that should be verified, not assumed. In short, SMCGL assume there to be no identification problem (in the sense of sims1980), and find no identification problem.

In this note, we detail the consequences of this high-standing omission on IF measures and SMCGL’s results. First, in simulations (using data generating processes proposed by Liang2014) of which we know the true causality structure, it is shown that IF will often conclude that $X_{t}$ is largely causing $Y_{t}$, when in fact the reverse is true. Second, we reconsider a key part of SMCGL’s study where the authors investigate the causal structure between different forcings and GMTA. Using an appropriate methodology that accounts for the correlated innovations, it is found that in most instances, the data by itself cannot back SMCGL’s claims. That is, unlike what the authors have put forward, it is not possible (within this framework) to claim that many forcings are causing GMTA’s increase as a direct implication of the data. In other words, from the data alone, it is impossible to discriminate between certain forcings mostly causing GMTA or the reverse. External assumptions based on physical knowledge could remedy that. Or different data. We show that in the case of CO₂, results in accord with the scientific consensus can be obtained when using concentration directly rather than its radiative content.

The alternative methodology that we use are structural Vector Autoregressions (SVAR) which are simultaneous dynamic systems of equations. They can characterize a linear dynamic system in discrete time. The methodology was introduced to macroeconomics by sims1980 and is now so in many fields, ranging from neuroscience (chen2011vector) to climate (VARCTIC). To document its reliability for this application, we report implied transient climate response (TCR) estimates of our alternative methodology. We find that those (i) largely depend on the necessary assumptions made about the simultaneous (with one year) impact of CO₂ forcing on GMTA but (ii) are in the range of recent estimates (OttoEtAl2013; MontamatStock2020) if one assumes simultaneous causality running from CO₂ forcing to GMTA.

Section 2 briefly review IFs and explain where the problematic assumption occurs. It also discusses relevant notions of Vector Autoregressions (VARs, sims1980) as a comprehensive framework to think about causality in multivariate time series systems. Section 3 display IFs possibly spurious behavior using simulated data where the true causality is known. Section 4 revisit the question of causality between different forcings and GMTA using appropriate tools. Section 5 concludes.

In this section, we review the basics of IFs as applied in empirical work, pin down the problematic assumption, explain why it is harmful through the lenses of a VAR, and propose an alternative measure of IF based on the VAR.

Liang2008 considers the data generating process (DGP)

where $\sigma_{i,i}$ is the variance of $i$, $\sigma_{ij,}$ the covariance between $i$ and $j$, and $\sigma_{i,\Delta i}$ is the covariance between $i$ and the $k^{th}$ difference of $i$ (we follow Liang2015 and set $k=1$). However, the validity of this appealing formula rests on the seemingly technical assumption of a diagonal $\mathbf{B}$. Our point is that this assumption is far from merely technical and very frequently wrong. Motivating the diagonal $\mathbf{B}$ among other assumptions, Liang2014 states

Since the dynamics is unknown, we first need to choose a model. As always, a linear model is the natural choice, at least at the initial stage of development.

The problem is that (i) of course, we must choose an empirical model, but we will try to avoid those that the data blatantly reject, (ii) $\mathbf{B}$ being diagonal has nothing to do with linearity, and (iii) there is no further potential "development" possible without this assumption, which, in effect, assumes the causal problem away. As a result, whenever the diagonal $\mathbf{B}$ is violated by the data, IFs – as currently used in empirical studies – provide spurious causality.

We now use a very popular framework, from macroeconomics, to think more clearly about $\mathbf{B}$. Liang2008’s flow of simplifications and assumptions makes his once sophisticated (1) collapse to that of a bivariate VAR with 1 lag

The seemingly innocuous assumption is much less so within a statistical framework relating assumptions directly to observable quantities. Indeed, the uncorrelatedness of $W$’s, which here translates to that of $\epsilon$, combined with $b_{21}=b_{12}=0$ have a very stark implication. Let the number of endogenous variables be $M$ and the number of lags $P$. Provided imposing $M=2$ and $P=1$ is reasonable (4) is valid if and only if regression residuals $\left[\begin{smallmatrix}u_{1t}\\ u_{2t}\end{smallmatrix}\right]=\left[\begin{smallmatrix}b_{11}&b_{12}\\ b_{21}&b_{22}\end{smallmatrix}\right]\left[\begin{smallmatrix}\varepsilon_{1t}\\ \varepsilon_{2t}\end{smallmatrix}\right]$ are not cross-correlated. This is easily testable: one needs to estimate equations 1 and 2 separately by least squares, collect the residuals and calculate their correlation $\rho_{u}$. If the latter is different from 0 (and this could be formally tested with a t-test), then Liang2008’s simple formula does not apply. While $\rho_{u}=0$ might be plausible in continuous time (or anything near it), this is a monumental stretch for data sampled at the yearly frequency (like in SMCGL). Fortunately, unlike true structural causality, the data can directly inform us on whether Liang2014’s formula is valid or not for specific pairs of time series.

This section showcases how IF can lead to a pretense of causal knowledge, with conclusions that are sometimes the exact opposite of the truth. We consider four data-generating processes (DGPs) where the true $P$ is one. The first two correspond to what is proposed in Liang2014. The last two are generic VAR positively autocorrelated processes with differing degrees of persistence. Following the notation above (equation (5)):

As highlighted in the previous section, IFs are calculated assuming $b_{ij}=0$. However, as will be reported and formally tested in section 4, this is frequently not the case for most time series, especially those sampled at low frequencies. Using a controlled simulation environment, we study how IFs behave for values of $\rho_{u}:=corr(u_{t}^{X_{1}},u_{t}^{X_{2}})\in[-1,1]$.⁴⁴4Variances of $u_{t}^{X_{1}}$ and $u_{t}^{X_{2}}$ are one. Note that IFs are invariant to $\rho_{u}\neq 0$ emerging from $b_{12}=0$ or $b_{21}=0$, which is obviously problematic given the causal content of $b_{ij}$.⁵⁵5It is important to note that while we consider cases where either $b_{12}=0$ or $b_{21}=0$, there is a continuum of possibilities between those. We do so for simplicity of exposition (it makes the problem dichotomous). Moreover, setting either $b_{12}=0$ or $b_{21}=0$ to 0 corresponds to a causal ordering which is by far the most common identification scheme used in practice – and what we will be using in section 4.

As our indicator of true underlying information flow, taking into account $\rho_{u}\neq 0$, we use the FEVD-based measure described earlier. Note that here, unlike the application to real data in section 4, we know what are the true $b_{ji}$’s. When $\rho_{u}\neq 0$, the correlation must be attributed to either $b_{12}$ or $b_{21}$, or a combination of both. In a bivariate setup, this amounts to setting $b_{ij}=0$, and attributing $\rho\neq 0$ to $b_{ji}$. Hence, it is possible to tell when standard IFs conclude falsehoods, because in simulations $b_{12}$ and $b_{21}$ are known.

In terms of notation, $\Upsilon^{i,j}_{i\rightarrow j}$ means the FEVD share at horizon $h=10$ with $i=1,2$ and $i\neq j$. The superscripts in $\Upsilon^{i,j}_{i\rightarrow j}$ determine the true ordering, hence $i$ ordered before $j$. The subscripts indicate that we plot the contribution of variable $i$ to the forecast error variance of variable $j$ at horizon $h=10$. The simulations have shown that $h=10$ is sufficient for convergence.

Figure 1 plots the absolute normalized information flows (NIF), $\tau_{i\rightarrow j}$ for DGP(1) through DGP(4), with varying $\rho_{u}$. A first observation is that $\tau_{i\rightarrow j}$ often varies significantly with $\rho_{u}$, even when the very formula underlying it assumes $\rho_{u}=0$. The relative strength of $\tau_{1\rightarrow 2}$ and $\tau_{2\rightarrow 1}$ can easily collapse to 0 or be much higher than what IFs report for $\rho_{u}=0$.

The fact that IFs vary with $\rho_{u}$ could give false hopes that they account for simultaneous relationships. We conduct a simple exercise to show it does not. An interesting empirical question is whether $i$ is causing $j$ more than $j$ is causing $i$, which is the appealing promise of IFs. This moves beyond testing for Granger Causality and aims at quantifying the flow of information. The shaded green region corresponds to the values of $\rho_{u}$ for which either the true underlying causality is that $X_{2}$ causes $X_{1}$ more than the reverse, irrespective of whether it emerges from either $b_{12}=0$ or $b_{12}=0$. For those values of $\rho_{u}$, when the true ordering is unknown, the qualitative conclusion about the sign of the net causality flow does not hinge on knowledge of $b_{ij}$. Analogously, the blue shade represents values of $\rho_{u}$ for which $X_{1}$ causing more $X_{2}$ is unanimous among $b_{ij}$ configurations. White regions are values of $\rho_{u}$ where conclusions about the sign of net causality flow cannot be determined from the data, i.e., the unknown $b_{ij}$ is necessary to settle. In other words, given the data available to the modeler, $X_{2}$ causing more $X_{1}$ and vice versa are both equally likely, and sorting it out decisively implies making a successful guess on the true value of $b_{ij}$.

If IFs were correctly calibrated, the green dotted line should only be above the blue in the shaded green regions, and the dotted blue line above the green one in blue regions. Figure 1 clearly demonstrates this not to be the case. For every DGP, there is a substantial range of values (white spaces) for which IFs clearly conclude $\tau_{i\rightarrow j}>\tau_{j\rightarrow i}$ when there is one pair of $(b_{ij},b_{ji})$ out of two for which the reverse is true. To address concerns about an horizon mismatch $h$ between IFs and FEVD, Figure 4 in Appendix B shows results for $h=2$. As it turns out, things worsen for IFs (the blue and green regions shrink) since the simultaneity problem is less diluted in dynamics at short horizons.

In sum, those simulations (largely inspired by Liang2015 DGPs) show two things. First, the tractable IFs (from formula (4)) are functions of $\rho_{u}$ even though they assume it to be zero. This compromises any statement on the strength of causal links. Second, for any DGP, there is an underlying pair $\left(b_{ij},b_{ji}\right)$ for which IFs’ conclusion about the net causal flow is the opposite of reality. This is not only an artifact of large $|\rho_{u}|$, as exemplified by DGP(4).

In this section, we first revisit SMCGL’s application of IFs to GMTA and forcings. Second, we look at what lies behind opaque FEVD measurements by reporting impulse response functions, and computing the transient climate response implied by simplistic bivariate VARs.

Global warming generated by man-made forcing is the prevalent generalization of the notion climate change. Numerous researchers have dedicated their works to this very relationship between anthropogenic forcing and constituents of global climate (Hansen2006; AndrewsEtAl2010; LiEtAl2013; notz2016observed) and the socio-economic consequences thereof (Nord14; ipcc2019). Despite overwhelming evidence for anthropogenic forcing being the main driver of global climate change (ipccGlobWarm2018), scientists have also observed that especially since the turn of the millennium, global temperature has plateaued despite ever-rising greenhouse gases and contrary to projections from key climate models.⁶⁶6Resolving this puzzle has led to reevaluate the role of oceans in the interplay of radiative forcing and the climatic response (Tollefson14; Marotzke2015).

IFs hinge exactly on these bivariate relationships, which are attractive in their clarity, but are certainly a stark oversimplification of a complex system. Nevertheless, to give empirical content to our critique of the methodology, we study the relationship between GMTA and 7 forcings from SMCGL using both IFs and our FEVD-based remedy. The sample of annual means ranges from 1850 - 2005.⁷⁷7Data on the Pacific Decadal Oscillation (PDO) ranges from 1900 - 2005. We follow SMCGL and take our data from therein referenced data providers. Table 1 summarizes the results.⁸⁸8The sample was restricted to 1850-2005 to match that of SMCGL. Table 5 reports results extending the sample to 2017.

In Table 1 we report estimated correlation between residuals of the bivariate VAR(1) implied by IFs. When forcing $P=1$, in 5 cases out of 7, the null that $\hat{\rho}_{u}^{i,GMTA}=0$ is rejected at least at the 10% level. When choosing $P$ with Bayesian Information Criterion (BIC), only $\hat{\rho}_{u}^{\text{Solar},GMTA}=0$ cannot be rejected. Hence, as repeatedly mentioned in the text, IFs assume something that can be, and is rejected by the data. Moreover, in the light of simulations carried earlier, the qualitative and quantitative insights from IFs are often spurious under such conditions.

Naturally, we concentrate on FEVD results which correctly account for $\rho_{u}\neq 0$. However, setting $\rho_{u}=0$ is only one of the empirical shortcomings of the empirical IF formula — it also sets $P$, the number of lags of each $X_{t}$, to be 1. Clearly, it is empirically plausible that $X_{i,t-2}$ or $X_{j,t-4}$ may have an impact on $X_{i,t}$ beyond what is channeled by single year lags ($X_{i,t-1}$ and $X_{j,t-1}$) . In other words, $P=1$ is extremely restrictive on climatic dynamics. When choosing $P$ with BIC, results align better with prior scientific knowledge. Nonetheless, for the sake of completeness, we report both results ($P=1$ and $P=P^{*}$, with $P^{*}$ being BIC’s choice).

With $P=P^{*}$, the fact that total forcing causes GMTA more than the reverse is without appeal. Nevertheless, the quantitative answer is, again, highly dependent on the ordering choice. After 15 years, total forcing anomalies are responsible for explaining between 28.0% and 47.4% of that of GMTA – depending on the preferred ordering. The net causal flow being higher from aerosol and solar to GMTA are also unanimous, but much smaller. Indecisive results are reported for Anthropogenic, Volcanic and PDO. Overall, Table 1 suggests that the data itself does not support the strong qualitative conclusions of SMCGL for their CO₂ measure and Anthropogenic.

When $P$ is forced to one, as in IFs, inconclusive results are reported for total forcing, aerosol, volcanic and PDO. $P=1$ specifications, irrespective of the ordering⁹⁹9There are other identification schemes which cannot be cast as ”orderings”. That is, there are rotations of $\bm{u}_{t}$ (even when $\rho_{u}=0$) giving different structural shocks. Hence, while the two orderings span a lot of possibilities (and those traditionally considered first in practice), they do represent the universe of rotations of $\bm{u}_{t}$ into $\bm{\epsilon}_{t}$., conclude that GMTA is causing more $CO_{2}$ and Anthropogenic than the reverse. Only solar forcing results are unanimous, and in line with what climatic common wisdom suggests. Note that it also agrees with results from original IFs, which is not surprising given that $\hat{\rho}_{Solar,GMTA}$ is in the close vicinity of 0. However, choosing a proper $P$ nearly doubles the share of forecast errors attributable to solar forcing.

This note is a cautionary tale about how seemingly innocuous simplifying assumptions can go wrong – especially when they are formulated without consulting the data. IFs, as proposed by SMCGL, is a concept that hinges on the assumption of zero correlation between the residuals of a bivariate VAR(1) process. In discrete time, especially with observations at lower frequencies, such an assumption is most often not justified. Both stylized simulations and an empirical application in the form of the transient climate response demonstrate that being negligent about cross-correlated residuals can lead to markedly different outcomes. Our results show that already in a bivariate system of CO₂ and GMTA, the resulting TCR depends on how one deals with correlated residuals.

FEVDs provide an alternative to IFs. Albeit being a decisive improvement over IFs, FEVDs are as good as their underlying statistical model. To avoid departing too much from the SMCGL’s framework, we only considered bivariate models. Climate systems obviously comprise numerous additional variables. wilson2010atmospheric considers four in a VAR setup, but there could be many more. Additionally, the dynamics were assumed to be linear and time invariant, an approximation that should be eventually tested. estrada2013statistically trend "breaks" model is one way to do it and they report results pointing in the same direction as ours. However, with machine learning tools becoming increasingly common use, more flexible alternatives could be used to yield further insights. Finally, we considered simplistic identification schemes which implied a causal ordering of variables. There is a plethora of more sophisticated schemes available (kilian2017svar) and those could be used in future work, especially when moving beyond bivariate systems. Another – simpler – avenue is the use of data sampled at higher frequencies (like daily) which, by construction, makes the simultaneity problem much less of a Damocles sword.