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.
54,324 characters · 11 sections · 47 citation commands
The Global Carbon Budget as a cointegrated system
} } }
\affil[1]{ \ Aarhus Center for Econometrics, Aarhus University} \affil[2]{ \ Center for Research in Energy, Aarhus University}
We present a cointegrated vector autoregressive (CVAR) model Johansen1995,Juselius2006 for the primary time series variables in the Global Carbon Budget (GCB) dataset provided by GCB2023. This model enables the statistical estimation of critical parameters of the global carbon cycle and the evaluation of their estimation uncertainty. The model integrates atmospheric carbon dioxide (CO$_2$) concentrations, anthropogenic emissions, and absorption by both the terrestrial biosphere (land sink) and the ocean (ocean sink). Central to the model is the global carbon budget equation: Changes in atmospheric concentrations are equal to the difference between anthropogenic emissions and uptake by the sinks. The model describes the sinks as functions of atmospheric CO$_2$ levels as well as the El Ni\ no-Southern Oscillation (ENSO) and Pacific Decadal Oscillation (PDO) cycles. The budget equation and the dependence of sinks on atmospheric concentrations jointly introduce simultaneity among the GCB variables. Emissions are modeled as a random walk with drift. All variables are trending, and the budget equation as well as the dependence of the sinks on concentrations constitute cointegrating relations.
The approach facilitates a data-driven examination of the global carbon cycle through a compact model incorporating both observational data and outputs from various large-scale Earth system models (ESMs). Historical GCB data are used for parameter estimation via maximum likelihood, and parameter uncertainty is assessed through statistical standard errors. Unlike ESMs or small-scale emulators, the CVAR approach quantifies uncertainty statistically. In contrast to earlier statistical studies of the GCB data set, the CVAR framework allows for formal hypothesis tests of the parametric restrictions motivated by the physical relations.
The Global Carbon Project\footnote{\url{https://www.globalcarbonproject.org}.} maintains an extensive database of time series variables describing the carbon cycle dynamics, aiding in understanding the transfer of anthropogenically emitted CO$_2$ to the atmosphere, oceans, and terrestrial biosphere. These data are updated annually and published in the Global Carbon Budget reports. Understanding carbon cycle dynamics is crucial for comprehending the overall climate system and climate change AR6chapter5.
The GCB data have been utilized for various statistical analyses. One research strand investigates whether the CO$_2$ absorption rate of the carbon sinks, measured by the airborne fraction or sink rate, is constant Raupach2008,Knorr2009,LeQuere2009,Gloor2010,Raupach2014,BHK2019,BHK2024. Other research uses the GCB data's residual, termed the budget imbalance, to verify the accuracy of reported CO$_2$ emissions by individual nations Peters2017,Bennedsen2021. Previous statistical analyses of the GCB data often limit model dimensions to univariate or bivariate settings and rarely consider all GCB variables together. Early studies, such as Enting1993 and Parkinson1998, evaluated parameter uncertainty in global carbon cycle models using statistical methodologies. BHK2023 specified a multivariate dynamic model of the GCB variables in a state space framework.
In this paper, we leverage the advantages of the CVAR framework that allows for specifying an unrestricted multivariate model that is agnostic about the physical relations of the variables and that nests the restricted, physically motivated model as a special case. Then we conduct a comprehensive (likelihood ratio) hypothesis test for whether the physical restrictions are supported by the data. Since we find that this is the case, we proceed with the restricted model. Earlier applications of the CVAR methodology to different aspects of climate research include Schmith2012, Pretis2020, and Castle2024.
The CVAR model can be used for both in-sample and out-of-sample analysis. In-sample, we demonstrate that the data from the Global Carbon Budget agree with the physical relationships. Out-of-sample, we explore projections of future paths of the global carbon cycle. We find that our projections align well with a common reduced-complexity climate model, in particular for a scenario with strong carbon cycle feedback effects.
The outline of the paper is as follows. In Section (ref), we explain the physical relations in the GCB data set and how they motivate a cointegration approach. We specify the restricted, physically motivated model. In Section (ref), we specify the unrestricted model by way of standard cointegration and specification tests. We then conduct a likelihood ratio test of the restricted against the unrestricted model, and we find that the data support the restricted, physically motivated model. We report the estimation results for the restricted model and discuss their physical meaning in the context of the system nature of the model. We also present residual diagnostics. In Section (ref), we apply the model to explore future projections of the paths of the system using scenarios from the Shared Socioeconomic Pathways (SSP) initiative Riahi2017. Section (ref) concludes. An appendix reports additional empirical results and conditional predictive evidence for the model, including validation and forecasting exercises.
Figure (ref) displays the time series data set from the Global Carbon Project studied in this paper. The GCB time series are annual from 1959 to 2022, measured in PgC (per year), and obtained from the global file of GCB2023.\footnote{Available at \url{https://globalcarbonbudgetdata.org}.} The top left panel in Figure (ref) shows the CO$_2$ uptake by the terrestrial biosphere (land sink), and the top right panel that by the oceans (ocean sink). The bottom left panel shows anthropogenic CO$_2$ emissions, calculated as the sum of fossil fuel and land-use change emissions minus the cement carbonation sink from the global file of GCB2023. The bottom right panel shows the “Keeling curve” of atmospheric CO$_2$ concentrations.
Figure (ref) shows that all time series are upwards trending. The source of trend in the land and ocean sink dynamics is atmospheric concentrations. BK1973 and gifford1993implications have proposed the most commonly used functions that relate atmospheric CO$_2$ concentrations and land uptake. BK1973 argue that the linear approximation to their function describes the relation well as long as concentrations do not deviate too far from pre-industrial concentrations. BHK2023 demonstrate on the GCB data set 1959--2020 that the linear approximation is accurate for both the land and ocean sinks. We therefore specify the sinks equations as
where $S^L_t$ and $S^O_t$ are land and ocean sinks, respectively, and $C_t$ is atmospheric concentrations. The deviation processes $X_{i,t} = \phi_i X_{i,t-1} + \varepsilon_{i,t}$, $i=1,2$, where $\varepsilon_{i,t}$ are mutually independent i.i.d.\ zero-mean variables, are specified as AR(1) time series. If $\phi_i = 0$, there is no autocorrelation in the deviations from the trend given by atmospheric concentrations.
BHK2023 show that anthropogenic emissions $E_t$ are well described by a random walk with drift $d>0$,
where $X_{3,t} = \phi_3 X_{3,t-1} + \varepsilon_{3,t}$ and $\varepsilon_{3,t}$ is an i.i.d.\ zero-mean variable. Atmospheric concentrations $C_t$ are given by the global carbon budget equation,
where $X_{4,t} = \phi_4 X_{4,t-1} + \varepsilon_{4,t}$ and $\varepsilon_{4,t}$ is an i.i.d.\ zero-mean variable. The budget equation states that changes in atmospheric concentrations, $\Delta C_t$, are given by anthropogenic emissions minus the combined land and ocean uptake plus an error, $X_{4,t}$, which corresponds to the budget imbalance variable in the global data file of GCB2023. The initial value of $C_t$ in 1959 is 672.87 PgC.
Figure (ref) suggests a possible COVID-19 effect in 2020, where, as a result of the global pandemic, economic activity and consequently emissions dropped substantially. We therefore include a COVID-19 dummy, $D2020_t = I_{\{t=2020\}}$, as an exogenous covariate in the emissions equation.
Large-scale climate variables have been found to explain part of the variation in the carbon sinks, and via those also some of the variation in concentrations feely1999influence,haverd2018new,BHK2023. These climate covariates are stationary, cyclical phenomena, while atmospheric concentrations are nonstationary. Concentrations are therefore by far the dominant explanatory variable in the sinks and central to the cointegrating relations. In this paper, we include two climate covariates, ENSO3.4 and PDO; see Figure (ref). ENSO3.4 is a standard measure of ENSO conditions based on sea-surface temperature anomalies in the central equatorial Pacific Trenberth1997; we use the monthly Ni\ no 3.4 series provided by the NOAA Physical Sciences Laboratory Nino34_PSL, computed from the HadISST dataset Rayner2003 on a 1981--2010 base period. The PDO is a pattern of decadal sea-surface temperature variability in the North Pacific Mantua1997; we use the monthly PDO index distributed by the NOAA National Centers for Environmental Information PDO_NCEI, derived from the ERSSTv5 dataset Huang2017. Since our analysis is at annual frequency, both indices are aggregated to calendar-year averages of the underlying monthly values.
For the land sink, El Ni\ no conditions are generally associated with reduced CO$_2$ uptake, whereas La Ni\ na conditions are associated with stronger uptake, consistent with the sensitivity of terrestrial carbon uptake to moisture and temperature conditions haverd2018new. For the ocean sink, La Ni\ na conditions are associated with reduced uptake because stronger upwelling brings carbon-rich waters to the surface, whereas El Ni\ no conditions weaken this mechanism feely1999influence. Hence, we expect the land sink to depend negatively on ENSO3.4 and the ocean sink positively. We include PDO to control for additional low-frequency climate variability beyond ENSO. Specifically, the PDO is a long-lived, ENSO-like mode of North Pacific climate variability, and the carbon-cycle literature attributes part of decadal ocean-sink variability to North Pacific processes. Since the net effects of PDO on global land and ocean uptake are less clear a priori, we do not have similarly strong prior expectations about the signs of its coefficients.
Because emissions are modeled as a random walk \hyperref[{E:emiss}]{\tagform@{\ref*{E:emiss}}}, and the deviation processes $X_{i,t}$, $i=1,\dots,4$, are all stationary, the model equations in \hyperref[{E:sinks}]{\tagform@{\ref*{E:sinks}}} and \hyperref[{E:GCB}]{\tagform@{\ref*{E:GCB}}} suggest that our system of $p=4$ variables is nonstationary I(1) and contains three linear cointegrating combinations that are stationary I(0). We model the variables collected in $Y_t=[S^L_t, S^O_t, E_t, C_t]'$ in the vector autoregressive (VAR) framework, and the analysis thus fits perfectly in the cointegrated VAR (CVAR) setup; see Johansen1995,Juselius2006. Within this systems framework, univariate properties are less relevant (if, e.g., one variable is stationary, then it will appear with a unit vector in the cointegration matrix). For this reason, we do not report univariate unit root or stationarity tests but refer to BHK2023.
The CVAR model for a $p$-dimensional vector time series $Y_t$ is given in vector error-correction model (VECM) form as
where the error term $U_t$ is assumed to be $p$-dimensional i.i.d.\ with mean zero and covariance matrix $\Sigma$. The long-run parameters $\alpha$ and $\beta$ are $p \times r$ matrices with $0 \leq r \leq p$. The rank $r$ is termed the cointegration rank, and the columns of $\beta$ constitute the $r$ cointegration vectors such that $\beta' Y_t$ are the stationary long-run equilibrium relations. The parameters in $\alpha$ are the adjustment parameters that represent the speed of adjustment towards equilibrium for each of the variables. To allow for a deterministic level (mean in the cointegrating relations) and a linear trend in the variables, we have included a so-called unrestricted constant term, $\mu$, and a restricted trend, $\alpha \rho t$. Finally, $Z_t$ allows inclusion of exogenous covariates, and the short-run dynamics are governed by the autoregressive parameters $\Gamma_1 , \ldots , \Gamma_k$ with lag-order denoted by $k$ (that is, there are $k+1$ lags in the VAR representation and equivalently $k$ lagged differences in the VECM representation).
The model specified in \hyperref[{E:sinks}]{\tagform@{\ref*{E:sinks}}} through \hyperref[{E:GCB}]{\tagform@{\ref*{E:GCB}}} can be written in structural VAR form as
where we have added $Z_t = [ENSO_t, PDO_t, D2020_t]'$ with the relevant coefficient matrix. Left-multiplying \hyperref[{E:VECM_w_SOI}]{\tagform@{\ref*{E:VECM_w_SOI}}} with the inverse of the leading matrix of contemporaneous relations, we obtain the VECM form \hyperref[{E:VECM_RF}]{\tagform@{\ref*{E:VECM_RF}}} with
where $c=1+b_1+b_2$, and the error vector is \[ U_t = \left[
\right] \left[
\right] = \left[
\right]. \] The model \hyperref[{E:VECM_w_SOI}]{\tagform@{\ref*{E:VECM_w_SOI}}} is equivalent to \hyperref[{E:VECM_RF}]{\tagform@{\ref*{E:VECM_RF}}} with coefficient matrices \hyperref[{E:VECM}]{\tagform@{\ref*{E:VECM}}}. The first row of $\beta'$ shows the cointegrating relation between $S^L_t$ and $C_t$, the second row the cointegrating relation between $S^O_t$ and $C_t$, and the third row the cointegrating relation of the carbon budget equation between $S^L_t$, $S^O_t$, and $E_t$. The adjustment coefficients in $\alpha$ show that $S^L_t$, $S^O_t$, and $C_t$ adjust to disequilibrium in all cointegrating relations, whereas $E_t$ does not adjust at all. This shows the special role of emissions as driver of the system, or, in other words, as the only source of deterministic and stochastic trend. It is apparent that the cointegrating rank is three, and the number of common stochastic trends therefore one, and given by the emissions.
The system \hyperref[{E:VECM_RF}]{\tagform@{\ref*{E:VECM_RF}}} implied by \hyperref[{E:VECM}]{\tagform@{\ref*{E:VECM}}} has $k=1$ lagged difference, and the coefficient $\Gamma_1$ of the first lag of differences, $\Delta Y_{t-1} = [\Delta S^L_{t-1}, \Delta S^O_{t-1}, \Delta E_{t-1}, \Delta C_{t-1}]'$, captures the short-run dynamics in stationary first differences of the system variables. It is apparent that autocorrelation in the dynamics of the deviations from the trend in the sinks $S^L_t$ and $S^O_t$, as captured by the autoregressive coefficients $\phi_1$ and $\phi_2$ in \hyperref[{E:sinks}]{\tagform@{\ref*{E:sinks}}}, does not play a role in the short-run dynamics of the cointegrated system. In fact, the first two columns of $\Gamma_1$ are zero. The reason can be understood from \hyperref[{E:VECM_w_SOI}]{\tagform@{\ref*{E:VECM_w_SOI}}}, where accounting for the contemporaneous relation $\Delta S^L_t-b_1\Delta C_t$ on the left-hand side removes all influence of the serial correlation parameter $\phi_1$ from the first differences and fully absorbs it in the dependence on the lagged level $S^L_{t-1}$. The same holds for $S^O_t$.
The global carbon budget equation, as a cointegrating relation, implies an influence of the climate covariates in $Z_t$ on $\Delta C_t$, by way of their influence on the sinks. This leads to a non-zero, but constrained, coefficient on $Z_t$ in the last row of $\Phi$. The error term $U_t$ is a linear combination of the structural errors $\varepsilon_t$, by way of the contemporaneous relations. Thus, even though the covariance matrix of $\varepsilon_t$ is diagonal, the covariance matrix of $U_t$ is not diagonal. However, the correlations are purely contemporaneous, and there is no serial correlation.
In this section, we first determine the lag order. Then we find a benchmark model that is unrestricted except in the cointegrating rank, which we determine using the Johansen trace tests. We then test the restrictions \hyperref[{E:VECM}]{\tagform@{\ref*{E:VECM}}} implied by the model \hyperref[{E:VECM_w_SOI}]{\tagform@{\ref*{E:VECM_w_SOI}}} against the unrestricted model, including exogeneity of emissions and the form of the cointegrating relations. Finally, we report estimation results and residual diagnostics for the restricted model \hyperref[{E:VECM_w_SOI}]{\tagform@{\ref*{E:VECM_w_SOI}}}.
We first determine the appropriate lag length. The model in \hyperref[{E:VECM_RF}]{\tagform@{\ref*{E:VECM_RF}}}--\hyperref[{E:VECM}]{\tagform@{\ref*{E:VECM}}} suggests $k=1$, but of course the unrestricted model \hyperref[{E:VECM_RF}]{\tagform@{\ref*{E:VECM_RF}}} could have any number of lags $k\geq 0$. Thus, we estimate the unrestricted model \hyperref[{E:VECM_RF}]{\tagform@{\ref*{E:VECM_RF}}} with full rank and no restrictions on the coefficient matrices for several choices of lag lengths. Because we want to allow for linear trends, we include a constant and a time trend in the unrestricted model. A summary of the results from unrestricted VAR models for $Y_t=[S^L_t, S^O_t, E_t, C_t]'$ are presented in Table (ref) with detailed single-equation residual diagnostics reported in Table (ref) in Appendix (ref).
The results in Table (ref) suggest that either $k=0$ or $k=1$ lags would be appropriate. Neither shows signs of misspecification. While the LR test for exclusion of lags suggests that $k=0$ lags could be sufficient, this non-rejection could be caused by the fact that there are many zero entries in the structural $\Gamma_1$ matrix; see \hyperref[{E:VECM}]{\tagform@{\ref*{E:VECM}}}. Furthermore, we prefer $k\geq 1$ such that the structural model is nested in the unrestricted model. Thus, we select $k=1$.
A priori, we expect three cointegrating relations to be present in $Y_t$ as given by \hyperref[{E:sinks}]{\tagform@{\ref*{E:sinks}}} and \hyperref[{E:GCB}]{\tagform@{\ref*{E:GCB}}}. Equivalently, emissions are the only source of stochastic trend in the system. Table (ref) confirms that the clear finding from our specification, without imposing any restrictions on the parameters, is a rank of three especially using bootstrap $P$-values for the trace test CRT2012. At this point, the statistics cannot tell whether emissions play any specific role.
Consequently, our chosen unrestricted benchmark model is the cointegrated VAR with one lagged difference, cointegration rank three, and including an unrestricted constant term and a restricted trend term, as well as the climate and dummy covariates $Z_t$ as exogenous explanatory variables. The VECM form of the benchmark model is thus given by \hyperref[{E:VECM_RF}]{\tagform@{\ref*{E:VECM_RF}}} with $p=4$, $k=1$, $r=3$, and $Z_t = [ENSO_t, PDO_t, D2020_t]'$, where $\alpha$ and $\beta$ are of dimension $p\times r = 4\times 3$ and $U_t$ is the $4\times 1$ reduced-form error vector.
Letting $\hat{U}_t$ denote the residuals, the maximized (quasi) log-likelihood function of the model is
where $T=63$ is the number of observations of $\Delta Y_t$, and \[ \hat{\Sigma} = \frac{1}{T-k} \sum_{t=k+1}^T \hat{U}_t\hat{U}_t' \] is the estimate of the covariance matrix of the reduced-form errors. Our estimated benchmark model has a numerically maximized log-likelihood of $-11.373$. Determining the degrees of freedom in the model \hyperref[{E:VECM_RF}]{\tagform@{\ref*{E:VECM_RF}}} accounts for the $p\times p$ matrix $\alpha \beta'$ of rank $r=3$ with $r(2p-r)=15$ free parameters. In addition, there is the intercept vector $\mu$ of dimension $p$, the trend vector $\rho$ of dimension $r$, and the short-run $p\times p$ coefficient matrix $\Gamma_1$, which together add 23 free parameters. Finally, the $p\times 3$ coefficient matrix on $Z_t$ adds 12 free parameters, for a total of 50 free parameters.
We first test for weak exogeneity and variable exclusion. Testing for weak exogeneity of each variable means testing whether the corresponding row in the $\alpha$ matrix is zero. In that case, the variable does not adjust to disequilibrium in any of the cointegrating relations. This could be either because it is not a relevant variable for the system or because it is a driver of the system. If the variable is not relevant for the cointegrating relations, then the corresponding row in the $\beta$ matrix is zero, i.e., variable exclusion from the cointegration space. In particular, if a variable cannot be excluded, and its $\alpha$ row is zero, then it is causing disequilibrium without adjusting to it, and it is therefore a driver of the system.
Table (ref) shows the results of LR tests for zero rows in $\beta$ (exclusion) and $\alpha$ (weak exogeneity). The tests for exclusion show that emissions $E_t$ cannot be excluded from the cointegrating relations in the unrestricted model \hyperref[{E:VECM_RF}]{\tagform@{\ref*{E:VECM_RF}}}, so emissions are important to explain the system dynamics. At the same time, the row in $\alpha$ corresponding to emissions is not distinguishable from zero (weak exogeneity), and therefore emissions do not adjust to disequilibrium. No other variable plays a similar role, so we conclude that emissions are the central driver of the system. Note that, this does not imply that emissions are exogenous with respect to the global economy, or to climate outcomes in a wider sense.
In Table (ref) in Appendix (ref) we present the coefficients in the common trends decomposition of the unrestricted model \hyperref[{E:VECM_RF}]{\tagform@{\ref*{E:VECM_RF}}}. The coefficient on emissions is normalized to one, and the remaining coefficients are much smaller in magnitude, which supports the interpretation of emissions as the common trend or driver of the system.
We also note that we cannot reject exclusion of atmospheric concentrations $C_t$ in the reduced-form unrestricted model. However, the inclusion of linear deterministic trends in the sinks in the unrestricted model apparently masks the necessity of $C_t$ in explaining the dynamics and the trending behavior in the sinks. This shows the importance of imposing parameter restrictions for capturing basic physical relations. Specifically, in the structural model \hyperref[{E:VECM_w_SOI}]{\tagform@{\ref*{E:VECM_w_SOI}}}, the intercept vector and coefficients $\Pi=\alpha\beta'$ and $\Gamma_1$ in \hyperref[{E:VECM}]{\tagform@{\ref*{E:VECM}}} are such that model \hyperref[{E:VECM_w_SOI}]{\tagform@{\ref*{E:VECM_w_SOI}}} does not allow for linear deterministic trends in the sinks, so that $C_t$ must be used to explain the trending behavior of the sinks.
The restricted, structural model \hyperref[{E:VECM_w_SOI}]{\tagform@{\ref*{E:VECM_w_SOI}}} has the parameter vector \[ \theta = [a_1, a_2, b_1, b_2, c_1, c_2, c_3, c_4, \delta , d, \phi_1, \phi_2, \phi_3, \phi_4 ] ', \] which is of dimension 14. The log-likelihood can be evaluated using \hyperref[{E:llh}]{\tagform@{\ref*{E:llh}}}, as in the unrestricted model. Model \hyperref[{E:VECM_w_SOI}]{\tagform@{\ref*{E:VECM_w_SOI}}} is restricted in the system matrices \hyperref[{E:VECM}]{\tagform@{\ref*{E:VECM}}} in the VECM, that is, in the way the residuals $\hat{U}_t$ are computed. The numerically maximized log-likelihood of model \hyperref[{E:VECM_w_SOI}]{\tagform@{\ref*{E:VECM_w_SOI}}} is $-27.967$. Testing the restricted specification \hyperref[{E:VECM_w_SOI}]{\tagform@{\ref*{E:VECM_w_SOI}}} against the unrestricted specification \hyperref[{E:VECM_RF}]{\textup{\tagform@{\ref*{E:VECM_RF}}}} therefore amounts to evaluating the LR test statistic \[ LR = -2(-27.967-(-11.373)) = 33.188, \] which is asymptotically $\chi^2$-distributed with $50-14=36$ degrees of freedom boswijk2004. The $P$-value is 0.60, and the restricted model cannot be rejected by the data on the sample period at any standard significance level. We therefore adopt the restricted, structural model \hyperref[{E:VECM_w_SOI}]{\textup{\tagform@{\ref*{E:VECM_w_SOI}}}} as our preferred model.
Table (ref) reports maximum likelihood estimates of the parameters of model \hyperref[{E:VECM_w_SOI}]{\tagform@{\ref*{E:VECM_w_SOI}}} together with standard errors. The sinks $S$ (for both $S=S^L$ and $S=S^O$) are measured in PgC and are linear in atmospheric concentrations $C$ in excess of pre-industrial levels, i.e., $S\approx b(C-C_{1750})$. The magnitude of the intercepts is, thus, approximately $a\approx -b C_{1750}$, where the 1750-level atmospheric concentrations are 593 PgC, or about 280 ppm.
The coefficient parameters of the sinks relative to atmospheric concentrations are very similar, and they compare with those found in BHK2023. The global carbon budget equation \hyperref[{E:GCB}]{\tagform@{\ref*{E:GCB}}} states that $C_t - C_{t-1}=E_t - S_t^L -S_t^O + X_t$, where $X_t=X_{4,t}-X_{1,t}-X_{2,t}$ is a stationary error. As a consequence of \hyperref[{E:sinks}]{\tagform@{\ref*{E:sinks}}}, \[ (1+b_1+b_2) C_t - C_{t-1} = E_t - (a_1+a_2) + X_t, \] which shows that the dynamics of atmospheric concentrations $C_t$ are determined by the stochastic trend in emissions $E_t$, which propagates in time according to the autoregressive lag polynomial $b(L)=1-(1+b_1+b_2)^{-1}L$. The sum $b_1+b_2$ keeps the dynamics of $C_t$ from a second unit root, in addition to the one in emissions. If the sinks saturate Canadell2007Saturation,LeQuere2007 and the coefficients $b_1$ and $b_2$ decrease, the dynamics of $C_t$ approach I(2).
The deviations of the sinks from their trends have different dynamics, as can be seen in Figure (ref). The trend deviations of the land sink are of higher variance with no visually discernible serial correlation, whereas the deviations of the ocean sink are tighter around the trend with more apparent patterns of serial correlation. This is supported by the estimates of the AR(1) parameters $\phi_1$ (just insignificant at 5%) and $\phi_2$ (significantly positive).
The serial correlation in first differences of emissions ($\phi_3$) is indistinguishable from zero, whereas the error in the budget equation has mild serial correlation ($\phi_4$). The drift $d$ in emissions is in line with estimates reported earlier BHK2023. The coefficients $c_1$ and $c_2$ on ENSO3.4 in the sinks have the expected sign and are significantly different from zero. The coefficients $c_3$ and $c_4$ on PDO are positive in both sink equations, but they are estimated less precisely, so we do not draw strong conclusions about their effects.
Figure (ref) shows the fit of model \hyperref[{E:VECM_w_SOI}]{\tagform@{\ref*{E:VECM_w_SOI}}} to the data of first differences of the system variables. It can be seen that the model captures variations in the sinks very well. A random walk with drift is a simple model that captures the main role of emissions as driver of the system and source of deterministic and stochastic trend in the cointegrating relations. The lower left panel of Figure (ref) shows that there is room for improvement in capturing short-run fluctuations of first-differenced emissions, but that is beyond the scope of this paper. In order to explain short-run fluctuations in emissions, macroeconomic data would have to be included in the analysis BHK2020.
Figure (ref) shows that first differences of atmospheric concentrations are upward trending, i.e., that the annual increase in atmospheric concentrations is rising. In the cointegration analysis of the unrestricted model \hyperref[{E:VECM_RF}]{\tagform@{\ref*{E:VECM_RF}}}, there is an unrestricted linear trend for $C_t$. In the restricted model \hyperref[{E:VECM_w_SOI}]{\tagform@{\ref*{E:VECM_w_SOI}}}, in contrast, $\Delta C_t$ is given by the budget equation \hyperref[{E:GCB}]{\tagform@{\ref*{E:GCB}}} as $\Delta C_t = E_t - S^L_t - S^O_t +X_{4,t}$. Thus, the unrestricted linear trend $dt$ in $E_t$ and its stochastic trend enter directly into the first differences of concentrations. In addition, the sinks have trending behavior as well, which similarly enters into the first differences of concentrations. The resulting trend in $\Delta C_t$, which is the trend in emissions minus the trends in the sinks, is not equal to zero, and first differences of atmospheric concentrations appear trend-stationary. Evidently, the parametric restrictions are crucial in capturing this important feature of the data.
Table (ref) shows residual diagnostics for model \hyperref[{E:VECM_w_SOI}]{\tagform@{\ref*{E:VECM_w_SOI}}}. There are no signs of misspecification. We note that without the COVID-19 dummy, both the emissions equation and the system tests for Gaussianity would have rejected strongly.
In this section, we apply the estimated restricted model to construct out-of-sample projections of the climate variables $C_t$, $S_t^L$, and $S_t^O$, conditional on a path for the exogenous emissions variable $E_t$, over the period 2023--2100. In Section (ref), we discuss the exogenous emissions data used for the projections, Section (ref) lays out the projection methodology, and Section (ref) presents the projection results.
Before turning to the scenario projections, we note that Appendix (ref) provides additional conditional predictive evidence for the restricted model. In a fixed-origin holdout exercise over 2008--2022, the model reproduces the land sink, ocean sink, and atmospheric concentrations reasonably well when conditioned on the realized emissions path, with the climate covariates contributing mainly to the short-run variation in the sink equations. In recursive pseudo-out-of-sample comparisons against less structured benchmark models, the restricted structural model performs particularly well for the land sink and atmospheric concentrations and remains competitive for the ocean sink. Taken together, these historical exercises support the use of the restricted model for the conditional forward-looking analysis conducted in this section.
We consider emissions trajectories implied by the Shared Socioeconomic Pathways Riahi2017, as generated by the climate model MAGICC MAGICC.\footnote{The MAGICC climate model can be run in a browser at \url{https://live.magicc.org/}.} These trajectories are hypothesized pathways of future anthropogenic CO$_2$ emissions extensively used in the sixth IPCC Assessment Report. SSP scenarios are labeled as “SSPX-Y”, where “X” specifies a certain socioeconomic “narrative” used in the construction of the scenario, while “Y” describes the level of radiative forcing (given in W/m$^2$) implied by the scenario in the year 2100.\footnote{The conversion from an emissions pathway to a forcing level in 2100 is done using climate models; see oneill2016 for details.} For instance, SSP1-1.9 refers to the SSP scenario in narrative 1, consistent with $1.9$ W/m$^2$ of radiative forcing in the year 2100 above pre-industrial levels. Here we do not focus on the various narratives, but only wish to select a number of emissions trajectories that represent a range of possible future emissions pathways. We therefore select five “canonical” SSP scenarios, SSP1-1.9, SSP4-3.4, SSP2-4.5, SSP3-7.0, and SSP5-8.5, as exogenous input into our projections oneill2016. These are shown in the left-most panels in Figure (ref). We refer to Riahi2017 for further details on the SSP scenarios.
Given an exogenous emissions scenario $E_t^*$ for $t = 2023, 2024, \ldots, 2100$, we form the sequence \[ d_t = E_{t}^\ast-E_{t-1}^\ast, \] where $E_{2022}^* = E_{2022} = 11.0980$ is set equal to the last in-sample value of the emissions from the GCB data set. The sequence $d_t$ is used as a time-varying drift parameter in the emissions equation of our restricted model and serves to ensure that the model-implied emissions are equal to those of the SSP scenario in the out-of-sample period. In particular, for $t \geq 2023$, we modify \hyperref[{E:emiss}]{\tagform@{\ref*{E:emiss}}} as follows
which, using that $d_t = E_{t}^\ast-E_{t-1}^\ast$, can also be written as
showing that using the time-varying drift $d_t$ and shutting down the disturbance term in the emissions equation \hyperref[{E:emiss}]{\tagform@{\ref*{E:emiss}}}, means that the model-implied emissions $E_t$ align with the input SSP trajectory $E_t^\ast$, $t=2023, 2024, \ldots, 2100$.
The sink specifications \hyperref[{E:sinks}]{\tagform@{\ref*{E:sinks}}} for the restricted model assume that sinks are linear in concentrations. While this assumption is appropriate for the in-sample period BHK2023, it may not be an accurate assumption for the out-of-sample period 2023--2100. The reason is that climate change may modify the future sinks/concentration relationships, a phenomenon known as the climate feedback effect Friedlingstein2015. Climate feedback effects mean that the sinks likely take up less CO$_2$ over 2023--2100 than the linear model \hyperref[{E:sinks}]{\tagform@{\ref*{E:sinks}}} would suggest. To capture these feedback effects on the sinks during the out-of-sample period, for $t \geq 2023$ we modify \hyperref[{E:sinks}]{\tagform@{\ref*{E:sinks}}} as
where, for $j = 1,2$ and $t \geq 2023$,
where $\widehat a_{1} = -6.1837$, $\widehat a_2 = -4.6213$, $\widehat b_{1} = 0.0113$, $\widehat b_2 = 0.0087$ are the estimates of $a_{\{1,2\}}, b_{\{1,2\}}$ obtained from the in-sample period 1959--2022, see Table (ref). The parameters $\gamma_{\{1,2\}} \in \mathbb{R}$ determine the degree of feedback in the land and ocean sinks, respectively. When $\gamma_{\{1,2\}} >0$, the activity of the sinks, as a response to the level of concentrations, weakens as $t$ increases. We note that it is possible that $S_t^{\{L,O\}} <0$ for some $t \geq 2023$, which would correspond to the carbon sink becoming a carbon source.
There appears to be no consensus in the climate literature on exactly how climate feedback effects modify the sink/concentration relationship in the future. Here, we parametrize the feedback effect parameters $\gamma_{\{1,2\}}$ by \[ \gamma_{\{1,2\}} = \gamma(p_{\{1,2\}}) = -\frac{\log(1-p_{\{1,2\}})}{28}, \] where $p_{\{1,2\}} \in [0,1]$ denote the relative degree of weakening of the sinks by mid-century, i.e., 28 years after 2022. For instance, $p_1 = 50\%$ would correspond to the case where the land sink activity in 2050 is half of what it would have been in the absence of feedback effects. We consider three different magnitudes of feedback effects specified by $(p_1,p_2) = (0,0)$, $(p_1,p_2) = (25\%,25\%)$, and $(p_1,p_2) = (50\%,50\%)$. We refer to the first specification as “no feedback”, the second as “low feedback”, and the third as “high feedback”. The no feedback specification corresponds to a situation where the sink parameters $a_{\{1,2\},t}, b_{\{1,2\},t}$ are fixed at the values estimated on the in-sample data, whereas the low (high) feedback specification corresponds to a situation where both the land and ocean sinks have weakened by 25% (50%) in 2050.
Using a given SSP scenario for emissions, and substituting \hyperref[{E:emiss2}]{\tagform@{\ref*{E:emiss2}}} and \hyperref[{E:sinks2}]{\tagform@{\ref*{E:sinks2}}} for \hyperref[{E:emiss}]{\tagform@{\ref*{E:emiss}}} and \hyperref[{E:sinks}]{\tagform@{\ref*{E:sinks}}}, respectively, we simulate from the estimated restricted model to generate out-of-sample projections of the climate variables over the period 2023--2100. We generate 100,000 trajectories from the model and report the $2.5\%$, $50\%$, and $97.5\%$ pointwise quantiles. These quantiles reflect sampling uncertainty, and not uncertainty about the emissions scenario, feedback parameters, or parameter estimation.
The first column of Figure (ref) shows the exogenous input emissions trajectories for the five SSP scenarios. The second to fourth columns of Figure (ref) present the projected climate variables for these scenarios, derived using the methodology outlined above. The projected CO$_2$ concentrations presented in the right-most column are of particular interest, as they are the main determinants of global temperature changes.
The emissions exhibit substantial variation across scenarios, and the differences directly influence the projected climate outputs. Additionally, the projections vary significantly depending on the degree of climate feedback incorporated. Notably, in the absence of feedback (blue projections), the sinks continue absorbing CO$_2$ according to \hyperref[{E:sinks}]{\tagform@{\ref*{E:sinks}}}, resulting in substantially lower CO$_2$ concentrations compared to cases with climate feedbacks included.
The magenta dashed lines in the last column of Figure (ref) represent the concentration pathway projected by the MAGICC climate model. These projections align closely with those from the high feedback specification of our restricted model. This alignment is significant for two reasons. First, despite the more complex equations used in MAGICC, our model produces comparable results for all five SSP scenarios studied here, underscoring its capability to replicate projections from detailed climate models. Second, this similarity arises only under the high feedback specification, where both sinks weaken by $p_1 = p_2 = 50\%$ in 2050. This suggests that consistency with the historical data record, on which our model is estimated, requires large climate feedback effects to materialize between 2023 and 2100.
For context, the third IPCC Assessment Report estimated mid-century feedback effects of 21%--43% for the land sink and 6%--25% for the ocean sink IPCC2001_3rd. However, more recent studies indicate that these estimates may be conservative, underestimating the potential magnitude of future feedback effects Friedlingstein2015.
In this paper, we have specified a cointegrated vector autoregressive (CVAR) model for the global variables in the Global Carbon Budget data set. In contrast to earlier comprehensive statistical models of the GCB data, the CVAR approach allows for formal hypothesis tests of the physically motivated functional forms of the land and ocean sinks, the emission dynamics, and the budget equation. The global carbon budget equation and the dependence of the sinks on atmospheric CO$_2$ concentrations imply simultaneous relations between the variables, and they constitute cointegrating relations at the same time. We specified both a restricted and physically motivated model as well as an unrestricted and physically agnostic model that nests the restricted model as a special case. A likelihood ratio test showed that the physically motivated, restricted model is supported by the data. We discussed the estimation results in light of the system nature of the variables. We then used the restricted model to explore future projections of the path of the global carbon cycle using SSP scenarios. Recent discussion in the literature Canadell2007Saturation,LeQuere2007 has suggested possible (partial) saturation in the land and ocean sinks in the near future. In the context of our modeling framework, this would constitute a nonlinear, time-dependent cointegrating relationship. Analyzing the GCB data in a model that allows for time-dependence and nonlinearity in the relation between the sinks and atmospheric concentrations is an interesting topic for future research.