# A Statistical Reduced Complexity Climate Model for
Probabilistic Analyses and Projections

Mikkel Bennedsen, Eric Hillebrand, Siem Jan Koopman  Department of Economics
and Business Economics and CoRE, Aarhus University, Fuglesangs Allé 4, 8210
Aarhus V, Denmark. E-mail<EMAIL_ADDRESS>Department of Economics and
Business Economics and CoRE, Aarhus University, Fuglesangs Allé 4, 8210 Aarhus
V, Denmark. E-mail<EMAIL_ADDRESS>Department of Econometrics, School
of Business and Economics, Vrije Universiteit Amsterdam, De Boelelaan 1105,
1081 HV Amsterdam, The Netherlands, E-mail<EMAIL_ADDRESS>
###### Abstract

We propose a new statistical reduced complexity climate model. The centerpiece
of the model consists of a set of physical equations for the global climate
system which we show how to cast in non-linear state space form. The
parameters in the model are estimated using the method of maximum likelihood
with the likelihood function being evaluated by the extended Kalman filter.
Our statistical framework is based on well-established methodology and is
computationally feasible. In an empirical analysis, we estimate the parameters
for a data set comprising the period $1959$–$2022$. A likelihood ratio test
sheds light on the most appropriate equation for converting the level of
atmospheric concentration of carbon dioxide into radiative forcing. Using the
estimated model, and different future paths of greenhouse gas emissions, we
project global mean surface temperature until the year $2100$. Our results
illustrate the potential of combining statistical modelling with physical
insights to arrive at rigorous statistical analyses of the climate system.

##  1 Introduction

Climate models play an important role in our understanding of the past,
present, and future climate. Most notably, the reports on the “physical
science basis” for global warming, compiled by the Intergovernmental Panel on
Climate Change (IPCC), rely heavily on output from large-scale Earth system
models (Solomon et al., 2007; Stocker and Qin, 2013; Masson-Delmotte and Zhai,
2021). However, since these large-scale models require enormous amounts of
computer power to run even a single instance, they are not suited for creating
large ensembles of climate system variables, which could be useful for, e.g.,
probabilistic analysis and/or for studying the effects of a large number of
different scenarios for future greenhouse gas emissions. In contrast to large-
scale Earth system models, reduced complexity climate models (RCMs) are
simpler climate models designed to run on desktop computers. The lower
computational requirements of RCMs make them eminently suitable for running
large ensembles and investigating the impact of a wide range of different
greenhouse gas scenarios, as well as for probabilistic analysis, which would
be infeasible using large scale models (Nicholls et al., 2020; Nicholls et
al., 2021). RCMs are by now recognized as important tools for understanding
the climate system. For instance, the “MAGICC” RCM (Meinshausen et al., 2011)
is instrumental in probabilistic analyses of the climate system (Meinshausen
et al., 2009) as well as in the construction of the influential Representative
Concentration Pathway (RCP) scenarios (Meinshausen et al., 2011). Both MAGICC
and the “FaIR” RCM (Leach et al., 2021) are extensively used in the 6’th
assessment report from the IPCC (AR6; Masson-Delmotte and Zhai, 2021, see,
e.g., Chapter 7).

We consider a statistically-based RCM where the entire chain from emissions to
temperatures is formulated as a statistical model with stochastic
representations of variables and their dependencies implied by the underlying
physics. This approach is in contrast to a physically-based RCM which is
formulated as a deterministic model using equations reflecting the underlying
physics. As a consequence, physically-based RCMs are not accounting for
stochastic variation and uncertainty in the modeling framework, and hence the
statistical properties of the calibrated/estimated model parameters are
unknown. Furthermore, the informative and insightful “probabilistic”
projections are not necessarily based on justified statistical foundations.
Most RCMs that model the entire chain from emissions to temperatures are
physically-based (Smith et al., 2024).111Notably, the most recent version of
the FaIR model (Smith et al., 2024) treats the energy balance module
statistically using the framework suggested in Cummins et al. (2020) and the
remaining modules deterministically, putting this model somewhere in-between
being statistically-based and physically-based.

In this paper, we propose and develop a fully statistical reduced complexity
climate model (Stat-RCM) that models the entire chain from emissions to
temperatures. We assume that the climate system consists of various latent
processes — such as the atmospheric concentration of carbon dioxide (CO2) and
the mean temperature at the Earth’s surface — and that historical data are
noisy observations of these variables. This framework allows us to cast the
climate model in a state space system, where the state equations represent the
latent “true” climate processes and the measurement equations express
observations as noisy measurements of these latent processes. We then specify
the transition equations of the state variables using functional forms which
are chosen to be compatible with the physics underlying the climate system.
The stochastic error processes in the state equation represent the difference
between the latent “true” climate system and our functional representation of
it.

Our proposed Stat-RCM treats the physical processes of the climate system at a
lower resolution compared to physically-based RCMs, such as MAGICC
(Meinshausen et al., 2011), FaIR (Leach et al., 2021), and Hector (Hartin et
al., 2015). The simpler physical representation of the climate system employed
by the Stat-RCM is necessitated by statistical identification, i.e., the
condition that the objective function, usually the log-likelihood, as a
function of the parameters of the model has a unique extremum. Identification
allows for the invocation of standard statistical results, such as a law of
large numbers and a central limit theorem. The former ensures that the
estimator of the parameter vector is consistent and the latter ensures that
the estimator is asymptotically normally distributed, which allows for
assessment of parameter uncertainty. We use our framework in this study to (i)
estimate the parameter vector of the model via the maximum likelihood method
using historical data from the period $1959$–$2022$; (ii) obtain standard
errors for the estimated parameters; (iii) investigate the finite sample
properties of the estimation method via Monte Carlo simulations; (iv) conduct
model selection on the parts of the model that are not informed by physical
first-principles; (v) extract estimates of the underlying “true” climate state
variables; (vi) perform statistical validation of the model by showing that it
can simulate climate variables that are compatible with the historical data
record; and (vii) construct probabilistic projections of the climate system,
conditional on a scenario for future greenhouse gas emissions. Our results
illustrate the potential for combining statistical modelling with physical
insight to arrive at rigorous probabilistic and statistical analyses of the
climate system.

The remainder of the paper is structured as follows. In Section 2, we
formulate the basic model equations of the Stat-RCM, and in Section 3 we show
how to cast these equations into a non-linear state space system. Section 4
presents the results from a Monte Carlo simulation study that aims to
investigate the finite sample properties of the proposed estimation procedure.
Section 5 shows and discusses the estimation results for the Stat-RCM based on
historical data from the period $1959$–$2022$. In Section 6 we adopt these
estimation results and use these to generate projections of the future
climate, conditional on given paths of greenhouse gases and other forcing
agents. Section 7 concludes. An appendix contains supplementary empirical
results and details on the implementation of the employed state space methods.

##  2 Stat-RCM: A new statistical reduced complexity climate model

An overview of Stat-RCM is given in Figure 1. The figure shows how
anthropogenic CO2 emissions from fossil fuels ($E^{FF}$) and land–use change
($E^{LUC}$) are affecting temperatures: first, total CO2 emissions
($E=E^{FF}+E^{LUC}$) are converted into atmospheric concentrations ($C$) in a
carbon cycle module; then atmospheric concentrations are converted into
radiative forcing ($F$) using a forcing equation; lastly, an energy balance
module converts radiative forcing into changes in surface temperature
($T^{m}$) and deep ocean temperature ($T^{d}$). CO2 emissions and the forcing
effects of non-CO2 greenhouse gases ($F^{Non-CO2}$) and various natural
phenomena such as volcanic eruptions and solar activity ($F^{Nat}$) are
incorporated into the model as covariates, i.e. these processes are “external”
and not determined inside the model.

The following subsections give further details on the three main modules in
Stat-RCM, namely the carbon cycle module, the forcing equation, and the energy
balance module. Similar models for the carbon cycle and energy balance modules
were considered separately in Bennedsen et al. (2023) and Bennedsen et al.
(2023), respectively. The Stat-RCM extends these models by connecting them
through the forcing equation, thereby modelling the entire chain from
emissions to temperatures, and by introducing climate feedbacks into the
carbon cycle module.

Figure 1: Diagram of Stat-RCM Forcing equation $C\rightarrow F$ Carbon cycle
module $E\rightarrow C$ External $E=E^{FF}+E^{LUC}$ Energy balance module
$F\rightarrow T$ Output $T^{m}$, $T^{d}$ External
$F^{\textrm{Ex}}=F^{\textrm{Non-CO2}}+F^{\textrm{Natural}}$

###  2.1 Carbon cycle module

The cornerstone of the carbon cycle module is the carbon budget equation,

$\displaystyle\frac{dC_{t}}{dt}=E^{FF}_{t}+E_{t}^{LUC}-S_{t}^{LND}-S_{t}^{OCN},$
(2.1)

where $C_{t}$ is atmospheric concentration of carbon dioxide (CO2) at time
$t$, while $E^{FF}_{t}$ and $E_{t}^{LUC}$ are emissions of CO2 at time $t$
from fossil fuel burning and land-use change, respectively, and, $S_{t}^{OCN}$
and $S_{t}^{LND}$ are fluxes in the ocean and terrestrial (land) sink,
respectively. Atmospheric concentrations, $C_{t}$, are given in gigatonnes of
carbon (GtC), while all terms on the right hand side of Equation (2.1) are
given in GtC _per year_. The carbon budget equation (2.1) captures the
principle that Earth’s carbon cycle is ‘closed’ in the sense that all CO2
emitted from human activities ($E^{FF}_{t}+E^{LUC}_{t}$) must be absorbed by
one of the three carbon sinks, namely the atmosphere ($dC_{t}/dt$), the
terrestrial biosphere ($S_{t}^{LND}$), and the oceanic biosphere
($S_{t}^{OCN}$); see Friedlingstein et al. (2022) for an extensive discussion
of the carbon budget.

The emissions time series $E_{t}^{FF}$ and $E_{t}^{LUC}$ are considered as
covariates (exogenous data), which “force” the carbon system. The stock time
series $C_{t}$ as well as the flux time series data $S_{t}^{LND}$ and
$S_{t}^{OCN}$ are included as part of the dynamic model system and interact
with the other variables in the model. As is common in reduced-form modelling
of the carbon cycle, we will assume that the land and ocean fluxes are related
to the level of atmospheric concentrations ($C$) and the state of the climate,
represented by the global surface temperature above pre-industrial levels
($T^{m}$). In particular, we specify the dynamics of the sinks as

$\displaystyle S_{t}^{OCN}$
$\displaystyle=g_{OCN}(C_{t},T_{t}^{m})-g_{OCN}(C_{1750},T_{1750}^{m})+\eta_{t}^{O},$
$\displaystyle S_{t}^{LND}$
$\displaystyle=g_{LND}(C_{t},T_{t}^{m})-g_{LND}(C_{1750},T_{1750}^{m})+\eta_{t}^{L},$

for non-linear functions $g_{OCN}(\cdot,\cdot)$ and $g_{LND}(\cdot,\cdot)$,
and where the terms $\eta_{t}^{O}$ and $\eta_{t}^{L}$ are zero-mean stochastic
processes, $C_{1750}$ and $T_{1750}^{m}$ denote the pre-industrial values of
$C$ and $T^{m}$, respectively, with default values given by
$C_{1750}=591.3060$ GtC and $T_{1750}^{m}=0^{\circ}$C. The terms
$g_{i}(C_{1750},T_{1750}^{m})$, for $i=OCN,LND$, are introduced to capture the
assumption that the climate system was in equilibrium in the pre-industrial
era, meaning that the net flux of CO2 between the atmosphere and the sinks
were zero on average, i.e.
$\mathbb{E}[S_{1750}^{OCN}]=\mathbb{E}[S_{1750}^{LND}]=0$, where
$\mathbb{E}[\cdot]$ denotes the expectation operator.

Different RCMs consider different functional forms for $g_{OCN}$ and
$g_{LND}$. Over the time period studied in the empirical section below
($1959$–$2022$), the functions $g_{OCN}$ and $g_{LND}$ can be taken to be
approximately linear in concentrations $C$ (Bennedsen et al., 2023). We follow
the MAGICC RCM (Meinshausen et al., 2011) approach and incorporate $T^{m}$ via
exponential functions,

$\displaystyle g_{OCN}(C_{t},T_{t}^{m})$
$\displaystyle=b_{1}C_{t}\exp(-c_{1}T_{t}^{m}),$ $\displaystyle
g_{LND}(C_{t},T_{t}^{m})$ $\displaystyle=b_{2}C_{t}\exp(-c_{2}T_{t}^{m}),$

where $b_{1},b_{2},c_{1},c_{2}\in\mathbb{R}$ are constant parameters. For
$b_{1},b_{2}>0$, the activity of the sinks depends positively on the level of
atmospheric concentrations $C$, and we provide evidence of this below. For the
land sink, this is due to the fertilization effect (e.g. Bacastow and Keeling,
1979), while for the ocean sink it is due to the fact that the uptake depends
on the difference in partial pressure of CO2 between the atmosphere and the
upper ocean (e.g. Joos et al., 1996). The dependence of the sink uptake on the
level of CO2 concentrations in the atmosphere, is often referred to as the
‘carbon cycle feedback’. Conversely, the dependence of the sink uptakes on the
state of the climate, represented by the global surface temperature $T^{m}$,
is referred to as the ‘climate feedback’ (e.g. Friedlingstein et al., 2003;
Fung et al., 2005; Friedlingstein, 2015).

###  2.2 Forcing equation

The forcing equation connects the level of CO2 concentrations (given in GtC)
to radiative forcing from CO2 (given in W/m2) above the pre-industrial level.
Here, it is specified as

$\displaystyle F_{t}^{CO2}=g_{F}(C_{t})-g_{F}(C_{1750})+\eta_{t}^{F},$ (2.2)

where $g_{F}(\cdot)$ is a non-linear function describing the link between CO2
concentrations and forcing, $\eta_{t}^{F}$ is a zero-mean stochastic process,
and the term $g_{F}(C_{1750})$ is introduced to reflect the assumption that
the climate system was in equilibrium in the pre-industrial era, i.e.
$\mathbb{E}[F_{1750}^{CO2}]=0$. In this study, we consider the following
functional form for $g_{F}()$

$\displaystyle
g_{F}(C_{t})=f_{1}\log\left(C_{t}+f_{2}C_{t}^{2}\right)+f_{3}\sqrt{C_{t}},$
(2.3)

where $f_{1},f_{2},f_{3}\in\mathbb{R}$ are unknown constant parameters. This
specification for $g_{F}()$ is general and it nests several functional forms
which have been used in the RCM literature previously. Setting $f_{1}=f_{2}=0$
will result in the specification used in MAGICC6 (Meinshausen et al., 2011);
setting $f_{2}=0$ will result in the specification used in FAIR (Leach et al.,
2021); while $f_{3}=0$ yields the form proposed in Hansen et al. (1998).
Below, we will exploit the statistical nature of the model to test which of
these functional forms are most appropriate when confronted with historical
data.

###  2.3 Energy balance module

For the energy balance module, we consider a two-box energy balance model
(Gregory, 2000)

$\displaystyle H_{m}\frac{dT_{t}^{m}}{dt}$ $\displaystyle=F_{t}-\lambda
T_{t}^{m}-\gamma(T_{t}^{m}-T_{t}^{d}),$ $\displaystyle
H_{d}\frac{dT_{t}^{d}}{dt}$ $\displaystyle=\gamma(T_{t}^{m}-T_{t}^{d}),$

where $H_{m}$ and $H_{d}$ are fixed unknown coefficients, $T_{t}^{m}$ and
$T_{t}^{d}$ are time-$t$ temperature anomalies (with respect to the pre-
industrial period) in the atmosphere upper-ocean mixed layer and the deep
ocean layer, respectively, $F_{t}$ is time-$t$ radiative forcing, and
$\lambda$ and $\gamma$ are fixed unknown coefficients. The term
$F_{t}=F_{t}^{CO2}+F_{t}^{Ex}$ is the total radiative forcing in the system
which we model as the sum of the forcing from CO2 ($F_{t}^{CO2}$ is modelled
within the dynamic system) and forcing from other sources ($F_{t}^{Ex}$
represents a set of covariates). The parameters $H_{m}$ and $H_{d}$ (given in
$W$ year $m^{-2}K^{-1}$) are the corresponding heat capacities for the two
layers, $\lambda$ (given in $W$ $m^{-2}K^{-1}$) is a climate feedback
parameter, and $\gamma$ (given in $W$ $m^{-2}K^{-1}$) is the coefficient of
heat exchange between the two layers.

In this study, we work with the two-box energy balance model as described
above. However, it is straightforward to expand our setting to a general model
with one box or multiple boxes. For instance, a three-box energy balance model
can be considered using related methods to those considered in Cummins et al.
(2020).

##  3 Non-linear state space representation of Stat-RCM

In this section, we show how to cast the climate model presented in Section 2
in a non-linear state space system. We do this by using an Euler
discretization of the relevant continuous-time equations. A state space system
consists of two sets of equations: $(i)$ state equations for the latent system
dynamic variables, and $(ii)$ measurement equations for the observed data
variables. In Section 3.1, we present the state equations which are obtained
from the physical equations above. In Section 3.2, we present the measurement
equations which specify the data observations as noisy measurements of the
latent system variables. In Section 3.3, we collect the state and measurement
equations into a state space system, which can then be estimated and analyzed
using statistical methods designed for such systems. We refer to Durbin and
Koopman (2012) for a textbook treatment on the statistical analysis of state
space systems.

In what follows, let $\Delta>0$ denote the length of a time step. In our
applications of the model presented below, we will set $\Delta=1$ year, but,
for generality, we specify the model for a generic time step.

###  3.1 State equations

The carbon cycle module equations are discussed in Section 2.1 and, after some
substitutions and minor manipulations, we have

$\displaystyle
C_{t+\Delta}=C_{t}+(E_{t+\Delta}-S_{t+\Delta}^{LND}-S_{t+\Delta}^{OCN})\Delta,\qquad$
$\displaystyle
S_{t+\Delta}^{OCN}=-b_{1}C_{1750}+b_{1}C_{t}\exp(-c_{1}T_{t}^{d})+\eta_{t+\Delta}^{O},$
$\displaystyle
S_{t+\Delta}^{LND}=-b_{2}C_{1750}+b_{2}C_{t}\exp(-c_{2}T_{t}^{d})+\eta_{t+\Delta}^{L},$

where $b_{1},b_{2},c_{1},c_{2}\in\mathbb{R}$ are fixed unknown coefficients,
$C_{1750}=591.3060$ GtC is the pre-industrial level of CO2 concentrations,
$E_{t}=E_{t}^{FF}+E_{t}^{LUC}$, and disturbances $\eta_{t+\Delta}^{O}$ and
$\eta_{t+\Delta}^{L}$ are mean-zero innovation terms. The forcing equation is
discussed in Section 2.2 and we consider the equations (2.2)–(2.3) as given by

$F_{t+\Delta}^{CO2}=g_{F}(C_{t})-g_{F}(C_{1750})+\eta_{t+\Delta}^{F},\qquad\qquad
g_{F}(C_{t})=f_{1}\log\left(C_{t}+f_{2}C_{t}^{2}\right)+f_{3}\sqrt{C_{t}},$

with mean-zero innovation $\eta_{t+\Delta}^{F}$ and
$f_{1},f_{2},f_{3}\in\mathbb{R}$ are fixed unknown coefficients. The energy
balance module equations are discussed in Section 2.3 and are given by

$\displaystyle T_{t+\Delta}^{m}$
$\displaystyle=\left(1-\frac{\gamma+\lambda}{H_{m}}\Delta\right)T^{m}_{t}+\frac{\gamma\Delta}{H_{m}}T^{d}_{t}+\frac{\Delta}{H_{m}}(F_{t}^{CO2}+F_{t}^{Ex})+\eta_{t+\Delta}^{m},$
$\displaystyle T^{d}_{t+\Delta}$
$\displaystyle=\frac{\gamma\Delta}{H_{d}}T^{m}_{t}+\left(1-\frac{\gamma\Delta}{H_{d}}\right)T^{d}_{t}+\eta_{t+\Delta}^{d},$

where $\gamma,\lambda,H_{m},H_{d}\in\mathbb{R}$ are fixed unknown
coefficients, $F_{t}^{Ex}=F_{t}^{Non-CO2}+F_{t}^{Nat}$ is forcing from non-CO2
sources (anthropogenic and natural), which are treated here as covariates, and
disturbances $\eta_{t+\Delta}^{m}$ and $\eta_{t+\Delta}^{d}$ are mean-zero
innovation terms.

###  3.2 Measurement equations

We assume that we have a number of observed variables available that we can
regard as, possibly noisy, measurements for the state variables. We will
denote most of these observed variables by the same name as their state
variable counterparts but we affix an asterisk to it. Many different observed
variables can be considered and the core selection typically depend on data
availability. In our study, we consider the following set of seven observed
variables with associated equations,

$\displaystyle C_{t}^{*}$ $\displaystyle=\mu_{C}+C_{t}+\epsilon_{t}^{C},\quad$
(CO2 concentrations) $\displaystyle S_{t}^{OCN,*}$
$\displaystyle=\mu_{O}+S_{t}^{OCN}+\epsilon_{t}^{OCN},\quad$ (CO2 ocean sink)
$\displaystyle S_{t}^{LND,*}$
$\displaystyle=\mu_{L}+S_{t}^{LND}+\epsilon_{t}^{LND},\quad$ (CO2 land sink)
$\displaystyle F_{t}^{CO2,*}$
$\displaystyle=\mu_{F}+F_{t}^{CO2}+\epsilon_{t}^{F},\quad$ (Forcing from CO2)
$\displaystyle T_{t}^{m,*}$
$\displaystyle=\mu_{m}+T^{m}_{t}+\epsilon_{t}^{m},\quad$ (Surface temperature)
$\displaystyle T_{t}^{d,*}$
$\displaystyle=\mu_{d}+T^{d}_{t}+\epsilon_{t}^{d},\quad$ (Deep ocean
temperature $0$–$2000$m) $\displaystyle O_{t}^{*}$
$\displaystyle=H_{d}\cdot\mu_{d}+H_{d}\cdot T^{d}_{t}+\epsilon_{t}^{O},\quad$
(Ocean heat content $0$–$2000$m)

where $H_{d}$ is the fixed unknown coefficient from the $T^{d}_{t}$ state
equation in Section 3.1, all $\epsilon^{\cdot}_{t}$ variables are mean-zero
error processes and the unknown fixed constants $\mu_{\cdot}\in\mathbb{R}$ are
introduced to ensure that all state variables, i.e. the variables without
asterisks, are benchmarked to the same base period. In our application below,
we will choose this to be the pre-industrial period, represented by the year
$1750$. This means that $\mu_{m}$, for instance, denotes the offset needed in
the observations $T^{m,*}$ to benchmark these to the pre-industrial period,
i.e. $T^{m,*}-\mu_{m}$ is the surface temperature anomaly relative to this
period. The last measurement equation for ocean heat content $O_{t}^{*}$ is
obtained from the relationship

$H_{d}\frac{dT^{d}_{t}}{dt}=\frac{dO_{t}}{dt},$

(e.g. Schwartz, 2007; Bennedsen et al., 2023) and the fact that the
observation time series $T_{t}^{d,*}$ and $O^{*}_{t}$ are anomalies assumed to
be benchmarked to the same historical period.

###  3.3 State space representation

Let
$y_{t}=(C_{t}^{*},S_{t}^{OCN,*},S_{t}^{LND,*},F_{t}^{CO2,*},T_{t}^{m,*},T_{t}^{d,*},O_{t}^{*})^{\prime}$
denote the $7\times 1$ vector of observations at time $t$. The resulting non-
linear state space model for the climate system described in Section 2 is

$y_{t}\ =\ \mu+Ax_{t}+\epsilon_{t},\qquad\qquad x_{t+\Delta}\ =\
B(x_{t})+W_{t}+R\eta_{t,\Delta},$

where the $7\times 1$ vector $\mu$ contains the intercepts in the measurement
equations of Section 3.2, the $7\times 6$ matrix $A$ enables the appropriate
selection of the element in the $6\times 1$ state vector
$x_{t}=(C_{t},S_{t}^{OCN},S_{t}^{LND},F_{t}^{CO2},T_{t}^{m},T_{t}^{d})^{\prime}$
for each variable in $y_{t}$, the $7\times 1$ vector $\epsilon_{t}$ contains
the measurement errors listed in Section 3.2, the $6\times 1$ non-linear
vector function $B(\cdot)$ captures the non-linear dynamic equations for the
latent state variables in Section 3.1, the $6\times 1$ vector $W_{t}$ consists
of zeros and composite covariates such as
$\Delta(F_{t}^{CO2}+F_{t}^{Ex})/H_{m}$, and the $6\times 5$ matrix $R$ selects
the appropriate element from the $5\times 1$ innovation vector
$\eta_{t,\Delta}=(\eta_{t}^{\prime
O},\eta_{t}^{L},\eta_{t}^{F},\eta_{t}^{m},\eta_{t}^{d})^{\prime}_{\Delta}$.
The state update equation comprises the dynamic processes of the physical
quantities in Section 3.1. The measurement error sequence $\epsilon_{t}$
captures the transitory deviations between observations and the corresponding
underlying state processes. We allow each element in $\epsilon_{t}$ to be a
serially correlated sequence modelled as a first-order stationary
autoregressive process, that is

$\epsilon_{t+\Delta}=\Phi\epsilon_{t}+\xi_{t,\Delta},\qquad\xi_{t,\Delta}\stackrel{{\scriptstyle
iid}}{{\sim}}N(0,\Delta\cdot P),$

where the autoregressive coefficient matrix $\Phi$ is diagonal, $N(0,V)$
refers to the multivariate normal distribution with a zero mean vector and
variance matrix $V$, and $iid$ means _independent and identically
distributed_. To impose the stationarity of $\epsilon_{t}$, we assume that
each diagonal element of $\Phi$ is smaller than unity in absolute value. Given
that variance matrix $P=Var(\xi_{t,\Delta})/\Delta$, stationarity implies that
$Var(\epsilon_{t})=(I-\Phi^{2})^{-1}P$. A detailed discussion on the dynamic
specification of $\epsilon_{t}$ is presented in Section 5.1. We further have
the innovation sequence

$\eta_{t,\Delta}\stackrel{{\scriptstyle iid}}{{\sim}}N(0,\Delta\cdot Q),$

with diagonal variance matrix $Q$, in recognition that the physical equations
in Section 3.1 are only approximations to the true climate system. It follows
that $R\eta_{t,\Delta}=x_{t+\Delta}-B(x_{t})-W_{t}$ can be regarded as the
error between the “true” climate system and our simplified non-linear version
of it. More details on the entities $\\{A,\,B(\cdot),\,R,\,\Phi,\,P,\,Q\\}$,
and how they are specified in the Stat-RCM, are provided in Appendix A. The
fixed unknown constants and coefficients are collected in the parameter vector
$\theta$ and are estimated by the method of maximum likelihood, where the
extended Kalman filter is used to evaluate the likelihood function. We refer
to Durbin and Koopman (2012, Chapter 10) for further details on this approach
to estimation.

##  4 Monte Carlo simulation study

We perform a Monte Carlo simulation study designed to gauge the finite sample
performance of the maximum likelihood estimation method in a controlled
setting. We use the historical CO2 emissions and forcing data over the period
$1959$–$2022$ as the covariates in the Stat-RCM model; see Section 5 for a
more in-depth discussion on these data. We then simulate the stochastic
processes $\epsilon_{t}$ and $\eta_{t}$ from their specifications provided in
Section 3.3, and use these to generate simulated paths for the state vector
$x_{t}$ and the observation vector $y_{t}$ over the period $1959$–$2022$,
resulting in $n=64$ observations. The model parameters used for simulation are
set equal to the estimates we have obtained from our empirical study in
Section 5. In particular, we have set $\mu_{C}=\mu_{L}=\mu_{O}=\mu_{F}=0$ and
$f_{2}=f_{3}=0$. The remaining parameters are estimated freely by numerically
maximizing the likelihood function obtained by the extended Kalman filter. We
repeat this simulation-estimation process $1000$ times to obtain $1000$ Monte
Carlo estimates of the parameters of the model.

Table 1 summarizes the results of the simulation experiment, where we focus on
the estimates of the physical parameters in the model,
$\tilde{\theta}=(b_{1},\,b_{2},\,c_{1},\,c_{2},\,f_{1},\,\gamma,\,\lambda,\,H_{m},\,H_{d})^{\prime}$,
and the constant offsets $\mu_{m}$ and $\mu_{d}$. Figure 2 presents histograms
of the $1000$ estimates for each of the parameters in $\widetilde{\theta}$,
where a vertical red line indicates the ‘true’ data-generating value of the
respective parameter. We can conclude that the parameters of the carbon cycle
module $\\{b_{1},\,b_{2},\,c_{1},\,c_{2}\\}$ and forcing equation
$\\{f_{1}\\}$ are estimated with a high degree of precision. The parameters
from the energy balance module $\\{\gamma,\,\lambda,\,H_{m},\,H_{d}\\}$, as
well as the constant offsets $\\{\mu_{m},\,\mu_{d}\\}$, appear more difficult
to estimate. In particular, the estimates of $H_{m}$ and $H_{d}$ appear to
show some small-sample bias. Overall, the finite sample properties of the
estimation procedure are adequate in this Monte Carlo study design and this
choice of parameter values.

Table 1: Monte Carlo simulation results

| $b_{1}$ | $b_{2}$ | $c_{1}$ | $c_{2}$ | $f_{1}$ | $\gamma$ | $\lambda$ | $H_{m}$ | $H_{d}$ | $\mu_{m}$ | $\mu_{d}$
---|---|---|---|---|---|---|---|---|---|---|---
True values: | $0.01$ | $0.02$ | $0.09$ | $0.09$ | $5.58$ | $1.44$ | $1.44$ | $8.97$ | $265.88$ | $0.28$ | $0.05$
MC Avg | $0.01$ | $0.02$ | $0.08$ | $0.10$ | $5.58$ | $1.24$ | $2.21$ | $8.23$ | $265.81$ | $0.39$ | $0.11$
(MC Std.) | $(0.00)$ | $(0.00)$ | $(0.03)$ | $(0.05)$ | $(0.01)$ | $(0.64)$ | $(1.25)$ | $(2.81)$ | $(0.25)$ | $(0.17)$ | $(15.78)$

Simulation results for physical parameters
$\tilde{\theta}=(b_{1},\,b_{2},\,c_{1},\,c_{2},\,f_{1},\,\gamma,\,\lambda,\,H_{m},\,H_{d})^{\prime}$
as well as the constant offsets $\mu_{m}$ and $\mu_{d}$. Number of Monte Carlo
replications is $1000$ and number of observations is $n=64$. The first row
reports the true parameters used for simulating the paths of the state and
measurement variables and the second and third row report the Monte Carlo
average and Monte Carlo standard deviation, respectively, calculated over the
$1000$ replications.

Figure 2: Histogram of the parameter estimates of
$\tilde{\theta}=(b_{1},\,b_{2},\,c_{1},\,c_{2},\,f_{1},\,\gamma,\,\lambda,\,H_{m},\,H_{d})^{\prime}$
obtained from $1000$ Monte Carlo replications, see Section 4. The vertical red
line indicates the true parameter used as data generating value in the
simulations.

##  5 The performance of Stat-RCM on the $1959$-$2022$ historical record

We consider a historical yearly time series data set from $1959$ to $2022$.
The data for the carbon cycle
$\\{C_{t}^{*},\,S_{t}^{OCN,*},\,S_{t}^{LND,*},\,E_{t}^{FF},\,E_{t}^{LUC}\\}$
are from the Global Carbon Project (Friedlingstein et al.,
2022)222https://www.globalcarbonproject.org.. The forcing data
$\\{F_{t}^{CO2,*},\,F_{t}^{Non-CO2},\,F_{t}^{Nat}\\}$ are the current best
estimates taken from the Sixth Assessment Report (AR6) of the IPCC (Smith et
al., 2021)333https://zenodo.org/record/5705391#.ZFIx_y9ByX2.. The global mean
surface temperature data ($T_{t}^{m,*}$) are from HadCRUT5 by the Met Office
Hadley Centre (Morice et al.,
2020)444https://www.metoffice.gov.uk/hadobs/hadcrut5/data/HadCRUT.5.0.2.0/download.html..
The deep ocean temperature data ($T_{t}^{d,*}$) and the ocean heat content
data ($O_{t}^{*}$) are from the Institute of Atmospheric Physics (IAP; Cheng
et al., 2017)555http://www.ocean.iap.ac.cn/, all data resources 1-4 are last
accessed on March 13, 2024..

The AR6 forcing data of Smith et al. (2021) end in $2019$ while our remaining
data run until $2022$. For the forcing variable $F_{t}^{CO2}$, which is a
model variable, it does not pose an issue since the Kalman filter and smoother
can easily handle missing data values. For the covariates, including the
forcing variables $F_{t}^{Non-CO2}$ and $F_{t}^{Nat}$, however, we need a
complete data record. We therefore impute the $2020$–$2022$ values of
$F_{t}^{Non-CO2}$ and $F_{t}^{Nat}$ using a linear trend, estimated from the
last five years of the data, i.e. $2015$–$2019$. The results obtained from the
Stat-RCM, reported below, are robust to alternative approaches to filling in
these missing values. This is not surprising since the forcing from CO2 is
much larger in magnitude than the forcing from non-CO2 greenhouse gases and
natural forcing.

The carbon cycle data from the Global Carbon Project and the forcing data for
the AR6 report are benchmarked to $1750$, meaning that we may set
$\mu_{C}=\mu_{L}=\mu_{O}=\mu_{F}=0$. We benchmark the surface temperature
data, $T_{t}^{m,*}$, to the base period $1850$–$1900$ and the deep ocean data,
$T_{t}^{d,*}$ and $O_{t}^{*}$, to $1940$. To ensure that the state variables
$T^{m}_{t}$ and $T_{t}^{d}$ are benchmarked to $1750$, we include the offset
parameters $\mu_{m}$ and $\mu_{d}$ in the relevant measurement equations and
estimate these freely.

###  5.1 The dynamic specification of the error structure

In Section 3.3, we have modelled the stochastic error variables in the
$7\times 1$ vector $\epsilon_{t}$ as autoregressive processes of order 1. This
type of process falls into the class of the autoregressive moving average
(ARMA) processes which is well-known in time series analysis (Brockwell and
Davis, 1996). The ARMA model can have different maximum lag lengths for the
autoregressive part, denoted by $p$, and for the moving average part, denoted
by $q$. We first have considered the ARMA specification, with $p=q=1$, for
each error term in $\epsilon_{t}$. In addition, we have included a correlation
coefficient, $\rho$, between the ARMA innovation term $\xi_{t}$ (see also
Section 3.3) associated with the measurement equations for $T^{d}_{t}$ and
$O_{t}$. These two variables are closely related and therefore this dependence
is introduced, see also Bennedsen et al. (2023). The initial estimates of the
ARMA parameters have revealed that parameters of the moving average part are
statistically insignificant at the $10\%$ level, while those of the
autoregressive part are significant. We therefore have opted for the
autoregressive part only and obtained the specification for $\xi_{t}$ as in
Section 3.3, also see Appendix A for further details. The estimate of the
dependence coefficient $\rho$ is statistically significant and therefore it is
kept in the specification for $\xi_{t}$.

###  5.2 Model selection: the forcing equation

The Beer-Lambert law suggests that the functional relationship between
atmospheric CO2 concentration and radiative forcing from CO2 is approximately
logarithmic (Lightfoot and Mamer, 2014). Various functional forms have been
suggested to improve upon this approximation (e.g. Hansen et al., 1988; Shi,
1992; Hansen et al., 1998). As discussed in Section 2.2, there appears to be
no consensus in the RCM literature on which of these functional forms is most
adequate for describing the link between atmospheric concentrations and
radiative forcings. Here, we will adopt statistical model selection methods to
determine which functional form is most appropriate for the historical data
under study. To the best of our knowledge, statistical model selection methods
have not been used to inform the construction of RCMs previously. Given the
statistical formulation of the Stat-RCM, model selection is relatively
straightforward to carry out.

In relation to the forcing equation (2.3), that is
$g_{F}(C_{t})=f_{1}\log\left(C_{t}+f_{2}C_{t}^{2}\right)+f_{3}\sqrt{C_{t}}$,
we are interested in testing the hypotheses

$\displaystyle H_{0}^{sqrt}:\ f_{1}$ $\displaystyle=f_{2}=0,$ $\displaystyle
H_{0}^{log}:\ f_{2}$ $\displaystyle=f_{3}=0,$ $\displaystyle
H_{0}^{log+sqrt}:\ f_{2}$ $\displaystyle=0,$ $\displaystyle H_{0}^{log2}:\
f_{3}$ $\displaystyle=0,$ $\displaystyle H_{0}^{Hansen98}:\ f_{1}$
$\displaystyle=5.04,\ f_{2}=0.00023507,\ f_{3}=0.$

The two hypotheses $H_{0}^{sqrt}$ and $H_{0}^{log}$ can be used to test
whether the forcing equation reduces to a square root function and a
logarithmic function, respectively. The hypothesis $H_{0}^{log}$ reduces the
forcing equation to the one used in MAGICC6 (Meinshausen et al., 2011, p.
1439). Under the hypothesis $H_{0}^{log+sqrt}$, the forcing equation is the
sum of a square root term and a logarithmic term, which is the forcing
equation used in the FaIR RCM (Leach et al., 2021, p. 3013). The forcing
equation implied by $H_{0}^{log2}$ is proposed in Hansen et al. (1998). The
hypothesis $H_{0}^{Hansen98}$ implies the same forcing equation, but where the
parameters are fixed to the specific values used in Hansen et al.
(1998).666Table 1 of Hansen et al. (1998) actually specifies the forcings
function $g_{F}(c_{t})=5.04\log(c_{t}+0.0005c_{t}^{2})$, where $c_{t}$ is
atmospheric concentrations given in ppm. Since $C_{t}=c_{t}\cdot
2.127$GtC/ppm, this translates into the parameters $f_{1}=5.04$ and
$f_{2}=0.0005/2.127=0.00023507$ in our framework. When estimating the
parameters from this forcing equation, i.e. when working under $H_{0}^{log2}$,
we find $\hat{f}_{1}=5.58$ and $\hat{f}_{2}=0.00000105$, indicating that the
forcing data used in this paper results in somewhat different parameter
estimates than what was used in Hansen et al. (1998). See also Ramaswami
(2001, Table 6.2., p. 358) for various suggestions for expressions of the
forcing equation, including some of those studied here.

Since the models implied by the hypotheses listed above are nested in the
general forcing equation (2.3), where $f_{1}$, $f_{2}$ and $f_{3}$ are
unrestricted coefficients, we can perform the various hypothesis tests using a
standard likelihood ratio test. Table 2 contains the maximized log-likelihood
values for the models with different forcing specifications, together with the
Bayesian Information Criteria (BIC; Schwarz, 1978), the numbers of estimated
parameters, and the $p$-values for the likelihood ratio tests. We find that
the model specification under $H_{0}^{log}:\ f_{2}=f_{3}=0$ above has the
lowest BIC value and hence has most support from the data. The two hypotheses
$H_{0}^{sqrt}$ and $H_{0}^{Hansen98}$ can formally be rejected at a $1\%$
significance level and $H_{0}^{log2}$ barely at a $5\%$ level. The rejection
of $H_{0}^{sqrt}$ indicates that a log-term is preferred to a square root term
when modelling the relationship between atmospheric CO2 concentrations and
radiative forcing, for this updated AR6 forcing data set. Finally, the
rejection of $H_{0}^{Hansen98}$ indicates that the parameter values used in
Hansen et al. (1998) do not fit this data set.

The reported results below are based on the model specification under
$H_{0}^{log}$ as it produces the lowest BIC value, and is not rejected by the
likelihood ratio test. The model specifications with a comparable fit as the
$H_{0}^{log}$ model, that is the model implied by $H_{0}^{log+sqrt}$ and
$H_{0}^{log2}$ and the unrestricted model, yield similar results as those from
the model with $H_{0}^{log}$.

Table 2: Likelihood ratio tests for the forcing equation

| Unrestricted | $H_{0}^{sqrt+log}$ | $H_{0}^{log2}$ | $H_{0}^{sqrt}$ | $H_{0}^{log}$ | $H_{0}^{Hansen98}$
---|---|---|---|---|---|---
Maximized log-likelihood: | 925.79 | 924.78 | 923.77 | 917.17 | 923.76 | 911.21
Bayesian Information Criterion: | -1714.34 | -1716.48 | -1714.45 | -1705.41 | -1718.60 | -1697.67
Number of parameters: | 33 | 32 | 32 | 31 | 31 | 30
$p$-value (Likelihood Ratio test): | - | 0.1552 | 0.0444 | 0.00018 | 0.1313 | 0

The Likelihood Ratio tests are against the “unrestricted” model with
coefficients $f_{1}$, $f_{2}$ and $f_{3}$ estimated freely.

###  5.3 Estimation results and assessment of model fit

Given the findings from the model selection exercise above, we work under the
hypothesis $H_{0}^{log}$, i.e. we set $f_{2}=f_{3}=0$ and estimate the
parameter $f_{1}$ in the forcing equations (2.2)–(2.3) together with the other
parameters. The estimation results from this model applied to the 1959–2022
data are presented in Table 3. For purposes of brevity, we only present the
estimates from the “physical” parameters in the model, i.e. the parameters
governing the main workings of the carbon cycle module,
$\\{b_{1},\,b_{2},\,c_{1},\,c_{2}\\}$, the forcing module, $f_{1}$, the energy
balance module, $\\{\gamma,\,\lambda,\,H_{m},\,H_{d}\\}$, as well as the
constant offsets from the measurement equations $\\{\mu_{m},\,\mu_{d}\\}$.
Estimates of the remaining parameters are presented in Appendix B.

The signs and magnitudes of the estimated carbon cycle feedback parameters
$\\{b_{1},\,b_{2}\\}$ are as anticipated and are found in previous studies
(e.g. Bennedsen et al., 2023). The climate feedback parameters
$\\{c_{1},\,c_{2}\\}$ are both estimated as
$\widehat{c}_{1},\widehat{c}_{2}\approx 0.09$, implying that a temperature
rise of $1^{\circ}C$, above pre-industrial levels, would correspond to a
climate feedback factor of $\exp(-T_{t}^{m}\cdot\widehat{c}_{i})=\exp(-1\cdot
0.09)\approx 0.9139$, for $i=1,2$. This is roughly a $9\%$ weakening of the
sinks compared to pre-industrial levels (notice that $T_{1750}^{m}=0$). A
similar magnitude of the climate feedback estimate obtained from historical
data has been reported previously (Zhang et al., 2021) but it is lower than
values found using large-scale climate models (Friedlingstein, 2015) The
reason for the low estimate of the climate feedback effect in the Stat-RCM is
probably due to the input data used, i.e. the historical data $1959$–$2022$,
where climate feedback effects are still rather weak. To obtain more precise
estimates of the climate feedback effect would require data where these
feedbacks are more manifest. This could for instance be achieved by fitting
the Stat-RCM to output from large-scale climate models run over longer periods
than the historical period, similar to what is done for other RCM emulators
(Nicholls et al., 2020). This avenue is left for future research.

In Table 3, we present the results for the forcing module and we find that the
estimated parameter $\widehat{f}_{1}$ is highly significant and has the
expected sign. The parameter estimates from the energy balance module align
well with previous studies (Cummins et al., 2020; Pretis, 2020; Bennedsen et
al., 2023). We notice that an estimate of the equilibrium climate sensitivity
(ECS) can be derived from our estimate of the climate feedback parameter
$\lambda$. In particular, the IPCC Sixth Assessment Report (AR6; Forster et
al., 2021) suggests the relationship $ECS=F_{2\times CO2}/\lambda$, where
$F_{2\times CO2}$ is the radiative forcing in response to a doubling of CO2
concentrations in the atmosphere. Using the best estimate of
$\widehat{F}_{2\times CO2}\approx 3.93$ ($\pm 0.47,\ 5\%$–$95\%$ CI) Wm-2 from
AR6 (updated), we find that $\widehat{ECS}=F_{2\times
CO2}/\widehat{\lambda}=3.93/1.42=2.78^{\circ}C$. This estimate aligns well
with the “likely” range of $2.5^{\circ}$–$4.0^{\circ}C$ obtained from
instrumental records and reported in AR6 (Forster et al., 2021). Under certain
assumptions, we can also obtain the (asymptotic) variance of $\widehat{ECS}$
which results in a standard error of $1.03$.777The standard error of
$\widehat{ECS}$ can be obtained as follows. Define
$\widehat{ECS}=g(\widehat{F}_{2\times
CO2},\widehat{\lambda})=\widehat{F}_{2\times CO2}/\widehat{\lambda}$, where
$g(x,y)$ is the function $g(x,y)=x/y$. The Jabobian of $g$ is $\partial
g=(1/\widehat{\lambda},-\widehat{F}_{2\times
CO2}/\widehat{\lambda}^{2})^{\prime}$. Let $\Sigma$ be the $2\times 2$
diagonal matrix with $Var(\widehat{F}_{2\times CO2})$ and
$Var(\widehat{\lambda})$ on the diagonal. An estimate of
$Var(\widehat{\lambda})$ can be obtained from the maximum likelihood method
discussed in Section 3.3. Assuming Gaussian distributions and a $90\%$
confidence interval of $(3.93\pm 0.47)$ as reported in Forster et al. (2021),
we obtain the estimate $\widehat{Var}(\widehat{F}_{2\times
CO2})=(0.47/1.6449)^{2}$. The variance of $\widehat{ECS}$ can then be
approximated via the ‘delta method’ as $\widehat{Var}(\widehat{ECS})=(\partial
g)^{\prime}\Sigma\partial g$. Hence, the Stat-RCM estimate and associated
standard error of equilibrium climate sensitivity is $\widehat{ECS}=2.78$ with
standard error of $1.03$.

For a correctly specified model, the standardized one-step ahead prediction
residuals will be a sequence of iid $N(0,1)$ variables (Durbin and Koopman,
2012). Diagnostics on these standardized prediction residuals are presented in
Table 4. There appears to be some autocorrelation left in the prediction
residuals at the longer lags for the land sink. The predictions residuals from
the forcing equation display some indications of non-Gaussianity and
heteroskedasticity due to outliers in the early $1970$s. Overall, however, the
diagnostics are very reasonable and indicate that the model is able to fit the
data well. The raw data, along with smoothed states as output from the
extended Kalman filter, are presented in Figure 3. The standardized one-step
ahead prediction residuals used for the residual diagnostics in Table 4 are
displayed in Figure 4.

Table 3: Stat-RCM parameter estimates

| $b_{1}$ | $b_{2}$ | $\delta_{1}$ | $\delta_{2}$ | $f_{1}$ | $\gamma$ | $\lambda$ | $C_{m}$ | $C_{d}$ | $\mu_{m}$ | $\mu_{d}$
---|---|---|---|---|---|---|---|---|---|---|---
Estimate: | $0.01$ | $0.02$ | $0.08$ | $0.09$ | $5.58$ | $1.46$ | $1.42$ | $9.37$ | $265.90$ | $0.30$ | $0.20$
Std. Err.: | $(0.00)$ | $(0.00)$ | $(0.02)$ | $(0.04)$ | $(0.01)$ | $(0.58)$ | $(0.51)$ | $(2.44)$ | $(0.41)$ | $(0.41)$ | $(0.36)$
$t$-stat: | $11.25$ | $7.83$ | $5.16$ | $2.53$ | $580.42$ | $2.51$ | $2.75$ | $3.84$ | $656.40$ | $0.74$ | $0.55$

Parameter estimates by the maximum likelihood method applied to the
$1959$–$2022$ data.

Table 4: Residual diagnostics

$y$ | Mean | Std. | Skew | Kurt | SC | JB | DW |  LB$(1)$ |  LB$(5)$ |  LB$(10)$ | ARCH
---|---|---|---|---|---|---|---|---|---|---|---
C | $-0.01$ | $1.08$ | $-0.28$ | $3.02$ | $0.03$ | $0.81$ | $1.93$ | $0.08$ | $1.72$ | $12.17$ | $0.01$
OCN | $-0.06$ | $1.00$ | $-0.18$ | $2.35$ | $0.14$ | $1.45$ | $1.71$ | $1.07$ | $6.64$ | $8.51$ | $0.18$
LND | $-0.07$ | $0.92$ | $0.17$ | $2.65$ | $-0.09$ | $0.62$ | $2.18$ | $0.64$ | $17.10$ | $47.91$ | $1.63$
Forc | $0.02$ | $1.01$ | $-0.55$ | $3.94$ | $0.05$ | $5.20$ | $1.89$ | $0.13$ | $2.06$ | $6.51$ | $5.29$
Temp | $-0.01$ | $1.00$ | $0.06$ | $2.40$ | $-0.01$ | $0.98$ | $2.01$ | $0.01$ | $7.19$ | $15.48$ | $0.01$
O-Temp | $0.04$ | $0.99$ | $-0.25$ | $2.62$ | $-0.18$ | $1.05$ | $2.37$ | $2.27$ | $10.30$ | $17.38$ | $2.88$
OHC | $0.04$ | $1.00$ | $-0.25$ | $2.63$ | $-0.18$ | $1.03$ | $2.37$ | $2.29$ | $10.39$ | $17.36$ | $2.86$

Diagnostics of standardized one-step ahead prediction residuals as output from
the extended Kalman filter. SC is the estimate of the serial correlation
coefficient $\phi$ in the regression $y_{t}=\phi y_{t-1}+\epsilon$; $JB$ is
the Jarque-Bera test statistic (Jarque and Bera, 1987): the null hypothesis of
Gaussianity is rejected if $JB$ is larger than the $95\%$ critical value of
$5.99$. $DW$ is the Durbin-Watson test statistic (Durbin and Watson, 1971):
$DW<2$ indicates positive serial correlation, $DW>2$ negative serial
correlation; $DW=2$ indicates no serial correlation. $LB(p)$ is the Ljung-Box
$Q$ test statistic (Ljung and Box, 1978) for autocorrelation, calculated using
$p$ lags. The $95\%$ critical value of the LB test are $3.84,11.07,18.31$ for
$p=1,5,10$, respectively. ARCH is the test statistic of the Engle (1988) test
for heteroskadasticity. The $95\%$ critical value of the ARCH test is
$3.8415$.

Figure 3: Smoothed states from the Stat-RCM (red dashed lines), estimated on
the historical data record $1959$–$2022$ (black lines), as output from the
extended Kalman smoother. For surface temperature, ocean temperature, and
ocean heat content, the estimated constant offsets
($\mu_{m},\mu_{d},H_{d}\cdot\mu_{d}$) have been added to the smoothed states
to make them comparable to the data. Figure 4: Standardized one-step ahead
prediction residuals from the Stat-RCM estimated on the historical data record
$1959$–$2022$.

###  5.4 State estimation

The Kalman filter methodology allows for the efficient estimation of the
latent variables in the $7\times 1$ state vector $x_{t}$. In particular, the
state estimate obtained from the Kalman smoother is an estimate of
$E[x_{t}|y_{1959:2022}]$ for $t=1959,1960,\ldots,2022$, where $x_{t}$ denotes
the state at time $t$ and $y_{1959:2022}$ denotes the full historical data
set. These smoothed states are presented in Figure 3.

The estimated state for the surface temperature $T_{t}^{m}$, i.e.
$\widehat{T}_{t}^{m}=\widehat{E}[T_{t}^{m}|y_{1959:2022}]$ for
$t=1959,1960,\ldots,2022$, can be seen as an estimate of the underlying long-
term temperature anomaly, i.e. an estimate where transitory effects, such as
measurement errors, captured in our model by the error process $\epsilon_{t}$,
have been filtered out. Thus, our estimate of the surface temperature state
can serve as a broad indicator of the overall state of the climate. In this
way, it can be used to assess how close we are to breaching international
temperature agreements such as the 2015 Paris Agreement of keeping “global
temperatures well below $2^{\circ}$C above pre-industrial times while pursuing
means to limit the increase to $1.5^{\circ}$C” (UNFCCC, 2015). Assessment and
detection of whether temperature targets have been breached is difficult, due
to the inherent noisiness and variability of temperature measurements (Betts
et al., 2023). In AR6, the IPCC itself proposed a way of detection such as
breach: _A breach is deemed detected if the average temperature over a
$20$-year horizon exceeds the target_ (e.g. IPCC, 2023). This means that there
will be a long lag, on the order of a decade, before a breach of a temperature
target will be detected. Such slow detection times run the risk of delaying
action that may be relevant in case of a breach. For this reason, several
methods have recently been suggested, which are able to monitor the underlying
global temperature level closer and offer faster detection times of possible
breaches. For instance, Betts et al. (2023) suggests to blend the past $10$
years of historical temperature data with projections of the next $10$ years
of temperature data, obtained from a climate model. While this method is able
to estimate a “$20$-year average” of the current temperature level without
delay, it has the downside that the credibility of the estimate depends on the
projected temperature, which in turn depends on the climate model used to
construct the projections. Alternative methods for keeping track of current
underlying global temperatures include Copernicus’ approach of calculating a
$30$-year linear trend into the
future,888https://cds.climate.copernicus.eu/apps/c3s/app-c3s-global-
temperature-trend-monitor?month:float=10&year:float=2023 NASA’s approach of
using lowess smoothing,999https://climate.nasa.gov/vital-signs/global-
temperature/ and the “Real-time Global Warming Index” of Haustein et al.
(2017), which estimates the amount of warming due to anthropogenic causes by
regressing temperature data on forcing time series.

In contrast, the estimated temperature state from the Stat-RCM can be seen as
a purely statistical and data-driven alternative to assess the current level
of underlying warming above pre-industrial levels. This estimate is only based
on historical data and it is compatible with the physical components of the
Stat-RCM, i.e. the carbon cycle module, the forcing module, and the energy
balance module. Based on the data used in this paper, the current best-
estimate of the underlying warming above the $1850$–$1900$ baseline (the
baseline usually used for defining temperature targets), obtained from the
Stat-RCM, is
$\widehat{T}_{2022}^{m}+\widehat{\mu}_{m}=\widehat{E}[T_{2022}^{m}|y_{1959:2022}]+\widehat{\mu}_{m}=1.27^{\circ}$C.
This number is substantially above the $2022$ temperature measurement,
$T_{2022}^{m,*}=1.16$, see Figure 3. This is mostly due to the $2022$ La Niña
conditions which have been filtered out of the estimate
$\widehat{T}_{2022}^{m}=\widehat{E}[T_{2022}^{m}|y_{1959:2022}]$, leading to a
more precise representation of the underlying warming trend, i.e. a trend free
from transitory deviations such as ENSO.

##  6 Scenario projections using Stat-RCM

We compute scenario projections of the future climate system from the
estimated Stat-RCM, conditional on paths for the covariate (external)
variables CO2 emissions, forcing from non-CO2 greenhouse gases, and forcing
from natural sources. We focus particularly on the probabilistic
projections.Section 6.1 lays out the projection methodology, which is based on
simulations from the estimated Stat-RCM. Section 6.2 contains a validation
exercise, where we use historical covariates from the period $1959$–$2022$ and
investigate whether the Stat-RCM is able to produce “projections” over
$1959$–$2022$ that agree with the historical data record. Section 6.3 presents
scenario projections for $2023$–$2100$ in a setting where CO2 emissions are
rapidly declining, a necessary condition for halting global warming and
achieving the goals as set out by international agreements, such as the Paris
Agreement (e.g. Sanderson et al., 2016; Luderer et al., 2018; Tokarska and
Gillett, 2018).

###  6.1 Epistemic and aleatoric uncertainty in climate projections

To gauge the degree of uncertainty introduced by the separate stochastic parts
of the Stat-RCM, we consider the following four different approaches of using
the model to project the climate system forward, conditional on paths of all
covariate variables.

  1. 1.

Deterministic run: Deterministic projection where the parameters of the Stat-
RCM, $\theta$, are set to their maximum likelihood estimates as reported in
Section 5, and based on the $1959$–$2022$ data; we have $\theta=\hat{\theta}$
and all error processes are set to zero, $\epsilon_{t}=\eta_{t}=0$ for all
$t$.

  2. 2.

Uncertainty coming from unknown parameters ($\theta$): Projection where the
physical parameters,
$\tilde{\theta}=(b_{1},b_{2},c_{1},c_{2},f_{1},\gamma,\lambda,H_{m},H_{d})^{\prime}$,
are randomly sampled from the asymptotic distribution, that is
$\tilde{\theta}\sim N(\hat{\theta},\widehat{\Sigma})$, where $\hat{\theta}$
and $\widehat{\Sigma}$ are the maximum likelihood estimates of the physical
parameters and their variance-covariance matrix, based on the $1959$–$2022$
data. The remaining parameters are set equal to their maximum likelihood
estimates. The error processes are set to zero, i.e. $\epsilon_{t}=\eta_{t}=0$
for all $t$.

  3. 3.

Uncertainty coming from unknown parameters and state innovations
($\theta,\eta$): Similar to the setup in 2, but now adding the stochastic
variation due to the innovation terms in the state equations $\eta$, and
setting the measurement error processes to zero i.e. $\epsilon_{t}=0$ for all
$t$.

  4. 4.

Uncertainty coming from unknown parameters, state innovations, and measurement
errors ($\theta,\eta,\epsilon$): Similar to the setup in 3, but now adding the
stochastic variation due to the measurement errors $\epsilon_{t}$.

In the first setup, the projection is deterministic, conditional on estimated
parameters and input paths of the covariates. The next three setups contain
stochastic variation to increasing degrees: First, by considering only
randomness arising from the uncertainty in parameters ($\theta$), then also
including uncertainty coming from the evolution of the state equations
($\theta$ and $\eta$), and finally also including the uncertainty coming from
other processes, such as measurement errors and ENSO effects ($\theta$, $\eta$
and $\epsilon$). We refer to the uncertainty about parameters as “epistemic”
uncertainty, in the sense that this uncertainty reflects our ignorance of the
data generating process behind the climate system, while uncertainty stemming
from the error processes $\eta$ and $\epsilon$ may be referred to as
“aleatoric”, in the sense that this uncertainty captures the internal random
variability of the climate system itself. To construct projections and
confidence bands from the model under the various uncertainty specifications,
we simulate $10^{5}$ different instances of the variables in the model using
given paths for the covariates and plot the pointwise quantiles of $2.5\%$,
$50\%$ and $97.5\%$.

###  6.2 Validation exercise

We adopt the model specification and its parameter estimates as reported in
Section 5 for the period $1959$–$2022$, to perform a validation exercise. We
take the covariates CO2 emissions ($E_{t}$), forcing from non-CO2 sources
($F_{t}^{Non-CO2}$), and forcing from natural sources ($F_{t}^{Nat}$) as
inputs and keep them fixed to their historical values. The projection results
are presented in Figure 5. The top left panel shows the covariates (input
variables) and the remaining panels show the quantiles of the simulated model
variables, along with the original observations (red) used for the estimation
of the model superimposed. Although the original observations are used for the
parameter estimates, we stress that they are not otherwise used as input for
this validation exercise. From Figure 5, we find that the model is well-
validated over the estimation sample $1959$–$2022$, in the sense that it can
produce simulations of the climate system with a high degree of similarity to
the observed historical data.

Figure 5: Validation exercise using covariate data and estimated parameters
from the historical data $1959$–$2022$. Top left panel: Covariate input for
the Stat-RCM. Remaining panels: Output of the Stat-RCM, including data (red
lines) for comparison. Black line: Path of deterministic run (Setup in 1).
Dark grey shading: $95\%$ pointwise quantile bands for simulations with
parameter uncertainty ($\theta$; Setup in 2). Dark blue shading: $95\%$
pointwise quantile bands for simulations with uncertainty coming from
parameters and innovations in the state equation ($\theta,\eta$; Setup in 3).
Light blue shading: $95\%$ pointwise quantile bands for simulations with
uncertainty coming from parameters, innovations in the state equation, and
transitory error terms ($\theta,\eta,\epsilon$; Setup in 4).

###  6.3 Projections to $2100$ in a strong mitigation scenario

Next, we use the Stat-RCM to project the variables until the year $2100$,
conditional on given future paths for the covariates of CO2 emissions, forcing
from non-CO2 greenhouse gases, and natural forcing. The path for CO2 emissions
is taken from the SSP119 scenario (Meinshausen et al., 2020), and forcings
from non-CO2 greenhouse gases consistent with SSP119 are obtained from the
MAGICC RCM.101010The SSP119 scenario can be run in MAGICC in a web browser via
the link
https://live.magicc.org/scenarios/bced417f-0f7f-4bb7-8359-792a0a8b0368/overview.
Here, the forcing from non-CO2 greenhouse gases can also be downloaded. Last
accessed May 12, 2023. The path for future natural forcing is set constant
equal to the last in-sample value. The values for the future paths of the
covariates are shown in the top-left panel of Figure 6, where a solid line
denotes historical data, from which the parameter estimates for Stat-RCM are
obtained, and dashed lines denote the future scenario used in the projection
of the model. We notice that the SSP119 scenario is very ambitious in the
sense that it implies rapid reductions in greenhouse gas emissions in the
short term, and substantial carbon dioxide removal from the atmosphere
(negative emissions) in the long term.

We assess the uncertainty using $10^{5}$ simulated trajectories of the model,
see also Section 6.1. Figure 6 shows the projection results for all the
climate state variables output by Stat-RCM. To visualize how the the different
layers of uncertainty impact the projections, Figure 7 contains the first $50$
simulated trajectories in the case of surface temperature, $T^{m}$. The left
panel of Figure 7 includes only uncertainty coming from parameters ($\theta$),
while the middle panel additionally includes uncertainty coming from the
innovations in the state equation ($\theta,\eta$), and the right panel adds
uncertainty from the measurement errors ($\theta,\eta,\epsilon$). From the
left panel, we learn that uncertainty coming from parameters ($\theta$) is
substantial. The resulting simulated trajectories are very smooth, due to the
absence of stochastic innovations in these simulations. From the middle panel,
we learn that adding uncertainty from the innovations in the state equation
($\theta,\eta$) increases uncertainty only slightly, which is to be expected
since $\widehat{\sigma}_{\eta,m}^{2}=Var(\eta_{t}^{m})\approx 0$ (see Table 5
in Appendix B). Conversely, from the right panel, we see that adding
uncertainty from the measurement errors ($\theta,\eta,\epsilon$) increases
uncertainty noticeably. These findings indicate that epistemic uncertainty is
more important than aleatoric uncertainty for projections of the underlying
surface temperature state variable $T_{t}^{m}$, while aleatoric uncertainty is
crucial in accounting for the variability in observations of surface
temperatures, $T_{t}^{m,*}$.

We may use the model to assess the probability of particular events happening
to the climate system in the future, such as surface temperatures exceeding
predefined thresholds, given a scenario for the covariate processes. For
instance, we might address the question “What is the probability of exceeding
$1.5^{\circ}$C at some point in $2023$–$2100$, given some future trajectory of
CO2 emissions and other covariate processes?” We can answer such questions
using simulations similar to those above, by reporting the fraction of times
the event in question happens in the simulations. Using this methodology for
the SSP119 setting discussed in this section, the $1.5^{\circ}$C threshold is
exceeded at some point in the period $2023$–$2100$ with a $90\%$ probability
when considering uncertainty from all parts of the model
($\theta,\eta,\epsilon$). When considering only uncertainty from uncertainty
from parameters ($\theta$) or uncertainty coming from parameters and
innovations in the state equation ($\theta,\eta$), however, the probability
drops to around $7\%$ and $8\%$, respectively. The $50$ simulated trajectories
shown in Figure 7 illustrate how these numbers come about: When considering
only uncertainty from parameters ($\theta$) or from parameters and innovations
in the state equation ($\theta,\eta$), the resulting trajectories are
reasonably smooth (Figure 7, left and middle panels), while if uncertainty
from measurement errors are included ($\theta,\eta,\epsilon$), the resulting
trajectories are very volatile (Figure 7, right panel). This volatility in the
temperature measurements, arising e.g. from ENSO effects and other natural
transitory processes, means that most trajectories will cross the
$1.5^{\circ}$C threshold (denoted by the horizontal red line in Figure 7) at
some point in the period $2023$–$2100$. At the same time, the underlying long-
term trend, absent the measurement error process $\epsilon$, may plausibly
stay below the threshold.

In summary, according to Stat-RCM, there is very high probability that the
$1.5^{\circ}$C threshold will be breached by temperature observations
($T_{t}^{m,*}$) at some point over the period $2023$–$2100$, but, if we filter
out the transitory deviations from the underlying temperature trend,
represented by $\epsilon$, then there is, in fact, a very good chance that the
underlying temperature trend ($T_{t}^{m}$) will stay below $1.5^{\circ}$C if
the (very ambitious) SSP119 scenario is followed.

Figure 6: Projections as based on covariates emissions and forcing data from
SSP119, $2023$–$2100$, and Stat-RCM parameters estimated using historical data
$1959$–$2022$. Top left panel: Covariates input to the Stat-RCM. Remaining
panels: Output of the Stat-RCM, including data for comparison. Black line:
Path of deterministic run (Setup in 1). Dark grey shading: $95\%$ pointwise
quantile bands for simulations with parameter uncertainty ($\theta$; Setup in
2). Dark blue shading: $95\%$ pointwise quantile bands for simulations with
uncertainty coming from parameters and innovations in the state equation
($\theta,\eta$; Setup in 3). Light blue shading: $95\%$ pointwise quantile
bands for simulations with uncertainty coming from parameters, innovations in
the state equation, and transitory error terms ($\theta,\eta,\epsilon$; Setup
in 4).

Figure 7: $50$ simulated trajectories of the surface temperature anomaly
$T^{m}$, benchmarked to the period $1850$–$1900$. Black line: Historical data
(1959–2022). Left panel: Simulations with parameter uncertainty ($\theta$;
Setup in 2). Middle panel: Simulations with uncertainty coming from parameters
and innovations in the state equation ($\theta,\eta$; Setup in 3). Right
panel: Simulations with uncertainty coming from parameters, innovations in the
state equation, and transitory error terms ($\theta,\eta,\epsilon$; Setup in
4). Red dashed line denotes $1.5^{\circ}$C.

##  7 Conclusion

In this study, we have proposed a new statistical reduced-complexity climate
model (Stat-RCM). We cast the Stat-RCM in a non-linear state space system that
facilitates estimation, filtering, and smoothing using standard statistical
methods (Durbin and Koopman, 2012). The model treats anthropogenic CO2
emissions, forcing from other greenhouse gases, and natural forcing as
covariates. The climate variables atmospheric concentrations, ocean sink, land
sink, forcing from CO2, surface temperature, and deep ocean temperature are
modelled as part of our climate model.

The stochastic formulation of Stat-RCM implies that tools from statistical
theory are available. A Monte Carlo study has shown that the proposed
estimation procedure, relying on maximizing the log-likelihood function as
computed by the extended Kalman filter, enjoys good finite sample properties.
The simulated data for this study have been given similar properties as those
in the historical data record $1959$–$2022$. When estimation is actually
carried out on the historical data, a statistical model selection procedure
has indicated that a forcing equation consisting of a single logarithmic term
is adequate for describing the relationship between atmospheric concentrations
of CO2 and the corresponding radiative forcing. The resulting model has been
validated using two different methods: first, by showing that the estimated
residuals are conform the theoretical expectations; second, by showing that
the estimated climate variables accurately reproduce the historical data over
the period 1959–2022. We regard $(i)$ Monte Carlo simulation studies for
assessing the quality of the estimators, $(ii)$ statistical model selection
procedures for informing model specification, and $(iii)$ statistical
validations for checking possible model misspecification, as important
features of empirical statistical modelling. However, to the best of our
knowledge, such statistical methods have not previously been applied to RCMs
that model the entire chain from emissions to temperatures, due to the
physics-based (deterministic) formulation of these models. Conversely, due to
the statistical nature of the Stat-RCM, these analyses can be performed using
well-established methodology.

By extracting the latent temperature state using the extended Kalman filter,
we have estimated the $2022$ long-term global temperature anomaly to be
$+1.27^{\circ}C$ with respect to the $1850$–$1900$ baseline. We also have used
the Stat-RCM to project the climate variables until $2100$, conditional on a
scenario for future CO2 emissions and radiative forcing from non-CO2 sources.
We have found that in the SSP119 scenario, where emissions are rapidly
declining, the $1.5^{\circ}$C Paris upper limit on the global surface
temperature anomaly, compared to a $1850$–$1900$ baseline, is not a forgone
conclusion. The Stat-RCM estimates a $90\%$ probability that the
$1.5^{\circ}$C threshold will be breached by a temperature measurement
$T_{t}^{m,*}$ in the period $2023$–$2100$ under the SSP119 scenario. However,
if we consider only the underlying temperature trend, represented by the
latent state $T_{t}^{m}$, the Stat-RCM estimates that there is only a $7\%$
probability that the threshold will be breached in the period $2023$–$2100$
under the SSP119 scenario. In the latter case, the uncertainty is epistemic
(i.e. stemming from unknown parameters), while, in the former case, aleatoric
uncertainty is also included (i.e. uncertainty from internal variability in
the climate system, including measurement errors and ENSO effects).

We have focused on estimating the Stat-RCM parameters from the historical data
$1959$–$2022$. The main application of RCMs, however, has been as “emulators”,
where the RCM is used to emulate the output of large-scale climate models
(Nicholls et al., 2020). In future work, we intend to employ the Stat-RCM as
an emulator by estimating it using output from large-scale climate models,
such as that those from the Coupled Model Intercomparison Project (CMIP;
Eyring et al., 2016). This will allow us to employ the statistical methodology
illustrated in this paper to output from CMIP models with the aim of
complementing the results from existing RCMs, increasingly used in the IPCC
reports, with those from the Stat-RCM.

## References

  * Bacastow and Keeling (1979) Bacastow, R. B. and C. D. Keeling (1979).  Models to predict future atmospheric CO2 concentrations.  In Workshop on the global effects of carbon dioxide from fossil fuels, pp. 72–90. US Department of Energy.
  * Bennedsen et al. (2023) Bennedsen, M., E. Hillebrand, and S. J. Koopman (2023).  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. (2023) Bennedsen, M., E. Hillebrand, and J. Z. Lykke (2023).  Global temperature projections from a statistical energy balance model using multiple sources of historical data.  Journal of Climate 36(19), 6817–6838.  Working paper available at https://arxiv.org/abs/2205.10269.
  * Betts et al. (2023) Betts, R. A., S. E. Belcher, L. Hermanson, A. K. Tank, J. A. Lowe, C. D. Jones, C. P. Morice, N. A. Rayner, A. A. Scaife, and P. A. Stott (2023).  Approaching $1.5^{\circ}$C: how will we know we’ve reached this crucial warming mark?  Nature 624, 33–35.
  * Brockwell and Davis (1996) Brockwell, P. J. and R. A. Davis (1996).  Introduction to time series and forecasting.  Springer New York.
  * Cheng et al. (2017) Cheng, L., K. E. Trenberth, J. Fasullo, T. Boyer, J. Abraham, and J. Zhu (2017).  Improved estimates of ocean heat content from 1960 to 2015.  Science Advances 3(3), e1601545.
  * Cummins et al. (2020) Cummins, D. P., D. B. Stephenson, and P. A. Stott (2020).  Optimal estimation of stochastic energy balance model parameters.  Journal of Climate 33(18), 7909–7926.
  * Durbin and Koopman (2012) Durbin, J. and S. J. Koopman (2012).  Time series analysis by state space methods.  Number 38. Oxford University Press.
  * Durbin and Watson (1971) Durbin, J. and G. S. Watson (1971).  Testing for serial correlation in least squares regression.  Biometrika 58(1), 1 – 19.
  * Engle (1988) Engle, R. (1988).  Autoregressive conditional heteroscedasticity with estimates of the variance of United Kingdom inflation.  Econometrica 96, 893–920.
  * Eyring et al. (2016) Eyring, V., S. Bony, G. A. Meehl, C. A. Senior, B. Stevens, R. J. Stouffer, and K. E. Taylor (2016).  Overview of the Coupled Model Intercomparison Project Phase 6 (CMIP6) experimental design and organization.  Geoscientific Model Development 9(5), 1937–1958.
  * Forster et al. (2021) Forster, P., T. Storelvmo, K. Armour, W. Collins, J.-L. Dufresne, D. Frame, D. Lunt, T. Mauritsen, M. Palmer, M. Watanabe, M. Wild, and H. Zhang (2021).  The Earth’s energy budget, climate feedbacks, and climate sensitivity.  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 Z. B. (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. Cambridge: Cambridge University Press. In Press.
  * Friedlingstein (2015) Friedlingstein, P. (2015).  Carbon cycle feedbacks and future climate change.  Philosophical Transactions of the Royal Society A 373, 20140421\.
  * Friedlingstein et al. (2003) Friedlingstein, P., J. L. Dufresne, P. M. Cox, and P. Rayner (2003).  How positive is the feedback between climate change and the carbon cycle?  Tellus B: Chemical and Physical Meteorology 55(2), 692–700.
  * Friedlingstein et al. (2022) Friedlingstein, P., M. O’Sullivan, and et al. (2022).  Global Carbon Budget 2022.  Earth System Science Data 14(11), 4811–4900.
  * Fung et al. (2005) Fung, I. Y., S. C. Doney, K. Lindsay, and J. John (2005).  Evolution of carbon sinks in a changing climate.  Proceedings of the National Academy of Sciences 102(32), 11201–11206.
  * Gregory (2000) Gregory, J. M. (2000).  Vertical heat transports in the ocean and their effect on time-dependent climate change.  Climate Dynamics 16(7), 501–515.
  * Hansen et al. (1988) Hansen, J., I. Fung, A. Lacis, D. Rind, S. Lebedeff, R. Ruedy, G. Russell, and P. Stone (1988).  Global climate changes as forecast by Goddard Institute for Space Studies three-dimensional model.  Journal of Geophysical Research: Atmospheres 93(D8), 9341–9364.
  * Hansen et al. (1998) Hansen, J. E., M. Sato, A. Lacis, R. Ruedy, I. Tegen, and E. Matthews (1998).  Climate forcings in the industrial era.  Proceedings of the National Academy of Sciences 95(22), 12753.
  * Hartin et al. (2015) Hartin, C. A., P. Patel, A. Schwarber, R. P. Link, and B. P. Bond-Lamberty (2015).  A simple object-oriented and open-source model for scientific and policy analyses of the global climate system – Hector v1.0.  Geoscientific Model Development 8(4), 939–955.
  * Haustein et al. (2017) Haustein, K., M. R. Allen, P. M. Forster, F. E. L. Otto, D. M. Mitchell, H. D. Matthews, and D. J. Frame (2017).  A real-time Global Warming Index.  Scientific Reports 7(1), 15417.
  * IPCC (2023) IPCC (2023).  Climate Change 2023 Synthesis Report - Summary for policy makers.  Technical report, Intergovernmental Panel on Climate Change (IPCC), https://www.ipcc.ch/report/ar6/syr/downloads/report/IPCC_AR6_SYR_SPM.pdf.
  * 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.
  * Joos et al. (1996) Joos, F., M. Bruno, R. Fink, U. Siegenthaler, T. Stocker, C. L. Quéré, and J. L. Sarmiento (1996).  An efficient and accurate representation of complex oceanic and biospheric models of anthropogenic carbon uptake.  Tellus 48B, 397–417.
  * Leach et al. (2021) Leach, N. J., S. Jenkins, Z. Nicholls, C. J. Smith, J. Lynch, M. Cain, T. Walsh, B. Wu, J. Tsutsui, and M. R. Allen (2021).  FaIRv2.0.0: a generalized impulse response model for climate uncertainty and future scenario exploration.  Geoscientific Model Development 14(5), 3007–3036.
  * Lightfoot and Mamer (2014) Lightfoot, H. D. and O. A. Mamer (2014).  Calculation of atmospheric radiative forcing (warming effect) of carbon dioxide at any concentration.  Energy & Environment 25(8), 1439–1454.
  * Ljung and Box (1978) Ljung, G. M. and G. E. P. Box (1978).  On a measure of lack of fit in time series models.  Biometrika 65(2), 297–303.
  * Luderer et al. (2018) Luderer, G., Z. Vrontisi, C. Bertram, O. Y. Edelenbosch, R. C. Pietzcker, J. Rogelj, H. S. De Boer, L. Drouet, J. Emmerling, O. Fricko, S. Fujimori, P. Havlík, G. Iyer, K. Keramidas, A. Kitous, M. Pehl, V. Krey, K. Riahi, B. Saveyn, M. Tavoni, D. P. Van Vuuren, and E. Kriegler (2018).  Residual fossil CO2 emissions in 1.5–2∘C pathways.  Nature Climate Change 8(7), 626–633.
  * Masson-Delmotte and Zhai (2021) Masson-Delmotte, V. and P. Zhai (2021).  IPCC Sixth Assessment Report. Climate Change 2021: The Physical Science Basis.  Cambridge University Press.
  * Meinshausen et al. (2009) Meinshausen, M., N. Meinshausen, W. Hare, S. C. B. Raper, K. Frieler, R. Knutti, D. J. Frame, and M. R. Allen (2009).  Greenhouse-gas emission targets for limiting global warming to 2∘C.  Nature 458(7242), 1158–1162.
  * Meinshausen et al. (2020) Meinshausen, M., Z. R. J. Nicholls, J. Lewis, M. J. Gidden, E. Vogel, M. Freund, U. Beyerle, C. Gessner, A. Nauels, N. Bauer, J. G. Canadell, J. S. Daniel, A. John, P. B. Krummel, G. Luderer, N. Meinshausen, S. A. Montzka, P. J. Rayner, S. Reimann, S. J. Smith, M. van den Berg, G. J. M. Velders, M. K. Vollmer, and R. H. J. Wang (2020).  The shared socio-economic pathway (SSP) greenhouse gas concentrations and their extensions to 2500.  Geoscientific Model Development 13(8), 3571–3605.
  * Meinshausen et al. (2011) Meinshausen, M., S. C. B. Raper, and T. M. L. Wigley (2011).  Emulating coupled atmosphere-ocean and carbon cycle models with a simpler model, MAGICC6 – Part 1: Model description and calibration.  Atmospheric Chemistry and Physics 11(4), 1417–1456.
  * Meinshausen et al. (2011) Meinshausen, M., S. J. Smith, K. Calvin, J. S. Daniel, M. L. T. Kainuma, J.-F. Lamarque, K. Matsumoto, S. A. Montzka, S. C. B. Raper, K. Riahi, A. Thomson, G. J. M. Velders, and D. P. P. van Vuuren (2011).  The RCP greenhouse gas concentrations and their extensions from 1765 to 2300.  Climatic Change 109(1), 213.
  * Morice et al. (2020) Morice, C. P., J. J. Kennedy, N. A. Rayner, J. Winn, E. Hogan, R. Killick, R. Dunn, T. Osborn, P. Jones, and I. Simpson (2020).  An updated assessment of near-surface temperature change from 1850: The HadCRUT5 dataset.  Journal of Geophysical Research: Atmospheres.
  * Nicholls et al. (2021) Nicholls, Z., M. Meinshausen, and et al. (2021).  Reduced complexity model intercomparison project phase 2: Synthesizing earth system knowledge for probabilistic climate projections.  Earth’s Future 9(6), 1–25.
  * Nicholls et al. (2020) Nicholls, Z. R. J., M. Meinshausen, J. Lewis, R. Gieseke, D. Dommenget, K. Dorheim, C.-S. Fan, J. S. Fuglestvedt, T. Gasser, U. Golüke, P. Goodwin, C. Hartin, A. P. Hope, E. Kriegler, N. J. Leach, D. Marchegiani, L. A. McBride, Y. Quilcaille, J. Rogelj, R. J. Salawitch, B. H. Samset, M. Sandstad, A. N. Shiklomanov, R. B. Skeie, C. J. Smith, S. Smith, K. Tanaka, J. Tsutsui, and Z. Xie (2020).  Reduced Complexity Model Intercomparison Project Phase 1: introduction and evaluation of global-mean temperature response.  Geoscientific Model Development 13(11), 5175–5190.
  * Pretis (2020) Pretis, F. (2020).  Econometric modelling of climate systems: The equivalence of energy balance models and cointegrated vector autoregressions.  Journal of Econometrics 214(1), 256–273.
  * Ramaswami (2001) Ramaswami, V. (2001).  Chapter 6: Radiative forcing of climate change.  Technical report, IPCC.
  * Sanderson et al. (2016) Sanderson, B. M., B. C. O’Neill, and C. Tebaldi (2016).  What would it take to achieve the Paris temperature targets?  Geophysical Research Letters 43(13), 7133–7142.
  * Schwartz (2007) Schwartz, S. E. (2007).  Heat capacity, time constant, and sensitivity of Earth’s climate system.  Journal of Geophysical Research 112(D24), 1–12.
  * Schwarz (1978) Schwarz, G. (1978).  Estimating the dimension of a model.  The Annals of Statistics 6(2), 461–464.
  * Shi (1992) Shi, G. (1992).  Radiative forcing and greenhouse effect due to the atmospheric trace gases.  Science in China (Series B) 35, 217–229.
  * Smith et al. (2024) Smith, C., D. P. Cummins, H.-B. Fredriksen, Z. Nicholls, M. Meinshausen, M. Allen, S. Jenkins, N. Leach, C. Mathison, and A.-I. Partanen (2024).  fair-calibrate v1.4.1: calibration, constraining and validation of the fair simple climate model for reliable future climate projections.  EGUsphere, 1–36.
  * Smith et al. (2021) Smith, C., Z. Nicholls, K. Armour, W. Collins, P. Forster, M. Meinshausen, M. Palmer, and M. Watanabe (2021).  The earth’s energy budget, climate feedbacks, and climate sensitivity supplementary material.  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 Z. B. (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. Cambridge: Cambridge University Press. In Press.
  * Solomon et al. (2007) Solomon, S., D. Qin, and M. Manning (2007).  IPCC Fourth Assessment Report. Climate Change 2007: The Physical Science Basis.  Cambridge University Press.
  * Stocker and Qin (2013) Stocker, T. F. and D. Qin (2013).  IPCC Fifth Assessment Report. Climate Change 2013: The Physical Science Basis.  Cambridge University Press.
  * Tokarska and Gillett (2018) Tokarska, K. B. and N. P. Gillett (2018).  Cumulative carbon emissions budgets consistent with 1.5 $\,{}^{\circ}$C global warming.  Nature Climate Change 8(4), 296–299.
  * UNFCCC (2015) UNFCCC (2015).  Adoption of the Paris Agreement, https://unfccc.int/resource/docs/2015/cop21/eng/l09r01.pdf.  Technical report, United Nations Framework Convention on Climate Change.
  * Zhang et al. (2021) Zhang, X., Y.-P. Wang, P. Rayner, P. Ciais, K. Huang, Y. Luo, S. Piao, Z. Wang, J. Xia, W. Zhao, X. Zheng, J. Tian, and Y. Zhang (2021).  A small climate-amplifying effect of climate-carbon cycle feedback.  Nature Communications 12, 2952.

##  Appendix A Details of the state space representation of Stat-RCM

Let
$y_{t}=(C_{t}^{*},S_{t}^{OCN,*},S_{t}^{LND,*},F_{t}^{CO2,*},T_{t}^{m,*},T_{t}^{d,*},O_{t}^{*})^{\prime}$
denote the $7\times 1$ vector of observations at time $t$ and $\Delta>0$ the
length of a time step between observations. In Section 3.3 we presented the
non-linear state space model of the climate system, given as follows

$y_{t}\ =\ \mu+Ax_{t}+\epsilon_{t},\qquad\qquad x_{t+\Delta}\ =\
B(x_{t})+W_{t}+R\eta_{t,\Delta},$

where
$x_{t}=(C_{t},S_{t}^{OCN},S_{t}^{LND},F_{t}^{CO2},T_{t}^{m},T_{t}^{d})^{\prime}$
denotes the $6\times 1$ latent state vector. At time $t$, there are thus $7$
observations, collected in $y_{t}$, for $6$ states, collected in $x_{t}$. The
$7\times 1$ vector
$\mu=(\mu_{C},\mu_{L},\mu_{O},\mu_{F},\mu_{m},\mu_{d},H_{d}\cdot\mu_{d})^{\prime}$
contains the intercepts in the measurement equations, and the $7\times 6$
matrix $A$ captures the relations between the observations $y_{t}$ and the
underlying state vector $x_{t}$ described in Section 3.2,

$\displaystyle A=\begin{pmatrix}1&0&0&0&0&0\\\ 0&1&0&0&0&0\\\ 0&0&1&0&0&0\\\
0&0&0&1&0&0\\\ 0&0&0&0&1&0\\\ 0&0&0&0&0&1\\\ 0&0&0&0&0&H_{d}\end{pmatrix}.$

The entries in the $7\times 1$ measurement error vector $\epsilon_{t}$ are
modelled as AR(1) processes,

$\epsilon_{t+\Delta}=\Phi\epsilon_{t}+\xi_{t,\Delta},$

with $\Phi$ being a $7\times 7$ diagonal matrix given by
$\Phi=\text{diag}(\phi_{1},\phi_{2},\ldots,\phi_{7})$, with
$\phi_{i}\in(-1,1)$ for all $i$, and $\xi_{t,\Delta}\stackrel{{\scriptstyle
iid}}{{\sim}}N(0,\Delta\cdot P)$ with $P$ being the $7\times 7$ matrix

$\displaystyle P=\begin{pmatrix}\sigma^{2}_{\epsilon,C}&0&0&0&0&0&0\\\
0&\sigma^{2}_{\epsilon,OCN}&0&0&0&0&0\\\
0&0&\sigma^{2}_{\epsilon,LND}&0&0&0&0\\\
0&0&0&\sigma^{2}_{\epsilon,F}&0&0&0\\\ 0&0&0&0&\sigma^{2}_{\epsilon,m}&0&0\\\
0&0&0&0&0&\sigma^{2}_{\epsilon,d}&\rho\sigma_{\epsilon,d}\sigma_{\epsilon,OHC}\\\
0&0&0&0&0&\rho\sigma_{\epsilon,d}\sigma_{\epsilon,OHC}&\sigma_{\epsilon,OHC}^{2}\end{pmatrix},$

where $\rho\in(-1,1)$ allows for correlation in the measurement errors of deep
ocean temperature $T_{t}^{d,*}$ and ocean heat content $O_{t}^{*}$, see
Bennedsen et al. (2023).

The non-linear $6\times 1$ mapping $B(\cdot)$ represents the state transition
equations as described in Section 3,

$\displaystyle B(x_{t})=\begin{pmatrix}C_{t}\\\
b_{1}C_{t}\exp(-c_{1}T_{t}^{d})\\\ b_{2}C_{t}\exp(-c_{2}T_{t}^{d})\\\
f_{1}\log(C_{t}+f_{2}C_{t}^{2})+f_{3}\sqrt{C_{t}}\\\
\left(1-\frac{\gamma+\lambda}{H_{m}}\Delta\right)T^{m}_{t}+\frac{\gamma\Delta}{H_{m}}T^{d}_{t}\\\
\frac{\gamma\Delta}{H_{d}}T^{m}_{t}+\left(1-\frac{\gamma\Delta}{H_{d}}\right)T^{d}_{t}\end{pmatrix},$

while the $6\times 1$ vector $W_{t}$ collects constants and covariates into
the state equation,

$\displaystyle W_{t}=\begin{pmatrix}C_{t}-\Delta(b_{1}+b_{2})C_{1750}+\Delta
E_{t+\Delta}\\\ -b_{1}C_{1750}\\\ -b_{2}C_{1750}\\\
-f_{1}\log(C_{1750}+f_{2}C_{1750}^{2})+f_{3}\sqrt{C_{1750}}\\\
\frac{\Delta}{H_{m}}(F_{t}^{CO2}+F_{t}^{Ex})\\\ 0\end{pmatrix}.$

The $5\times 1$ innovation sequence $\eta_{t,\Delta}\stackrel{{\scriptstyle
iid}}{{\sim}}N(0,\Delta\cdot Q)$ is modelled with a $5\times 5$ diagonal
covariance matrix
$Q=\text{diag}(\sigma_{\eta,OCN}^{2},\sigma_{\eta,LND}^{2},\sigma_{\eta,F}^{2},\sigma_{\eta,m}^{2},\sigma_{\eta,d}^{2})$
and the $6\times 5$ matrix $R$ completes the specification for the state
variables,

$\displaystyle R=\begin{pmatrix}-\Delta&-\Delta&0&0&0\\\ 1&0&0&0&0\\\
0&1&0&0&0\\\ 0&0&1&0&0\\\ 0&0&0&1&0\\\ 0&0&0&0&1\end{pmatrix}.$

Letting $\phi=(\phi_{1},\ldots,\phi_{7})^{\prime}$,
$\sigma_{\epsilon}^{2}=(\sigma_{\epsilon,C}^{2},\sigma_{\epsilon,OCN}^{2},\sigma_{\epsilon,LND}^{2},\sigma_{\epsilon,F}^{2},\sigma_{\epsilon,m}^{2},\sigma_{\epsilon,d}^{2},\sigma_{\epsilon,O}^{2})^{\prime}$,
$\sigma_{\eta}^{2}=(\sigma_{\eta,OCN}^{2},\sigma_{\eta,LND}^{2},\sigma_{\eta,F}^{2},\sigma_{\eta,m}^{2},\sigma_{\eta,d}^{2})$
we arrive at the $37\times 1$ vector $\theta$ of unknown parameters, given by

$\displaystyle\theta=(b_{1},b_{2},c_{1},c_{2},f_{1},f_{2},f_{3},\gamma,\lambda,H_{m},H_{d},\mu_{C},\mu_{L},\mu_{O},\mu_{F},\mu_{m},\mu_{d},\rho,\phi^{\prime},\sigma_{\epsilon}^{\prime},\sigma_{\eta}^{\prime})^{\prime}.$

Note that we in our implementation of the Stat-RCM impose several restrictions
on the elements in $\theta$, so that the resulting dimension of the parameter
vector is smaller than $37$. Our preferred specification of the Stat-RCM sets
$\mu_{C}=\mu_{L}=\mu_{O}=\mu_{F}=0$ and $f_{2}=f_{3}=0$, resulting in $\theta$
being a $31\times 1$ vector, see Section 5 of the main paper.

In our implementation, we estimate the parameter vector $\theta$ by the method
of maximum likelihood, using the extended Kalman filter (Durbin and Koopman,
2012, Chapter 10). This procedure requires the Jacobian of the transition
function $B(\cdot)$, which is given by

$\displaystyle\frac{\partial B(x_{t})}{\partial
x_{t}^{\prime}}=\begin{pmatrix}1&0&0&0&0&0\\\
b_{1}\exp(-c_{1}T_{t}^{d})&0&0&0&0&-b_{1}c_{1}C_{t}\exp(c_{1}T_{t}^{d})\\\
b_{2}\exp(-c_{2}T_{t}^{d})&0&0&0&0&-b_{2}c_{2}C_{t}\exp(c_{2}T_{t}^{d})\\\
\frac{f_{1}}{C_{t}+f_{2}C_{t}^{2}}\left(1+2f_{2}C_{t}\right)+\frac{f_{3}}{2\sqrt{C_{t}}}&0&0&0&0&0\\\
0&0&0&0&\left(1-\frac{\gamma+\lambda}{H_{m}}\Delta\right)&\frac{\gamma\Delta}{H_{m}}\\\
0&0&0&0&\frac{\gamma\Delta}{H_{d}}&\left(1-\frac{\gamma\Delta}{H_{d}}\right)\end{pmatrix}.$

##  Appendix B Additional parameter estimates from the Stat-RCM

Section 5.3 of the main paper presents the estimates of the physical
parameters of the Stat-RCM when applied to the historical data over the period
$1959$–$2019$. The estimated parameters relating to the state ($\eta_{\cdot}$)
and measurement ($\epsilon_{\cdot}$) error processes are shown in Tables 5–6.
We see that some estimates of the variances of the state disturbances
$\eta_{\cdot}$ are estimated to be close to zero, in particular
$\widehat{\sigma}_{\eta,F}^{2},\widehat{\sigma}_{\eta,m}^{2}\approx 0$. This
indicates that for these state variables – corresponding to forcing from CO2,
$F_{t}^{CO2}$, and surface temperature, $T_{t}^{m}$ – the internal state
dynamics as specified by the three modules of the Stat-RCM are sufficient to
capture the evolution of the states. For the remaining state variables, we
have
$\widehat{\sigma}_{\eta,OCN}^{2},\widehat{\sigma}_{\eta,LND}^{2},\widehat{\sigma}_{\eta,d}^{2}>0$,
indicating that for these states, an additional disturbance term must be
included to capture the evolution of the states of the sinks $S_{t}^{OCN}$ and
$S_{t}^{LND}$ and the deep ocean temperature $T_{t}^{d}$.

It is also worth remarking that $\widehat{\rho}=0.87$, highlighting the close
connection between the data on deep ocean temperature ($T_{t}^{d,*}$) and
ocean heat content ($O_{t}^{*}$). We refer to Bennedsen et al. (2023) for a
more in-depth discussion on the relation between these two time series.

Table 5: Stat-RCM parameter estimates II

| $\sigma_{\eta,OCN}^{2}$ | $\sigma_{\eta,LND}^{2}$ | $\sigma_{\eta,F}^{2}$ | $\sigma_{\eta,m}^{2}$ | $\sigma_{\eta,d}^{2}$ | $\sigma_{\epsilon,C}^{2}$ | $\sigma_{\epsilon,OCN}^{2}$ | $\sigma_{\epsilon,LND}^{2}$ | $\sigma_{\epsilon,F}^{2}$ | $\sigma_{\epsilon,m}^{2}$ | $\sigma_{\epsilon,d}^{2}$ | $\sigma_{\epsilon,O}^{2}$ | $\rho$
---|---|---|---|---|---|---|---|---|---|---|---|---|---
Estimate: | $0.00$ | $0.70$ | $0.00$ | $0.00$ | $0.01$ | $0.45$ | $0.12$ | $0.42$ | $0.01$ | $0.09$ | $0.00$ | $0.04$ | $0.87$
Std. Err.: | $(0.03)$ | $(0.08)$ | $(0.00)$ | $(0.02)$ | $(0.00)$ | $(0.06)$ | $(0.01)$ | $(0.07)$ | $(0.00)$ | $(0.01)$ | $(0.00)$ | $(0.02)$ | $(0.21)$
$t$-stat: | $0.03$ | $9.17$ | $0.00$ | $0.01$ | $11.11$ | $7.20$ | $11.01$ | $5.73$ | $8.97$ | $11.07$ | $1.18$ | $2.71$ | $4.16$

Parameter estimates by the maximum likelihood method applied to the
$1959$–$2022$ data.

Table 6: Stat-RCM: Parameter estimates III

| $\phi_{1}$ | $\phi_{2}$ | $\phi_{3}$ | $\phi_{4}$ | $\phi_{5}$ | $\phi_{6}$ | $\phi_{7}$
---|---|---|---|---|---|---|---
Estimate: | $0.78$ | $0.56$ | $0.39$ | $0.57$ | $0.15$ | $0.89$ | $0.89$
Std. Err.: | $(0.08)$ | $(0.11)$ | $(0.21)$ | $(0.15)$ | $(0.13)$ | $(0.07)$ | $(0.07)$
$t$-stat (Asym): | $10.00$ | $4.92$ | $1.90$ | $3.80$ | $1.16$ | $12.12$ | $12.11$

Parameter estimates by the maximum likelihood method applied to the
$1959$–$2022$ data.