# A New Approach to the CO2 Airborne Fraction:
Enhancing Statistical Precision and Tackling Zero Emissions

Mikkel Bennedsen, Eric Hillebrand, Siem Jan Koopman

###### Abstract

The global fraction of anthropogenically emitted CO2 that stays in the
atmosphere, the CO2 airborne fraction, has been fluctuating around a constant
value over the recent period of 1959–2021. The consensus estimate of the
airborne fraction is around $44\%$, indicating that approximately $44\%$ of
emitted CO2 stays in the atmosphere and the remaining $56\%$ is absorbed by
the oceanic and terrestrial biospheres. In this study, we show that the
conventional estimator of the airborne fraction suffers from a number of
statistical deficiencies, such as non-existence of moments and a non-Gaussian
limiting distribution. We formulate a new estimator of the airborne fraction
that does not suffer from these deficiencies. We show that the new estimator
is superconsistent, has a Gaussian limiting distribution, and is able to
reduce estimation uncertainty substantially. Our empirical analysis leads to
an estimate of the airborne fraction over 1959–2021 of $47.6\%$ ($\pm 1.1\%$;
$1\sigma$). Hence, our empirical results imply a higher, and better
constrained, airborne fraction than the current consensus. Our theoretical
results indicate that a conventional analysis of the airborne fraction may be
problematic, especially in future low emissions scenarios such as those
implied by “net zero” international agreements. Our newly proposed estimator
avoids these issues and thus provides a promising approach to future studies
of the airborne fraction, especially in low emission scenarios.

The amount of anthropogenically emitted CO2 that stays in the atmosphere, the
so-called _airborne fraction_ (AF), is an important quantity in the study of
the carbon-cycle of the Earth system (Bacastow and Keeling, 1973; Siegenthaler
and Oeschger, 1978; Gloor et al., 2010). Whether the AF has been increasing,
decreasing, or remained constant over the period from 1959 to today, where
atmospheric measurements of CO2 concentrations have been available, has been
investigated and debated in the literature. Earlier studies found evidence of
an increasing AF (e.g. Canadell et al., 2007; Raupach et al., 2008; Le Quéré
et al., 2009), but these studies have been critiqued on methodological grounds
(e.g. Knorr, 2009; Ballantyne et al., 2015). Later studies suggest that the AF
has remained constant around $44\%$ and this now appears to be the consensus
view (Raupach et al., 2014; Bennedsen et al., 2019; Canadell et al., 2021). A
recent study argues for a declining trend in the AF (van Marle et al., 2022),
but this finding can also be attributed to a structural break in the data
around $1990$ (Bennedsen et al., 2023b). Raupach (2013) shows that the AF, as
a ratio of states in a so-called “LinExp” system, is given by the constant
$r\,/\,(r+k)$ where $r$ is the exponential growth rate of carbon emissions and
$k$ is the linear absorption rate of the carbon sinks. Statistical estimates
of these parameters yield a point estimate of the AF of
$0.016\,/\,(0.016+0.02)\approx 0.44$ (Bennedsen et al., 2023a).

Previous studies have analyzed the AF as the ratio of yearly changes in
atmospheric CO2 concentrations ($G_{t}$, numerator) and anthropogenic CO2
emissions ($E_{t}$, denominator) (e.g. Canadell et al., 2007; Raupach et al.,
2008; Le Quéré et al., 2009; Knorr, 2009; Ballantyne et al., 2012; Raupach et
al., 2014; Ballantyne et al., 2015; Keenan et al., 2016; Bennedsen et al.,
2019; van Marle et al., 2022). Following this literature, the conventional
estimator of the AF has been defined as the sample average of the yearly ratio
$G_{t}/E_{t}$. In this paper, we show that this estimator suffers from a
number of statistical deficiencies due to its definition as the ratio of two
stochastic processes. One problem is that the denominator in the ratio,
$E_{t}$, may have a positive probability density at zero, implying that the
estimator does not possess any moments. This lack of moments implies, for
instance, that the mean and variance of the estimator do not exist. Although
this issue may be of limited concern during the period studied in this paper
(1959–2021), it becomes important in future scenarios where CO2 emissions
decrease, such as those consistent with “net-zero” CO2 emissions, a committed
goal of the international community (Rockström et al., 2017; Riahi et al.,
2022). For example, issues with the conventional analysis of the AF when
emissions decrease are recently encountered in Pressburger et al. (2023),
where the future AF implied by output from a climate model is studied. Another
problem is that the time series of changes in atmospheric concentrations
$G_{t}$ and emissions $E_{t}$ are trending, which implies that a central limit
theorem does not hold for the estimator. Therefore, its limiting distribution
may be non-Gaussian, unless a separate assumption of Gaussianity is imposed.

We show that the time series of changes in atmospheric CO2 concentrations,
$G_{t}$, cointegrates with the time series of anthropogenic CO2 emissions,
$E_{t}$, over the period 1959–2021. Cointegration means that, even though
$G_{t}$ and $E_{t}$ are individually non-stationary time series, the linear
regression equation $G_{t}=\alpha E_{t}+u_{t}$, with constant parameter
$\alpha$, has a stationary error $u_{t}$. In other words, the non-stationary
processes $G_{t}$ and $E_{t}$ are cointegrated when a linear combination
exists that yields a stationary process. The parameter $\alpha$ denotes the
fraction of emissions that is added to the atmosphere in year $t$, i.e.,
$\alpha$ is the AF. We show that estimating the AF $\alpha$ using an ordinary
least squares regression where yearly changes in atmospheric CO2
concentrations $G_{t}$ are regressed on yearly anthropogenic CO2 emissions
$E_{t}$ yields an estimator with superior statistical properties compared to
the conventional AF estimator based on the ratio $G_{t}\,/\,E_{t}$. The new
estimator does not suffer from the defects described above, and, under mild
assumptions, converges to the true value $\alpha$ at the fast rate of
$T^{3/2}$, where $T$ is the sample size. A central limit theorem applies, i.e.
the estimator has a Gaussian distribution as the sample size increases to
infinity for a wide range of distributions of the regression residual $u_{t}$.

We apply the new and the conventional estimators of the AF to yearly data from
the Global Carbon Project (Friedlingstein et al., 2022) over the period
1959–2021, and we find that the new estimator improves precision, as measured
by the standard error of the estimator, by approximately $10\%$ compared to
the conventional estimator. Our estimate of the AF over the period 1959–2021
is $45.3\%$ with an associated standard error of $1.5\%$. Hence, this estimate
is not significantly different from the consensus estimate of $44\%$. To
further reduce estimation uncertainty, we also consider the inclusion of
additional data on the El Niño-Southern Oscillation (ENSO) and volcanic
activity in the analysis. In this case, our proposed estimator improves
precision by approximately $16\%$, compared to the conventional estimator of
the AF. Our best estimate of the AF over the period 1959–2021 is $47.6\%$ with
an associated standard error of $1.1\%$, which leads to a $95\%$ confidence
interval of $[45.5\%,49.8\%]$ for the AF. The consensus estimate of $44\%$
(Canadell et al., 2021, p. 676) falls just outside this confidence interval.
Our empirical results thus imply a higher, and better constrained, AF than the
current consensus. The higher value of the estimate is mainly driven by the
inclusion of ENSO and volcanic activity; the tighter confidence interval is
mainly driven by our new estimation approach.

## Atmospheric changes, emissions, and cointegration

Figures 1a) and 1b) show the yearly changes in atmospheric CO2 ($G_{t}$) and
the yearly CO2 emissions from anthropogenic sources ($E_{t}$), respectively.
Figure 1c) shows the ratio of these two variables, $G_{t}/E_{t}$. Data are
obtained from the Global Carbon Project and cover the period $1959$–$2021$; in
the empirical section we provide further details on the data. The most
conspicuous feature of the two data series $G_{t}$ and $E_{t}$ are that both
time series are trending upwards. Trending behavior is indicative of time
series being non-stationary. A simple least-squares statistical analysis where
the non-stationarity of a time series is not accounted for yields invalid
inference and should be avoided (Granger and Newbold, 1974). However, the
notion of cointegration (see, for example, chapter 19 in Hamilton (1994) for a
textbook treatment) allows us to remain working with the trending time series
$G_{t}$ and $E_{t}$ while still obtaining valid statistical inference.
Cointegration methods have been applied in earlier climate studies, see, e.g.,
Schmith et al. (2012). The time series $G_{t}$ and $E_{t}$ are said to be
cointegrated when both $G_{t}$ and $E_{t}$ are non-stationary, and the error
term $u_{t}$ in the regression equation $G_{t}=\alpha E_{t}+u_{t}$ is a
stationary process. We adopt the Dickey-Fuller test of Dickey and Fuller
(1979) to determine whether a time series is non-stationary. The test
statistic is for the null hypothesis of a unit root, that is, of having a unit
value for the autoregressive dependence of $G_{t}$ on its lagged value
$G_{t-1}$ (and $E_{t}$ on $E_{t-1}$). The results from this test strongly
suggest that both time series are non-stationary (Appendix B.1). The null
hypothesis of no cointegration (that is, $u_{t}$ is non-stationary) can be
tested formally using the Engle-Granger statistic (Engle and Granger, 1987),
which is a Dickey-Fuller test on the residuals in the regression $G_{t}=\alpha
E_{t}+u_{t}$, adjusted for the fact that these residuals are not observed but
must be estimated. The null hypothesis of a unit root in $u_{t}$ is firmly
rejected (Appendix B.1), and we can therefore conclude that the two time
series are cointegrated.

The cointegration analysis also provides evidence that, over the period
studied here (1959–2021), the AF parameter $\alpha$ is constant over time,
i.e., the AF does not have trending behaviour. If the parameter $\alpha$ was
changing in a specific direction, this would introduce a trend in the
residuals $u_{t}$. The result from the Engle-Granger test suggests that such a
trend is not present. This is confirmed graphically in Figure 1d) and is in
line with recent studies (Raupach et al., 2014; Bennedsen et al., 2019;
Canadell et al., 2021; Bennedsen et al., 2023b). A Jarque-Bera test (Jarque
and Bera, 1987) of the null of normality applied to the estimated residuals
for $u_{t}$ results in a $p$-value of $25.3\%$, implying that we cannot reject
that $u_{t}$ is a Gaussian process.

Figure 1: _a) Atmospheric concentration changes ( $G_{t}$) data used in the
study over the period $1959$–$2019$. b): Emissions ($E_{t}$) data used in the
study over the period $1959$–$2019$. c): The ratio of atmospheric changes to
emissions ($G_{t}/E_{t}$). d): Estimated residuals
$\hat{u}_{t}=G_{t}-\hat{\alpha}E_{t}$ from the model $G_{t}=\alpha
E_{t}+u_{t}$, where $\alpha$ is estimated via least-squares regression._

## A new estimator of the airborne fraction

The conventional estimator of the AF can be regarded as the result of a
regression based on the statistical model given by

$\displaystyle\frac{G_{t}}{E_{t}}=\alpha+u^{(1)}_{t},$ (1)

where $G_{t}$ are the yearly changes in atmospheric concentrations of CO2,
$E_{t}$ are yearly CO2 emissions, the constant parameter $\alpha$ is the AF,
and $u^{(1)}_{t}$ is the disturbance modelled as a zero-mean error process,
for $t=1,2,\ldots,T$, with $T$ denoting the number of yearly observations in
the sample. The disturbance $u^{(1)}_{t}$ captures deviations of the data
${G_{t}}/{E_{t}}$ from the constant value $\alpha$ due to measurement errors
and internal variability of the climate system. For the statistical model (1),
it is straightforward to estimate the AF parameter $\alpha$ using the sample
average of the data ${G_{t}}/{E_{t}}$, yielding the estimator

$\hat{\alpha}_{1}=\frac{1}{T}\sum_{t=1}^{T}\frac{G_{t}}{E_{t}}.$

The model in equation (1) expresses that, on average, the fraction $\alpha$ of
emissions $E_{t}$ is absorbed in the atmosphere, resulting in atmospheric
concentrations increasing with the amount $G_{t}$. An alternative way to
express this association between $G_{t}$ and $E_{t}$ is through the model
formulation

$\displaystyle G_{t}=\alpha E_{t}+u^{(2)}_{t},$ (2)

for $t=1,2,\ldots,T$, where the disturbance $u^{(2)}_{t}$ is also a zero-mean
error process. The relationship between the disturbances in equations (1) and
(2) is given by $u^{(2)}_{t}=E_{t}\cdot u^{(1)}_{t}$. Our earlier
cointegration analysis suggests that $u_{t}^{(2)}$ is a stationary process. In
this formulation and under these assumptions, the parameter $\alpha$ can be
estimated straightforwardly using a simple least-squares calculation, yielding
the estimator

$\hat{\alpha}_{2}=\left(\sum_{t=1}^{T}E_{t}^{2}\right)^{-1}\sum_{t=1}^{T}E_{t}G_{t}.$

### Statistical properties of the two estimators

We show in Appendix A that the statistical properties of the estimators
$\hat{\alpha}_{1}$ and $\hat{\alpha}_{2}$ depend on the dynamics of CO2
emissions $E_{t}$. In Figure 2a) we present the differences of yearly CO2
emissions, $\Delta E_{t}=E_{t}-E_{t-1}$, and in Figure 2b) we present the
autocorrelation function of $\Delta E_{t}$. Figure 2a) suggests that the model
$\Delta E_{t}=b+\xi_{t}$, where $b$ is a constant and $\xi_{t}$ is a zero-mean
error process, is appropriate for $\Delta E_{t}$, while Figure 2b) suggests
that $\xi_{t}$ is actually a white noise process. This implies the model
$E_{t}=E_{t-1}+b+\xi_{t}$ for CO2 emissions, which can be written as
$E_{t}=E_{0}+bt+x_{t}$, where $x_{t}=\sum_{\tau=1}^{t}\xi_{\tau}$ is a random
walk. We can estimate $b$ straightforwardly to find $\hat{b}=0.1070$ (standard
error $0.0239$). This estimate is close to the corresponding estimate found in
Bennedsen et al. (2023a). A Jarque-Bera test (Jarque and Bera, 1987) of the
null of normality applied to the estimated residuals $\hat{\xi}_{t}=\Delta
E_{t}-\hat{b}$ results in a $p$-value of $40.0\%$, meaning that we cannot
reject that $\hat{\xi}_{t}$ is a Gaussian process. It emerges from Figure 2a)
that the variability of the differences in CO2 emissions increases in the
early 1990s. This is most likely due to increased volatility of estimates of
emissions from land-use and land-cover change, starting in the early 1990s
(see Appendix B.2 and also van Marle et al., 2022). In Appendix B.3 we re-
analyze the subsample 1992–2021 as a robustness check, and we find results
very similar to those reported below.

Figure 2: _Yearly changes in emissions ( $\Delta E_{t}=E_{t}-E_{t-1}$) over
the period $1960$–$2019$. b) Autocorrelation of $\Delta E_{t}$. _

Under the assumption that $E_{t}=E_{0}+bt+x_{t}$, we can show that (Appendix
A, Proposition 1)

$\displaystyle T^{3/2}\cdot(\hat{\alpha}_{2}-\alpha)\stackrel{{\scriptstyle
d}}{{\rightarrow}}N(0,\sigma_{2}^{2}),$

as $T\rightarrow\infty$, where $\stackrel{{\scriptstyle d}}{{\rightarrow}}$
denotes _convergence in distribution_ , $N$ is the normal distribution, and
the variance $\sigma_{2}^{2}$ can be estimated using standard methods. It
might appear surprising that the asymptotic distribution of $\hat{\alpha}_{2}$
is Gaussian, since it is well-documented that non-Gaussian limit distributions
often appear in the context of cointegration (e.g. Phillips and Park, 1988).
However, the Gaussian limit obtained here is a result of the properties of
$E_{t}$, namely that it contains the time trend $E_{0}+bt$. This time trend
dominates the stochastic trend in $E_{t}$, which is the random walk $x_{t}$,
and permits using standard least-squares inference in the case of
$\hat{\alpha}_{2}$ (Park, 1992). The estimator $\hat{\alpha}_{2}$ is
“superconsistent”, in the sense that it converges faster to $\alpha$ than the
rate $T^{1/2}$ commonly encountered in central limit theorem statements (see,
e.g., Hamilton, 1994, chapter 16, for a textbook treatment). Intuitively, the
trending behavior of the series results in a high signal-to-noise ratio that
leads to fast convergence of coefficient estimates.

For the estimator $\hat{\alpha}_{1}$ we can show (Appendix A, Proposition 2)
that, for arbitrarily small $\gamma>0$,

$\displaystyle
T^{1-\gamma}\cdot(\hat{\alpha}_{1}-\alpha)\stackrel{{\scriptstyle
p}}{{\rightarrow}}0,$

as $T\rightarrow\infty$, where $\stackrel{{\scriptstyle p}}{{\rightarrow}}$
denotes _convergence in probability_. This results implies that
$\hat{\alpha}_{1}$ is also a superconsistent estimator of $\alpha$, albeit at
a slower rate than $\hat{\alpha}_{2}$. To the best of our knowledge, even
though the estimator $\hat{\alpha}_{1}$ has often been used to assess the AF
previously, this is the first time it is noted that $\hat{\alpha}_{1}$ is
superconsistent. This somewhat surprising result is obtained from the
assumption that $G_{t}=\alpha E_{t}+u_{t}$, yielding
$G_{t}/E_{t}=\alpha+u_{t}/E_{t}$, with $u_{t}$ being a stationary error
process and $E_{t}$ a trending process, for which we found evidence in the
cointegration analysis above. Further details are provided in Appendix A.

If the process $E_{t}$ has a positive probability density at zero, then
$\hat{\alpha}_{1}$ does not have any moments (Appendix A, Proposition 3). In
particular, the expectation $E(\hat{\alpha}_{1})$ and variance
$Var(\hat{\alpha}_{1})$ do not exist. This follows from the fact that
$\hat{\alpha}_{1}$ is defined using the ratio $G_{t}/E_{t}$, and if values
$E_{t}\approx 0$ are possible, then the resulting ratio $G_{t}/E_{t}$ is not
integrable. Given the model expression $E_{t}=E_{0}+bt+x_{t}$ for CO2
emissions, $E_{t}$ has a positive probability density at zero if, e.g.,
$x_{t}$ is a random walk of Gaussian increments, as we have found evidence of
above. On our sample, however, $E_{0}+bt$ is much larger in magnitude than
$x_{t}$. If we impose the assumption that $x_{t}=0$, then the moments of
$\hat{\alpha}_{1}$ do exist. In particular, in this case, $\hat{\alpha}_{1}$
is an unbiased estimator of $\alpha$, i.e. $E(\hat{\alpha}_{1})=\alpha$, and
the variance $Var(\hat{\alpha}_{1})$ has a simple form that can be easily
estimated. However, we can show that even when $x_{t}=0$, then, in general, a
central limit theorem does _not_ hold for $\hat{\alpha}_{1}$ (Appendix A,
Proposition 4(i)). That is, without further assumptions, we cannot conclude
that $\hat{\alpha}_{1}$ is Gaussian, even asymptotically. If one is willing to
further assume that both $x_{t}=0$ _and_ $u_{t}$ is Gaussian, then we can
conclude that (Appendix A, Proposition 4(ii))

$\displaystyle T\cdot(\hat{\alpha}_{1}-\alpha)\stackrel{{\scriptstyle
d}}{{\rightarrow}}N(0,\sigma_{1}^{2}),$

as $T\rightarrow\infty$, where the variance $\sigma_{1}^{2}$ can be estimated
using standard methods.

### An assessment of the two estimators

The discussion above shows that the construction of $\hat{\alpha}_{1}$, based
on the ratio $G_{t}/E_{t}$, is statistically problematic, especially when
$E_{t}$ has a positive probability density at zero. While this may not be an
issue on the sample studied here (1959–2021), where the constant and drift
$E_{0}+bt$ dominate the randomness of $x_{t}$, it does become a problem if
$E_{t}$ approaches zero. That is, in low(er) emission regimes, an estimator of
the AF based on the ratio $G_{t}/E_{t}$ has poor statistical properties, such
as non-existence of moments. In any case, $\hat{\alpha}_{1}$ is
superconsistent, converging at the rate $T$, meaning that it may, in practice,
provide a precise estimate of $\alpha$. However, our newly proposed estimator
$\hat{\alpha}_{2}$ converges at the even faster rate $T^{3/2}$, and it does
not have the statistical difficulties of $\hat{\alpha}_{1}$. By avoiding the
ratio $G_{t}/E_{t}$, and instead leveraging the cointegrating relation between
$G_{t}$ and $E_{t}$, a more precise and reliable estimator is obtained.

Although the estimator $\hat{\alpha}_{2}$ is _asymptotically_ (i.e. as $T$
increases) more precise than $\hat{\alpha}_{1}$, it is an empirical question
which estimator is more precise in smaller samples. In the next section, we
use historical data to estimate both $Var(\hat{\alpha}_{1})$ and
$Var(\hat{\alpha}_{2})$ to assess which variance is larger on our data sample.
We note again that $Var(\hat{\alpha}_{1})$ is only well-defined under the
strong assumption that $E_{t}$ does not have positive probability density at
zero, e.g. if $x_{t}=0$. Further, in the next section, we construct confidence
intervals for our estimators, based on the Gaussian distribution. The
confidence interval for $\hat{\alpha}_{1}$ is valid only under the additional
assumption that $u_{t}$ is Gaussian, whereas the validity of the confidence
interval for $\hat{\alpha}_{2}$ can be established by asymptotic arguments and
does therefore not depend on any specific distributional assumption on
$u_{t}$.

## Estimating the airborne fraction

We use data on yearly changes in atmospheric CO2 ($G_{t}$), yearly CO2
emissions from fossil fuels ($E_{t}^{FF}$), and yearly CO2 emissions from
land-use and land cover change ($E_{t}^{LULCC}$). Total anthropogenic CO2
emissions are then $E_{t}=E_{t}^{FF}+E_{t}^{LULCC}$. The data series are
measured in gigatonnes of carbon per year (GtC/yr), obtained from the Global
Carbon Project111Data from the Global Carbon Project can be found here:
https://www.icos-cp.eu/science-and-impact/global-carbon-budget/2022, last
accessed April 18, 2023.  and presented in Figure 1.

The estimates $\hat{\alpha}_{1}$ and $\hat{\alpha}_{2}$ are obtained from
least-squares regressions applied to models (1) and (2), respectively, using
the Global Carbon Project data. To account for possible serial correlation and
heteroskedasticity in the error terms $u$, we obtain estimates of the standard
errors using the heteroskedasticity and autocorrelation consistent (HAC)
estimator of Newey and West (1987). The results are displayed in the first two
columns of Table 1. The estimates largely agree on the magnitude of the AF,
$\hat{\alpha}_{1}=44.41\%$ and $\hat{\alpha}_{2}=45.33\%$. However, the
estimated standard deviations of the error terms in the two models are very
different, namely
$\widehat{SD}(u_{t}^{(1)})=\sqrt{\widehat{Var}(u_{t}^{(1)})}=0.13$ and
$\widehat{SD}(u_{t}^{(2)})=\sqrt{\widehat{Var}(u_{t}^{(1)})}=0.92$. The
discrepancy is due to the very different nature of the two models, in
particular, we have $u_{t}^{(1)}=u_{t}^{(2)}/E_{t}$, where $E_{t}>>1$ in the
sample 1959–2021. The standard error of $\hat{\alpha}_{2}$ is $10\%$ lower
than the standard error of $\hat{\alpha}_{1}$, indicating that the faster
convergence rate of this estimator outweighs the fact that the error process
$u_{t}^{(2)}$ in (2) has a larger variance than the error process
$u_{t}^{(1)}$ in (1).

Table 1: Least-squares regression output of the linear models in Equations (1)–(4) for the Global Carbon Project data for the period 1959–2021. “$\hat{\alpha}$” denotes the estimate of the AF $\alpha$, “SE($\hat{\alpha}$)” denotes the estimated standard error of $\hat{\alpha}$, and “Relative SE” denotes the ratio of the standard error to the standard error from $\hat{\alpha}_{1}$, corresponding to the left-most column. Hence, the left-most column has a Relative SE of one by definition, while a Relative SE below one in the other columns indicate a lower standard error for the estimator in question. “$CI_{95\%}(\alpha)$” denotes a $95\%$ confidence interval for the AF $\alpha$, based on the Gaussian distribution. “$\widehat{SD}(u_{t})$” denotes the estimated standard deviation of the error term $u_{t}$ and $R^{2}$ is the coefficient of determination. |  (1) |  (2) |  (3) |  (4)
---|---|---|---|---
$\hat{\alpha}$ | 0.4441 | 0.4533 | 0.4785 | 0.4763
SE($\hat{\alpha}$) | 0.0163 | 0.0147 | 0.0129 | 0.0109
Relative SE | 1.0000 | 0.8986 | 0.7881 | 0.6712
$CI_{95\%}(\alpha)$ | $[0.4121,0.4761]$ | $[0.4246,0.4820]$ | $[0.4533,0.5037]$ | $[0.4548,0.4977]$
$\widehat{SD}(u_{t})$ | 0.1284 | 0.9153 | 0.0894 | 0.6327
$R^{2}$ | 0 | 0.5841 | 0.5310 | 0.8077

To reduce the variance of the error processes $u$, and thus achieve more
precise estimates of the AF $\alpha$, we introduce exogenous variables in the
models (1)–(2). In particular, we control for the effects of ENSO and volcanic
activity through two yearly indices, denoted by $ENSO_{t}$ and $VAI_{t}$,
respectively, a common practice in the literature (e.g. Raupach et al., 2008;
van Marle et al., 2022).222ENSO data is the Niño 3 SST Index from the National
Oceanic and Atmospheric Administration (NOAA), available at
https://psl.noaa.gov/gcos_wgsp/Timeseries/Data/nino3.long.anom.data, last
accessed April 18, 2023, and data from volcanic activity is from Ammann et al.
(2003). The ENSO data contain a slight trend, which was removed by de-trending
before using the data in the analysis of this paper, to avoid having the trend
influence the estimates of the AF. The ENSO and VAI data series used here are
shown in Figure 3 in Appendix B.2. Hence, we consider the models

$\displaystyle\frac{G_{t}}{E_{t}}$ $\displaystyle=\ \alpha\ +\
\tilde{\gamma}_{1}ENSO_{t}+\tilde{\gamma}_{2}VAI_{t}+u^{(3)}_{t},$ (3)
$\displaystyle G_{t}$ $\displaystyle=\alpha
E_{t}+\gamma_{1}ENSO_{t}+\gamma_{2}VAI_{t}+u^{(4)}_{t},$ (4)

for $t=1,2,\ldots,T$, where $\tilde{\gamma}_{i},\,\gamma_{i}$, $i=1,2$, are
regression coefficients, and the disturbances $u^{(j)}_{t}$ follow a zero-mean
error process for $j=3,4$. For both models, the parameters can be estimated
using least-squares regression. Let $\hat{\alpha}_{3}$ and $\hat{\alpha}_{4}$
denote the least-squares estimators of $\alpha$ from models (3) and (4),
respectively. The estimation results are presented in columns three and four
of Table 1. The results show that controlling for the effects of ENSO and
volcanic activity increases the estimate of the AF considerably, resulting in
$\hat{\alpha}_{3}=47.85\%$ and $\hat{\alpha}_{4}=47.63\%$. The estimates of
the standard deviations of the error terms $(\widehat{SD}(u_{t}^{(3)})=0.09$
and $\widehat{SD}(u_{t}^{(4)})=0.63)$ decrease substantially compared to
models (1)–(2), indicating that the explanatory variables ENSO and VAI help in
explaining the variation in the data, as also indicated by the $R^{2}$ values
reported in the bottom row of Table 1. The decreased variance of the error
terms means that our estimates of the constant AF $\alpha$ become more
precise. The standard error on $\hat{\alpha}_{4}$ is approximately $16\%$
lower than the estimator $\hat{\alpha}_{3}$, and approximately $33\%$ lower
than when using the conventional estimator $\hat{\alpha}_{1}$. Our preferred
estimator $\hat{\alpha}_{4}$ results in an estimate of the AF of $47.6\%$
($\pm 1.1\%$; $1\sigma$) with an associated $95\%$ confidence interval of
$[45.5\%,49.8\%]$.

## Discussion

In our empirical study, the findings suggest a slightly higher AF than the
$44\%$ that have so far been accepted as the consensus estimate (Canadell et
al., 2021, p. 676). Our best estimate of the AF over 1959-2021, using the most
recent edition of the Global Carbon Project data set and controlling for ENSO
and volcanic activity, is $47.6\%$ ($\pm 1.1\%$; $1\sigma$) with an associated
$95\%$ confidence interval of $[45.5\%,49.8\%]$. Excluding ENSO and volcanic
activity from the analysis results in an estimate more in line with the
consensus, namely $45.3\%$ ($\pm 1.5\%$; $1\sigma$). An analysis of two
alternative data sets yield slightly higher estimates of the AF than those
from the Global Carbon Project (Appendix B.2). Restricting the analysis to the
more recent 1992–2021 sample, we find similar results to those reported here,
namely an AF of approximately $47\%$ ($\pm 1.1\%$; $1\sigma$), regardless of
the data set used (Appendix B.3).

The theoretical assessment of the conventional estimators $\hat{\alpha}_{1}$
and $\hat{\alpha}_{3}$ can be summarized by three statistical features that
all arise from the definition as the mean of the ratio $G_{t}/E_{t}$. First,
the moments of $\hat{\alpha}_{1}$ and $\hat{\alpha}_{3}$ do not exist,
including their mean and variance, unless it is explicitly assumed that the
probability density of $E_{t}$ vanishes at zero. This last assumption is not
sustainable in future scenarios because many international agreements and
global policies are committed to “net zero” emission goals. Second, a central
limit theorem does not hold for $\hat{\alpha}_{1}$ and $\hat{\alpha}_{3}$.
Consequently, confidence intervals and $p$-values for various test statistics
based on the Gaussian distribution are not valid, unless the more restrictive
assumption is made that $x_{t}=0$ and the error term $u_{t}$ in the model
$G_{t}/E_{t}=\alpha+u_{t}$ is Gaussian. Third, the estimators
$\hat{\alpha}_{1}$ and $\hat{\alpha}_{3}$ both converge to the AF $\alpha$ at
superconsistent rate $T$, i.e., faster than $T^{1/2}$. However, the
convergence rate is considerably lower when compared to our newly proposed AF
estimators $\hat{\alpha}_{2}$ and $\hat{\alpha}_{4}$ that converge at the rate
$T^{3/2}$. Our newly proposed estimators can therefore be expected to be more
precise than the conventional estimators $\hat{\alpha}_{1}$ and
$\hat{\alpha}_{3}$, which is what we find on the historical data. In a
simulation study we find similar improvements of $\hat{\alpha}_{2}$ compared
to $\hat{\alpha}_{1}$ and show that these improvements continue as the samples
size is increased (Appendix C.1). In addition, we show that $\hat{\alpha}_{2}$
and $\hat{\alpha}_{4}$ follow a central limit theorem, which implies that the
construction of confidence intervals and $p$-values based on the Gaussian
distribution is valid.

The possible lack of moments for $\hat{\alpha}_{1}$ and $\hat{\alpha}_{3}$ is
of minor importance for empirical findings obtained from analyzing the
historical sample 1959–2021, since emissions $E_{t}$ have been much larger
than zero during this period (Figure 1b). Also, the lack of a central limit
theorem is likely of limited importance on the historical sample period, since
our cointegration analysis shows that the null of Gaussian error terms $u_{t}$
in the regression equation $G_{t}=\alpha E_{t}+u_{t}$ cannot be rejected on
the basis of the Jarque-Bera test. Therefore, previous studies of the AF on
the historical period relying on the ratio $G_{t}/E_{t}$ (such as Canadell et
al., 2007; Raupach et al., 2008; Le Quéré et al., 2009; Knorr, 2009;
Ballantyne et al., 2012; Raupach et al., 2014; Ballantyne et al., 2015; Keenan
et al., 2016; Bennedsen et al., 2019; van Marle et al., 2022, among others)
are most likely not influenced by the statistical issues pointed out in this
paper to any considerable degree. Therefore, when applied to the historical
record 1959–2021, the main advantage of our new estimator compared to the
conventional estimator is increased precision, resulting in tighter confidence
bounds on the historical AF. However, the conventional approach of analyzing
the AF using the ratio $G_{t}/E_{t}$ results in inefficient, or even invalid,
statistical analyses when applied to scenarios where global CO2 emissions
decrease to “net zero” (e.g. Rockström et al., 2017; Riahi et al., 2022). For
completeness, we illustrate these issues in a simulation study in Appendix
C.2. Such situations could arise both in future analyses of observational data
(at least if international “net zero” commitments are kept) and in
contemporary analyses of output from climate models where low emission
scenarios are considered. Our newly proposed AF estimator avoids these issues
and thus provides a promising approach to future studies of the AF, especially
in low emission scenarios.

## Acknowledgements

We thank Morten Ø. Nielsen for helpful discussions regarding convergence of
stochastic processes. M.B. gratefully acknowledges funding from the
Independent Research Fund Denmark under grant no. 0219-00001B.

## References

  * Ammann et al. (2003) Ammann, C. M., G. A. Meehl, W. M. Washington, and C. S. Zender (2003).  A monthly and latitudinally varying volcanic forcing dataset in simulations of 20th century climate.  Geophysical Research Letters 30(12).
  * Bacastow and Keeling (1973) Bacastow, R. and C. D. Keeling (1973).  Atmospheric carbon dioxide and radiocarbon in the natural cycle: Ii. changes from a. d. 1700 to 2070 as deduced from a geochemical model.  In Carbon and the biosphere conference proceedings; Upton, New York, USA.
  * Ballantyne et al. (2012) Ballantyne, A. P., C. B. Alden, J. B. Miller, P. P. Tans, and J. W. C. White (2012).  Increase in observed net carbon dioxide uptake by land and oceans during the past 50 years.  Nature 488(7409), 70–72.
  * Ballantyne et al. (2015) Ballantyne, A. P., R. Andres, R. Houghton, B. D. Stocker, R. Wanninkhof, W. Anderegg, L. A. Cooper, M. DeGrandpre, P. P. Tans, J. B. Miller, C. Alden, and J. W. C. White (2015).  Audit of the global carbon budget: estimate errors and their impact on uptake uncertainty.  Biogeosciences 12(8), 2565–2584.
  * Bennedsen et al. (2019) Bennedsen, M., E. Hillebrand, and S. J. Koopman (2019).  Trend analysis of the airborne fraction and sink rate of anthropogenically released CO2.  Biogeosciences 16(18), 3651–3663.
  * Bennedsen et al. (2023a) Bennedsen, M., E. Hillebrand, and S. J. Koopman (2023a).  A multivariate dynamic statistical model of the global carbon budget 1959–2020.  Journal of the Royal Statistical Society Series A: Statistics in Society 186(1), 20–42.
  * Bennedsen et al. (2023b) Bennedsen, M., E. Hillebrand, and S. J. Koopman (2023b).  On the evidence of a trend in the CO2 airborne fraction.  Nature 616(7956), E1–E3.
  * Canadell et al. (2021) Canadell, J., P. Monteiro, M. Costa, L. C. da Cunha, P. Cox, A. Eliseev, S. Henson, M. Ishii, S. Jaccard, C. Koven, A. Lohila, P. Patra, S. Piao, J. Rogelj, S. Syampungani, S. Zaehle, and K. Zickfeld (2021).  Global carbon and other biogeochemical cycles and feedbacks.  In V. Masson-Delmotte, P. Zhai, A. Pirani, S. Connors, C. Péan, S. Berger, N. Caud, Y. Chen, L. Goldfarb, M. Gomis, M. Huang, K. Leitzell, E. Lonnoy, J. Matthews, T. Maycock, T. Waterfield, O. Yelekçi, R. Yu, and B. Zhou (Eds.), Climate Change 2021: The Physical Science Basis. Contribution of Working Group I to the Sixth Assessment Report of the Intergovernmental Panel on Climate Change, pp. 673–816. Cambridge University Press, Cambridge, United Kingdom and New York, NY, USA.
  * Canadell et al. (2007) Canadell, J. G., C. Le Quéré, M. R. Raupach, C. B. Field, E. T. Buitenhuis, P. Ciais, T. J. Conway, N. P. Gillett, R. A. Houghton, and G. Marland (2007).  Contributions to accelerating atmospheric CO2 growth from economic activity, carbon intensity, and efficiency of natural sinks.  Proceedings of the National Academy of Sciences 104(47), 18866–18870.
  * Davidson and de Jong (1997) Davidson, J. and R. M. de Jong (1997).  Strong laws of large numbers for dependent heterogeneous processes: a synthesis of recent and new results.  Econometric Reviews 16(3), 251–279.
  * de Jong (1995) de Jong, R. M. (1995).  Laws of large numbers for dependent heterogeneous processes.  Econometric Theory 11(2), 347–358.
  * Dickey and Fuller (1979) Dickey, D. A. and W. A. Fuller (1979).  Distribution of the estimators for autoregressive time series with a unit root.  Journal of the American Statistical Association 74(366a), 427–431.
  * Engle and Granger (1987) Engle, R. F. and C. W. J. Granger (1987).  Co-integration and error correction: Representation, estimation, and testing.  Econometrica 55(2), 251–276.
  * Fazekas and Klesov (2001) Fazekas, I. and O. Klesov (2001).  A general approach to the strong law of large numbers.  Theory of Probability & Its Applications 45(3), 436–449.
  * Friedlingstein et al. (2022) Friedlingstein, P., M. O’Sullivan, M. W. Jones, R. M. Andrew, L. Gregor, J. Hauck, C. Le Quéré, I. T. Luijkx, A. Olsen, G. P. Peters, W. Peters, J. Pongratz, C. Schwingshackl, S. Sitch, J. G. Canadell, P. Ciais, R. B. Jackson, S. R. Alin, R. Alkama, A. Arneth, V. K. Arora, N. R. Bates, M. Becker, N. Bellouin, H. C. Bittig, L. Bopp, F. Chevallier, L. P. Chini, M. Cronin, W. Evans, S. Falk, R. A. Feely, T. Gasser, M. Gehlen, T. Gkritzalis, L. Gloege, G. Grassi, N. Gruber, O. Gürses, I. Harris, M. Hefner, R. A. Houghton, G. C. Hurtt, Y. Iida, T. Ilyina, A. K. Jain, A. Jersild, K. Kadono, E. Kato, D. Kennedy, K. Klein Goldewijk, J. Knauer, J. I. Korsbakken, P. Landschützer, N. Lefèvre, K. Lindsay, J. Liu, Z. Liu, G. Marland, N. Mayot, M. J. McGrath, N. Metzl, N. M. Monacci, D. R. Munro, S.-I. Nakaoka, Y. Niwa, K. O’Brien, T. Ono, P. I. Palmer, N. Pan, D. Pierrot, K. Pocock, B. Poulter, L. Resplandy, E. Robertson, C. Rödenbeck, C. Rodriguez, T. M. Rosan, J. Schwinger, R. Séférian, J. D. Shutler, I. Skjelvan, T. Steinhoff, Q. Sun, A. J. Sutton, C. Sweeney, S. Takao, T. Tanhua, P. P. Tans, X. Tian, H. Tian, B. Tilbrook, H. Tsujino, F. Tubiello, G. R. van der Werf, A. P. Walker, R. Wanninkhof, C. Whitehead, A. Willstrand Wranne, R. Wright, W. Yuan, C. Yue, X. Yue, S. Zaehle, J. Zeng, and B. Zheng (2022).  Global Carbon Budget 2022.  Earth System Science Data 14(11), 4811–4900.
  * Gloor et al. (2010) Gloor, M., J. L. Sarmienti, and N. Gruber (2010).  What can be learned about carbon cycle climate feedbacks from the CO2 airborne fraction?  Atmospheric Chemistry and Physics 10, 7739 – 7751.
  * Granger and Newbold (1974) Granger, C. W. J. and P. Newbold (1974).  Spurious regression in econometrics.  Journal of Econometrics 2, 111 – 120.
  * Hamilton (1994) Hamilton, J. D. (1994).  Time Series Analysis.  Princeton University Press.
  * Hill (1999) Hill, M. R. (1999).  Carbon dioxide emissions from the Russian Federation – problems and choices.  Energy & Environment 10(1), 51–78.
  * Houghton and Nassikas (2017) Houghton, R. A. and A. A. Nassikas (2017).  Global and regional fluxes of carbon from land use and land cover change 1850-2015.  Global Biogeochemical Cycles 31, 456 – 472.
  * Jarque and Bera (1987) Jarque, C. M. and A. K. Bera (1987).  A test for normality of observations and regression residuals.  International Statistical Review 2, 163–172.
  * Karlin and Taylor (1975) Karlin, S. and H. M. Taylor (1975).  A First Course in Stochastic Processes (Second Edition ed.).  Academic Press.
  * Keenan et al. (2016) Keenan, T. F., I. C. Prentice, J. G. Canadell, C. A. Williams, H. Wang, M. Raupach, and G. J. Collatz (2016).  Recent pause in the growth rate of atmospheric CO2 due to enhanced terrestrial carbon uptake.  Nature communications 7(1), 1–10.
  * Knorr (2009) Knorr, W. (2009).  Is the airborne fraction of anthropogenic CO2 emissions increasing?  Geophysical Research Letters 36.
  * Kuczmaszewska (2005) Kuczmaszewska, A. (2005).  The strong law of large numbers for dependent random variables.  Statistics & Probability Letters 73(3), 305–314.
  * Le Quéré et al. (2009) Le Quéré, C., M. R. Raupach, J. G. Canadell, G. Marland, L. Bopp, P. Ciais, T. J. Conway, S. C. Doney, R. A. Feely, P. Foster, P. Friedlingstein, K. Gurney, R. A. Houghton, J. I. House, C. Huntingford, P. E. Levy, M. R. Lomas, J. Majkut, N. Metzl, J. P. Ometto, G. P. Peters, I. C. Prentice, J. T. Randerson, S. W. Running, J. L. Sarmiento, U. Schuster, S. Sitch, T. Takahashi, N. Viovy, G. R. van der Werf, and F. I. Woodward (2009).  Trends in the sources and sinks of carbon dioxide.  Nature Geoscience 2, 831 – 836.
  * Newey and West (1987) Newey, W. K. and K. D. West (1987).  A simple, positive semi-definite, heteroskedasticity and autocorrelationconsistent covariance matrix.  Econometrica 55(394), 703–708.
  * Park (1992) Park, J. Y. (1992).  Canonical cointegrating regressions.  Econometrica 60(1), 119–143.
  * Phillips and Durlauf (1986) Phillips, P. C. B. and S. N. Durlauf (1986).  Multiple time series regression with integrated processes.  The Review of Economic Studies 53(4), 473–495.
  * Phillips and Park (1988) Phillips, P. C. B. and J. Y. Park (1988).  Asymptotic equivalence of ordinary least squares and generalized least squares in regressions with integrated regressors.  Journal of the American Statistical Association 83(401), 111–115.
  * Piegorsch and Casella (1985) Piegorsch, W. W. and G. Casella (1985).  The existence of the first negative moment.  The American Statistician 39(1), 60–62.
  * Pressburger et al. (2023) Pressburger, L., K. Dorheim, T. F. Keenan, H. McJeon, S. J. Smith, and B. Bond-Lamberty (2023).  Quantifying airborne fraction trends and the destination of anthropogenic CO2 by tracking carbon flows in a simple climate model.  Environmental Research Letters 18(5).
  * Raupach (2013) Raupach, M. R. (2013).  The exponential eigenmodes of the carbon-climate sytem, and their implications for ratios of responses to forcings.  Earth System Dynamics 4, 31 – 49.
  * Raupach et al. (2008) Raupach, M. R., J. G. Canadell, and C. L. Quéré (2008).  Anthropogenic and biophysical contributions to increasing atmospheric CO2 growth rate and airborne fraction.  Biogeosciences 5, 1601 – 1613.
  * Raupach et al. (2014) Raupach, M. R., M. Gloor, J. L. Sarmiento, J. G. Canadell, T. L. Frölicher, T. Gasser, R. A. Houghton, C. Le Quéré, and C. M. Trudinger (2014).  The declining uptake rate of atmospheric CO2 by land and ocean sinks.  Biogeosciences 11(13), 3453–3475.
  * Riahi et al. (2022) Riahi, K., R. Schaeffer, J. Arango, K. Calvin, C. Guivarch, T. Hasegawa, K. Jiang, E. Kriegler, R. Matthews, G. Peters, A. Rao, S. Robertson, A. Sebbit, J. Steinberger, M. Tavoni, and D. van Vuuren (2022).  Mitigation pathways compatible with long-term goals.  In P. Shukla, J. Skea, R. Slade, A. A. Khourdajie, R. van Diemen, D. McCollum, M. Pathak, S. Some, P. Vyas, R. Fradera, P. Shukla, J. Skea, R. Slade, A. A. Khourdajie, R. van Diemen, D. McCollum, M. Pathak, S. Some, P. Vyas, R. Fradera, M. Belkacemi, A. Hasija, G. Lisboa, S. Luz, and J. Malley (Eds.), IPCC, 2022: Climate Change 2022: Mitigation of Climate Change. Contribution of Working Group III to the Sixth Assessment Report of the Intergovernmental Panel on Climate Change. Cambridge University Press, Cambridge, UK and New York.
  * Rockström et al. (2017) Rockström, J., O. Gaffney, J. Rogelj, M. Meinshausen, N. Nakicenovic, and H. J. Schellnhuber (2017).  A roadmap for rapid decarbonization.  Science 355(6331), 1269–1271.
  * Roe and Baker (2007) Roe, G. H. and M. B. Baker (2007).  Why is climate sensitivity so unpredictable?  Science 318(5850), 629–632.
  * Schmith et al. (2012) Schmith, T., S. Johansen, and P. Thejll (2012).  Statistical analysis of global surface temperature and sea level using cointegration methods.  Journal of Climate 25(22), 7822–7833.
  * Siegenthaler and Oeschger (1978) Siegenthaler, U. and H. Oeschger (1978).  Predicting future atmospheric carbon dioxide levels.  Science 199(4327), 388–395.
  * van Marle et al. (2022) van Marle, M. J. E., D. van Wees, R. A. Houghton, R. D. Field, J. Verbesselt, and G. R. van der Werf (2022).  New land-use-change emissions suggest a declining CO2 airborne fraction.  Nature 603, 450–454.

##  Appendix A Theoretical analysis of the two estimators

Consider the cointegrated system

$\displaystyle y_{t}$ $\displaystyle=\alpha z_{t}+u_{t},$ $\displaystyle
z_{t}$ $\displaystyle=z_{0}+bt+x_{t},$ $\displaystyle x_{t}$
$\displaystyle=\sum_{\tau=1}^{t}\xi_{t},$

where $z_{0},b\in\mathbb{R}$, $u_{t}$ and $\xi_{t}$ are covariance-stationary
error terms, and $t=1,2,\ldots,T$. In the context of the analysis of the AF,
we would have $y_{t}=G_{t}$ and $z_{t}=E_{t}$. Recall the two estimators of
$\alpha$:

$\displaystyle\hat{\alpha}_{1}$
$\displaystyle=\frac{1}{T}\sum_{t=1}^{T}\frac{y_{t}}{z_{t}},$
$\displaystyle\hat{\alpha}_{2}$
$\displaystyle=\left(\sum_{t=1}^{T}z_{t}^{2}\right)^{-1}\sum_{t=1}^{T}z_{t}y_{t}.$

The second estimator $\hat{\alpha}_{2}$ has been thoroughly studied in the
literature on cointegration. If $b=0$, i.e. if there is no linear trend, then
the stochastic trend $x_{t}$ in $z_{t}$ implies that
$T(\hat{\alpha}_{2}-\alpha)$ has a non-Gaussian limiting distribution
(Phillips and Park, 1988). If $b\neq 0$, i.e. if there is a linear trend
present in $z_{t}$, then the linear trend dominates the stochastic trend,
which implies that $T^{3/2}(\hat{\alpha}_{2}-\alpha)$ has a Gaussian limit.
These results are well-known, see e.g. Park (1992), but for convenience we
state the result relevant for the specific setting in this paper.

######  Proposition 1 (Park (1992)).

Let $b\neq 0$ and suppose that $\nu_{t}=(u_{t},\xi_{t})^{\prime}$ is
stationary and ergodic with zero-mean. Define
$B_{n}(r):=n^{-1/2}\sum_{i=1}^{\lfloor nr\rfloor}\nu_{i}$ for $r\in[0,1]$.
Suppose that $B_{n}\stackrel{{\scriptstyle d}}{{\rightarrow}}B$ such that $B$
is a Brownian motion with $Var(B_{1})$ positive definite. Then

$\displaystyle T^{3/2}(\hat{\alpha}_{2}-\alpha)\stackrel{{\scriptstyle
d}}{{\rightarrow}}N(0,\sigma_{2}^{2}),$

as $T\rightarrow\infty$, where $\sigma^{2}_{2}>0$.

###### Proof.

The assumptions made on $\nu_{t}$ ensure that Assumption 2.1 in Park (1992) is
satisfied. Since $b\neq 0$, Assumption 2.2 in Park (1992) also holds. The
result now follows from Lemma 3.1 and Corollary 3.2 in Park (1992). ∎

######  Remark 1.

When $u_{t}$ is independent and identically distributed (iid), usual least-
squares inference is valid, meaning that $Var(\hat{\alpha}_{2})$ can be
estimated by the usual least-squares estimator

$\displaystyle\widehat{Var}(\hat{\alpha}_{2})=\left(\sum_{t=1}^{T}z_{t}^{2}\right)^{-1}T^{-1}\sum_{t=1}^{T}\left(y_{t}-\hat{\alpha}_{2}z_{t}\right)^{2}.$

If $u_{t}$ may contain autocorrelation and heteroskedasticity, we can instead
estimate $Var(\hat{\alpha}_{2})$ using a heteroskedasticity and
autocorrelation consistent (HAC) estimator, which is what we do in this paper
(Newey and West, 1987).

######  Remark 2.

The conditions imposed on $\nu_{t}$ in Proposition 1 are rather weak and are
satisfied for a large class of models. For instance, they are clearly
satisfied if $u_{t}$ and $\xi_{t}$ are independent stationary and ergodic
sequences. We refer to Phillips and Durlauf (1986) and references therein for
various specific conditions ensuring the condition on $\nu_{t}$ holds.

The properties of the estimator $\hat{\alpha}_{1}$ are less well known.
Indeed, to the best of our knowledge, the statistical properties of
$\hat{\alpha}_{1}$ have not been studied previously. To this end, we write

$\displaystyle\hat{\alpha}_{1}=\frac{1}{T}\sum_{t=1}^{T}\frac{y_{t}}{z_{t}}=\alpha+\frac{1}{T}\sum_{t=1}^{T}\frac{u_{t}}{z_{0}+bt+x_{t}}=\alpha+T^{-1}\sum_{t=1}^{T}U_{t},$

where $U_{t}=\frac{u_{t}}{z_{0}+bt+x_{t}}$. We first show that
$\hat{\alpha}_{1}$ is a superconsistent estimator of $\alpha$.

######  Proposition 2.

Let $b\neq 0$ and assume that $u_{t}$ is iid with $E(u_{1})=0$ and
$Var(u_{1})<\infty$, and that $x_{t}$ obeys a strong law of large numbers,
i.e. $t^{-1}x_{t}\stackrel{{\scriptstyle a.s}}{{\rightarrow}}0$ as
$t\rightarrow\infty$. Then, for all $\gamma>0$, it holds that

$\displaystyle T^{1-\gamma}(\hat{\alpha}_{1}-\alpha)\stackrel{{\scriptstyle
p}}{{\rightarrow}}0,$

as $T\rightarrow\infty$.

###### Proof.

Let $\epsilon,\delta>0$ be given. Suppose, for ease of notation, that $b>0$.
The proof when $b<0$ follows similar arguments. We want to show that

$\displaystyle
P(T^{1-\gamma}|\hat{\alpha}_{1}-\alpha|\geq\epsilon)\leq\delta,$ (5)

for all $T$ sufficiently large.

By assumption $x_{t}$ obeys a strong law of large numbers. By the continuous
mapping theorem, the same holds for $|x_{t}|$, meaning that $P(A)=1$, where

$\displaystyle
A=\\{\omega\in\Omega:\lim_{t\rightarrow\infty}t^{-1}|x_{t}(\omega)|=0\\}.$

Let $T\geq 1$ and define the random variable

$\displaystyle\tilde{x}_{T}:=\sup_{t\geq T}t^{-1}|x_{t}|.$

For all $\omega\in A$ we have
$\lim_{t\rightarrow\infty}t^{-1}|x_{t}(\omega)|=0$ and thus
$\limsup_{t\rightarrow\infty}t^{-1}|x_{t}(\omega)|=\lim_{t\rightarrow\infty}t^{-1}|x_{t}(\omega)|=0$.
By the definition of limes superior, we get

$\displaystyle\lim_{T\rightarrow\infty}\tilde{x}_{T}(\omega):=\lim_{T\rightarrow\infty}\sup_{t\geq
T}t^{-1}|x_{t}(\omega)|=\limsup_{t\rightarrow\infty}t^{-1}|x_{t}(\omega)|=\lim_{t\rightarrow\infty}t^{-1}|x_{t}(\omega)|=0.$

Since $P(A)=1$, we conclude that $\tilde{x}_{T}\stackrel{{\scriptstyle
a.s.}}{{\rightarrow}}0$ as $T\rightarrow\infty$, which implies that

$\displaystyle\sup_{t\geq T}t^{-1}|x_{t}|=\tilde{x}_{T}\stackrel{{\scriptstyle
p}}{{\rightarrow}}0$ (6)

as $T\rightarrow\infty$. Let $\nu\in(0,1)$ and $T\geq 1$, and define the set

$\displaystyle A_{T}=\left\\{\sup_{t\geq T}t^{-1}|x_{t}|\leq b\nu/2\right\\}.$

The convergence in (6) implies that there exists a $T_{1}\geq 1$ such that
$P(A_{T})\geq 1-\delta/2$ for all $T\geq T_{1}$. Likewise, since
$z_{0}\in\mathbb{R}$ is a constant, there exists a $T_{2}\geq 1$ such that
$\frac{|z_{0}|}{T}\leq b\nu/2$. Set $T_{3}=\max(T_{1},T_{2})$, and write

$\displaystyle T^{1-\gamma}(\hat{\alpha}_{1}-\alpha)$
$\displaystyle=T^{-\gamma}\sum_{t=1}^{T}\frac{u_{t}}{z_{0}+bt+x_{t}}$
$\displaystyle=T^{-\gamma}\sum_{t=1}^{T_{3}}\frac{u_{t}}{z_{0}+bt+x_{t}}+T^{-\gamma}\sum_{t=T_{3}+1}^{T}\frac{u_{t}}{z_{0}+bt+x_{t}}$
$\displaystyle=T^{-\gamma}B_{1,T_{3}}+T^{-\gamma}B_{T_{3}+1,T},$

where $B_{s,r}:=\sum_{t=s}^{r}\frac{u_{t}}{z_{0}+bt+x_{t}}$. Since $T_{3}$ is
fixed there are only finitely many terms in $B_{1,T_{3}}$, and hence
$T^{-\gamma}B_{1,T_{3}}\stackrel{{\scriptstyle p}}{{\rightarrow}}0$ as
$T\rightarrow\infty$. For $B_{T_{3}+1,T}$, we write

$\displaystyle T^{-\gamma}|B_{T_{3}+1,T}|$
$\displaystyle=T^{-\gamma}\left|\sum_{t=T_{3}+1}^{T}\frac{u_{t}}{z_{0}+bt+x_{t}}\right|$
$\displaystyle=T^{-\gamma}\left|\sum_{t=T_{3}+1}^{T}\frac{u_{t}}{bt}\frac{1}{1+z_{0}b^{-1}t^{-1}+x_{t}b^{-1}t^{-1}}\right|$
$\displaystyle\leq
T^{-\gamma}\sum_{t=T_{3}+1}^{T}\frac{|u_{t}|}{bt}\frac{1}{|1+z_{0}b^{-1}t^{-1}+x_{t}b^{-1}t^{-1}|}.$

Assuming that we are on $A_{T_{3}}$, we therefore have

$\displaystyle T^{-\gamma}|B_{T_{3}+1,T}|\leq
T^{-\gamma}\sum_{t=T_{3}+1}^{T}\frac{|u_{t}|}{bt}\frac{1}{|1-\nu|}.$

This last term converges to zero in $L^{2}$. To see this, note that

$\displaystyle
Var\left(T^{-\gamma/2}\sum_{t=T_{3}+1}^{T}\frac{|u_{t}|}{bt}\frac{1}{|1-\nu|}\right)$
$\displaystyle=T^{-\gamma}\sum_{t=T_{3}+1}^{T}\frac{Var(u_{t})}{b^{2}t^{2}}\frac{1}{(1-\nu)^{2}}$
$\displaystyle\sim
T^{-\gamma}Var(u_{1})\sum_{t=T_{3}+1}^{\infty}\frac{1}{b^{2}t^{2}}\frac{1}{(1-\nu)^{2}},$

where the sum converges. This shows that
$T^{-\gamma/2}\sum_{t=T_{3}+1}^{T}\frac{|u_{t}|}{bt}\frac{1}{|1-\nu|}$
converges in $L^{2}$ to a constant and thus that
$T^{-\gamma}|B_{T_{3}+1,T}|\stackrel{{\scriptstyle L^{2}}}{{\rightarrow}}0$ as
$T\rightarrow\infty$. We conclude that
$T^{-\gamma}B_{T_{3}+1,T}\stackrel{{\scriptstyle p}}{{\rightarrow}}0$ as
$T\rightarrow\infty$. Since this holds on $A_{T_{3}}$ and $P(A_{T_{3}})\geq
1-\delta/2$, then, by choosing $T\geq T_{3}$ large enough, (5) holds, which
completes the proof.

∎

######  Remark 3.

It is straightforward to relax the assumption that $u_{t}$ is iid, made in
Proposition 2. In practice, we use a HAC estimator to estimate the variance of
our estimators of $\hat{\alpha}_{2}$, to control for possible autocorrelation
and heteroskedasticity in $u_{t}$.

######  Remark 4.

The assumption in Proposition 2 that $x_{t}=\sum_{i=1}^{t}\xi_{i}$ obeys a
strong law of large numbers is satisfied under quite mild conditions. For
instance, it is satisfied if $\xi_{t}$ is iid with finite variance
(Kolmogorov’s Law of Large Numbers), $\xi_{t}$ is stationary and ergodic (the
Ergodic Theorem, e.g. Karlin and Taylor, 1975, Theorem 9.5.5, p. 487),
$\xi_{t}$ is a mixingale (e.g. de Jong, 1995; Davidson and de Jong, 1997), or
under certain mixing conditions on $\xi_{t}$ (e.g. Fazekas and Klesov, 2001;
Kuczmaszewska, 2005).

Although $\hat{\alpha}_{1}$ is a superconsistent estimator of $\alpha$, under
what are arguably mild conditions, it is the case that $U_{t}$, and
consequently $\hat{\alpha}_{1}$, does not possess any moments. This is, for
example, the case if the denominator $z_{0}+bt+x_{t}$ is a continuous random
variable with positive probability density at zero. In such circumstances,
$E(\hat{\alpha}_{1})$ and $V(\hat{\alpha}_{1})$ are not well-defined and
talking about whether the estimator is, e.g., unbiased is meaningless.
Although the results of non-existence of moments is well-known (e.g. Piegorsch
and Casella, 1985), it may not be immediately obvious. We therefore state it
as a proposition and provide a straightforward proof.

######  Proposition 3.

Let $X$ be a random variable with continuous probability density function
$f(x)$ such that $f(0)>0$. Then $E(X^{-\beta})$ does not exist for any
$\beta\geq 1$.

###### Proof.

Note that because $f$ is continuous and has $f(0)>0$, then there exists an
$\epsilon>0$ such that $f(x)>0$ for all $x\in[0,\epsilon]$. (If $X$ cannot
take positive values, then a similar arguments holds, but for
$x\in[-\epsilon,0]$.) Let $A:=\inf_{x\in[0,\epsilon]}f(x)>0$ and write

$\displaystyle
E(|X|^{-\beta})=\int_{-\infty}^{\infty}|x|^{-\beta}f(x)dx\geq\int_{0}^{\epsilon}|x|^{-\beta}f(x)dx\geq
A\int_{0}^{\epsilon}x^{-\beta}dx=\infty$

when $\beta\geq 1$. ∎

######  Remark 5.

Statistical issues with studying a ratio with a stochastic variable in the
denominator have been encountered before in climate science. Roe and Baker
(2007) discuss the “long tails” of the probability density function (PDF) for
the climate sensitivity, i.e. the temperature response to a doubling of
atmospheric CO2 over its preindustrial value. They formulate the problem as
$\Delta T=\lambda\cdot\Delta R_{f}$, where $\Delta T$ is the climate
sensitivity and $\Delta R_{f}$ is the change in downward radiative flux,
resulting from a doubling of CO2 concentrations. The unknown parameter
$\lambda$ is given by $\lambda_{0}/(1-f)$, where $\lambda_{0}$ is a “reference
climate sensitivity” and $f$ is the “total feedback factor”. It is clear that
if $f$ is a random variable such that the PDF of $f$ is continuous PDF with
positive density at $1$, then we are in a situation akin to the one studied in
this paper and in Proposition 3. In fact, Roe and Baker (2007) rely on a
number of published studies on the climate sensitivity and feedbacks in the
climate system and argue that such a situation is likely. In their
calculations they assume that the distribution of $f$ is Gaussian, which
satisfies the conditions of a continuous PDF with positive density at $1$.
This leads to very long tails in the PDF of the climate sensitivity $\Delta T$
and, as our Proposition 3 shows, non-existence of moments.

In the setting of this paper, we find that the deterministic term $z_{0}+bt$
is quite large compared to the stochastic term $x_{t}$, implying that any
probability density of the denominator $z_{0}+bt+x_{t}$ at zero is likely very
small. Hence, we proceed to study the estimator $\hat{\alpha}_{1}$ under the
assumption that $x_{t}=0$, i.e. that $z_{t}=z_{0}+bt$ is a deterministic
linear time trend. This avoids the problem of dividing by a random variable
and thus of having a positive probability density at zero for the term in the
denominator. Assume, for simplicity, that $u_{t}$ is an iid sequence. In this
case, we find that

$\displaystyle E(\hat{\alpha}_{1})=\alpha,$

and

$\displaystyle
V(\hat{\alpha}_{1})=T^{-2}\sum_{t=1}^{T}E(U_{t}^{2})=T^{-2}\sigma_{u}^{2}\sum_{t=1}^{T}\frac{1}{(z_{0}+bt)^{2}}\sim
T^{-2}\sigma_{u}^{2}c,$

as $T\rightarrow\infty$, where $c=\sum_{t=1}^{\infty}\frac{1}{(z_{0}+bt)^{2}}$
is a finite constant. Hence, in this case, we see that $\hat{\alpha}_{1}$ is
an unbiased estimator of $\alpha$, converging at the rate $T^{-1}$. However,
the estimator does not obey a Central Limit Theorem. To see this, write

$\displaystyle
T(\hat{\alpha}_{1}-\alpha)=\sum_{t=1}^{T}\frac{u_{t}}{z_{0}+bt},$ (7)

and observe that $\frac{u_{t}}{z_{0}+bt}\stackrel{{\scriptstyle
p}}{{\rightarrow}}0$ as $t\rightarrow\infty$, meaning that the first terms in
the sum in (7) dominate the entire sum, thus making central limit theorem-type
results invalid. That is, we cannot conclude that $\hat{\alpha}_{1}$ is
(asymptotically) Gaussian unless we impose an additional assumption that the
error term $u_{t}$ is itself Gaussian. For convenience, we again state these
results formally.

######  Proposition 4.

Let $b\neq 0$ and assume that $u_{t}$ is iid with $E(u_{1})=0$ and
$Var(u_{1})<\infty$ and consider the sequence of random variables for
$C_{T}=\sum_{t=1}^{T}\frac{u_{t}}{z_{0}+bt}$, $T=1,2,\ldots$. The following
holds:

  1. (i)

A Central Limit Theorem does not hold in general for $C_{T}$, i.e. we cannot
conclude that $C_{T}$ has a Gaussian limit as $T\rightarrow\infty$.

  2. (ii)

If $u_{t}\sim N(0,\sigma_{u}^{2})$ with $\sigma_{u}^{2}>0$, then

$\displaystyle C_{T}\sim N(0,\sigma_{1,T}^{2}),$

where $\sigma_{1,T}^{2}=\sigma_{u}^{2}\sum_{t=1}^{T}(z_{0}+bt)^{-2}$, and, as
$T\rightarrow\infty$,

$\displaystyle C_{T}\stackrel{{\scriptstyle
d}}{{\rightarrow}}N(0,\sigma_{1}^{2}),$

where $\sigma_{1}^{2}=\sigma_{u}^{2}\sum_{t=1}^{\infty}(z_{0}+bt)^{-2}$.

###### Proof.

(i) To see that a Central Limit Theorem cannot hold in general, suppose that
$u_{t}$ is such that $E(u_{t}^{3})=s\neq 0$. Then

$\displaystyle
E(C_{T}^{3})=\sum_{t=1}^{T}\frac{E(u_{t}^{3})}{(z_{0}+bt)^{3}}=s\sum_{t=1}^{T}\frac{1}{(z_{0}+bt)^{3}}\rightarrow
sc_{3}\neq 0,$

as $T\rightarrow\infty$, where
$c_{3}=\sum_{t=1}^{\infty}\frac{1}{(z_{0}+bt)^{3}}<\infty$ and $c_{3}\neq 0$.
This shows that the asymptotic skewness of $C_{T}$ is not zero, implying that
the limiting distribution of $C_{T}$ is not Gaussian.

(ii) That $C_{T}\sim N(0,\sigma_{1,T}^{2})$ follows trivially by the
assumption that $u_{t}$ is iid Gaussian. The same holds for the second result.
Formally, we may deduce it using characteristic functions. For
$\theta\in\mathbb{R}$, we have

$\displaystyle E\left(e^{i\theta
C_{T}}\right)=\prod_{t=1}^{T}E\left(e^{i\frac{\theta}{z_{0}+bt}u_{t}}\right)=\prod_{t=1}^{T}e^{-\frac{1}{2}\frac{\theta^{2}\sigma_{u}^{2}}{(z_{0}+bt)^{2}}}=e^{-\frac{1}{2}\theta^{2}\sigma_{u}^{2}\sum_{t=1}^{T}\frac{1}{(z_{0}+bt)^{2}}}\rightarrow
e^{-\frac{1}{2}\theta^{2}\sigma_{1}^{2}},$

as $T\rightarrow\infty$, which concludes the proof. ∎

In the case of the conventional estimator $\hat{\alpha}_{1}$, Proposition 4
shows that if we want to base inference on the normal distribution, we have to
make distributional assumptions on the error term $u_{t}$. Indeed, if we are
willing to assume that $u_{t}\sim N(0,\sigma_{u}^{2})$, then we get

$\displaystyle T(\hat{\alpha}_{1}-\alpha)\sim N(0,\sigma_{1,T}^{2}),$

and, as $T\rightarrow\infty$,

$\displaystyle T(\hat{\alpha}_{1}-\alpha)\stackrel{{\scriptstyle
d}}{{\rightarrow}}N(0,\sigma_{1}^{2}),$

where $\sigma_{1,T}$ and $\sigma_{1}^{2}$ are as in Proposition 4(ii). Note,
if $u_{t}$ is iid, then $\sigma_{1,T}^{2}$ and $\sigma_{1}^{2}$ can be
estimated straightforwardly from the sample variance of $y_{t}/z_{t}$. If
$u_{t}$ may contain autocorrelation and heteroskedasticity, we can instead
estimate $\sigma_{1,T}^{2}$ and $\sigma_{1}^{2}$ using a HAC estimator, which
is what we do in this paper.

##  Appendix B Further empirical results

###  B.1 Unit root and cointegration tests

Let $y_{t}$, $t=1,2,\ldots,T$, be a time series and $L\geq 0$ an integer. We
consider the following null model for $y_{t}$,

$\displaystyle y_{t}=a+y_{t-1}+\sum_{l=1}^{L}\beta_{l}\Delta
y_{t-l}+\epsilon_{t},$ (8)

where $a,\beta_{1},\ldots,\beta_{L}\in\mathbb{R}$ and $\epsilon_{t}$ is an
error term. We use the convention that $\sum_{l=1}^{L}=0$ if $L=0$. The
alternative model for $y_{t}$ is

$\displaystyle y_{t}=b+c\cdot t+\phi y_{t-1}+\sum_{l=1}^{L}\beta_{l}\Delta
y_{t-l}+\epsilon_{t},$ (9)

where $b,c\in\mathbb{R}$ and $\phi\in(-1,1)$.

We use the Augmented Dickey-Fuller test (ADF; Dickey and Fuller, 1979) to test
the null ($H_{0}$) that a data series is generated by the model in (8) against
the alternative ($H_{1}$) that the data series is generated by (9). Depending
on restrictions made on the three parameters $a,b,c$, this leads to three
different tests:

  1. 1.

AR: $a=b=c=0$. This is the test for $y_{t}$ having a unit root ($H_{0}$)
against the alternative that $y_{t}$ is an AR process.

  2. 2.

ARD: $a=c=0$. This is the test for $y_{t}$ having a unit root ($H_{0}$)
against the alternative that $y_{t}$ is an AR process with drift.

  3. 3.

TS: This is the test for $y_{t}$ having a unit root and a drift ($H_{0}$)
against the alternative that $y_{t}$ is an AR process around a linear trend
(i.e. that $y_{t}$ is “trend-stationary”).

We implement the ADF test using the built-in function adftest in the MATLAB
Econometrics Toolbox. The $p$-values for testing $H_{0}$ against $H_{1}$ for
these three tests and the two data series $y_{t}=G_{t}$ and $y_{t}=E_{t}$ are
shown in the top six rows of Table 2. For $y_{t}=E_{t}$ the $p$-values are
very large regardless of the particular test conducted (AR, ARD, or TS) and
the number of lags $L$ used. Hence, there is strong evidence that $E_{t}$
contains a stochastic trend. For $y_{t}=G_{t}$ the $p$-values are very large
for the AR-test, regardless of the number of lags $L$ used. For the ARD-test,
the $p$-values are below $10\%$ for $L=0$ and $L=1$ and very large for $L\geq
2$. For the TS-test, the $p$-values are below $10\%$ regardless of the number
of lags $L$ used. Thus, the AR-test and ARD-test provide evidence that the
time series of atmospheric CO2 changes $G_{t}$ contains a stochastic trend,
while the TS-test indicates that the trend might be modelled using a
deterministic trend $b+c\cdot t$. In summary, the ADF tests provide very
strong evidence for the presence of trends in $E_{t}$ and $G_{t}$ and suggests
that both series are non-stationary.

We next employ the Engle-Granger test (Engle and Granger, 1987) to test
whether $G_{t}$ and $E_{t}$ are cointegrated. Formally, we consider the model

$\displaystyle G_{t}=\alpha E_{t}+u_{t}$ (10)

and test whether the error term $u_{t}$ is stationary. The parameter $\alpha$
is estimated by least-squares regression and the estimated residuals are
computed,

$\displaystyle\widehat{u}_{t}=G_{t}-\widehat{\alpha}E_{t},$

where $\widehat{\alpha}$ is the least-squares estimate of $\alpha$ obtained
from (10). In a second step, we investigate whether the error terms $u_{t}$
are stationary by performing the AR-version of the ADF test on the estimated
residuals $\widehat{u}_{t}$. The Engle-Granger test adjusts the critical
values of the ADF test to account for the fact that the error term $u_{t}$ is
unobserved and that the estimated residuals $\widehat{u}_{t}$ are used in
their place as data in the ADF test. We implement the Engle-Granger test using
the built-in function egcitest in the MATLAB Econometrics Toolbox. The
$p$-values for the Engle-Granger test are shown in the bottom row in Table 2.
The null hypothesis of a unit root in $u_{t}$ is firmly rejected for all lags
$L$ and we can therefore conclude that the two time series are cointegrated.

###  B.2 Analysis of alternative data sets

We use data on yearly changes in atmospheric CO2 ($G_{t}$), yearly CO2
emissions from fossil fuels ($E_{t}^{FF}$), and yearly CO2 emissions from
land-use and land cover change ($E_{t}^{LULCC}$). Total anthropogenic CO2
emissions are then $E_{t}=E_{t}^{FF}+E_{t}^{LULCC}$. The uncertainty in
measurements of the AF stems in large part from uncertainties in the magnitude
of LULCC emissions (van Marle et al., 2022). We therefore follow van Marle et
al. (2022) and estimate the AF using three different data sets for LULCC
emissions (GCP, H&N, vMa), while data for atmospheric CO2 changes and CO2
emissions from fossil fuels are obtained from the most recent edition of the
Global Carbon Budget data set. (Friedlingstein et al., 2022).333Data from the
Global Carbon Project can be found at https://www.icos-cp.eu/science-and-
impact/global-carbon-budget/2022, last accessed April 18, 2023. Data for the
various LULCC time series can be found at
https://doi.org/10.5281/zenodo.5617953, last accessed April 18, 2023. The GCP
LULCC data are from the Global Carbon Project Friedlingstein et al. (2022),
the H&N LULCC data are from Houghton and Nassikas (2017), and the vMa LULCC
data are from van Marle et al. (2022). We thus have three different data sets,
only differing in the LULCC time series. The data are available for the period
$1959$–$2021$ (Figure 3).

Estimation results for $\alpha$ using the alternative data sets are reported
in Table 3. For ease of comparison, we also include the results presented in
Table 1 using the GCP data set. The table shows that the overall estimates of
$\alpha$ vary across the three data sets, with the GCP data generally
resulting in lower estimates of the AF than the H&N and vMa data. The
reduction in uncertainty obtained by using $\hat{\alpha}_{2}$ instead of
$\hat{\alpha}_{1}$ is similar across the three data sets, however. In Appendix
B.3 we re-analyze a later sub-sample (1992–2021) and find that in this case,
the estimates of $\alpha$ from the three data sets largely agree. In fact,
they are very similar to the estimate of $\alpha$ obtained from the GCP data
set using the entire sample 1959–2021.

###  B.3 Analysis of the subsample 1992–2021

In this section, we analyze the AF over the recent subsample 1992-2021. This
is motivated by the contradictory findings of trends in the AF mentioned in
the introduction, the possibility of a structural break in the data around
$1990$ (Bennedsen et al., 2023b), substantial disagreements in land-use and
land-cover change emissions data before $1992$ (van Marle et al., 2022, see
also Figure 3), and the presumption of better data quality in the latter part
of the sample (e.g. Hill, 1999).

The estimation results are shown in Table 4. The results show that the
estimates of $\alpha$ are very similar across the three different data sets
for the sub-sample 1992–2021. This is to be expected since the different data
sets agree much more closely on the latter part of the record than on the
earlier (Figure 3). The estimates of $\alpha$ obtained from the sub-sample
1992–2021 are very close to those obtained using the GCP data on the entire
1959–2021 sample.

We find that our new estimator $\hat{\alpha}_{2}$ improves the precision,
measured by the standard error of the estimator, by approximately $18\%$
compared to the conventionally used estimator $\hat{\alpha}_{1}$. Including
data on ENSO and volcanic activity in the analysis, our proposed estimator
$\hat{\alpha}_{4}$ improves the precision with approximately $47\%$, compared
to the conventional estimator of the AF $\hat{\alpha}_{1}$. Using data from
the Global Carbon Project (e.g. Friedlingstein et al., 2022), our best
estimate of the AF is $46.8\%$ with an associated standard error of $1.1\%$,
resulting in a $95\%$ confidence interval of $[44.7\%,48.8\%]$ for the AF. The
improvement in precision of our new estimator is even more pronounced in this
subsample.

##  Appendix C Simulation studies

###  C.1 Performance of the estimators in conditions similar to the data
sample $1959$–$2021$

We simulate data from the models in (1)–(2). To facilitate comparison, we set
the data-generating value of $\alpha$ to $0.4441$, which is the estimated
value obtained using the model in (1) and the GCP data set, see Table 1. The
variance of the error terms are set to $Var(\epsilon_{t}^{(1)})=0.1284^{2}$
and $Var(u_{t}^{(2)})=0.9153^{2}$, again corresponding to the estimates from
the GCP data (Table 1). Using these values, we simulate data
$y_{t}=G_{t}/E_{t}$ from model (1) directly. To simulate from model (2), we
first need to simulate an emissions trajectory. To this end, we specify
emissions as a random walk plus drift, i.e. $E_{t}=E_{t-1}+b+\xi_{t}$, where
$b$ is drift parameter (see the main paper as well as Bennedsen et al., 2023a,
for a justification of this model).

Using the GCP data for $1959$–$2021$, we estimate $b$ and $Var(\xi_{t})$ from
the model $\Delta E_{t}:=E_{t}-E_{t-1}=b+\xi_{t}$. We also set $E_{0}=4.3433$,
which is the first observation of $E_{t}$ from the GCP data. Then, we simulate
emissions time series using the dynamic equation
$E_{t}=E_{t-1}+\hat{b}+\xi_{t}$, $t=1,2,\ldots,T$, where the variance of
$\xi_{t}$ is set to the value $\widehat{Var}(\xi_{t})$ found on the GCP data.

We simulate $10^{6}$ instances of $G_{t}$ and $E_{t}$ for $t=1,2,\ldots,T$.
For each instance, we estimate $\alpha$ using the estimators
$\hat{\alpha}_{1}$ and $\hat{\alpha}_{2}$, and calculate the squared errors of
the resulting estimates. The root mean squared errors (RMSE) of the two
estimators are obtained by averaging the squared errors over the $10^{6}$
simulations. We run these simulation studies for various sample sizes, namely
$T=63,64,\ldots,142$. $T=63$ corresponds to the sample size relevant for the
data in the main paper (1959–2021), while $T=142$ corresponds to data until
$2100$ (1959–2100). The RMSEs of the two estimators, as functions of sample
size $T$, are shown in the left plot of Figure 4. The relative RMSE, i.e.
RMSE($\hat{\alpha}_{2}$)$/$RMSE($\hat{\alpha}_{1}$), as a function of $T$, is
shown in the right plot of Figure 4.

From the left plot of Figure 4, we see that the RMSE of both estimators are
decreasing in $T$. The decrease is much more pronounced for the estimator
$\hat{\alpha}_{2}$ compared to the estimator $\hat{\alpha}_{1}$. This is to be
expected, since we saw above that $\hat{\alpha}_{2}$ converges to $\alpha$ at
the rate $T^{3/2}$, while $\hat{\alpha}_{1}$ converges at the slower rate
$T^{1-\gamma}$, $\gamma>0$. The difference between the performance of the two
estimators is illustrated in the right plot of Figure 4, where the relative
RMSE is plotted. For $T=63$, corresponding to the left-most observation, the
RMSE of $\hat{\alpha}_{2}$ is approximately $10\%$ lower than the RMSE of
$\hat{\alpha}_{1}$. These simulation results align with the empirical
estimates, where it was found that $\hat{\alpha}_{2}$ was approximately $10\%$
more precise than $\hat{\alpha}_{1}$ (Tables 1 and 3). The relative
performance continues to diverge as the sample size gets larger; for $T=142$,
the reduction in RMSE is on the order of $50\%$.

###  C.2 Performance of the estimators in low emission scenarios

In this section, we perform a stylized simulation experiment to highlight the
issues that may arise when basing analyses of the AF on the ratio
$G_{t}/E_{t}$, and how these issues are not present when analyzing $G_{t}$ and
$E_{t}$ in a system, as we have proposed to do instead.

We consider a setup where emissions are linearly increasing in $T_{0}=63$
periods (representing $1959$–$2021$), followed by an exponential decrease for
$78$ periods (representing $2022$–$2100$). That is, we consider the following
model for emissions

$\displaystyle E_{t}=E_{t}^{*}+\xi_{t},$

where

$\displaystyle E_{t}^{*}=\begin{cases}E_{0}+b\cdot t&t=1,2,\ldots,T_{0},\\\
E_{T_{0}}\cdot(1+g)^{t-T_{0}}&t=T_{0}+1,T_{0}+1,\ldots,T,\end{cases}$

where $T=63+78=141$ (representing $1959$–$2100$). For a simulated emissions
path $E_{t}$, we simulate a corresponding path for the atmospheric CO2 changes
$G_{t}$ via

$\displaystyle G_{t}=\alpha\cdot E_{t}+u_{t},$

for $t=1,2,\ldots,T$. Note that this stylized model assumes that the AF
$\alpha$ is constant over the entire period $1,2,\ldots,T$ (representing
$1959$–$2100$). While this assumption is likely unrealistic in practice, since
we have theoretical reasons to think that the AF will decrease as CO2
emissions decrease (e.g. Raupach, 2013), it will allow us to illustrate the
issues with using the ratio $G_{t}/E_{t}$ in the context of the present paper.
We leave investigations of how our framework can be applied to situations
where the AF is time-varying, e.g. by analysing the system
$G_{t}=\alpha_{t}E_{t}+u_{t}$ where $\alpha_{t}$ is modelled using a time-
varying parameter model, for future research.

As above, we set $Var(\xi_{t})$ and $b$ equal to the estimated values found on
the GCP data. The exponential growth rate of emissions $g$ is set to
$g=-8.10\%$, ensuring that emissions $E_{t}$ is halved in $2030$, compared to
the $2010$ level, which has been found to be a necessary condition for
achieving the goals of the Paris Agreement (e.g. Rockström et al., 2017; Riahi
et al., 2022). Thus, the simulated data used in this section represent a
combination of the historical record ($1959$–$2021$) followed by an ambitious
effort to reduce CO2 emissions in line with international agreements
($2022$–$2100$). The AF $\alpha$ is set to the same value used in the previous
section, i.e. $\alpha=0.4441$, which is the AF found via $\hat{\alpha}_{1}$
applied to the GCP data (Table 1). Lastly, we set $Var(u_{t})=0.9153^{2}$
corresponding to the estimate of $Var(u_{t}^{(2)})$ from the GCP data (Table
1). Typical simulated trajectories of $E_{t}$ and $G_{t}$ resulting from this
setup are shown in Figure 5, along with the corresponding GCP data series for
$E_{t}$ and $G_{t}$ over $1959$–$2021$ for comparison.

Given simulated paths for emissions $E_{t}$ and atmospheric changes $G_{t}$,
we may estimate $\alpha$ using both the conventional estimator (1) and our
newly proposed estimator (2). Figure 6 shows the estimated moments
$E(|\hat{\alpha}-\alpha|)$ and $E(|\hat{\alpha}-\alpha|^{2})$, as a function
of the sample size (represented by the last year in the estimation period),
for both of these estimators using $10^{6}$ simulations of $(G_{t},E_{t})$,
$t=1,2,\ldots,T$. We consider sample sizes from $63$ (corresponding to data
over the period $1959$–$2021$) to $141$ (corresponding to data over the period
$1959$–$2100$). Figure 6 clearly illustrates how the non-existence of the
moments of $\hat{\alpha}_{1}$ becomes problematic as we enter into a low-
emission regime.444Note that because the moments do not exist in the case of
$\hat{\alpha}_{1}$, then the estimates of $E(|\hat{\alpha}-\alpha|)$ and
$E(|\hat{\alpha}-\alpha|^{2})$ based on $\hat{\alpha}_{1}$ do not converge as
the number of simulated instances is increased. Hence, different initial seeds
in the simulation exercise can lead to wildly different estimation results for
$\hat{\alpha}_{1}$, even though $10^{6}$ simulated instances are used in the
estimations. In contrast, since the moments of $\hat{\alpha}_{2}$ do exist,
the estimates of $E(|\hat{\alpha}-\alpha|)$ and $E(|\hat{\alpha}-\alpha|^{2})$
based on $\hat{\alpha}_{2}$ are very stable across initial seeds in the
simulation exercise. Figure 7 shows the RMSE of the two estimators, similar to
Figure 4 presented in the previous simulation study.

Table 2: $p$-values for the Augmented Dickey-Fuller (ADF) test for a unit root and the Engle-Granger test for cointegration between $G_{t}$ and $E_{t}$ on the GCP data (1959–2021), see Appendix B.1 for details. First six rows are ADF tests and last row is Engle-Granger test. “AR”: The null model is a unit root model while the alternative model is an AR model. “ARD”: The null model is a unit root model while the alternative model is an AR model with drift. “TS”: The null model is a unit root model with a constant (i.e. a drift in levels) while the alternative model is an AR model with a time trend. “Engle-Granger”: The null hypothesis is that the error term $u_{t}$ in $G_{t}=\alpha E_{t}+u_{t}$ has a unit root, i.e. that $G_{t}$ and $E_{t}$ are not cointegrated. The value “$L$” denotes the number of lags included into both the null and alternative models, see equations (8)–(9). | $L=0$ | $L=1$ | $L=2$ | $L=3$ | $L=4$ | $L=5$
---|---|---|---|---|---|---
$y_{t}=G_{t}$ (AR) | 0.2616 | 0.5123 | 0.6987 | 0.8266 | 0.8606 | 0.8702
$y_{t}=G_{t}$ (ARD) | 0.0045 | 0.0905 | 0.3604 | 0.4954 | 0.5074 | 0.4485
$y_{t}=G_{t}$ (TS) | 0.0010 | 0.0010 | 0.0036 | 0.0212 | 0.0238 | 0.0072
$y_{t}=E_{t}$ (AR) | 0.9990 | 0.9990 | 0.9990 | 0.9984 | 0.9915 | 0.9967
$y_{t}=E_{t}$ (ARD) | 0.9454 | 0.9365 | 0.9105 | 0.8696 | 0.8280 | 0.8897
$y_{t}=E_{t}$ (TS) | 0.1418 | 0.3016 | 0.4662 | 0.4300 | 0.3126 | 0.5520
Engle-Granger | 0.0010 | 0.0010 | 0.0023 | 0.0156 | 0.0181 | 0.0041

Figure 3: Data used in the study over the period $1959$–$2019$. The dashed
vertical black line denotes the year $1992$, the starting point for the
analysis in Appendix B.3.

Table 3: Least-squares regression output of the linear models in Equations (1)–(4) for the period 1959–2021 for the three data sets GCP, H&N, and vMa. “$\hat{\alpha}$” denotes the estimate of the AF $\alpha$, “SE($\hat{\alpha}$)” denotes the estimated standard error of $\hat{\alpha}$, and “Relative SE” denotes the ratio of the standard error to the standard error from $\hat{\alpha}_{1}$, corresponding to the left-most column. Hence, the left-most column has a Relative SE of one by definition, while a Relative SE in the other columns indicate a lower standard error for the estimator in question. “ $CI_{95\%}(\alpha)$ ” denotes a $95\%$ confidence interval for the AF $\alpha$, based on the Gaussian distribution. “$\widehat{SD}(u_{t})$” denotes the estimated standard deviation of the error term $u_{t}$ and $R^{2}$ is the coefficient of determination. “$\hat{\gamma}_{i}$”, $i=1,2$, denotes the estimates of the $\gamma_{i}$ parameter from Equations (3)–(4) and “SD($\hat{\gamma}_{i}$)” is their estimated standard deviation. Data: GCP |  (1) |  (2) |  (3) |  (4)
---|---|---|---|---
$\hat{\alpha}$ | 0.4441 | 0.4533 | 0.4785 | 0.4763
SE($\hat{\alpha}$) | 0.0163 | 0.0147 | 0.0129 | 0.0109
Relative SE | 1.0000 | 0.8986 | 0.7881 | 0.6712
$CI_{95\%}(\alpha)$ | $[0.4121,0.4761]$ | $[0.4246,0.4820]$ | $[0.4533,0.5037]$ | $[0.4548,0.4977]$
$\widehat{SD}(u_{t})$ | 0.1284 | 0.9153 | 0.0894 | 0.6327
$R^{2}$ | 0 | 0.5841 | 0.5310 | 0.8077
$\hat{\gamma}_{1}$ (ENSO) | - | - | 0.1410 | 1.0185
SD($\hat{\gamma}_{1}$) | - | - | 0.0175 | 0.1186
$\hat{\gamma}_{2}$ (VAI) | - | - | -2.3219 | -15.2890
SD($\hat{\gamma}_{2}$) | - | - | 0.2950 | 2.3676
Data: H&N |  (1) |  (2) |  (3) |  (4)
$\hat{\alpha}$ | 0.4693 | 0.4672 | 0.5032 | 0.4880
SE($\hat{\alpha}$) | 0.0173 | 0.0153 | 0.0148 | 0.0117
Relative SE | 1.0000 | 0.8841 | 0.8549 | 0.6789
$CI_{95\%}(\alpha)$ | $[0.4354,0.5032]$ | $[0.4373,0.4972]$ | $[0.4743,0.5322]$ | $[0.4650,0.5110]$
$\widehat{SD}(u_{t})$ | 0.1366 | 0.9264 | 0.0992 | 0.6790
$R^{2}$ | 0 | 0.5740 | 0.4895 | 0.7785
$\hat{\gamma}_{1}$ (ENSO) | - | - | 0.1465 | 0.9948
SD($\hat{\gamma}_{1}$) | - | - | 0.0211 | 0.1328
$\hat{\gamma}_{2}$ (VAI) | - | - | -2.2923 | -13.7762
SD($\hat{\gamma}_{2}$) | - | - | 0.4194 | 3.4531
Data: vMa |  (1) |  (2) |  (3) |  (4)
$\hat{\alpha}$ | 0.4854 | 0.4725 | 0.5175 | 0.4893
SE($\hat{\alpha}$) | 0.0183 | 0.0153 | 0.0179 | 0.0132
Relative SE | 1.0000 | 0.8331 | 0.9747 | 0.7186
$CI_{95\%}(\alpha)$ | $[0.4495,0.5214]$ | $[0.4426,0.5025]$ | $[0.4825,0.5525]$ | $[0.4635,0.5151]$
$\widehat{SD}(u_{t})$ | 0.1425 | 0.9246 | 0.1093 | 0.7192
$R^{2}$ | 0 | 0.5755 | 0.4302 | 0.7515
$\hat{\gamma}_{1}$ (ENSO) | - | - | 0.1452 | 0.9455
SD($\hat{\gamma}_{1}$) | - | - | 0.0211 | 0.1301
$\hat{\gamma}_{2}$ (VAI) | - | - | -2.1667 | -11.6429
SD($\hat{\gamma}_{2}$) | - | - | 0.5137 | 3.8757

Table 4: Least-squares regression output of the linear models in Equations (1)–(4) for the period 1992–2021 for the three data sets GCP, H&N, and vMa. “$\hat{\alpha}$” denotes the estimate of the AF $\alpha$, “SE($\hat{\alpha}$)” denotes the estimated standard error of $\hat{\alpha}$, calculated using the HAC estimator of Newey and West (1987), and “Relative SE” denotes the ratio of the standard error to the standard error from $\hat{\alpha}_{1}$, corresponding to the left-most column. Hence, the left-most column has a Relative SE of one by definition, while a Relative SE in the other columns indicate a lower standard error for the estimator in question. “ $CI_{95\%}(\alpha)$ ” denotes a $95\%$ confidence interval for the AF $\alpha$, based on the Gaussian distribution. “$\widehat{SD}(u_{t})$” denotes the estimated standard deviation of the error term $u_{t}$. “$\hat{\gamma}_{i}$”, $i=1,2$, denotes the estimates of the $\gamma_{i}$ parameter from Equations (3)–(4) and “SD($\hat{\gamma}_{i}$)” is their estimated standard deviation. Data: GCP |  (1) |  (2) |  (3) |  (4)
---|---|---|---|---
$\hat{\alpha}$ | 0.4507 | 0.4550 | 0.4694 | 0.4675
SE($\hat{\alpha}$) | 0.0199 | 0.0164 | 0.0128 | 0.0106
Relative SE | 1.0000 | 0.8224 | 0.6446 | 0.5324
$CI_{95\%}(\alpha)$ | $[0.4117,0.4898]$ | $[0.4229,0.4871]$ | $[0.4442,0.4946]$ | $[0.4467,0.4883]$
$\widehat{SD}(u_{t})$ | 0.1100 | 0.9458 | 0.0677 | 0.5956
$R^{2}$ | 0 | 0.3597 | 0.6471 | 0.7636
$\hat{\gamma}_{1}$ (ENSO) | - | - | 0.1183 | 1.0532
SD($\hat{\gamma}_{1}$) | - | - | 0.0179 | 0.1200
$\hat{\gamma}_{2}$ (VAI) | - | - | -2.2707 | -17.6307
SD($\hat{\gamma}_{2}$) | - | - | 0.1832 | 1.2699
Data: H&N |  (1) |  (2) |  (3) |  (4)
$\hat{\alpha}$ | 0.4548 | 0.4594 | 0.4739 | 0.4724
SE($\hat{\alpha}$) | 0.0201 | 0.0164 | 0.0129 | 0.0105
Relative SE | 1.0000 | 0.8184 | 0.6410 | 0.5244
$CI_{95\%}(\alpha)$ | $[0.4155,0.4941]$ | $[0.4273,0.4916]$ | $[0.4487,0.4991]$ | $[0.4518,0.4930]$
$\widehat{SD}(u_{t})$ | 0.1112 | 0.9463 | 0.0679 | 0.5907
$R^{2}$ | 0 | 0.3591 | 0.6529 | 0.7674
$\hat{\gamma}_{1}$ (ENSO) | - | - | 0.1196 | 1.0528
SD($\hat{\gamma}_{1}$) | - | - | 0.0177 | 0.1128
$\hat{\gamma}_{2}$ (VAI) | - | - | -2.3166 | -17.9201
SD($\hat{\gamma}_{2}$) | - | - | 0.1756 | 1.1641
Data: vMa |  (1) |  (2) |  (3) |  (4)
$\hat{\alpha}$ | 0.4535 | 0.4584 | 0.4709 | 0.4688
SE($\hat{\alpha}$) | 0.0192 | 0.0157 | 0.0132 | 0.0107
Relative SE | 1.0000 | 0.8151 | 0.6857 | 0.5583
$CI_{95\%}(\alpha)$ | $[0.4158,0.4912]$ | $[0.4277,0.4892]$ | $[0.4450,0.4967]$ | $[0.4478,0.4899]$
$\widehat{SD}(u_{t})$ | 0.1073 | 0.9064 | 0.0679 | 0.5886
$R^{2}$ | 0 | 0.4120 | 0.6272 | 0.7691
$\hat{\gamma}_{1}$ (ENSO) | - | - | 0.1159 | 1.0197
SD($\hat{\gamma}_{1}$) | - | - | 0.0211 | 0.1301
$\hat{\gamma}_{2}$ (VAI) | - | - | -2.1210 | -15.5501
SD($\hat{\gamma}_{2}$) | - | - | 0.2068 | 1.3084

Figure 4: Simulation study of the properties of the estimators
$\hat{\alpha}_{1}$ and $\hat{\alpha}_{2}$. See Appendix C.1 for details.

Figure 5: A typical simulated instance of an emissions trajectory (a) and a
corresponding trajectory for atmospheric CO2 changes (b) in a low-emission
regime (blue lines), see Appendix C.2 for details. For comparison, we show the
historical emissions (a) and atmospheric CO2 changes (b) from the GCP data
over $1959$–$2021$ (black lines).

Figure 6: Estimated moments of the two estimators $\hat{\alpha}_{1}$ and
$\hat{\alpha}_{2}$ using $10^{6}$ simulations in a low-emission regime, see
Appendix C.2 for details. The left-most data point ($2021$ on the x-axis) in
both panels corresponds to a sample size of $63$, while the right-most data
point ($2100$ on the x-axis) corresponds to a sample size of $141$.

Figure 7: Simulation study of the properties of the estimators
$\hat{\alpha}_{1}$ and $\hat{\alpha}_{2}$ in a low-emission regime. See
Appendix C.2 for details.