Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
65,691 characters · 10 sections · 67 citation commands
On Spurious Causality, CO$_2$, and Global Temperature
\qquad \qquad Maximilian G\"obel}
\newgeometry{left=1.8cm, right = 1.8cm, bottom = 2.0cm, top = 2.0cm} \thispagestyle{empty}
Keywords: information flows, vector autoregressions, global warming, climate econometrics
\restoregeometry
\pagenumbering{arabic}
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 \href{https://www.nature.com/articles/srep21691/metrics}{popular} article, \textcolor{PineGreen}{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.\footnote{In 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$_2$, 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$_2$ forcing on GMTA but (ii) are in the range of recent estimates OttoEtAl2013,MontamatStock2020 if one assumes simultaneous causality running from CO$_2$ forcing to GMTA.
Section (ref) 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 (ref) display IFs possibly spurious behavior using simulated data where the true causality is known. Section (ref) revisit the question of causality between different forcings and GMTA using appropriate tools. Section (ref) 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 $\textbf{X} = \left(X_1,X_2\right) \in \mathbb{R}^2$ are the state variables and $\textbf{F} = \left(F_1,F_2\right)$. $\textbf{W} = \left(W_1,W_2\right)$ is a standard Wiener process, with $\Delta \textbf{W} = \frac{d\textbf{W}}{dt}$ and $E\left(\Delta W_i\right) = 0$ and $E\left(\Delta W_i\right)^2 = \Delta t$. $\textbf{B}$ is a $2 \times 2$ matrix and its entries $b_{ij}$ govern how perturbations instantaneously impact the system. At this point, the only non-innocuous assumption is that of a bivariate system. The information flow from $X_2$ to $X_1$, $T_{2 \rightarrow 1}$, is then defined as
where $\frac{dH_1}{dt}$ is the evolution of the marginal entropy of $X_1$ and $\frac{dH_{1\bcancel{2}}}{dt}$ denotes $\frac{dH_1}{dt}$ where the spillovers from $X_2$ are excluded. After lengthy derivations and reasonable assumptions (like that $W$'s, soon to be called structural shocks, are uncorrelated) it is obtained that
where $\rho$ is the joint probability density of variables $X_i$ and $X_j$, and $\rho_i$ denotes the marginal density of series $X_i$. Given that the model in (ref) is not readily identified from the data, this is not yet operational. After further derivations and assuming $\mathbf{B} = \left[
\right]$, Liang2014 transforms (ref) into a workable formula made of empirical moments
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
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 (ref) 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 (ref) is valid if and only if regression residuals $\left[
\right] = \left[
\right] \left[
\right]$ are \textit{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 \cite{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 exposition follows closely VARCTIC. In time series analysis, the "identification" problem originates from simultaneity in the data. We can learn whether $$X_{t-1} \rightarrow Y_{t} \quad \text{or} \quad Y_{t-1} \rightarrow X_{t}$$ is more plausible. This is predictive causality in the sense of granger1969. However, the data itself cannot discriminate between $$X_{t} \rightarrow Y_{t} \quad \text{and} \quad X_{t} \leftarrow Y_{t}.$$ In words, a correlation between $X_t$ and $Y_t$ can be generated by two different causal structures. Liang2008's solution is to assume such relationships do not exist --- yet, they do. Within a VAR, the problem boils down to the need for identifying $C$ in
where $\boldsymbol{y}_t$ is an $M$ by 1 vector -- meaning the dynamic system incorporates $M$ variables. $\Psi_p$'s parameterizes how each of these variables is predicted by its own lags and lags of the $M-1$ remaining variables. $P$ is the number of lags being included. The matrix $C$ characterizes how the $M$ different variables interact contemporaneously --- e.g., how total forcing affects GMTA within the same year (a time unit $t$ in SMCGL). Finally, the structural disturbances are mutually uncorrelated with mean zero: \[ \boldsymbol{\varepsilon}_t = \left[ \varepsilon_{1,t},\enskip ... \enskip , \varepsilon_{M,t} \right] ~ \sim ~ \, N \left (0 , ~ I_M \right ). \] Equation ((ref)) is the so-called structural form of the VAR, which cannot be estimated because $C$ is not identified by the data. SMCGL uses formula (ref) which implicitly assumes a constrained version of (ref) with $M=2$, $P=1$, and, most importantly, $C$ being a diagonal matrix. The validity of their analysis hinges upon those constraints not being rejected by the data. In section (ref), we find that the data disagrees with at least two of them.
Equation (ref) is a structural model which can be used to answer causal questions directly. However, the elements of $C$ are not plain regression coefficients, and cannot be estimated as such --- they would be biased. It does not mean that they do not exist. The implications of their existence can best be understood by looking at an estimable "reduced-form" VAR
where $\boldsymbol{c} = C^{-1} \Psi_0$ and $\Phi_p = C^{-1} \Psi_p$ are both regression coefficients obtained by running least squares separately on each equation. $\boldsymbol{u}_t$ are now regression residuals \[ \boldsymbol{u}_t = \left[ u_{1,t},\enskip ... \enskip , u_{M,t} \right] ~ \sim ~ \, N \left (0 , ~ \Sigma_u \right ) \] which will be cross-correlated if the true $C$ is not diagonal. As mentioned earlier, Liang2014's simplifying assumptions translate into $\Sigma_u=C^{-1'}C^{-1}$ being diagonal which is often at odds with the data. All the parameters of (ref) can be estimated with traditional methods, but the model is not "structural" and cannot be used for causal inference. Structural VARs, which aim at uncovering "structural" causality (instead of predictive causality à la granger1969) acknowledge $\Sigma_u$ being non-diagonal and provide ways to obtain $C$. As a byproduct, they can procure valid measures of information flow.
The raw material of causal measures are $\boldsymbol{\varepsilon}_t$, the structural disturbances entering the systems. However, those, like structural causality, are not directly extractable from the data: we only have $\boldsymbol{u}_t$ and translating those back to $\boldsymbol{\varepsilon}_t$ necessitates $C$. The latter is not directly attainable, but can be retrieved using the mapping $\Sigma_u=C^{-1'}C^{-1}$. In words, this means the covariance matrix of regression residuals from (ref) can be used as raw material to retrieve the "structural" $C$. Mechanically, the identification problem emerges because there are many $C$'s satisfying $\hat{\Sigma}_u=C^{-1'}C^{-1}$ --- numerous causal structures deliver the same residuals' cross-correlations.
The strategy we opt for is the traditional Choleski decomposition of $\hat{\Sigma}_u$. This is one of many identification strategies for the VAR kilian2017svar. But among the catalog of methodologies, the Choleski decomposition is certainly popular (if not the most popular) in applied work and is simple to implement. Mechanically, it provides a lower-triangular matrix $C$, satisfying $\hat{\Sigma}_u=\mathbf{B}'\mathbf{B}$ where $\mathbf{B}$ is the same as from equation (ref), but with dimensions $M \times M$. Its purpose is to transform cross-correlated regressions residuals $u_t$ (equation (ref)) into uncorrelated structural shocks $\boldsymbol{\varepsilon}_t$ (equation (ref)). This is done by reversing the relationship $\boldsymbol{u}_t=C\boldsymbol{\varepsilon}_t$.
Uncorrelatedness is essential to study how GMTA responds to a given forcing, keeping everything else constant. Such a causal claim would be impossible when considering an impulse from correlated residuals $u_t$ as those always co-move. A Choleski decompostion of $\Sigma_u$ is one way to transform the observed $u_t$ into the unobserved fundamental shocks $\boldsymbol{\varepsilon}_t$. The assumption underlying such an approach to orthogonalization is a causal ordering of shocks. The ordering restricts how variables interact with each other within the same year, conditional on the previous state of the system. In our bivariate setup, ordering forcing $j$ after GMTA implies that forcing cannot impact GMTA within the same year. Ordering GMTA after a given forcing implies that GMTA cannot impact forcing $j$ within the same year. SMCGL implicitly assumes both restrictions at the same time, which results in a rejected overidentified model. In contrast, when the model is just identified (when only one restriction is imposed), this choice cannot be validated by the data itself as it does not alter the likelihood.
When revisiting SMCGL's empirical work, we consider both orderings and document how sensitive conclusions are to that necessary choice.\footnote{Of course, there are identification schemes outside of the family of "orderings" obtained by Choleski decomposition, but those are beyond the scope of this paper -- and unnecessary to make our main point.} There are cases where the sign of net IF (a qualitative notion) between $j$ and GMTA depends on the ordering choice, and cases where it does not. For instance, we will find that whatever is assumed about $C$ in a proper VAR system, total forcing appears to be causing GMTA much more than the reverse. In other cases, like $CO_2$-induced radiative forcing, this cannot be simply ruled out by the data.
In a complete multivariate system like a VAR, the errors of the $h$-steps ahead forecast $y_{t+h,m}$ can be related back to structural shocks -- that is, the anomalies driving the dynamics of the system. For instance, we can compute the share of the forecast error variance of GMTA 10 years from now that is attributable to $CO_2$ anomalies. Intuitively, if $CO_2$ is causing GMTA to increase, its exogenous impulses should be an important driver of GMTA's variance rather than GMTA anomalies themselves -- a high information flow/transfer between the two variables. Accordingly, VAR forecast errors are $$ \boldsymbol{u}_{t+h} = \boldsymbol{y}_{t+h}-\hat{\boldsymbol{y}}_{t+h}=\sum_{h'=0}^{h-1} \Theta_{h'} \boldsymbol{\varepsilon}_{t+h-h'} $$ where $\Theta_{h}$ is a function of $\Phi_p$'s and $C$. See kilian2017svar for further details. The forecast error variance decomposition (FEVD) of the whole system at horizon $h$ can be analytically calculated using the entries of matrix $\Theta_{h'}$. Precisely, $$
$$ An information flow "share" of $i$ for $j$ can be characterized by the share of forecast error variance of variable $j$ attributable to structural shocks of $i$.\footnote{For a discussion on how to think about "shocks" in a physical system, see \cite{VARCTIC}.} While those measures can be assessed for any horizon $h$ (which can contain useful information), we focus on the cumulative sum since it is a measure of total flow. For non-stationary VARs (like those of the empirical section) we use a horizon of $h=15$ years. In the case of stationary VAR process, FEVD measures quickly converge to their long-run value as $h$ increases. Moreover, in that case, our FEVD-based IF measures can be more directly motivated from the Wold representation of a VAR process \citep{kilian2017svar}. Finally, note that whenever the $b_{12}=b_{21}=0$ assumption is approximately true, the FEVD approach and IFs give qualitatively similar answers.
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 (ref)):
As highlighted in the previous section, IFs are calculated assuming $b_{ij}=0$. However, as will be reported and formally tested in section (ref), 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]$.\footnote{Variances 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}$.\footnote{It 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 (ref).}
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 (ref), 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 (ref) 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 (ref) 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 (ref) in Appendix (ref) 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 (ref)) 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.\footnote{Resolving 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.\footnote{Data on the Pacific Decadal Oscillation (PDO) ranges from 1900 - 2005.} We follow SMCGL and take our data from therein referenced data providers. Table (ref) summarizes the results.\footnote{The sample was restricted to 1850-2005 to match that of SMCGL. Table (ref) reports results extending the sample to 2017.}
In Table (ref) 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 (ref) suggests that the data itself does not support the strong qualitative conclusions of SMCGL for their CO$_2$ 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\footnote{There are other identification schemes which cannot be cast as "orderings". That is, there are rotations of $\boldsymbol{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 $\boldsymbol{u}_t$ into $\boldsymbol{\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.
Results for CO$_2$ are rather surprising. Irrespective of the ordering, GMTA is reported to cause CO$_2$ more than the reverse, a finding contradicting SMCGL's results and common wisdom.\footnote{Hansen2006 states that global warming started to accelerate not prior to the 1970s. Only after 1975 did the global temperature increase by approximately 0.2$^\circ$C per decade.} However, SMCGL's CO$_2$ measure exhibit ill dynamic behavior that cannot possibly be that of a natural quantity. This becomes obvious when \href{https://drive.google.com/file/d/1yhaJi92dvY_Lax0H5BFIzJeNc517Gn44/view?usp=sharing}{plotting} their CO$_2$-ERF measure in first difference: it evolves according to a series of dichotomous jumps. This CO$_2$-ERF series, in which ERF stands for effective radiative forcing, is a concept presented in MyhreEtAl2013. The series itself is based on EtminanEtAl2016.
Given the strange jumping behavior of the CO$_2$-ERF series, and to make our calculations comparable to other findings in the literature BrunsEtAl2020,MontamatStock2020, we derive RF from CO$_2$ concentration as follows: we use the well established Meinshausen2017 data set on annual global means of CO$_2$ concentration, measured in parts per million (ppm). We follow MyhreEtAl2013 and transform the increase in CO$_2$ concentration in year $t$ measured in ppm relative to the concentration in a given base year, CO$_{2,base}$, into radiative forcing, $RF^{\text{CO}_2}_t$ -- measured in W/m$^\text{2}$ -- as follows: $RF^{\text{CO}_2}_t = 5.35 \times \text{ln}\left(\text{CO}_{2,t} / \text{CO}_{2,base} \right)$. We use 1850 as base year following BrunsEtAl2020.
As it turns out, considering this less contentious CO$_2$ series does not resolve the apparently counterintuitive finding that GMTA explains a larger portion of the forecast error variance of CO$_2$ than vice versa. Such a finding has also been reported using a different methodology in koutsoyiannis2020atmospheric. We explore a last avenue, that of using annual CO$_2$ emissions rather than $RF^{\text{CO}_2}_t$. This last attempt is successful in reconciling the FEVD approach with the traditional wisdom that CO$_2$ is causing GMTA "more" than the reverse. This finding is independent of the ordering choice.
In sum, based on this particular time series evidence, the causal link between certain forcings and GMTA remains disputable. What is less disputable are the effects of total forcing, CO$_2$ emissions, and solar forcing which all explain an important share of GMTA anomalies independently of arbitrary ordering preferences. Nonetheless, we see such analyses as rather primitive and potentially misleading. For instance, $X$ causing ”more” $Y$ does not mean that the reverse causality is not quantitatively important, or climatologically relevant. We now turn to a more promising way to extract meaning out of selected bivariate VARs.
To open the black box of those rather abstract measurements, we report in Figure (ref) impulse response functions (IRF) for the bivariate models of CO$_2$-GMTA and total forcing-GMTA. Since sims1980, the dominant approach for studying the properties of the VAR around its deterministic path has been IRFs to structural shocks. Their dynamic effect can be analyzed as that of a randomly assigned treatment because those have been transformed to be uncorrelated, which provides the "keeping everything else constant" interpretation.
The IRF of a variable $i$ to a one standard deviation shock of $\varepsilon_{j,t}$ is defined as
Thus, it is the expected difference, $h$ months after "impact", between a bivariate system that responded to an unexpected CO$_2$ increase, and the same system where no such increase occurred. In a linear VAR, the above takes a closed-form solution in terms of the matrices from (ref). Models are now estimated with Bayesian methods, optimizing the hyperparameters of a standard Minnesota prior, and choosing the number of lags as reported in the gray-shaded rows of Table (ref). The primary motivation is to obtain valid inference even in the presence of nonstationarity. Point estimates are nearly identical to that of OLS. For details on such choices in the context of climate data and a more thorough (yet introductory) treatment of IRFs, see VARCTIC.
In Figure (ref) we show the effect of an unexpected increase in annual emissions on GMTA, while Figure (ref) shows the response of GMTA radiative forcing generated by an unexpected rise in cumulative emissions.
Qualitatively, the impact of annual emissions and cumulative CO$_2$-induced forcing shocks on global temperature is vastly similar. In accord with findings in VARCTIC for the effect of CO$_2$ on Arctic sea ice extent, the impact of total forcing and CO$_2$ shocks is highly durable. This time, it is on GMTA rather than sea ice extent. As reflected in Figures (ref) and (ref), this result is qualitatively independent of the ordering choice. In both cases (and for both forcing variables), the effect of forcing takes about two years to completely settle in. However, it is clear that the reported short-run impact strongly depends on the identification assumptions, which SMCGL completely abstract from.
Whether the slightly negative short-run response of GMTA in the first panel of Figure (ref) favors the ordering $\left\lbrace \text{\textit{GMTA}}, \,\text{\textit{CO$_2$}} \right\rbrace$ over $\left\lbrace \text{\textit{CO$_2$}}, \, \text{\textit{GMTA}} \right\rbrace$ is debatable: ForsterEtAl2020 find the reduction in global CO$_2$ emissions during the COVID-19 pandemic to have resulted in a short-run rise of global temperature. The key mechanism is a decline in the cooling-effect of aerosols as a result of less SO$_2$ emissions. The authors project a rise in global temperature over the first 24 months following the pandemic-induced reduction in global nitrogen oxide (NO$_x$) emissions.\footnote{Especially NO$_2$ is found to be well-correlated with CO$_2$ emissions ForsterEtAl2020.}
Lastly, a robustness check. Our variables are clearly nonstationary. While this does not pose a problem for Bayesian inference, it is nevertheless natural to wonder if results would be significantly altered by including a time trend. Accordingly, the second column reports the same IRFs, but for VAR specifications augmented with trends. Those show that adding such an exogenous explanatory variable does not change the dynamics of a GMTA response to an unexpected shock to CO$_2$. The addition of a trend to the bivariate model of GMTA and total forcing, allows GMTA to slowly revert to a lower impact -- which is nevertheless highly persistent.
The transient climate response (TCR) is a frequently used metric to measure the impact of rising atmospheric CO$_2$ concentration on temperature. It is not only an indication of the trajectory of ongoing climate change, but also serves as a benchmark to evaluate the results of climate model projections PhillipsLeirvikStorelvmo2020. The TCR is defined as the increase in temperature, between $h_0$ and $h_T$, under the assumption that CO$_2$ increases annually by 1%. $h_T$ is defined as that point in time when -- due to the steady annual increase of 1% -- CO$_2$ concentration is twice as high as at date $h_0$ Pretis2020,MontamatStock2020. Such a doubling of CO$_2$ would occur approximately after 70 years OttoEtAl2013. Following the transformation of ppm to W/m$^\text{2}$ as suggested by MyhreEtAl2013 and MontamatStock2020, a doubling of CO$_2$ under an annual increase in concentration of 1% would generate a radiative forcing of $5.35 \times \text{ln}\left(2\right) \approx 3.7 \text{ W/m}^2$.\footnote{{An annual increase of 1% in atmospheric CO$_2$ concentration results in a doubling of CO$_2$ after approximately 70 years, which is described more formally as: $h \times ln\left(\frac{1.01}{1}\right) = ln\left(2\right) \; ,$ for $h \approx 70$ MontamatStock2020. }}
Typical estimates for TCR fall within a range of 1$^{\circ}$C-2.5$^{\circ}$C with a 66% probability, as summarized in the IPCC 5th Assessment Report BindoffEtAl2013. More recent estimates are well aligned with this range. BrunsEtAl2020 reports a point estimate of TCR ranging from 1.17$^{\circ}$C to 1.85$^{\circ}$C, depending on the type of data and model specification. Pretis2020 embeds a two-component energy balance model into a cointegrated vector autoregressive model. His estimates vary across model specification and range from 1.24$^{\circ}$C to 1.38$^{\circ}$C. PhillipsLeirvikStorelvmo2020 report a global transient climate sensitivity of 2.05$^{\circ}$C. The IV regression in MontamatStock2020 allows for a differentiation of TCR measurements across different horizons, normalized to giving 70-year-horizon estimates. Their point estimates range in the neighborhood of 1.5$^{\circ}$C within a 95% confidence interval of roughly 0.9$^{\circ}$C to 2.1$^{\circ}$C.
Despite being much more simplistic than the models of the aforementioned works on TCR, the bivariate VARs do also allow to estimate the impact of a doubling of CO$_2$ on temperature. Here we make use of the concept of IRFs, as presented in Equation (ref). In particular, we estimate the impact on temperature when $RF^{\text{CO}_2}$ increases by one standard deviation of its reduced-form residuals, $\sigma_{\varepsilon_{j}}$, where $j = RF^{\text{CO}_2}$, from a bivariate VAR of $RF^{\text{CO}_2}$ and GMTA. Recalling the definition of TCR, a doubling of CO$_2$ concentration, which is achieved by an annual increase of 1% in atmospheric CO$_2$ concentration, generates an additional radiative forcing of approximately $5.35 \times \text{ln}\left(2\right) \approx 3.7 \text{ W/m}^2$. In our case, the shock $\sigma_{\varepsilon_{j}}$ to CO$_2$ is a one-time event at horizon $h = 0$, but its effects are distributed over horizons $h = 1, 2, 3, ..., H$. This allows us to measure the cumulative increase, $\Xi_i$, in $i \in \left\lbrace RF^{\text{CO}_2}, GMTA \right\rbrace$ generated by $\sigma_{\varepsilon_{j}}$, where $j = RF^{\text{CO}_2}$, at any horizon $h$: \[ \Xi_{j,h} = \sum^h_{s=0} IRF(j\rightarrow i, s) \enskip , \] where $IRF(j\rightarrow i, s)$ is defined as in Equation (ref). Adapting the formula of OttoEtAl2013, we estimate TCR$_h$ as the increase in GMTA at horizon $h$ as follows:
where $\Xi_{GMTA,h}$ is the cumulative increase in global temperature at horizon $h$, resulting from the shock $\sigma_{\varepsilon_{j}}$, where $j = RF^{\text{CO}_2}$, at horizon $h = 0$. Likewise, $\Xi_{RF^{\text{CO}_2},h}$ is the cumulative increase in radiative forcing of CO$_2$ at horizon $h$, resulting from the shock $\sigma_{\varepsilon_{j}}$, where $j = RF^{\text{CO}_2}$, at horizon $h = 0$.
In Table (ref) we present median point estimates of TCR$_h$ for $h=20$ and $h=70$ (as in MontamatStock2020) from bivariate VARs of $RF^{\text{CO}_2}$ and GMTA -- with and without an exogenous time trend. Thus, Table (ref) reports the TCR corresponding to the model specifications in Figure (ref) (a) and (b). That is, we use Bayesian estimation techniques and deploy a Minnesota prior on our parameter estimates. Our VAR has four lags and the estimation is based on annual observations between 1850 and 2005.
The main message of Table (ref) is twofold: first, the ordering of the variables heavily influences the final results for $\text{TCR}_{20}$, demonstrating the importance of respecting the possibility of cross-correlated residuals. {Mechanically, the discordance brought up by the ordering choice vanishes at much longer horizons, and, as a result, $\text{TCR}_{70}$ estimates are largely similar. However, at that horizon, it is the choice of whether or not to include a trend that can alter results significantly.} Second, even though the TCR$_h$ point estimates of the trend models are rather located at the upper bound of the IPCC range of 1$^{\circ}$C-2.5$^{\circ}$C, a simplistic bivariate VAR model including a constant and a time trend as additional exogenous regressors, is capable of providing a reasonable approximation of the rise in global mean temperature, triggered by a doubling of atmospheric CO$_2$ concentration.
TCR estimates can be helpful in choosing which ordering is most plausible. For instance, only by ordering CO$_2$ first do we get $\text{TCR}_{20}$ to fall within the IPCC range. IRFs can also help sort things out. Ordering CO$_2$ second leads to a surprisingly lasting negative effect of CO$_2$ shocks on GMTA. An increasingly popular approach to VAR identification in macroeconomics is to use sign restrictions, where implausible IRF draws (based on economic theory) are tossed out Uhlig2005. This dispenses the researcher from formulating a likely contentious causal ordering of variables. In a climate application, one could identify the VAR by rejecting specifications generating implausible $IRF(\text{CO$_2$}\rightarrow \text{GMTA}, h)$ or TCRs. Applying this sort of reasoning leads us to favor -- to nobody's surprise -- the specification where simultaneous causality runs from CO$_2$ to GMTA. However, it is important to stress that this choice is obtained from prior knowledge on what is deemed reasonable and what is not, rather than our two time series.
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$_2$ 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.