Next Article in Journal
Forecasting FOMC Forecasts
Next Article in Special Issue
Climate Finance: Mapping Air Pollution and Finance Market in Time Series
Previous Article in Journal
Prais–Winsten Algorithm for Regression with Second or Higher Order Autoregressive Errors
Previous Article in Special Issue
Monitoring Cointegrating Polynomial Regressions: Theory and Application to the Environmental Kuznets Curves for Carbon and Sulfur Dioxide Emissions
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

On Spurious Causality, CO2, and Global Temperature

1
Département des Sciences Économiques, Université du Québec à Montréal, Montréal, QC H2L 2C4, Canada
2
Lisbon School of Economics and Management, Universidade de Lisboa, 1200-781 Lisboa, Portugal
*
Author to whom correspondence should be addressed.
Submission received: 30 July 2021 / Revised: 23 August 2021 / Accepted: 26 August 2021 / Published: 7 September 2021
(This article belongs to the Collection Econometric Analysis of Climate Change)

Abstract

:
Stips et al. (2016) use information flows (Liang (2008, 2014)) to establish causality from various forcings to global temperature. We show that the formulas being used hinge on a simplifying assumption that is nearly always rejected by the data. We propose the well-known forecast error variance decomposition based on a Vector Autoregression as an adequate measure of information flow, and find that most results in Stips et al. (2016) cannot be corroborated. Then, we discuss which modeling choices (e.g., the choice of CO 2 series and assumptions about simultaneous relationships) may help in extracting credible estimates of causal flows and the transient climate response simply by looking at the joint dynamics of two climatic time series.

1. Introduction

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 are 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 (often) by 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 are available. To what extent do different forcings cause global mean surface temperature anomalies (GMTA) Stips et al. (2016)? 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 Liang (2008, 2014, 2015, 2016). The authors make grand claims about the technique being able to extract “true” and “rigorous” causality, which clashes with the usual 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 the Liang (2008) formulas that are used throughout assume that a certain statistical quantity is zero when it is not. This can be tested and it is rejected almost all the time. More precisely, the Liang (2014) 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 the last year’s values of both. That is unlikely to happen, and the data will testify to that. The assumption is presented as “linearity” (Liang 2014), 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 (Sims 1980). The problem is that for discretely-sampled time series, the assumption is almost always rejected by lower-frequency data.1 In short, SMCGL assume there to be no identification problem (in the sense of Sims (1980)), 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 Liang (2014)) 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 themselves 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 accordance with the scientific consensus can be obtained when using concentration directly rather than its radiative content.
The alternative methodology that we use is 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 Sims (1980) and has since then entered various fields, ranging from neuroscience (Chen et al. 2011) to climate (Goulet Coulombe Göbel 2021). 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 (Montamat and Stock 2020; Otto et al. 2013) if one assumes simultaneous causality running from CO 2 forcing to GMTA.
The remainder of this article is structured as follows: Section 2 briefly reviews IFs and explains where the problematic assumption occurs. It also discusses relevant notions of Vector Autoregressions (VARs, Sims 1980) as a comprehensive framework to think about causality in multivariate time series systems. Section 3 displays IFs possibly spurious behavior using simulated data where the true causality is known. Section 4 revisits the question of causality between different forcings and GMTA using appropriate tools. Section 5 concludes the paper.

2. Information Flows, Non-Innocuous Assumptions and VARs

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. Liang (2008) initiate the derivation of IF by considering the following data generating process (DGP):
d X = F X , t d t + B x , t d W ,
where X = X 1 , X 2 R 2 are the state variables and F = F 1 , F 2 . W = W 1 , W 2 is a standard Wiener process, with Δ W = d W d t and E Δ W i = 0 and E Δ W i 2 = Δ t . B is a 2 × 2 matrix and its entries b i j 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 1 , is then defined as
T 2 1 = d H 1 d t d H 1 2 d t
where d H 1 d t is the evolution of the marginal entropy of X 1 and d H 1 2 d t denotes d H 1 d t where the spillovers from X 2 are excluded. After lengthy derivations and reasonable assumptions, such as Euler–Bernstein approximations, Taylor-series expansions or that the components of W —which are soon to be called structural shocks—are uncorrelated, it is obtained that
T j i = E 1 ρ i F i ρ i X i + 1 2 E 1 ρ i 2 g i i ρ i X i 2 ,
where ρ is the joint probability density of variables X i and X j , and ρ i denotes the marginal density of series X i . Given that the model in (1) is not readily identified from the data, this is not yet operational. After further derivations and assuming B = b i i 0 0 b j j , Liang (2014) transforms (3) into a workable formula made of empirical moments
T j i = σ i , i σ i , j σ j , Δ i σ i , j 2 σ i , Δ i σ i , i 2 σ j , j σ i , i σ i , j 2
where σ i , i is the variance of i, σ i j , the covariance between i and j, and σ i , Δ i is the covariance between i and the k t h difference of i (we follow Liang (2015) and set k = 1 ). However, the validity of this appealing formula rests on the seemingly technical assumption of a diagonal B . Our point is that this assumption is far from merely technical and very frequently wrong. Motivating the diagonal B among other assumptions, Liang (2014) states
Since the dynamics are 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) 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 B is violated by the data, IFs—as currently used in empirical studies—provide spurious causality.
We now use a very popular framework from macroeconometrics to think more clearly about B . Liang (2008)’s flow of simplifications and assumptions make his once sophisticated (1) collapse to that of a bivariate VAR with one lag
X 1 , t X 2 , t = c 1 c 2 c + a 11 1 a 12 1 a 21 1 a 22 1 A X 1 , t 1 X 2 , t 1 + b 11 b 12 b 21 b 22 B ε 1 t ε 2 t .
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 ϵ , 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 u 1 t u 2 t = b 11 b 12 b 21 b 22 ε 1 t ε 2 t 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 ρ u . If the latter is different from zero (and this could be formally tested with a t-test), then Liang (2008)’s simple formula does not apply. While ρ 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 the true structural causality, the data can directly inform us on whether Liang (2014)’s formula is valid or not for specific pairs of time series.

2.1. Acknowledging the Identification Problem: Vector Autoregressions

This exposition closely follows Goulet Coulombe and Göbel (2021). In time series analysis, the "identification" problem originates from simultaneity in the data. We can learn whether
X t 1 Y t or Y t 1 X t
is more plausible. This is predictive causality in the sense of Granger (1969). However, the data themselves cannot discriminate between
X t Y t and X t Y t .
In essence, a correlation between X t and Y t can be generated by two different causal structures. Liang (2008)’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
C y t = Ψ 0 + p = 1 P Ψ p y t p + ε t ,
where y t is an M by one vector—meaning the dynamic system incorporates M variables. Ψ 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:
ε t = ε 1 , t , , ε M , t N 0 , I M .
Equation (6) is the so-called structural form of the VAR, which cannot be estimated because C is not identified by the data. SMCGL uses formula (4) that implicitly assumes a constrained version of (6) 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 4, we find that the data disagrees with at least two of them.
Equation (6) is a structural model that 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
y t = c + p = 1 P Φ p y t p + u t ,
where c = C 1 Ψ 0 and Φ p = C 1 Ψ p are both regression coefficients obtained by running least squares separately on each equation. u t are now regression residuals
u t = u 1 , t , , u M , t N 0 , Σ u
which will be cross-correlated if the true C is not diagonal. As mentioned earlier, Liang (2014)’s simplifying assumptions translate into Σ u = C 1 C 1 being diagonal, which is often at odds with the data. All of the parameters of (7) 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 Granger 1969) acknowledge Σ u being non-diagonal and provide ways to obtain C. As a by-product, they can procure valid measures of information flow.
The raw material of causal measures are ε t , the structural disturbances entering the systems. However, those like structural causality are not directly extractable from the data: we only have u t and translating those back to ε t necessitates C. The latter is not directly attainable, but can be retrieved using the mapping Σ u = C 1 C 1 . In words, this means that the covariance matrix of regression residuals from (7) can be used as raw material to retrieve the “structural” C. Mechanically, the identification problem emerges because there are many C’s satisfying Σ ^ 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 Σ ^ u . This is one of many identification strategies for the VAR (Kilian and Lütkepohl 2017). However, 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 Σ ^ u = B B where B is the same as from Equation (5), but with dimensions M × M . Its purpose is to transform cross-correlated regression residuals u t (Equation (7)) into uncorrelated structural shocks ε t (Equation (6)). This is done by reversing the relationship u t = C ε 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 decomposition of Σ u is one way to transform the observed u t into the unobserved fundamental shocks ε 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 over-identified 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.2 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.

2.2. An Adequate Measure of IF Based on the VAR

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
u t + h = y t + h y ^ t + h = h = 0 h 1 Θ h ε t + h h
where Θ h is a function of Φ p ’s and C. See Kilian and Lütkepohl (2017) 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 Θ h . Precisely,
MSPE h = E ( u t + h u t + h ) = h = 0 h 1 Θ h Θ h .
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.3 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 (Kilian and Lütkepohl 2017). Finally, note that whenever the b 12 = b 21 = 0 assumption is approximately true, the FEVD approach and IFs give qualitatively similar answers.

3. Simulations

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 Liang (2014). The last two are generic VAR positively-autocorrelated processes with differing degrees of persistence. Following the notation above (Equation (5)):
D G P ( 1 ) : = { A = 0.5 0.5 0 0.6 , c = 0.1 0.7
D G P ( 2 ) : = { A = 0.5 0.9 0.2 0.5 , c = 0 0
D G P ( 3 ) : = { A = 0.5 0.2 0.5 0.25 , c = 0 0
D G P ( 4 ) : = { A = 0.25 0.1 0.2 0.1 , c = 0 0
As highlighted in the previous section, IFs are calculated assuming b i j = 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 ρ u : = c o r r ( u t X 1 , u t X 2 ) [ 1 , 1 ] .4 Note that IFs are invariant to ρ u 0 emerging from b 12 = 0 or b 21 = 0 , which is obviously problematic given the causal content of b i j .5
As our indicator of true underlying information flow, taking into account ρ u 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 j i ’s. When ρ u 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 i j = 0 and attributing ρ 0 to b j i . 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, Υ i j i , j means the FEVD share at horizon h = 10 with i = 1 , 2 and i j . The superscripts in Υ i j i , 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), τ i j for DGP(1) through DGP(4), with varying ρ u . A first observation is that τ i j often varies significantly with ρ u , even when the very formula underlying it assumes ρ u = 0 . The relative strength of τ 1 2 and τ 2 1 can easily collapse to 0 or be much higher than what IFs report for ρ u = 0 .
The fact that IFs vary with ρ u could give false hopes that they account for simultaneous relationships. We conduct a simple exercise to show they do 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 ρ 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 ρ 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 i j . Analogously, the blue shade represents values of ρ u for which X 1 causing more X 2 is unanimous among b i j configurations. White regions are values of ρ u where conclusions about the sign of net causality flow cannot be determined from the data, i.e., the unknown b i j 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 i j .
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 τ i j > τ j i when there is one pair of ( b i j , b j i ) out of two for which the reverse is true. To address concerns about a horizon mismatch h between IFs and FEVD, Figure A1 (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 Liang (2015) DGPs) show two things. First, the tractable IFs (from formula (4)) are functions of ρ 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 b i j , b j i for which the IFs’ conclusion about the net causal flow is the opposite of reality. This is not only an artifact of large | ρ u | , as exemplified by DGP(4).

4. GMTA and Radiative Forcing Revisited

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 (Andrews et al. 2010; Hansen et al. 2006; Li et al. 2013; Notz and Stroeve 2016) and the socio-economic consequences thereof (Nordhaus 2014; Shukla et al. 2019). Despite overwhelming evidence for anthropogenic forcing being the main driver of global climate change (Masson-Delmotte et al. 2018), 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.6
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 seven forcings from SMCGL using both IFs and our FEVD-based remedy. The sample of annual means ranges from 1850 to 2005.7 We follow SMCGL and take our data from therein referenced data providers. Table 1 summarizes the results.8
In Table 1, we report the estimated correlation between residuals of the bivariate VAR(1) implied by IFs. When forcing P = 1 , in five cases out of seven, the null that ρ ^ u i , G M T A = 0 is rejected at least at the 10% level. When choosing P with Bayesian Information Criterion (BIC), only ρ ^ u Solar , G M T A = 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 out earlier, the qualitative and quantitative insights from IFs are often spurious under such conditions.
Naturally, we concentrate on FEVD results, which correctly account for ρ u 0 . However, setting ρ 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 as one. 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 is also unanimous, but much smaller. Indecisive results are reported for Anthropogenic, Volcanic, and PDO. Overall, Table 1 suggests that the data themselves do 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 ordering9, it can be concluded 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 ρ ^ S o l a r , G M T A is in the close vicinity of zero. However, choosing a proper P nearly doubles the share of forecast errors attributable to solar forcing.

4.1. What’s Up with CO 2 ?

Results for CO 2 are rather surprising. Irrespective of the ordering, GMTA is reported to cause more CO 2 than the reverse, a finding contradicting SMCGL’s results and common wisdom.10 However, SMCGL’s CO 2 measure exhibits ill dynamic behavior that cannot possibly be that of a natural quantity. This becomes obvious when plotting their CO 2 -ERF measure in first difference: it evolves according to a series of dichotomous jumps.11 This CO 2 -ERF series, in which ERF stands for effective radiative forcing, is a concept presented in Myhre et al. (2013). The series itself is based on Etminan et al. (2016).
Given the strange jumping behavior of the CO 2 -ERF series, and to make our calculations comparable to other findings in the literature (Bruns et al. 2020; Montamat and Stock 2020), we derive RF from CO 2 concentration as follows: we use the well-established Meinshausen et al. (2017) data set on annual global means of CO 2 concentration, measured in parts per million (ppm). We follow Myhre et al. (2013) 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 , b a s e , into radiative forcing, R F t CO 2 —measured in W/m 2 —as follows: R F t CO 2 = 5.35 × ln CO 2 , t / CO 2 , b a s e . We use 1850 as the base year following Bruns et al. (2020).
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 Koutsoyiannis and Kundzewicz (2020). We explore a last avenue, that of using annual CO 2 emissions rather than R F t CO 2 . 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 are 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.

4.2. Impulse Response Functions

To open the black box of those rather abstract measurements, we report in Figure 3 impulse response functions (IRF) for the bivariate models of CO 2 -GMTA and total forcing-GMTA. Since Sims (1980), 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 ε j , t is defined as
I R F ( j i , h ) = E ( y i , t + h | y t , ε t , j = σ ε j ) E ( y i , t + h | y t , ε t , j = 0 ) .
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 (6). 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 1. The primary motivation is to obtain a valid inference even in the presence of non-stationarity. 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 Goulet Coulombe and Göbel (2021).
In Figure 2, we show the effect of an unexpected increase in annual emissions on GMTA, while Figure 3 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 accordance with findings in Goulet Coulombe and Göbel (2021) 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 Figure 2 and Figure 3, 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 2 favors the ordering GMTA , CO 2 over CO 2 , GMTA is debatable: Forster et al. (2020) 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.12
Lastly, we conduct a robustness check. Our variables are clearly non-stationary. While this does not pose a problem for Bayesian inference, it is nevertheless natural to wonder if the 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.

4.3. Are VAR Estimates Quantitatively Reasonable?

The transient climate response (TCR) is a frequently used metric of 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 (Phillips et al. 2020). 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 (Montamat and Stock 2020; Pretis 2020). Such a doubling of CO 2 would occur approximately after 70 years (Otto et al. 2013). Following the transformation of ppm to W/m 2 as suggested by Myhre et al. (2013) and Montamat and Stock (2020), a doubling of CO 2 under an annual increase in concentration of 1% would generate a radiative forcing of 5.35 × ln 2 3.7 W / m 2 .13
Typical estimates for TCR fall within a range of 1 ° C–2.5 ° C with a 66% probability, as summarized in the IPCC 5th Assessment Report (Bindoff et al. 2013). More recent estimates are well aligned with this range. Bruns et al. (2020) report a point estimate of TCR ranging from 1.17 ° C to 1.85 ° C, depending on the type of data and model specification. Pretis (2020) embeds a two-component energy balance model into a cointegrated vector autoregressive model. His estimates vary across model specification and range from 1.24 ° C to 1.38 ° C. Phillips et al. (2020) report a global transient climate sensitivity of 2.05 ° C. The IV regression in Montamat and Stock (2020) 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 ° C within a 95% confidence interval of roughly 0.9 ° C to 2.1 ° C.
Despite being much more simplistic than the models of the aforementioned works on TCR, bivariate VARs also allow for the estimation of the impact of a doubling of CO 2 on temperature. Here we make use of the concept of IRFs, as presented in Equation (9). In particular, we estimate the impact on temperature when R F CO 2 increases by one standard deviation of its reduced-form residuals, σ ε j , where j = R F CO 2 , from a bivariate VAR of R F 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 × ln 2 3.7 W / m 2 . In our case, the shock σ ε 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, Ξ i , in i R F CO 2 , G M T A generated by σ ε j , where j = R F CO 2 , at any horizon h:
Ξ j , h = s = 0 h I R F ( j i , s ) ,
where I R F ( j i , s ) is defined as in Equation (9). Adapting the formula of Otto et al. (2013), we estimate TCR h as the increase in GMTA at horizon h as follows:
T C R h = Ξ G M T A , h × 5.35 × ln 2 Ξ R F CO 2 , h ,
where Ξ G M T A , h is the cumulative increase in global temperature at horizon h, resulting from the shock σ ε j , where j = R F CO 2 , at horizon h = 0 . Likewise, Ξ R F CO 2 , h is the cumulative increase in radiative forcing of CO 2 at horizon h, resulting from the shock σ ε j , where j = R F CO 2 , at horizon h = 0 .
In Table 2, we present median point estimates of TCR h for h = 20 and h = 70 (as in Montamat and Stock (2020)) from bivariate VARs of R F CO 2 and GMTA—with and without an exogenous time trend. Thus, Table 2 reports the TCR corresponding to the model specifications in Figure 3a,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. Results for 1850–2017 are available in Table A3 (Appendix C).
The main message of Table 2 is two-fold: first, the ordering of the variables heavily influences the final results for 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, 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 ° C–2.5 ° 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 achieve TCR 20 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 (Uhlig 2005). 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 I R F ( CO 2 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.

5. Conclusions

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, are 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.14 Nevertheless, provided reasonable modeling choices, the practical implications of our VAR results are in line with previous works (Richardson et al. 2016; Schurer et al. 2018), though when including an exogenous time-trend, he estimates range at the upper bound of estimates found in Otto et al. (2013), Montamat and Stock (2020), or Forster et al. (2021). That is, for a doubling of atmospheric CO 2 concentration our estimates suggest GMTA to increase by 1.99 ° C (2.17 ° C with trend) after 20 years and 2.06 ° C (2.58 ° C with trend) after 70 years.
FEVDs provide an alternative to IFs. Although 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 of numerous additional variables. Granville Tunnicliffe (2015) 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. Estrada et al. (2013)’s trend “breaks” model is one way to do it and they report results pointing in the same direction as ours. However, with use of machine learning tools becoming increasingly common, more flexible alternatives could be used to yield further insights—like those developed for macroeconomic time series in Goulet Coulombe (2020). Finally, we considered simplistic identification schemes that implied a causal ordering of variables. There is a plethora of more sophisticated schemes available (Kilian and Lütkepohl 2017) 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 (daily for example) which, by construction, makes the simultaneity problem much less of a Damocles sword.

Author Contributions

All authors contributed equally. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Data are either available from public sources, as indicated in Table A1, or from the authors’ websites.

Acknowledgments

We are grateful to Boyuan Zhang for many helpful comments and suggestions.

Conflicts of Interest

The authors declare no conflict of interest.

Appendix A. Data Sources

For the empirical part of Section 4, we refer to the variables presented in Table 1 of SMCGL. Due to the sparse description of the data sources, our data set reduces to the annual global mean surface temperature anomalies (GMTA) and seven other representatives of radiative forcing. For five variables we could match the reported correlation coefficients as well as the IF values of SMCGL, Table 1. For total forcing and solar these measures do not coincide. For the analysis in Section 4.2 and Section 4.3, we use the CO 2 series based on Meinshausen et al. (2017). We transform ppm into W/m 2 as described in Section 4.
Table A1. List of Variables.
Table A1. List of Variables.
AbbreviationDescriptionData Source
Total Forcingannual; 1850–2005KNMI Climate Explorer
Anthropogenicannual; 1850–2005KNMI Climate Explorer
CO 2 -ERF (W/m 2 )annual; 1850–2005KNMI Climate Explorer
Aerosolannual; 1850–2005KNMI Climate Explorer
Solarannual; 1850–2005KNMI Climate Explorer
Volcanicannual; 1850–2005KNMI Climate Explorer
PDOannual averages of monthly observations; 1900–2005KNMI Climate Explorer
GMTAannual; global; 1900–2005HadCRUT4
CO 2 (Million Tonnes/Year)annual; global production-based emissions; 1850–2005Our World in Data—not in use
CO 2 (ppm)annual; global
1850–2014 (Meinshausen et al. 2017):
2015–2017 (NOAA-ESRL):

IAC ETH Zürich
NOAA-ESRL

Appendix B. Additional Simulation Results

In Figure A1 below, we compare NIF and FEVD h = 2 . One might be concerned that results in Figure 1 suffer from a horizon mismatch between NIF and FEVD h = 10 . Results suggest that reducing h worsens NIFs’ problems by shrinking the “safe” regions. This is intuitive; the effect of assumptions on simultaneous relationships becomes milder as we move further from h = 0 .
Figure A1. Ranking of NIFs ( τ i j ) vs. Ranking of FEVDs ( Υ i j i , j ); Different Levels of Correlation; Horizon h = 2 . Notes: Color specification: the light blue area shows the region, in which FEVDs suggest Υ 1 2 i , j > Υ 2 1 i , j , regardless of the ordering. Similarly, the light green area shows the region, in which FEVDs suggest Υ 2 1 i , j > Υ 1 2 i , j —independent of the ordering. The white/non-colored areas shows those regions for which the ordering of X 1 and X 2 does matter. An unambiguous determination of Υ 1 2 i , j > Υ 2 1 i , j , respectively, Υ 2 1 i , j > Υ 1 2 i , j , is not possible.
Figure A1. Ranking of NIFs ( τ i j ) vs. Ranking of FEVDs ( Υ i j i , j ); Different Levels of Correlation; Horizon h = 2 . Notes: Color specification: the light blue area shows the region, in which FEVDs suggest Υ 1 2 i , j > Υ 2 1 i , j , regardless of the ordering. Similarly, the light green area shows the region, in which FEVDs suggest Υ 2 1 i , j > Υ 1 2 i , j —independent of the ordering. The white/non-colored areas shows those regions for which the ordering of X 1 and X 2 does matter. An unambiguous determination of Υ 1 2 i , j > Υ 2 1 i , j , respectively, Υ 2 1 i , j > Υ 1 2 i , j , is not possible.
Econometrics 09 00033 g0a1

Appendix C. Additional Empirical Results

In this section, we show results for the sample period to 1850–2017.
Table A2. Empirical Results for the Bivariate Relationship Between Various Forcings and GMTA Sample Period: 1850–2017.
Table A2. Empirical Results for the Bivariate Relationship Between Various Forcings and GMTA Sample Period: 1850–2017.
CorrelationNormalized IF
(IF × 100)
FEVD
Lags
(P)
Correlation of
Residuals ( ρ u )
Ordering: i, GMTAOrdering: GMTA, i
i→ GMTAGMTA → i i GMTA GMTA i i GMTA GMTA i
Total Forcing0.8236.829.540.23 ***48.015.628.329.3
(17.2)(15.3)10.29 ***55.714.430.036.6
Anthropogenic0.9143.7−18.94−0.16 **6.55.04.913.5
(39.9)(−0.6)1−0.15 **4.25.12.413.7
CO 2 —ERF (W/m 2 )
SMCGL
0.9143.5−13.04−0.105.759.25.315.4
(39.0)(−0.3)1−0.112.13.11.28.0
Aerosol−0.8437.9−45.74−0.17 **3.92.96.10.0
(19.4)(−1.3)1−0.002.033.32.033.1
Solar31.47.02.180.054.60.93.11.4
(1.1)(0.6)10.066.71.04.42.0
Volcanic0.111.3−0.340.18 **10.10.52.72.4
(0.2)(−0.2)10.20 ***7.10.30.63.7
PDO0.15−2.3−0.440.4 ***31.00.83.913.7
(−0.3)(−0.3)10.38 ***9.50.30.713.5
CO 2 (Mt/yr)0.8942.01.02−0.127.91.98.80.3
(31.0)(0.00)1−0.105.40.05.40.7
CO 2 (W/m 2 )0.9143.4−13.420.18 **5.016.16.16.2
(38.3)(−0.3)10.061.53.41.01.6
Notes: i corresponds to the type of radiative forcing, listed in the left most column. The second column ("Correlation") documents the correlation between GMTA and variable i. FEVD values are taken at horizon h = 15 , which translates into the contribution of variable i in the variance of the forecast error of variable j a decade and a half after the in-sample end date. Numbers in bold underline the highest absolute causal flow among a (i,GMTA) pair for a given measure. “*”, “**”, and “***” means that the null of the residuals cross-correlation of residuals is rejected at the 10%,5%, and 1% level, respectively. Sample period: 1850–2017.
Table A3. Transient Climate Response Sample Period: 1850–2017.
Table A3. Transient Climate Response Sample Period: 1850–2017.
Ordering Without Trend With Trend
TCR 20 TCR 70 TCR 20 TCR 70
CO 2 , GMTA 1.46 ° C1.94 ° C 1.76 ° C2.35 ° C
GMTA, CO 2 0.58 ° C1.79 ° C 0.97 ° C2.22 ° C

Notes

1
In Tawia Hagan et al. (2019), IFs 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.
2
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.
3
For a discussion on how to think about “shocks” in a physical system, see Goulet Coulombe and Göbel (2021).
4
Variances of u t X 1 and u t X 2 are one.
5
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 zero corresponds to a causal ordering that is by far the most common identification scheme used in practice, which happens to be what we will be using in Section 4.
6
Resolving this puzzle has led to re-evaluating the role of oceans in the interplay of radiative forcing and the climatic response (Marotzke and Forster 2015; Tollefson 2014).
7
Data on the Pacific Decadal Oscillation (PDO) range from 1900 to 2005.
8
The sample was restricted to 1850–2005 to match that of SMCGL. Table A2 reports results extending the sample to 2017.
9
There are other identification schemes that cannot be cast as “orderings”. That is, there are rotations of u t (even when ρ 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 u t into ϵ t .
10
Hansen et al. (2006) states that global warming did not start to accelerate prior to the 1970s. Only after 1975 did the global temperature increase by approximately 0.2 ° C per decade.
11
12
Especially NO 2 is found to be well-correlated with CO 2 emissions (Forster et al. 2020).
13
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 × l n 1.01 1 = l n 2 , for h 70 (Montamat and Stock 2020).
14
With the proliferation of new methods to extract information flows between time series (e.g., fractal regressions (Kristoufek and Ferreira 2018) or multiscale transfer entropy (Zhao et al. 2018)), there is much research to be done on dealing with the simultaneity problem within those heterogeneous frameworks.

References

  1. Andrews, T., P. M. Forster, O. Boucher, N. Bellouin, and A. Jones. 2010. Precipitation, Radiative Forcing and Global Temperature Change. Geophysical Research Letters 37. [Google Scholar] [CrossRef] [Green Version]
  2. Bindoff, N. L., P. A. Stott, P. A. AchutaRao, M. R. Allen, N. Gillett, D. Gutzler, K. Hansingo, G. Hegerl, Y. Hu, S. Jain, and et al. 2013. Chapter 10—Detection and Attribution of Climate Change: From Global to Regional. In Climate Change 2013: The Physical Science Basis. IPCC Working Group I Contribution to AR5. Cambridge: Cambridge University Press. [Google Scholar]
  3. Bruns, S. B., Z. Csereklyei, and D. I. Stern. 2020. A Multicointegration Model of Global Climate Change. Journal of Econometrics 214: 175–97, Annals Issue: Econometric Models of Climate Change. [Google Scholar] [CrossRef]
  4. Chen, Gang, Daniel R. Glen, Ziad S. Saad, J. Paul Hamilton, Moriah E. Thomason, Ian H. Gotlib, and Robert W. Cox. 2011. Vector autoregression, structural equation modeling, and their synthesis in neuroimaging data analysis. Computers in Biology and Medicine 41: 1142–55. [Google Scholar] [CrossRef] [Green Version]
  5. Estrada, Francisco, Pierre Perron, and Benjamín Martínez-López. 2013. Statistically derived contributions of diverse human influences to twentieth-century temperature changes. Nature Geoscience 6: 1050–5. [Google Scholar] [CrossRef] [Green Version]
  6. Etminan, M., G. Myhre, E. J. Highwood, and K. P. Shine. 2016. Radiative Forcing of Carbon Dioxide, Methane, and Nitrous Oxide: A significant Revision of the Methane Radiative Forcing. Geophysical Research Letters 43: 12614–23. [Google Scholar] [CrossRef]
  7. Forster, P. M., H. I. Forster, M. J. Evans, M. J. Gidden, C. D. Jones, C. A. Keller, R. D. Lamboll, C. L. Quéré, J. Rogelj, D. Rosen, and et al. 2020. Current and Future Global Climate Impacts Resulting from COVID-19. Nature Climate Change 10: 913–19. [Google Scholar] [CrossRef]
  8. Forster, P., T. Storelvmo, K. Armour, W. Collins, J. L. Dufresne, D. Frame, D. J. Lunt, T. Mauritsen, M. D. Palmer, M. Watanabe, and et al. 2021. The Earth’s Energy Budget, Climate Feedbacks, and Climate Sensitivity. In Climate Change 2021: The Physical Science Basis. Contribution of Working Group I to the Sixth Assessment Report of the Intergovernmental Panel on Climate Change. Edited by V. Masson-Delmotte, P. Zhai, A. Pirani, S. Connors, C. Péan, S. Berger, Y. Caud, N. Chen, L. Goldfarb, M. Gomis and et al. Cambridge: Cambridge University Press. [Google Scholar]
  9. Goulet Coulombe, Philippe. 2020. The Macroeconomy as a Random Forest. Available online: https://ssrn.com/abstract=3633110 (accessed on 31 August 2021).
  10. Goulet Coulombe, Philippe, and Maximilian Göbel. 2021. Arctic amplification of anthropogenic forcing: A vector autoregressive analysis. Journal of Climate. Forthcoming. [Google Scholar]
  11. Granger, Clive W. J. 1969. Investigating causal relations by econometric models and cross-spectral methods. Econometrica: Journal of the Econometric Society 37: 424–38. [Google Scholar] [CrossRef]
  12. Granville Tunnicliffe, Wilson. 2015. Atmospheric co2 and global temperatures: The strength and nature of their dependence. In Working Paper. Lancaster: Lancaster University. [Google Scholar]
  13. Hansen, J., M. Sato, R. Ruedy, K. Lo, D. W. Lea, and M. Medina-Elizade. 2006. Global Temperature Change. Proceedings of the National Academy of Sciences 103: 14288–93. [Google Scholar] [CrossRef] [Green Version]
  14. Kilian, Lutz, and Helmut Lütkepohl. 2017. Structural Vector Autoregressive Analysis. Cambridge: Cambridge University Press. [Google Scholar]
  15. Koutsoyiannis, Demetris, and Zbigniew W Kundzewicz. 2020. Atmospheric temperature and co2: Hen-or-egg causality? Sci 2: 83. [Google Scholar] [CrossRef]
  16. Kristoufek, Ladislav, and Paulo Ferreira. 2018. Capital asset pricing model in portugal: Evidence from fractal regressions. Portuguese Economic Journal 17: 173–83. [Google Scholar] [CrossRef] [Green Version]
  17. Li, C., D. Notz, S. Tietsche, and J. Marotzke. 2013. The Transient versus the Equilibrium Response of Sea Ice to Global Warming. Journal of Climate 26: 5624–36. [Google Scholar] [CrossRef] [Green Version]
  18. Liang, X. S. 2008. Information Flow within Stochastic Dynamical System. Phys. Rev. E 78: 031113. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  19. Liang, X. S. 2014. Unraveling the Cause-Effect Relation between Time Series. Physical Review E 90: 052150. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  20. Liang, X. S. 2015. Normalizing the Causality between Time Series. Physical Review E 92: 022126. [Google Scholar] [CrossRef] [Green Version]
  21. Liang, X. S. 2016. Information Flow and Causality as rigorous Notions ab initio. Physical Review E 94: 052201. [Google Scholar] [CrossRef] [Green Version]
  22. Marotzke, J., and P. M. Forster. 2015. Forcing, Feedback and Internal Variability in Global Temperature Trends. Nature 517: 565–70. [Google Scholar] [CrossRef] [Green Version]
  23. Masson-Delmotte, V., P. Zhai, H.-O. Pörtner, D. Roberts, J. Skea, P. R. Shukla, A. Pirani, W. Moufouma-Okia, C. Péan, R. Pidcock, and et al. 2018. IPCC, 2018: Global Warming of 1.5 °C. An IPCC Special Report on the Impacts of Global Warming of 1.5 °C above Pre-Industrial Levels and related Global Greenhouse Gas Emission Pathways, in the Context of Strengthening the Global Response to the Threat of Climate Change. Sustainable Development, and Efforts to Eradicate Poverty. Available online: https://www.ipcc.ch/2018/10/08/summary-for-policymakers-of-ipcc-special-report-on-global-warming-of-1-5c-approved-by-governments/ (accessed on 31 August 2021).
  24. Meinshausen, M., E. Vogel, A. Nauels, K. Lorbacher, N. Meinshausen, D. M. Etheridge, P. J. Fraser, S. A. Montzka, P. J. Rayner, C. M. Trudinger, and et al. 2017. Historical greenhouse gas concentrations for climate modelling (cmip6). Geoscientific Model Development 10: 2057–116. [Google Scholar] [CrossRef] [Green Version]
  25. Montamat, G., and J. H. Stock. 2020. Quasi-experimental estimates of the transient climate response using observational data. Climatic Change 160: 361–71. [Google Scholar] [CrossRef]
  26. Myhre, G., D. Shindell, F.-M. Bréon, W. Collins, J. Fuglestvedt, J. Huang, D. Koch, J.-F. Lamarque, D. Lee, B. Mendoza, and et al. 2013. Anthropogenic and Natural Radiative Forcing. In Climate Change 2013: The Physical Science Basis. Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change. Edited by T. Stocker, D. Qin, G.-K. Plattner, M. Tignor, S. K. Allen, J. Doschung, A. Nauels, Y. Xia, V. Bex and P. Midgley. Cambridge: Cambridge University Press, pp. 659–740. [Google Scholar] [CrossRef]
  27. Nordhaus, W. 2014. Estimates of the social cost of carbon: Concepts and results from the dice-2013r model and alternative approaches. Journal of the Association of Environmental and Resource Economists 1: 273–312. [Google Scholar] [CrossRef]
  28. Notz, D., and J. Stroeve. 2016. Observed arctic sea-ice loss directly follows anthropogenic co2 emission. Science 354: 747–50. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  29. Otto, A., F. E. L. Otto, O. Boucher, J. Church, G. Hegerl, P. M. Forster, N. P. Gillett, J. Gregory, G. C. Johnson, R. Knutti, and et al. 2013. Energy Budget Constraints on Climate Response. Nature Geoscience 6: 415–6. [Google Scholar] [CrossRef]
  30. Phillips, P. C. B., T. Leirvik, and T. Storelvmo. 2020. Econometric Estimates of Earth’s Transient Climate Sensitivity. Journal of Econometrics 214: 6–32. [Google Scholar] [CrossRef]
  31. Pretis, Felix. 2020. Econometric Modelling of Climate Systems: The Equivalence of Energy Balance Models and Cointegrated Vector Autoregressions. Journal of Econometrics 214: 256–73. [Google Scholar] [CrossRef]
  32. Richardson, M., K. Cowtan, E. Hawkins, and M. B. Stolpe. 2016. Reconciled climate response estimates from climate models and the energy budget of earth. Nature Climate Change 6: 931–5. [Google Scholar] [CrossRef] [Green Version]
  33. Schurer, A., G. Hegerl, A. Ribes, D. Polson, C. Morice, and S. Tett. 2018. Estimating the Transient Climate Response from Observed Warming. Journal of Climate 31: 8645–63. [Google Scholar] [CrossRef]
  34. Shukla, P. R., J. Skea, E. Calvo Buendia, V. Masson-Delmotte, H.-O. Pörtner, D. C. Roberts, P. Zhai, R. Slade, S. Connors, R. van Diemen, and et al. 2019. Climate Change and Land: An IPCC Special Report on Climate Change, Desertification, Land Degradation, Sustainable Land Management, Food Security, and Greenhouse Gas Fluxes in Terrestrial Ecosystems. Available online: https://spiral.imperial.ac.uk/bitstream/10044/1/76618/2/SRCCL-Full-Report-Compiled-191128.pdf (accessed on 6 September 2021).
  35. Sims, Christopher A. 1980. Macroeconomics and reality. Econometrica: Journal of the Econometric Society 48: 1–48. [Google Scholar] [CrossRef] [Green Version]
  36. Stips, A., D. Macias, C. Coughlan, E. Garcia-Gorriz, and X. S. Liang. 2016. On the Causal Structure between CO2 and Global Temperature. Scientific Reports 6: 21691. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  37. Tawia Hagan, Daniel Fiifi, Guojie Wang, X San Liang, and Han AJ Dolman. 2019. A time-varying causality formalism based on the liang–kleeman information flow for analyzing directed interactions in nonstationary climate systems. Journal of Climate 32: 7521–37. [Google Scholar] [CrossRef]
  38. Tollefson, J. 2014. Climate Change: The Case of the Missing Heat. Nature 505: 276–8. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  39. Uhlig, H. 2005. What are the effects of monetary policy on output? results from an agnostic identification procedure. Journal of Monetary Economics 52: 381–419. [Google Scholar] [CrossRef] [Green Version]
  40. Zhao, Xiaojun, Yupeng Sun, Xuemei Li, and Pengjian Shang. 2018. Multiscale transfer entropy: Measuring information transfer on multiple time scales. Communications in Nonlinear Science and Numerical Simulation 62: 202–12. [Google Scholar] [CrossRef]
Figure 1. Ranking of NIFs ( τ i j ) vs. Ranking of FEVDs ( Υ i j i , j ); Different Levels of Correlation; Horizon h = 10 . Notes: Color specification: the light blue area shows the region in which FEVDs suggest Υ 1 2 i , j > Υ 2 1 i , j , regardless of the ordering. Similarly, the light green area shows the region in which FEVDs suggest Υ 2 1 i , j > Υ 1 2 i , j —independent of the ordering. The white/non-colored areas show those regions for which the ordering of X 1 and X 2 does matter. An unambiguous determination of Υ 1 2 i , j > Υ 2 1 i , j , respectively, Υ 2 1 i , j > Υ 1 2 i , j , is not possible.
Figure 1. Ranking of NIFs ( τ i j ) vs. Ranking of FEVDs ( Υ i j i , j ); Different Levels of Correlation; Horizon h = 10 . Notes: Color specification: the light blue area shows the region in which FEVDs suggest Υ 1 2 i , j > Υ 2 1 i , j , regardless of the ordering. Similarly, the light green area shows the region in which FEVDs suggest Υ 2 1 i , j > Υ 1 2 i , j —independent of the ordering. The white/non-colored areas show those regions for which the ordering of X 1 and X 2 does matter. An unambiguous determination of Υ 1 2 i , j > Υ 2 1 i , j , respectively, Υ 2 1 i , j > Υ 1 2 i , j , is not possible.
Econometrics 09 00033 g001
Figure 2. IRFs: Annual Emissions and Global Temperature. Notes: We show impulse response functions from bivariate Vector Autoregressions of GMTA and annual CO 2 emissions. The right column includes a time trend as an additional exogenous regressor. The solid line is the median of 10,000 draws from the posterior distribution. Hyperparameters were optimized. The shaded area is the 68% credible region. Lags are those reported in Table 1 under the gray shaded rows.
Figure 2. IRFs: Annual Emissions and Global Temperature. Notes: We show impulse response functions from bivariate Vector Autoregressions of GMTA and annual CO 2 emissions. The right column includes a time trend as an additional exogenous regressor. The solid line is the median of 10,000 draws from the posterior distribution. Hyperparameters were optimized. The shaded area is the 68% credible region. Lags are those reported in Table 1 under the gray shaded rows.
Econometrics 09 00033 g002
Figure 3. IRFs: Cumulative Emissions and Global Temperature. Notes: We show impulse response functions from bivariate Vector Autoregressions of GMTA and CO 2 , and Total Forcing, respectively. The right column includes a time trend as an additional exogenous regressor. The solid line is the median of 10,000 draws from the posterior distribution. Hyperparameters were optimized. The shaded area is the 68% credible region. Lags are those reported in Table 1 under the gray shaded rows.
Figure 3. IRFs: Cumulative Emissions and Global Temperature. Notes: We show impulse response functions from bivariate Vector Autoregressions of GMTA and CO 2 , and Total Forcing, respectively. The right column includes a time trend as an additional exogenous regressor. The solid line is the median of 10,000 draws from the posterior distribution. Hyperparameters were optimized. The shaded area is the 68% credible region. Lags are those reported in Table 1 under the gray shaded rows.
Econometrics 09 00033 g003
Table 1. Empirical Results for the Bivariate Relationship—Between Various Forcings and GMTA.
Table 1. Empirical Results for the Bivariate Relationship—Between Various Forcings and GMTA.
CorrelationNormalized IF (IF × 100)FEVD
LagsCorrelation ofOrdering: i, GMTAOrdering: GMTA, i
i→ GMTAGMTA → i(P)Residuals ( ρ u ) i GMTA GMTA i i GMTA GMTA i
Total Forcing0.7330.620.840.23 ***47.413.028.025.4
(15.3)(11.1)10.29 ***51.49.627.628.7
Anthropogenic0.8639.8−20.04−0.19 **6.53.73.813.4
(35.7)(−0.6)1−0.19 **5.05.82.217.1
CO 2 -ERF (W/m 2 ) SMCGL0.8639.6−15.24−0.14 *6.58.45.617.4
(35.1)(−0.4)1−0.15 *2.84.71.112.8
Aerosol−0.8235.9−24.54−0.19 **2.90.62.11.6
(24.3)(−0.4)1−0.103.54.01.81.2
Solar0.4913.56.780.058.51.66.62.5
(3.8)(2.3)10.0816.64.212.36.8
Volcanic0.090.9−0.540.18 **10.90.83.13.6
(0.2)(−0.4)10.20 **7.11.40.63.7
PDO0.17−1.2−0.640.35 ***31.10.96.310.3
(−0.2)(−0.5)10.34 ***9.10.50.210.7
CO 2 (Mt/yr)0.8237.1−4.32−0.108.92.110.70.6
(27.0)(−0.0)1−0.054.20.04.40.4
CO 2 (W/m 2 )0.8639.5−14.040.23 ***5.216.82.94.7
(34.5)(−0.3)10.071.64.10.91.8
Notes: i corresponds to the type of radiative forcing, listed in the left most column. The second column (“Correlation”) documents the correlation between GMTA and variable i. FEVD values are taken at horizon h = 15 , which translates into the contribution of variable i in the variance of the forecast error of variable j a decade and a half after the in-sample end date. Numbers in bold underline the highest absolute causal flow among a (i,GMTA) pair for a given measure. “*”, “**”, and “***’ means that the null of the residuals cross-correlation of residuals is rejected at the 10%, 5%, and 1% level, respectively.
Table 2. Transient Climate Response.
Table 2. Transient Climate Response.
Ordering Without Trend With Trend
TCR 20 TCR 70 TCR 20 TCR 70
CO 2 , GMTA 1.99 ° C2.06 ° C 2.17 ° C2.58 ° C
GMTA, CO 2 0.57 ° C1.82 ° C 0.85 ° C2.39 ° C
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Goulet Coulombe, P.; Göbel, M. On Spurious Causality, CO2, and Global Temperature. Econometrics 2021, 9, 33. https://0-doi-org.brum.beds.ac.uk/10.3390/econometrics9030033

AMA Style

Goulet Coulombe P, Göbel M. On Spurious Causality, CO2, and Global Temperature. Econometrics. 2021; 9(3):33. https://0-doi-org.brum.beds.ac.uk/10.3390/econometrics9030033

Chicago/Turabian Style

Goulet Coulombe, Philippe, and Maximilian Göbel. 2021. "On Spurious Causality, CO2, and Global Temperature" Econometrics 9, no. 3: 33. https://0-doi-org.brum.beds.ac.uk/10.3390/econometrics9030033

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop