# Continuous-time state-space methods for $\delta^{18}$O and $\delta^{13}$C

Mikkel Bennedsen(a), Eric Hillebrand(a∗), Siem Jan Koopman(b), Kathrine
Larsen(a)

Authors in alphabetical order
(a)Department of Economics and Business Economics, Aarhus University, Denmark.
(a∗)Corresponding author<EMAIL_ADDRESS>
(c)Department of Econometrics, Vrije Universiteit Amsterdam, The Netherlands

###### Abstract

Time series analysis of $\delta^{18}$O and $\delta^{13}$C measurements from
benthic foraminifera for purposes of paleoclimatology is challenging. The time
series reach back tens of millions of years, they are relatively sparse in the
early record and relatively dense in the later, the time stamps of the
observations are not evenly spaced, and there are instances of multiple
different observations at the same time stamp. The time series appear non-
stationary over most of the historical record with clearly visible temporary
trends of varying directions. In this paper, we propose a continuous-time
state-space framework to analyze the time series. State space models are
uniquely suited for this purpose, since they can accommodate all the
challenging features mentioned above. We specify univariate models and joint
bivariate models for the two time series of $\delta^{18}$O and $\delta^{13}$C.
The models are estimated using maximum likelihood by way of the Kalman filter
recursions. The suite of models we consider has an interpretation as an
application of the Butterworth filter. We propose model specifications that
take the origin of the data from different studies into account and that allow
for a partition of the total period into sub-periods reflecting different
climate states. The models can be used, for example, to impute evenly time-
stamped values by way of Kalman filtering. They can also be used, in future
work, to analyze the relation to proxies for CO2 concentrations.

##  1 Introduction

Measurements of the deviation in the ratio of the oxygen isotopes 18O and 16O
in ocean sediment cores containing fossil benthic foraminifera allow for
obtaining a time series record covering the Cenozoic, i.e., the last 66
million years (Westerhold et al., 2020; Zachos et al., 2001). Since the ratio
of the isotopes depends on prevailing temperatures and the amount of water
locked in ice sheets at the time of formation of the shell, these measurements
can serve as a proxy for global temperatures. Similarly, deviations in the
ratio of the carbon isotopes 13C and 12C are used as a proxy for the abundance
of vegetation.

In this paper, we propose a class of state space models for the time series of
$\delta^{18}$O and $\delta^{13}$C. This class of models addresses the central
challenges that the time series pose: non-stationarity with trends in changing
directions, irregular spacing of the time stamps of the observations, multiple
observations at certain time stamps, different studies of origin of the data,
different climate states covered by the data, and covariation in the two
series, so that one can inform the estimation of a model for the other. The
models are specified and estimated by maximum likelihood in the time domain,
and they have a direct interpretation as Butterworth filters in the frequency
domain (Gómez, 2001; Harvey and Trimbur, 2003). Their theory is well
developed, see, for example, Harvey (1989) and Durbin and Koopman (2012) for
textbook treatments, and they can be estimated in relatively short time on
laptop computers.

The data set employed in this paper is provided in the study by Westerhold et
al. (2020), where $\delta^{18}$O and $\delta^{13}$C measurements from 34
different studies are collated and harmonized in a single data file. This file
contains 24321 time-stamped lines. The time stamps range from 67.10113 million
years ago (MYA) to 0.000564 MYA. In this paper, we ignore uncertainty in the
time stamps. We order the observations such that the series start at -67.10113
million years and end at -0.000564 million years. Of the 24321 time stamps,
23722 are unique, meaning that there are 599 time stamps with more than one
observation on either or both of the series.

Figure 1: Left panel: $\delta^{18}$O data, middle panel: $\delta^{13}$C data,
right panel: Time series of differences in the time stamp. Note that the
$y$-axis for $\delta^{18}$O is reversed, following common practice.

We use the correlation-corrected observations (column L “benthic d18O VPDB
Corr” and column K “benthic d13C VPDB Corr”) from the data file, for which
there are four time stamps for which neither $\delta^{18}$O nor $\delta^{13}$C
are available (unique observation numbers 19043, 19054, 19057, 21089; data
file observation numbers 2823, 4895, 4898, 4909). The maximum number of
observations at any time stamp is four, often, but not always, originating
from different studies.

Figure 1 shows the time series of $\delta^{18}$O (left panel) and
$\delta^{13}$C (middle panel). The right panel shows the time series of
differences in the time stamps of the observations in million years. The
smallest difference that occurs is a single year (unique obs. no. 22130, data
file obs. no. 1778), the biggest is 115366 years (unique obs. no. 7287, data
file obs. no. 16932). The plot in the right panel of Figure 1 shows that the
data set is denser in the later part of the record and sparser in the earlier,
as to be expected.

We show that a data set as complex as this can be analyzed by means of state
space methods tailored to the continuous-time nature of the time stamps.
Specifically, we employ so-called unobserved component models, and we specify
the time series under study as consisting of a latent component we aim to
estimate plus measurement noise (Harvey and Koopman, 2000, 2009). Multiple
alternative measurements can occur. The Kalman filtering and smoothing
recursions allow for the extraction of this unobserved component.

We consider several alternative models that differ in the integration order of
the unobserved component. In the simplest case, the random-walk-plus-noise
model, the unobserved component is a random walk, that is, integrated of order
one. This stochastic trend is a flexible structure that can accommodate the
changing trend directions in the series and, at the same time, allow for
specifying a non-stationary time series model. With increasing orders of
integration of the latent component, the degree of smoothing of the extracted
component increases. We describe the relation of the integration orders to
orders of a Butterworth filter applied in the frequency domain, and we show
the gain functions implied by the maximum likelihood estimation of our time-
domain models (Gómez, 1999, 2001; Harvey and Trimbur, 2003).

The separate specification of measurement and transition equations in a state
space model, with their own respective error processes, allows for reflecting
the different sources of the data (34 different studies, or, alternatively, 23
different benthic foraminifera species), by way of specifying different
measurement equation variances. The variances of the random variables that
drive the transition equation are differentiated according to the six
different climate states identified in Westerhold et al. (2020). We consider
bivariate models for $\delta^{18}$O and $\delta^{13}$C jointly, such that we
can exploit their covariation to inform the estimation of our models.

The rest of the paper is organized as follows: In Section 2, we apply random-
walk-plus-noise models to the data set, both without and with differentiating
variances. We consider univariate and bivariate specifications. In Section 3,
we describe the relation of unobserved component models with the Butterworth
filter. In Section 4, we explore higher integration orders of the unobserved
component, which correspond to higher orders of Butterworth filters. Section 6
concludes.

##  2 Random walk plus noise models

We propose the class of unobserved component models for $\delta^{18}$O and
$\delta^{13}$C, since these models can accommodate non-stationarity, unevenly
spaced observations, missing data, and multiple observations for the same
object of interest. Unobserved component models can be written in univariate
form as

$y_{t}=\mu_{t}+\psi_{t},$

where $\mu_{t}$ is the “unobserved component” of interest that is observed
with a random error $\psi_{t}$. The time index $t$ is an element of an ordered
set that can be evenly or unevenly spaced. We can think of $\mu_{t}$ as the
“signal” and $\psi_{t}$ as the “noise”.

The simplest specification in this class that we are going to consider for the
two time series studied here is the random walk plus noise model

$\displaystyle y_{t}$ $\displaystyle=\mu_{t}+\varepsilon_{t},$
$\displaystyle\varepsilon_{t}\stackrel{{\scriptstyle
i.i.d.}}{{\sim}}\mathsf{N}(0,\sigma_{\varepsilon}^{2}),$
$\displaystyle\mu_{t+\Delta t}$ $\displaystyle=\mu_{t}+\eta_{t},$
$\displaystyle\eta_{t}\stackrel{{\scriptstyle
i.i.d.}}{{\sim}}\mathsf{N}(0,\sigma_{\eta}^{2}\Delta t),$ (1)

$t=\sum_{s<t}\Delta s$. In this case $\psi_{t}=\varepsilon_{t}$ is simply an
identically and independently distributed normal random variable. The “signal”
$\mu_{t}$ is a random walk and thus non-stationary. The time increment $\Delta
t$ is varying with time. As an example, the first three and the last three
observations in the $\delta^{18}$O and $\delta^{13}$C series are shown in
Table 1, where the time increments are reported in the third column.

Table 1: First three and last three observations of the $\delta^{18}$O and $\delta^{13}$C time series with time stamps and time increments. The data are from the accompanying data file of Westerhold et al. (2020), file name aba6853_tables_s8_s34.xlsx, tab “Table S33”, column L “benthic d18O VPDB Corr” and column K “benthic d13C VPDB Corr”. The data are sorted according to the time stamp in column D “age_tuned”. Obs | Obs (data file) | time (MYA) | $\Delta t$ |  $\delta^{18}$O |  $\delta^{13}$C
---|---|---|---|---|---
1 | 24321 | 67.101133 | NA | 0.800 | 1.376
2 | 24320 | 67.098975 | 0.002158 | 0.800 | 1.343
3 | 24319 | 67.096818 | 0.002157 | 0.800 | 1.320
⋮ |  | ⋮ | ⋮ | ⋮ | ⋮
23720 | 3 | 0.003099 | 0.001032 | 3.56 | -0.07
23721 | 2 | 0.002689 | 0.000410 | 3.57 | 0.43
23722 | 1 | 0.000564 | 0.002125 | 3.52 | -0.80

The variance of a random walk and its continuous-time counterpart Brownian
motion scales with time (e.g. Karatzas and Shreve, 1991, p. 47). The variance
of an increment of Brownian motion $\sigma B_{t}$ is given by

$\mathds{E}[\sigma(B_{t}-B_{s})]^{2}=\sigma^{2}(t-s),\,0\leq s<t,$

and this motivates our definition of the variance of $\eta_{t}$. Loosely
speaking, the random variable $\eta_{t}$ captures the variation occurring in
the process $\mu_{t}$ in the time interval $[t,t+\Delta t]$, and so if
$\mu_{t}$ is to model a random walk, its variance on this increment must be
$\sigma_{\eta}^{2}\Delta t$.

Model (1) can be estimated by maximum likelihood derived from the distribution
assumptions. The evaluation of the likelihood is a byproduct of the Kalman
filter recursions, see Durbin and Koopman (2012), in particular Ch. 7, for a
textbook treatment. The first equation in model (1) constitutes a measurement
equation and the second a transition equation of a state space model. The
state variable $\mu_{t}$ is unobserved and inferred from observations $y_{t}$.
The parameters to be estimated by maximum likelihood are
$\sigma_{\varepsilon}^{2}$ and $\sigma_{\eta}^{2}$. Estimates of the
unobserved state $\mu_{t}$ can be obtained in predicted form
($\mathds{E}[\mu_{t+\Delta t}|\\{y_{s},\;s\leq t\\}]$), filtered form
($\mathds{E}[\mu_{t+\Delta t}|\\{y_{s},\;s\leq t+\Delta t\\}]$), and smoothed
form ($\mathds{E}[\mu_{t+\Delta t}|\\{y_{s},\;s\leq\sum\Delta t\\}]$) from the
Kalman filtering and smoothing recursions.

We implement the model with diffuse initialization for the non-stationary
state process $\mu_{t}$ in the KFAS package in R (Helske, 2017, 2019) and in
OxMetrics/SsfPack (Doornik, 2021).

Figure 2: Top left panel: $\delta^{18}$O data (blue) and smoothed state
$\mathds{E}[\mu_{t+\Delta t}|\\{y_{s},\;s\leq\sum\Delta t\\}]$ (red) according
to model (1) from Kalman smoothing recursions. Top right panel: same for
$\delta^{13}$C. Note that the $y$-axis for $\delta^{18}$O is reversed,
following common practice. Bottom left panel: standardized residuals for
$\delta^{18}$O, bottom right panel: standardized residuals for $\delta^{13}$C.

Figure 2 shows the data for $\delta^{18}$O and $\delta^{13}$C together with
the smoothed values of $\mu_{t}$ from the Kalman smoothing recursions in the
top row. The bottom row shows the standardized prediction residuals from the
maximum likelihood estimation (Durbin and Koopman, 2012, p. 38, eq. (2.65)).
Clearly, these are not i.i.d. with constant variance, and while the standard
errors obtained in the estimation give some sense of the dispersion of the
parameter estimate, they cannot be interpreted in the sense of a central limit
theorem. With this caveat, the maximum likelihood estimates of model (1) are
reported in Table 2. The Bayes information criterion is a measure of the fit
that penalizes for the number of parameters. This becomes important when we
increase the complexity of the model in the following.

Table 2: Maximum likelihood parameter estimates for model (1). Standard errors in parentheses. BIC: Bayes information criterion. |  $\delta^{18}$O |  $\delta^{13}$C
---|---|---
$\hat{\sigma}_{\eta}^{2}$ | 1.8364 (0.0480) | 1.2135 (0.0403)
$\hat{\sigma}_{\varepsilon}^{2}$ | 0.0205 (0.0002) | 0.0340 (0.0004)
log-likelihood | 7218.57 | 2992.10
BIC | -14416.99 | -5964.05

###  2.1 Differentiating measurement equation variances according to studies
and species

The data set prepared in Westerhold et al. (2020) originates from 34 different
sources, according to the data file (column “benthic Source” in tab “Table
S33” of the file aba6853_tables_s8_s34.xlsx). It seems plausible that the
measurement error variance is at least in part dependent on the source the
data was taken from. This dependence is implemented in model (1) by modifying
the measurement equation to

$y_{i,t}=\mu_{t}+\varepsilon_{i,t},\quad\varepsilon_{i,t}\stackrel{{\scriptstyle
i.i.d.}}{{\sim}}\mathsf{N}(0,\sigma_{\varepsilon,i}^{2}),\quad
i=1,\ldots,34,\\\ $ (2)

where $i$ indicates the study of origin of a given data point. This
generalization of the measurement equation extends the parameter vector to be
estimated by maximum likelihood to 35 coordinates (the 34
$\sigma_{\varepsilon,i}^{2}$ plus $\sigma_{\eta}^{2}$).

Figure 3: Top left panel: $\delta^{18}$O data (blue) and smoothed state
$\mathds{E}[\mu_{t+\Delta t}|\\{y_{s},\;s\leq\sum\Delta t\\}]$ (red) according
to model (2) from Kalman smoothing recursions. Top right panel: same for
$\delta^{13}$C. Note that the $y$-axis for $\delta^{18}$O is reversed,
following common practice. Bottom left panel: standardized residuals for
$\delta^{18}$O, bottom right panel: standardized residuals for $\delta^{13}$C.

Figure 3, top row, shows the data and the smoothed values of $\mu_{t}$ from
the Kalman recursions. The bottom row shows the standardized prediction
residuals, which appear to have less time variation in the variance in
comparison with Figure 2. Table 3 reports the maximum likelihood estimates.
The fit as measured by the log-likelihood improves substantially compared with
model (1). The same holds for the BIC, despite the much higher number of
parameters.

A pattern is apparent that the measurement equation variance decreases with
year of publication of the study. A linear regression of the maximum
likelihood estimates of the study variances for the $\delta^{18}$O series on
the year of publication and a constant yields a highly significantly negative
slope estimate of about -0.001 per year. The adjusted coefficient of
determination of the regression is 22%, and the F-test of the regression has a
$p$-value of 0.002. The study coded as “McCarren et al. 2008” in the data file
is a bit of an outlier. Excluding this study and its estimated variance from
the regression yields an even more highly significant slope estimate of about
the same size and an adjusted coefficient of determination of 35%. The same
regression for the estimated variances from the $\delta^{13}$C series yields a
highly significantly negative slope estimate of -0.004, an adjusted
coefficient of determination of 18%, and an F-test of the regression with
$p$-value 0.007. Excluding the “McCarren et al. 2008” study yields a highly
significant slope estimate of -0.0035, an adjusted coefficient of
determination of 45%, and an F-test of the regression with $p$-value 0.00001.

Table 3: Maximum likelihood estimates of model (2). Standard errors in parentheses. “Source” refers to the entry in column “benthic Source” in tab “Table S33” of the data file aba6853_tables_s8_s34.xlsx of Westerhold et al. (2020). The identifier “McCarren et al. 2008 et al. 2008” was merged with “McCarren et al. 2008” and “Bickert et al.1997” with “Bickert et al. 1997”. The study Westerhold et al. 2020 is coded in the data file as “this study”. BIC: Bayes information criterion. Source | No. obs. |  $\delta^{18}$O | No. obs. |  $\delta^{13}$C
---|---|---|---|---
$\sigma_{\varepsilon,i}^{2}$ |  |  |  |
“Barnet et al. 2017” | 548 | 0.0015 (0.0002) | 548 | 0.0037 (0.0004)
“Hull et al. 2020” | 29 | 0.0136 (0.0039) | 29 | 0.0136 (0.0040)
“Barnet et al. 2019” | 879 | 0.0053 (0.0004) | 879 | 0.0030 (0.0003)
“Littler et al. 2014” | 2164 | 0.0095 (0.0004) | 2164 | 0.0066 (0.0003)
“McCarren et al. 2008” | 47 | 0.1064 (0.0314) | 47 | 0.3525 (0.0969)
“Lauretano et al. 2015” | 636 | 0.0343 (0.0025) | 636 | 0.0207 (0.0014)
“Stap et al. 2010” | 372 | 0.0305 (0.0025) | 372 | 0.0274 (0.0023)
“Thomas et al. 2018” | 86 | 0.0159 (0.0028) | 86 | 0.0203 (0.0035)
“Lauretano et al. 2016; 2018” | 541 | 0.0072 (0.0007) | 541 | 0.0113 (0.0010)
“Sexton et al. 2011” | 301 | 0.0124 (0.0016) | 301 | 0.0164 (0.0019)
Westerhold et al. 2020 | 1723 | 0.0072 (0.0004) | 1723 | 0.0047 (0.0003)
“Westerhold et al. 2018 LLTM” | 30 | 0.0132 (0.0045) | 30 | 0.0575 (0.0182)
“Boscolo Galazzo et al. 2014” | 34 | 0.0227 (0.0065) | 34 | 0.0140 (0.0041)
“Lear et al. 2004” | 182 | 0.0072 (0.0013) | 182 | 0.0069 (0.0012)
“Coxall et al. 2005” | 831 | 0.0089 (0.0006) | 831 | 0.0091 (0.0006)
“Coxall & Wilson, 2011” | 413 | 0.0195 (0.0017) | 413 | 0.0115 (0.0010)
“Riesselmann et al. 2007” | 43 | 0.0161 (0.0046) | 43 | 0.0354 (0.0092)
“Paelike2006 PI Herrle” | 258 | 0.0211 (0.0023) | 258 | 0.0135 (0.0016)
“Wade and Paelike 2004” | 152 | 0.0432 (0.0058) | 152 | 0.0326 (0.0047)
“Liebrand et al. 2011;2016” | 3251 | 0.0114 (0.0004) | 3251 | 0.0107 (0.0003)
“Holbourn et al. 2015” | 892 | 0.0207 (0.0014) | 892 | 0.0107 (0.0007)
“Tian et al. 2014” | 61 | 0.0461 (0.0096) | 61 | 0.0283 (0.0061)
“Holbourn et al. 2014” | 1370 | 0.0144 (0.0007) | 1370 | 0.0086 (0.0005)
“Holbourn et al. 2007;2013;2018” | 1987 | 0.0064 (0.0003) | 1987 | 0.0076 (0.0003)
“Drury et al. 2017” | 1129 | 0.0084 (0.0005) | 1129 | 0.0183 (0.0009)
“Tian et al. 2018” | 551 | 0.0110 (0.0009) | 551 | 0.0276 (0.0019)
UniCambridge Hodell / Westerhold et al. 2020 | 547 | 0.0152 (0.0011) | 547 | 0.0109 (0.0008)
“Bell et al. 2014” | 260 | 0.0072 (0.0008) | 260 | 0.0169 (0.0017)
MARUM Bremen / Westerhold et al. 2020 | 45 | 0.0092 (0.0021) | 45 | 0.0129 (0.0029)
“Franz & Tiedemann2002” | 1304 | 0.0648 (0.0027) | 1050 | 0.1204 (0.0057)
“Tiedemann & Franz 1997” | 330 | 0.0353 (0.0032) | 264 | 0.0595 (0.0062)
“Billups et al. 1998” | 538 | 0.0313 (0.0022) | 538 | 0.0614 (0.0041)
“Bickert et al. 1997” | 2562 | 0.0611 (0.0022) | 2562 | 0.1561 (0.0047)
“deMenocal et al. 1997” | 163 | 0.0265 (0.0035) | 163 | 0.1398 (0.0169)
$\sigma_{\eta}^{2}$ |  | 1.6235 (0.0506) |  | 1.1547 (0.0368)
log-likelihood | 24259 | 10036.92 | 23939 | 9750.60
BIC |  | -19719.83 |  | -19150.07

An alternative specification that takes the different origins of the data into
account differentiates the variances in the measurement equation according to
benthic foraminifera species. There are 23 distinct species listed in column I
(“benthic Species”) of the data file accompanying Westerhold et al. (2020).
Therefore, we can modify the measurement equation as follows:

$y_{i,t}=\mu_{t}+\varepsilon_{i,t},\quad\varepsilon_{i,t}\stackrel{{\scriptstyle
i.i.d.}}{{\sim}}\mathsf{N}(0,\sigma_{\varepsilon,i}^{2}),\quad
i=1,\ldots,23,\\\ $ (3)

where $i$ now indicates the different species of benthic foraminifera. In the
appendix, we report graphs of the smoothed signal and maximum likelihood
estimates. They show that judging both by the maximized log-likelihood value
and by visual inspection of the standardized prediction residuals, the
specification that differentiates by studies is superior to this one.

###  2.2 Differentiating transition equation variances according to climate
states

Westerhold et al. (2020) identify six different climate states in the records
of $\delta^{18}$O and $\delta^{13}$C with the help of recurrence analysis. The
periods are shown in Table 4.

Table 4: Climate states identified in Westerhold et al. (2020).  Million Years Ago | Climate State | Obs (unique time stamp) | Obs (data file)
---|---|---|---
67.10113 to 56 | Warmhouse 2 | 1–2690 | 21561–24321
56 to 47 | Hothouse | 2691–5693 | 18531–21560
47 to 34 | Warmhouse 1 | 5694–7471 | 16744–18530
34 to 13.9 | Coolhouse 1 | 7472–13933 | 10067–16743
13.9 to 3.3 | Coolhouse 2 | 13934–20179 | 3762–10066
3.3 to 0.000564 | Icehouse | 20180–23722 | 1–3761

We employ these periods to differentiate the variances driving the unobserved
component $\mu_{t}$ according to climate states. For example, the data plots
show that the variance of the observations seem to increase towards the later
part of the sample. This is particularly conspicuous for the “icehouse” regime
in the $\delta^{13}$C record. A priori, it is not clear whether such a
differentiation of the variances should be placed into the measurement or the
transition equation. We opt for the transition equation, for two reasons: (1)
we already differentiate for originating studies in the measurement equation
and want to avoid having to interact these two differentiating dimensions
there, and (2) the periods refer to climate states, and we aim for the
unobserved component $\mu_{t}$ to have a physical interpretation as a climate
state variable that is observed in $y_{t}$ with measurement noise. Therefore,
it seems reasonable to let $\mu_{t}$ undergo the variance shifts according to
climate states.

The transition equation in model (1) then becomes the following.

$\mu_{t+\Delta t}=\mu_{t}+\eta_{t},\quad\eta_{t}\stackrel{{\scriptstyle
i.i.d.}}{{\sim}}\mathsf{N}(0,\sigma_{\eta,j}^{2}\Delta t),\quad
j\in\\{1,\ldots,6\\},$ (4)

where $j\in\\{1,\ldots,6\\}$ refers to the six climate state regimes shown in
Table 4. The measurement equation (2) is employed, so that the model also
differentiates the studies of origin in the measurement variances.

Table 5: Maximum likelihood estimates of model (4). Standard errors in parentheses. “Source” refers to the entry in column “benthic Source” in tab “Table S33” of the data file aba6853_tables_s8_s34.xlsx of Westerhold et al. (2020). The identifier “McCarren et al. 2008 et al. 2008” was merged with “McCarren et al. 2008” and “Bickert et al.1997” with “Bickert et al. 1997”. The study Westerhold et al. 2020 is coded in the data file as “this study”. BIC: Bayes information criterion. Source | No. obs. |  $\delta^{18}$O | No. obs. |  $\delta^{13}$C
---|---|---|---|---
$\sigma_{\varepsilon,i}^{2}$ |  |  |  |
“Barnet et al. 2017” | 548 | 0.0031 (0.0003) | 548 | 0.0049 (0.0004)
“Hull et al. 2020” | 29 | 0.0127 (0.0035) | 29 | 0.0126 (0.0037)
“Barnet et al. 2019” | 879 | 0.0080 (0.0005) | 879 | 0.0037 (0.0003)
“Littler et al. 2014” | 2164 | 0.0117 (0.0005) | 2164 | 0.0068 (0.0003)
“McCarren et al. 2008” | 47 | 0.1707 (0.0518) | 47 | 0.1702 (0.0640)
“Lauretano et al. 2015” | 636 | 0.0315 (0.0023) | 636 | 0.0200 (0.0014)
“Stap et al. 2010” | 372 | 0.0308 (0.0024) | 372 | 0.0261 (0.0023)
“Thomas et al. 2018” | 86 | 0.0170 (0.0029) | 86 | 0.0189 (0.0033)
“Lauretano et al. 2016; 2018” | 541 | 0.0089 (0.0008) | 541 | 0.0086 (0.0009)
“Sexton et al. 2011” | 301 | 0.0154 (0.0018) | 301 | 0.0124 (0.0018)
Westerhold et al. 2020 | 1723 | 0.0116 (0.0005) | 1723 | 0.0066 (0.0004)
“Westerhold et al. 2018 LLTM” | 30 | 0.0194 (0.0061) | 30 | 0.0758 (0.0227)
“Boscolo Galazzo et al. 2014” | 34 | 0.0243 (0.0063) | 34 | 0.0141 (0.0039)
“Lear et al. 2004” | 182 | 0.0095 (0.0014) | 182 | 0.0081 (0.0013)
“Coxall et al. 2005” | 831 | 0.0094 (0.0006) | 831 | 0.0097 (0.0006)
“Coxall & Wilson, 2011” | 413 | 0.0196 (0.0016) | 413 | 0.0122 (0.0010)
“Riesselmann et al. 2007” | 43 | 0.0188 (0.0043) | 43 | 0.0379 (0.0093)
“Paelike2006 PI Herrle” | 258 | 0.0212 (0.0023) | 258 | 0.0140 (0.0016)
“Wade and Paelike 2004” | 152 | 0.0433 (0.0058) | 152 | 0.0329 (0.0047)
“Liebrand et al. 2011;2016” | 3251 | 0.0116 (0.0004) | 3251 | 0.0112 (0.0004)
“Holbourn et al. 2015” | 892 | 0.0215 (0.0015) | 892 | 0.0119 (0.0008)
“Tian et al. 2014” | 61 | 0.0469 (0.0097) | 61 | 0.0295 (0.0062)
“Holbourn et al. 2014” | 1370 | 0.0147 (0.0008) | 1370 | 0.0091 (0.0005)
“Holbourn et al. 2007;2013;2018” | 1987 | 0.0063 (0.0003) | 1987 | 0.0068 (0.0003)
“Drury et al. 2017” | 1129 | 0.0083 (0.0005) | 1129 | 0.0173 (0.0009)
“Tian et al. 2018” | 551 | 0.0110 (0.0009) | 551 | 0.0269 (0.0019)
UniCambridge Hodell / Westerhold et al. 2020 | 547 | 0.0151 (0.0011) | 547 | 0.0105 (0.0008)
“Bell et al. 2014” | 260 | 0.0072 (0.0008) | 260 | 0.0163 (0.0016)
MARUM Bremen / Westerhold et al. 2020 | 45 | 0.0092 (0.0021) | 45 | 0.0127 (0.0029)
“Franz & Tiedemann2002” | 1304 | 0.0620 (0.0027) | 1050 | 0.1137 (0.0056)
“Tiedemann & Franz 1997” | 330 | 0.0295 (0.0030) | 264 | 0.0562 (0.0066)
“Billups et al. 1998” | 538 | 0.0313 (0.0022) | 538 | 0.0602 (0.0041)
“Bickert et al. 1997” | 2562 | 0.0395 (0.0014) | 2562 | 0.1280 (0.0041)
“deMenocal et al. 1997” | 163 | 0.0146 (0.0027) | 163 | 0.1020 (0.0135)
$\sigma_{\eta,j}^{2}$ |  |  |  |
67.10113 – 56 MYA (“Warmhouse 2”) | 2761 | 0.3365 (0.0342) | 2761 | 0.6426 (0.0474)
56 – 47 MYA (“Hothouse”) | 3030 | 0.9574 (0.0879) | 3030 | 1.9851 (0.1335)
47 – 34 MYA (“Warmhouse 1”) | 1786 | 0.1608 (0.0242) | 1784 | 0.4026 (0.0438)
34 – 13.9 MYA (“Coolhouse 1”) | 6669 | 1.4655 (0.0859) | 6671 | 0.8437 (0.0569)
13.9 – 3.3 MYA (“Coolhouse 2”) | 6282 | 1.7187 (0.0980) | 6139 | 1.6038 (0.0991)
3.3 – 0.000564 MYA (“Icehouse”) | 3731 | 11.5945 (0.7145) | 3554 | 8.3269 (0.8329)
log-likelihood | 24259 | 10825.05 | 23939 | 10038.34
BIC |  | -21247.13 |  | -19673.73

The maximum likelihood estimates of the resulting model with climate state
variances in the transition equation and study variances in the measurement
equation are reported in Table 5. It can be seen in comparison with Table 3
that the study variances in the measurement equation are essentially
unaffected by the addition of climate-state dependent variances in the
transition equation. (The regression of the study variances on the year of
publication for $\delta^{18}$O is only significantly negative if the “McCarren
et al. 2008” study is excluded. For the $\delta^{13}$C series it remains
highly significantly negative throughout.) The climate-state dependent
variances vary strongly across the different periods. By far the highest
variance is estimated in the “icehouse” period (3.3 to 6e-4 MYA),
corresponding to the visual impression of the late record mentioned above. The
BIC indicates an additional improvement compared with model (2).

###  2.3 Bivariate random walk plus noise models

The comovement in $\delta^{18}$O and $\delta^{13}$C allows for an improvement
of the estimation of the unobserved $\mu$-components in $\delta^{18}$O and
$\delta^{13}$C by modelling the two time series jointly and estimating the
covariance of the error terms $\eta$.

To conceptualize the time dependence in the covariances of increments of
$\mu^{\delta^{18}O}_{t}$ and $\mu^{\delta^{13}C}_{t}$, consider two
independent Brownian motions $B_{1}(t)$ and $B_{3}(t)$. Let

$B_{2}(t)=\rho B_{1}(t)+\sqrt{1-\rho^{2}}B_{3}(t).$

Here, $B_{1}$ stands for the $\mu$-component in $\delta^{18}$O and $B_{2}$ for
the $\mu$-component in $\delta^{13}$C. We have for $s<t$ that

$\displaystyle\,\mbox{cov}[B_{1}(t)-B_{1}(s),B_{2}(t)-B_{2}(s)]=\mathds{E}[(B_{1}(t)-B_{1}(s))(B_{2}(t)-B_{2}(s))]$
$\displaystyle\quad=\mathds{E}[(B_{1}(t)-B_{1}(s))(\rho(B_{1}(t)-B_{1}(s))+\sqrt{1-\rho^{2}}(B_{3}(t)-B_{3}(s)))]$
$\displaystyle\quad=\rho\mathds{E}[(B_{1}(t)-B_{1}(s))^{2}]=\rho(t-s),$

because of the independence of $B_{1}$ and $B_{3}$. More generally, for
$s<t<u<v$,

$\displaystyle\,\mbox{cov}[B_{1}(v)-B_{1}(t),B_{2}(u)-B_{2}(s)]=\rho\mathds{E}[(B_{1}(v)-B_{1}(t))(B_{1}(u)-B_{1}(s))]$
$\displaystyle\quad=\rho\mathds{E}\left[\left((B_{1}(v)-B_{1}(u))+(B_{1}(u)-B_{1}(t))\right)\left((B_{1}(u)-B_{1}(t))+(B_{1}(t)-B_{1}(s))\right)\right]$
$\displaystyle\quad=\rho\mathds{E}\left[(B_{1}(u)-B_{1}(t))^{2}\right]=\rho(u-t),$

because

$\displaystyle\mathds{E}[(B_{1}(v)-B_{1}(u))(B_{1}(u)-B_{1}(t))]=\mathds{E}[(B_{1}(v)-B_{1}(u))(B_{1}(t)-B_{1}(s))]$
$\displaystyle\quad=\mathds{E}[(B_{1}(u)-B_{1}(t))(B_{1}(t)-B_{1}(s))]=0.$

Thus, the covariance of two increments of $\mu^{\delta^{18}O}$ and
$\mu^{\delta^{13}C}$ is modelled as

$\mathds{E}[(\mu^{\delta^{18}O}_{v}-\mu^{\delta^{18}O}_{t})(\mu^{\delta^{13}C}_{u}-\mu^{\delta^{13}C}_{s})]=\rho\;\text{length}((t,v)\cap(s,u)).$

In all cases in the data, either $s=t$ or $u=v$, and therefore
$\text{length}((t,v)\cap(s,u))=\min\\{v-t,u-s\\}$.

A simple bivariate random walk plus noise model for the two series that
extends equation (1) consists of the measurement equations

$\displaystyle y^{\delta^{18}O}_{t}$
$\displaystyle=\mu^{\delta^{18}O}_{t}+\varepsilon^{\delta^{18}O}_{t},$
$\displaystyle\varepsilon^{\delta^{18}O}_{t}\stackrel{{\scriptstyle
i.i.d.}}{{\sim}}\mathsf{N}(0,\sigma_{\varepsilon,\delta^{18}O}^{2}),$
$\displaystyle y_{t}^{\delta^{13}C}$
$\displaystyle=\mu_{t}^{\delta^{13}C}+\varepsilon_{t}^{\delta^{13}C},$
$\displaystyle\varepsilon_{t}^{\delta^{13}C}\stackrel{{\scriptstyle
i.i.d.}}{{\sim}}\mathsf{N}(0,\sigma_{\varepsilon,\delta^{13}C}^{2}),$

where $\varepsilon_{t}^{\delta^{18}O}$ and $\varepsilon_{t}^{\delta^{13}C}$
are independent. The interesting interaction of the bivariate specification
plays out in the covariance matrix of the transition equation, where we allow
for correlation of the transition error from the unobserved component in
$\delta^{18}$O with the one in $\delta^{13}$C:

$\displaystyle\mu_{t+\Delta t}^{\delta^{18}O}$
$\displaystyle=\mu_{t}^{\delta^{18}O}+\eta_{t}^{\delta^{18}O},$
$\displaystyle\mu_{t+\Delta t}^{\delta^{13}C}$
$\displaystyle=\mu_{t}^{\delta^{13}C}+\eta_{t}^{\delta^{13}C}.$

Collect the state disturbances in the random vector
$\eta_{t}=(\eta_{t}^{\delta^{18}O},\eta_{t}^{\delta^{13}C})^{T}$, and assume

$\eta_{t}\stackrel{{\scriptstyle i.i.d.}}{{\sim}}\mathsf{N}(0,Q_{t}\Delta
t),\qquad
Q_{t}=\left[\begin{array}[]{cc}\sigma^{2}_{\eta,\delta^{18}O}&\rho\sigma_{\eta,\delta^{18}O}\sigma_{\eta,\delta^{13}C}\\\
\rho\sigma_{\eta,\delta^{18}O}\sigma_{\eta,\delta^{13}C}&\sigma^{2}_{\eta,\delta^{13}C}\end{array}\right].$
(5)

The estimation results of the model are reported in Table 6. The results are
qualitatively very close to Table 2. The correlation parameter is estimated at
-0.17, and even though the corresponding standard error is much smaller in
magnitude, the log-likelihood is close to the sum of the two log-likelihoods
from the separate estimation. A plot of the smoothed values and residuals can
visually not be distinguished from Figure 2, and so it is omitted here.

Table 6: Maximum likelihood parameter estimates for the simplest bivariate random walk plus noise model. Standard errors in parentheses. BIC: Bayes information criterion. |  $\delta^{18}$O |  $\delta^{13}$C
---|---|---
$\hat{\sigma}_{\eta}^{2}$ | 1.8154 (0.0475) | 1.1991 (0.0396)
$\hat{\sigma}_{\varepsilon}^{2}$ | 0.0205 (0.0002) | 0.0341 (0.0004)
$\hat{\rho}$ | -0.1700 (0.0209)
log-likelihood | 10229.85
BIC | -20409.34

We generalize the bivariate approach to distinguish measurement equation
variances according to studies and transition equation variances according to
climate states in the periods in Table 4. That is, there are six different
covariance matrices of the type (5) for the different climate states and 34
study variances for $\delta^{18}$O and $\delta^{13}$C, respectively, for a
total of 86 parameters. The system matrices are spelled out in the appendix A.

Table 7: Maximum likelihood estimates of model (5). Standard errors in parentheses. “Source” refers to the entry in column “benthic Source” in tab “Table S33” of the data file aba6853_tables_s8_s34.xlsx of Westerhold et al. (2020). The identifier “McCarren et al. 2008 et al. 2008” was merged with “McCarren et al. 2008” and “Bickert et al.1997” with “Bickert et al. 1997”. The study Westerhold et al. 2020 is coded in the data file as “this study”. BIC: Bayes information criterion. Source | No. obs. |  $\delta^{18}$O | No. obs. |  $\delta^{13}$C
---|---|---|---|---
$\sigma_{\varepsilon,i}^{2}$ |  |  |  |
“Barnet et al. 2017” | 548 | 0.0026 (0.0002) | 548 | 0.0043 (0.0004)
“Hull et al. 2020” | 29 | 0.0137 (0.0038) | 29 | 0.0139 (0.0040)
“Barnet et al. 2019” | 879 | 0.0066 (0.0004) | 879 | 0.0029 (0.0002)
“Littler et al. 2014” | 2164 | 0.0096 (0.0004) | 2164 | 0.0054 (0.0003)
“McCarren et al. 2008” | 47 | 0.0611 (0.0139) | 47 | 0.0783 (0.0195)
“Lauretano et al. 2015” | 636 | 0.0298 (0.0022) | 636 | 0.0184 (0.0013)
“Stap et al. 2010” | 372 | 0.0289 (0.0023) | 372 | 0.0257 (0.0022)
“Thomas et al. 2018” | 86 | 0.0134 (0.0023) | 86 | 0.0148 (0.0027)
“Lauretano et al. 2016; 2018” | 541 | 0.0053 (0.0005) | 541 | 0.0048 (0.0006)
“Sexton et al. 2011” | 301 | 0.0075 (0.0009) | 301 | 0.0045 (0.0009)
Westerhold et al. 2020 | 1723 | 0.0076 (0.0004) | 1723 | 0.0038 (0.0003)
“Westerhold et al. 2018 LLTM” | 30 | 0.0128 (0.0039) | 30 | 0.0499 (0.0146)
“Boscolo Galazzo et al. 2014” | 34 | 0.0242 (0.0064) | 34 | 0.0146 (0.0041)
“Lear et al. 2004” | 182 | 0.0083 (0.0013) | 182 | 0.0080 (0.0013)
“Coxall et al. 2005” | 831 | 0.0090 (0.0007) | 831 | 0.0096 (0.0006)
“Coxall & Wilson, 2011” | 413 | 0.0180 (0.0016) | 413 | 0.0114 (0.0010)
“Riesselmann et al. 2007” | 43 | 0.0171 (0.0043) | 43 | 0.0334 (0.0085)
“Paelike2006 PI Herrle” | 258 | 0.0210 (0.0023) | 258 | 0.0133 (0.0015)
“Wade and Paelike 2004” | 152 | 0.0422 (0.0057) | 152 | 0.0318 (0.0046)
“Liebrand et al. 2011;2016” | 3251 | 0.0109 (0.0004) | 3251 | 0.0107 (0.0004)
“Holbourn et al. 2015” | 892 | 0.0180 (0.0013) | 892 | 0.0100 (0.0007)
“Tian et al. 2014” | 61 | 0.0431 (0.0090) | 61 | 0.0264 (0.0057)
“Holbourn et al. 2014” | 1370 | 0.0134 (0.0007) | 1370 | 0.0082 (0.0006)
“Holbourn et al. 2007;2013;2018” | 1987 | 0.0061 (0.0003) | 1987 | 0.0066 (0.0003)
“Drury et al. 2017” | 1129 | 0.0083 (0.0005) | 1129 | 0.0174 (0.0009)
“Tian et al. 2018” | 551 | 0.0109 (0.0009) | 551 | 0.0270 (0.0019)
UniCambridge Hodell / Westerhold et al. 2020 | 547 | 0.0150 (0.0011) | 547 | 0.0105 (0.0008)
“Bell et al. 2014” | 260 | 0.0072 (0.0008) | 260 | 0.0162 (0.0016)
MARUM Bremen / Westerhold et al. 2020 | 45 | 0.0092 (0.0021) | 45 | 0.0128 (0.0029)
“Franz & Tiedemann2002” | 1304 | 0.0615 (0.0027) | 1050 | 0.1132 (0.0055)
“Tiedemann & Franz 1997” | 330 | 0.0307 (0.0031) | 264 | 0.0591 (0.0066)
“Billups et al. 1998” | 538 | 0.0312 (0.0022) | 538 | 0.0602 (0.0041)
“Bickert et al. 1997” | 2562 | 0.0395 (0.0014) | 2562 | 0.1309 (0.0041)
“deMenocal et al. 1997” | 163 | 0.0148 (0.0027) | 163 | 0.1019 (0.0131)
$\sigma_{\eta,j}^{2}$ |  |  |  |
67.10113 – 56 MYA (“Warmhouse 2”) | 2761 | 0.5924 (0.0522) | 2761 | 0.9224 (0.0635)
56 – 47 MYA (“Hothouse”) | 3030 | 1.9261 (0.1187) | 3030 | 3.3403 (0.1757)
47 – 34 MYA (“Warmhouse 1”) | 1786 | 0.6984 (0.0724) | 1784 | 0.9711 (0.0829)
34 – 13.9 MYA (“Coolhouse 1”) | 6669 | 1.8470 (0.1086) | 6671 | 1.0720 (0.0722)
13.9 – 3.3 MYA (“Coolhouse 2”) | 6282 | 1.7873 (0.1020) | 6139 | 1.6644 (0.1027)
3.3 – 0.000564 MYA (“Icehouse”) | 3731 | 11.5363 (0.7112) | 3554 | 7.0819 (0.6603)
$\rho$ |  |  |  |
67.10113 – 56 MYA (“Warmhouse 2”) | 0.8212 (0.0185)
56 – 47 MYA (“Hothouse”) | 0.9580 (0.0051)
47 – 34 MYA (“Warmhouse 1”) | 0.9053 (0.0122)
34 – 13.9 MYA (“Coolhouse 1”) | 0.6043 (0.0253)
13.9 – 3.3 MYA (“Coolhouse 2”) | 0.1278 (0.0365)
3.3 – 0.000564 MYA (“Icehouse”) | -0.8052 (0.0322)
log-likelihood | 22120.73
BIC | -43375.08

The correlation of the state disturbances presents itself at a much higher
resolution in Table 7 compared to Table 6. There emerges a pattern that high
positive correlation in the early record covering warm- and hothouses weakens
in the fourth and fifth climate state covering coolhouses and finally reverses
sign in the final, icehouse state. Positive correlation in the innovations of
two random walks means that there is a tendency that the series move together
in local upward or downward runs, and they move in opposite directions in the
presence of negative correlation. This pattern was found earlier in Turner
(2014), for example. Correlation is substantial in all periods except the
fifth, and thus joint modeling of $\delta^{18}$O and $\delta^{13}$C is
advantageous in the sense that the estimation of the unobserved states
$\mu_{t}^{(\delta^{18}\text{O},\delta^{13}\text{C})}$ can benefit from the
larger information set. This is reflected in the log-likelihood of the joint
model reported in Table 7, which is substantially higher than the sum of the
two separate log-likelihoods in Table 6 (20863.39), at the cost of additional
parameters. It is also apparent from a comparison of the estimated variances
in the two tables: Almost all measurement equation variances for the different
studies are reduced in the joint model compared to the separate estimation,
and in the few instances where the variance increased, the difference is not
large. (The regression of the study variances on year of publication of the
study yields highly significantly negative slope estimates for both proxy
series.) On the other hand, the variances of the transition equation for the
different climate states are higher for the first five periods in the joint
model compared to the separate models, meaning that more variation in the data
is captured by the unobserved component of interest, the “signal”, and less
variation remains in the “noise.” The BIC is more highly negative than the sum
of the BICs in Table 5. For these reasons, the joint model (5) reported in
Table 7 is the preferred model of this section on random-walk-plus-noise
models. Figures 9 to 14 in appendix C show the smoothed values of the
unobserved components $\mu_{t}^{(\delta^{18}\text{O},\delta^{13}\text{C})}$,
together with the data, and the standardized prediction residuals for the six
different climate states.

##  3 Butterworth filtering in the time domain

The models proposed in the previous section naturally lend themselves to an
interpretation where a signal $\mu_{t}$ is extracted from observations $y_{t}$
and distinguished from noise $\psi_{t}$. In this section, we present the well-
known connection to signal extraction in the frequency domain and show how the
random-walk-plus-noise models of the previous section and the integrated-
random-walk-plus-noise models of the next section are time-domain counterparts
of the Butterworth filter.

Consider a time series that consists of an autoregressive component
$\mu_{t}=\phi_{1}\mu_{t-1}+\ldots+\phi_{p}\mu_{t-p}+\eta_{t}$, where
$\eta_{t}$ is white noise with mean zero and variance $\sigma_{\eta}^{2}$,
plus a white noise contamination $\varepsilon_{t}$ with mean zero and variance
$\sigma_{\varepsilon}^{2}$ that is independent of $\eta_{t}$:

$y_{t}=\Phi^{-1}(L)\eta_{t}+\varepsilon_{t},\\\ $ (6)

where $\Phi(L)=1-\phi_{1}L-\ldots\phi_{p}L^{p}$ is the lag polynomial
corresponding to the autoregressive component with roots outside the unit
circle. For a zero-mean stationary time series such as this, define the
autocovariance-generating function

$g_{y}(z)=\sum_{j=-\infty}^{\infty}\gamma(j)z^{j},\;z\in\mathds{C},$

with $\gamma(j)=\mathds{E}(y_{t}y_{t-j})$. Since $\eta_{t},\varepsilon_{t}$
independent, we have that

$g_{y}(z)=g_{\mu}(z)+g_{\varepsilon}(z),$

where

$g_{\mu}(z)=\frac{\sigma_{\eta}^{2}}{|\Phi(z)|^{2}}=\frac{\sigma_{\eta}^{2}}{\Phi(z)\Phi(z^{-1})},$

and $g_{\varepsilon}(z)=\sigma_{\varepsilon}^{2}$. The convention
$|\Phi(z)|^{2}=\Phi(z)\Phi(z^{-1})$ is motivated by the fact that we consider
$z\in\mathds{C}$ only on the unit circle, i.e., $z=\exp(-i\lambda)$,
$\lambda\in\mathds{R}$, and all $\gamma(j)\in\mathds{R}$, thus
$|g(z)|^{2}=g(z)\overline{g(z)}=g(z)g(\bar{z})$, cf. Whittle (1983), p. 12.

The Wiener-Kolmogorov filter for the doubly-infinite sample for model (6) is
given by

$\hat{\mu}_{t}=\sum_{j=-\infty}^{\infty}b_{j}y_{t-j}=B(L)y_{t},$

with response function

$B(z)=\frac{g_{\mu}(z)}{g_{\mu}(z)+g_{\varepsilon}(z)}=\frac{1}{1+q^{-1}|\Phi(z)|^{2}},$

where $q=\sigma_{\eta}^{2}/\sigma_{\varepsilon}^{2}$ is the signal-to-noise
ratio. For the AR(1) case $\Phi(z)=1-\phi z$, and $\phi\to 1$, we get the
random walk plus noise model

$y_{t}=\sum_{j=-\infty}^{t}\eta_{j}+\varepsilon_{t}=(1-L)^{-1}\eta_{t}+\varepsilon_{t}.$

Bell (1984) showed that the Wiener-Kolmogorov filter is asymptotically valid
in the nonstationary case under standard assumptions, see also Gómez (1999).
The random walk plus noise model therefore has response function

$B(z)=\frac{1}{1+q_{1}^{-1}|1-z|^{2}}.$

The signal-to-noise ratios $q_{m}$ depend on the integration order $m$ in a
manner that is explained below. The integrated random walk plus noise model

$y_{t}=\sum_{s=-\infty}^{t}\sum_{j=-\infty}^{s}\eta_{j}+\varepsilon_{t}=(1-L)^{-2}\eta_{t}+\varepsilon_{t}.$

has response function

$B(z)=\frac{1}{1+q_{2}^{-1}|1-z|^{4}}.$

Integrating the random walk $m$ times yields

$y_{t}=(1-L)^{-m}\eta_{t}+\varepsilon_{t}$

with response function

$B(z)=\frac{1}{1+q_{m}^{-1}|1-z|^{2m}}.$

For $z=e^{-i\lambda}$, we then obtain as the gain function $G(\lambda)$ of the
$m$-th order integrated random walk plus noise model

$G(\lambda)=|B(e^{-i\lambda})|=\frac{1}{1+q_{m}^{-1}[(1-e^{-i\lambda})(1-e^{i\lambda})]^{m}}.$

Since the filter is symmetric, and the response function is non-negative, the
gain function is equal to the response function. With
$(1-e^{-i\lambda})(1-e^{i\lambda})=4\sin^{2}(\lambda/2)$, we have that

$G(\lambda)=\frac{1}{1+q_{m}^{-1}2^{2m}\sin^{2m}(\lambda/2)},$ (7)

and, if $q_{m}\leq 2^{2m}$, letting

$\lambda_{h}:=2\sin^{-1}\left(q_{m}^{\frac{1}{2m}}/2\right),$ (8)

we can write

$G(\lambda)=\frac{1}{1+\left(\frac{\sin\frac{\lambda}{2}}{\sin\frac{\lambda_{h}}{2}}\right)^{2m}}.$

This is the gain function of the Butterworth filter of order $m$ with cutoff
frequency $\lambda_{h}/2$, see Harvey and Koopman (2000); Harvey and Trimbur
(2003). Estimating the $m$-th order integrated random walk plus noise model by
maximum likelihood and obtaining variance estimates $\hat{\sigma}_{\eta}^{2}$
and $\hat{\sigma}_{\varepsilon}^{2}$ thus implies an estimated cutoff
frequency

$\frac{\hat{\lambda}_{h}}{2}=\sin^{-1}\left(\hat{q}_{m}^{\frac{1}{2m}}/2\right).$

The higher the order $m$ of the filter, the smoother the filtered signal
$\hat{\mu}_{t}$, the smaller the estimated variance $\hat{\sigma}_{\eta}^{2}$
relative to $\hat{\sigma}_{\varepsilon}^{2}$ and the lower the estimated
signal-to-noise ratio $\hat{q}_{m}$. In this sense, $q_{m}$ depends on $m$.

##  4 Integrated random walk plus noise models

###  4.1 Specifications

The insights of the last section into the connection of integrated random walk
models ($m>1$) with the Butterworth filter suggest considering models
$y_{t}=\mu_{t}+\varepsilon_{t}$, $\varepsilon\stackrel{{\scriptstyle
i.i.d.}}{{\sim}}\mathsf{N}(0,\sigma_{\varepsilon}^{2})$, where, for example,
$m=2$ and

$\mu_{t_{\nu}}=\sum_{i_{1}=1}^{\nu-2}\sum_{i_{2}=1}^{i_{1}}\eta_{t_{i_{2}}},\quad\eta_{t_{\nu}}\stackrel{{\scriptstyle
i.i.d.}}{{\sim}}\mathsf{N}(0,\sigma_{\eta,2}^{2}\Delta t_{\nu}).$

These models smooth the data to a higher degree than the random walk plus
noise model ($m=1$), since the filter weighs lower frequencies more with
increasing $m$. Let

$\mu_{t_{\nu}}^{(1)}=\sum_{i_{1}=1}^{\nu-1}\eta_{t_{i_{1}}}=\sum_{i_{1}=1}^{\nu-2}\eta_{t_{i_{1}}}+\eta_{t_{\nu-1}}=\mu_{t_{\nu-1}}^{(1)}+\eta_{t_{\nu-1}},$

then

$\mu_{t_{\nu}}=\sum_{i_{1}=1}^{\nu-3}\sum_{i_{2}=1}^{i_{1}}\eta_{t_{i_{2}}}+\sum_{i_{1}=1}^{\nu-2}\eta_{t_{i_{1}}}=\mu_{t_{\nu-1}}+\mu_{t_{\nu-1}}^{(1)}.$

Then, the transition equation of a corresponding state space model is
specified as

$\left[\begin{array}[]{c}\mu_{t+\Delta t}\\\ \mu_{t+\Delta
t}^{(1)}\end{array}\right]=\left[\begin{array}[]{cc}1&1\\\
0&1\end{array}\right]\left[\begin{array}[]{c}\mu_{t}\\\
\mu_{t}^{(1)}\end{array}\right]+\left[\begin{array}[]{c}0\\\
\eta_{t}\end{array}\right],\;\eta_{t}\stackrel{{\scriptstyle
i.i.d.}}{{\sim}}\mathsf{N}(0,\sigma_{\eta,2}^{2}\Delta t).$ (9)

The measurement equation is as before

$y_{t}=\mu_{t}+\varepsilon_{t},\quad\varepsilon\stackrel{{\scriptstyle
i.i.d.}}{{\sim}}\mathsf{N}(0,\sigma_{\varepsilon}^{2}),$

and the distinctions according to studies in the measurement equation and
according to periods in the transition equation can be introduced as outlined
in Section 2. Appendix A spells out the system matrices.

For $m=3$, let

$\mu_{t_{\nu}}^{(1)}=\sum_{i_{1}=1}^{\nu-1}\eta_{t_{i_{1}}}$

as above, and

$\mu_{t_{\nu}}^{(2)}=\sum_{i_{1}=1}^{\nu-2}\sum_{i_{2}=1}^{i_{1}}\eta_{t_{i_{2}}}=\sum_{i_{1}=1}^{\nu-3}\sum_{i_{2}=1}^{i_{1}}\eta_{t_{i_{2}}}+\sum_{i_{1}=1}^{\nu-2}\eta_{t_{i_{1}}}=\mu_{t_{\nu-1}}^{(2)}+\mu_{t_{\nu-1}}^{(1)}.$

Then, let

$\mu_{t_{\nu}}=\sum_{i_{1}=1}^{\nu-3}\sum_{i_{2}=1}^{i_{1}}\sum_{i_{3}=1}^{i_{2}}\eta_{t_{i_{3}}}=\sum_{i_{1}=1}^{\nu-4}\sum_{i_{2}=1}^{i_{1}}\sum_{i_{3}=1}^{i_{2}}\eta_{t_{i_{3}}}+\sum_{i_{1}=1}^{\nu-3}\sum_{i_{2}=1}^{i_{1}}\eta_{t_{i_{2}}}=\mu_{t_{\nu-1}}+\mu_{t_{\nu-1}}^{(2)}.$

The transition equation of the corresponding state space model is specified as

$\left[\begin{array}[]{c}\mu_{t+\Delta t}\\\ \mu_{t+\Delta t}^{(2)}\\\
\mu_{t+\Delta t}^{(1)}\end{array}\right]=\left[\begin{array}[]{ccc}1&1&0\\\
0&1&1\\\ 0&0&1\end{array}\right]\left[\begin{array}[]{c}\mu_{t}\\\
\mu_{t}^{(2)}\\\ \mu_{t}^{(1)}\end{array}\right]+\left[\begin{array}[]{c}0\\\
0\\\ \eta_{t}\end{array}\right],\;\eta_{t}\stackrel{{\scriptstyle
i.i.d.}}{{\sim}}\mathsf{N}(0,\sigma_{\eta,3}^{2}\Delta t).$ (10)

For general $m$, let

$\mu_{t_{\nu}}=\sum_{i_{1}=1}^{\nu-m}\sum_{i_{2}=1}^{i_{1}}\cdots\sum_{i_{m}=1}^{i_{m-1}}\eta_{t_{i_{m}}},$

and for $j=1,\ldots,m-1$, define

$\mu_{t_{\nu}}^{(j)}=\sum_{i_{1}=1}^{\nu-j}\sum_{i_{2}=1}^{i_{1}}\cdots\sum_{i_{j}=1}^{i_{j-1}}\eta_{t_{i_{j}}}=\mu_{t_{\nu-1}}^{(j)}+\mu_{t_{\nu-1}}^{(j-1)}$

and thus

$\mu_{t_{\nu}}=\mu_{t_{\nu-1}}+\mu_{t_{\nu-1}}^{(m-1)}.$

The transition equation of the state space model is

$\left[\begin{array}[]{c}\mu_{t+\Delta t}\\\ \mu_{t+\Delta t}^{(m-1)}\\\
\vdots\\\ \mu_{t+\Delta t}^{(2)}\\\ \mu_{t+\Delta
t}^{(1)}\end{array}\right]=\left[\begin{array}[]{ccccccc}1&1&0&0&\cdots&0&0\\\
0&1&1&0&\cdots&0&0\\\ \vdots&\vdots&\vdots&\vdots&\ddots&\vdots&\vdots\\\
0&0&0&0&\cdots&1&1\\\
0&0&0&0&\cdots&0&1\end{array}\right]\left[\begin{array}[]{c}\mu_{t}\\\
\mu_{t}^{(m-1)}\\\ \vdots\\\ \mu_{t}^{(2)}\\\
\mu_{t}^{(1)}\end{array}\right]+\left[\begin{array}[]{c}0\\\ 0\\\ \vdots\\\
0\\\ \eta_{t}\end{array}\right],\;\eta_{t}\stackrel{{\scriptstyle
i.i.d.}}{{\sim}}\mathsf{N}(0,\sigma_{\eta,m}^{2}\Delta t).$ (11)

###  4.2 Applications

We pursue this class of models in the application to $\delta^{18}$O and
$\delta^{13}$C in two ways: First, in order to study the gain functions and
the extent of smoothing that the different orders of integration imply, we
estimate a single measurement equation variance and a single transition
equation variance for $m=1,\ldots,6$. Second, we differentiate measurement
equation variances according to studies and transition equation variances
according to climate states.

In the constant-variances case, the signal-to-noise ratios are computed as

$\hat{q}_{m}=\frac{\hat{\sigma}_{\eta,m}^{2}\overline{\Delta
t}}{\hat{\sigma}_{\varepsilon}^{2}},$ (12)

where $\overline{\Delta t}$ is the average time increment per observation from
the entire sample. Loosely speaking, the transition equation variance
$\sigma_{\eta,m}^{2}$ has continuous-time dimension, whereas the measurement
equation variance $\sigma_{\varepsilon}^{2}$ has per-observation dimension.
Multiplying the transition equation variance by the average time increment per
observation makes them dimension-congruent.

Table 8 reports the signal-to-noise ratios computed from Equation (12) and the
estimated $\sigma_{\eta,m}^{2}$ and $\sigma_{\varepsilon}^{2}$. (The estimated
$\sigma_{\varepsilon}^{2}$ also vary with $m$, since greater smoothing in the
estimated state leaves more data variance to be captured by the measurement
equation, but this is not very important for our analysis, and so we do not
index $\sigma_{\varepsilon}^{2}$ with $m$ in order to keep notation as simple
as possible.) The corresponding cutoff frequencies $\lambda_{h}$ are computed
from Equation (8). Note that the BICs cannot be compared across $m$, since $m$
reflects a design choice of the degree of smoothing, and stronger smoothing
implies a deterioration of fit.

Table 8: Estimated signal-to-noise ratios and cutoff frequencies from model (11). BIC: Bayes information criterion. |  $\delta^{18}$O |  $\delta^{13}$C
---|---|---
| $\hat{q}_{m}$ | $\hat{\lambda}_{h}$ | log-likelihood | BIC | $\hat{q}_{m}$ | $\hat{\lambda}_{h}$ | log-likelihood | BIC
$m=1$ | 0.2537 | 0.5092 | 7218.57 | -14416.99 | 0.1009 | 0.3191 | 2992.10 | -5964.05
$m=2$ | 0.0095 | 0.0974 | 6237.73 | -12455.31 | 0.0020 | 0.0445 | 2204.29 | -4388.43
$m=3$ | 0.0001 | 0.0122 | 4812.72 | -9605.40 | 9.1488e-6 | 0.0030 | 1167.39 | -2314.64
$m=4$ | 2.8447e-7 | 5.3336e-4 | 3545.53 | -7070.91 | 2.3601e-8 | 1.5363e-4 | 352.74 | -685.33
$m=5$ | 9.6940e-10 | 3.1135e-5 | 2641.94 | -5263.72 | 9.4511e-12 | 3.0743e-6 | -332.63 | 685.42
$m=6$ | 3.0777e-12 | 1.7543e-6 | 1820.51 | -3620.88 | 3.5973e-15 | 5.9978e-8 | -836.30 | 1692.74

The order of integration $m$ in model (11), equal to the order of the
corresponding Butterworth filter, is a design choice by the researcher that
reflects the degree of desired smoothing. Given this design choice, the
estimation maximizes the log-likelihood of model (11) with respect to
$\sigma_{\eta,m}^{2}$ and $\sigma_{\varepsilon}^{2}$ by way of the Kalman
filter recursions (Durbin and Koopman, 2012). Therefore, conditional on the
data and $m$, the signal-to-noise ratios and the implied gain functions from
Equation (7) are optimal in this sense. Table 8 shows that the estimated
optimal signal-to-noise ratios decline steeply in $m$, which is expected,
since a smoother state process $\mu$ captures less variation in the data, has
a smaller variance $\sigma_{\eta,m}^{2}$ of the underlying white-noise driver
and more variation in the data is captured by $\sigma_{\varepsilon}^{2}$. The
maximized log-likelihood decreases sharply in $m$. This also reflects the
higher degree of smoothing of the estimated latent process $\mu$ with
increasing $m$, which makes it increasingly unlikely that the residuals are
draws from a normal distribution.

Figure 4: Left panel: Estimated gain functions for $\delta^{18}$O for $m=1,\ldots,6$. Right panel: Same for $\delta^{13}$C. Table 9: Maximum likelihood estimates of model (11) for $\delta^{18}$O. Standard errors in parentheses. “Source” refers to the entry in column “benthic Source” in tab “Table S33” of the data file aba6853_tables_s8_s34.xlsx of Westerhold et al. (2020). The identifier “McCarren et al. 2008 et al. 2008” was merged with “McCarren et al. 2008” and “Bickert et al.1997” with “Bickert et al. 1997”. The study Westerhold et al. 2020 is coded in the data file as “this study”. BIC: Bayes information criterion.  Source | $m=2$ | $m=4$ | $m=6$ | $m=8$
---|---|---|---|---
$\sigma_{\varepsilon,i}^{2}$ |  $\delta^{18}$O
“Barnet et al. 2017” | 0.0046 (0.0003) | 0.0053 (0.0003) | 0.0056 (0.0004) | 0.0057 (0.0004)
“Hull et al. 2020” | 0.0146 (0.0039) | 0.0158 (0.0042) | 0.0159 (0.0043) | 0.0161 (0.0043)
“Barnet et al. 2019” | 0.0109 (0.0006) | 0.0134 (0.0007) | 0.0144 (0.0007) | 0.0152 (0.0008)
“Littler et al. 2014” | 0.0140 (0.0005) | 0.0158 (0.0005) | 0.0183 (0.0007) | 0.0204 (0.0007)
“McCarren et al. 2008” | 0.3357 (0.0727) | 0.3843 (0.0811) | 0.4849 (0.1011) | 0.5193 (0.1079)
“Lauretano et al. 2015” | 0.0299 (0.0022) | 0.0314 (0.0022) | 0.0353 (0.0024) | 0.0367 (0.0024)
“Stap et al. 2010” | 0.0312 (0.0024) | 0.0328 (0.0025) | 0.0338 (0.0025) | 0.0341 (0.0026)
“Thomas et al. 2018” | 0.0189 (0.0031) | 0.0229 (0.0037) | 0.0304 (0.0048) | 0.0293 (0.0045)
“Lauretano et al. 2016; 2018” | 0.0161 (0.0011) | 0.0207 (0.0013) | 0.0246 (0.0015) | 0.0258 (0.0016)
“Sexton et al. 2011” | 0.0230 (0.0020) | 0.0244 (0.0020) | 0.0270 (0.0022) | 0.0284 (0.0024)
Westerhold et al. 2020 | 0.0141 (0.0005) | 0.0153 (0.0005) | 0.0159 (0.0006) | 0.0162 (0.0006)
“Westerhold et al. 2018 LLTM” | 0.0267 (0.0073) | 0.0269 (0.0073) | 0.0266 (0.0072) | 0.0268 (0.0072)
“Boscolo Galazzo et al. 2014” | 0.0286 (0.0072) | 0.0337 (0.0084) | 0.0350 (0.0087) | 0.0362 (0.0090)
“Lear et al. 2004” | 0.0163 (0.0019) | 0.0236 (0.0026) | 0.0247 (0.0027) | 0.0255 (0.0027)
“Coxall et al. 2005” | 0.0137 (0.0008) | 0.0183 (0.0009) | 0.0188 (0.0009) | 0.0192 (0.0010)
“Coxall & Wilson, 2011” | 0.0216 (0.0017) | 0.0250 (0.0018) | 0.0256 (0.0018) | 0.0261 (0.0019)
“Riesselmann et al. 2007” | 0.0192 (0.0042) | 0.0205 (0.0046) | 0.0220 (0.0051) | 0.0227 (0.0055)
“Paelike2006 PI Herrle” | 0.0248 (0.0024) | 0.0305 (0.0027) | 0.0309 (0.0028) | 0.0315 (0.0028)
“Wade and Paelike 2004” | 0.0496 (0.0060) | 0.0605 (0.0070) | 0.0610 (0.0071) | 0.0614 (0.0071)
“Liebrand et al. 2011;2016” | 0.0147 (0.0004) | 0.0229 (0.0006) | 0.0234 (0.0006) | 0.0236 (0.0006)
“Holbourn et al. 2015” | 0.0324 (0.0017) | 0.0502 (0.0025) | 0.0522 (0.0026) | 0.0539 (0.0027)
“Tian et al. 2014” | 0.0573 (0.0110) | 0.0988 (0.0181) | 0.1009 (0.0184) | 0.1014 (0.0185)
“Holbourn et al. 2014” | 0.0208 (0.0010) | 0.0402 (0.0016) | 0.0402 (0.0016) | 0.0414 (0.0016)
“Holbourn et al. 2007;2013;2018” | 0.0103 (0.0004) | 0.0205 (0.0007) | 0.0208 (0.0007) | 0.0210 (0.0007)
“Drury et al. 2017” | 0.0101 (0.0005) | 0.0162 (0.0007) | 0.0163 (0.0007) | 0.0165 (0.0007)
“Tian et al. 2018” | 0.0127 (0.0009) | 0.0169 (0.0010) | 0.0169 (0.0010) | 0.0173 (0.0011)
UniCambridge Hodell / Westerhold et al. 2020 | 0.0176 (0.0012) | 0.0260 (0.0016) | 0.0261 (0.0016) | 0.0263 (0.0016)
“Bell et al. 2014” | 0.0081 (0.0008) | 0.0147 (0.0013) | 0.0149 (0.0013) | 0.0154 (0.0014)
MARUM Bremen / Westerhold et al. 2020 | 0.0093 (0.0020) | 0.0155 (0.0033) | 0.0155 (0.0033) | 0.0158 (0.0033)
“Franz & Tiedemann2002” | 0.0643 (0.0027) | 0.0690 (0.0028) | 0.0702 (0.0028) | 0.0727 (0.0029)
“Tiedemann & Franz 1997” | 0.0310 (0.0029) | 0.0374 (0.0032) | 0.0377 (0.0032) | 0.0386 (0.0032)
“Billups et al. 1998” | 0.0346 (0.0023) | 0.0408 (0.0025) | 0.0411 (0.0025) | 0.0422 (0.0026)
“Bickert et al. 1997” | 0.0462 (0.0015) | 0.0543 (0.0017) | 0.0597 (0.0019) | 0.0645 (0.0024)
“deMenocal et al. 1997” | 0.0233 (0.0031) | 0.0273 (0.0034) | 0.0288 (0.0035) | 0.0298 (0.0036)
log-likelihood | 9265.22 | 7234.54 | 6072.71 | 5033.15
BIC | -18127.47 | -14066.12 | -11742.45 | -9663.34
Table 10: Maximum likelihood estimates of model (11) for $\delta^{13}$C. Standard errors in parentheses. “Source” refers to the entry in column “benthic Source” in tab “Table S33” of the data file aba6853_tables_s8_s34.xlsx of Westerhold et al. (2020). The identifier “McCarren et al. 2008 et al. 2008” was merged with “McCarren et al. 2008” and “Bickert et al.1997” with “Bickert et al. 1997”. The study Westerhold et al. 2020 is coded in the data file as “this study”. BIC: Bayes information criterion. Source | $m=2$ | $m=4$ | $m=6$ | $m=8$
---|---|---|---|---
$\sigma_{\varepsilon,i}^{2}$ |  $\delta^{13}$C
“Barnet et al. 2017” | 0.0073 (0.0005) | 0.0083 (0.0005) | 0.0097 (0.0006) | 0.0102 (0.0006)
“Hull et al. 2020” | 0.0117 (0.0033) | 0.0113 (0.0031) | 0.0120 (0.0032) | 0.0127 (0.0034)
“Barnet et al. 2019” | 0.0056 (0.0003) | 0.0079 (0.0005) | 0.0117 (0.0006) | 0.0138 (0.0007)
“Littler et al. 2014” | 0.0092 (0.0003) | 0.0117 (0.0004) | 0.0163 (0.0005) | 0.0176 (0.0006)
“McCarren et al. 2008” | 0.6930 (0.1494) | 0.8323 (0.1740) | 0.9165 (0.1908) | 0.9725 (0.2020)
“Lauretano et al. 2015” | 0.0203 (0.0014) | 0.0225 (0.0014) | 0.0245 (0.0016) | 0.0265 (0.0016)
“Stap et al. 2010” | 0.0244 (0.0020) | 0.0255 (0.0021) | 0.0268 (0.0021) | 0.0278 (0.0022)
“Thomas et al. 2018” | 0.0200 (0.0032) | 0.0217 (0.0035) | 0.0248 (0.0040) | 0.0278 (0.0046)
“Lauretano et al. 2016; 2018” | 0.0188 (0.0014) | 0.0320 (0.0021) | 0.0366 (0.0023) | 0.0378 (0.0024)
“Sexton et al. 2011” | 0.0249 (0.0023) | 0.0316 (0.0027) | 0.0325 (0.0027) | 0.0331 (0.0028)
Westerhold et al. 2020 | 0.0109 (0.0004) | 0.0130 (0.0005) | 0.0147 (0.0005) | 0.0174 (0.0006)
“Westerhold et al. 2018 LLTM” | 0.0962 (0.0259) | 0.0932 (0.0245) | 0.0920 (0.0238) | 0.1108 (0.0293)
“Boscolo Galazzo et al. 2014” | 0.0158 (0.0041) | 0.0173 (0.0044) | 0.0203 (0.0051) | 0.0276 (0.0068)
“Lear et al. 2004” | 0.0135 (0.0016) | 0.0179 (0.0020) | 0.0184 (0.0020) | 0.0184 (0.0020)
“Coxall et al. 2005” | 0.0131 (0.0007) | 0.0167 (0.0008) | 0.0171 (0.0009) | 0.0172 (0.0009)
“Coxall & Wilson, 2011” | 0.0140 (0.0011) | 0.0158 (0.0011) | 0.0162 (0.0012) | 0.0163 (0.0012)
“Riesselmann et al. 2007” | 0.0380 (0.0088) | 0.0435 (0.0111) | 0.0994 (0.0330) | 0.1777 (0.0441)
“Paelike2006 PI Herrle” | 0.0177 (0.0017) | 0.0222 (0.0020) | 0.0223 (0.0020) | 0.0222 (0.0020)
“Wade and Paelike 2004” | 0.0360 (0.0047) | 0.0396 (0.0050) | 0.0414 (0.0053) | 0.0431 (0.0055)
“Liebrand et al. 2011;2016” | 0.0138 (0.0004) | 0.0178 (0.0005) | 0.0183 (0.0005) | 0.0185 (0.0005)
“Holbourn et al. 2015” | 0.0184 (0.0010) | 0.0263 (0.0013) | 0.0279 (0.0014) | 0.0287 (0.0014)
“Tian et al. 2014” | 0.0358 (0.0069) | 0.0426 (0.0079) | 0.0429 (0.0079) | 0.0434 (0.0080)
“Holbourn et al. 2014” | 0.0142 (0.0006) | 0.0184 (0.0007) | 0.0196 (0.0008) | 0.0200 (0.0008)
“Holbourn et al. 2007;2013;2018” | 0.0112 (0.0004) | 0.0159 (0.0006) | 0.0200 (0.0007) | 0.0231 (0.0009)
“Drury et al. 2017” | 0.0201 (0.0010) | 0.0299 (0.0013) | 0.0326 (0.0014) | 0.0332 (0.0014)
“Tian et al. 2018” | 0.0290 (0.0019) | 0.0373 (0.0023) | 0.0413 (0.0025) | 0.0425 (0.0026)
UniCambridge Hodell / Westerhold et al. 2020 | 0.0124 (0.0008) | 0.0157 (0.0010) | 0.0197 (0.0012) | 0.0206 (0.0013)
“Bell et al. 2014” | 0.0180 (0.0017) | 0.0256 (0.0023) | 0.0270 (0.0024) | 0.0272 (0.0024)
MARUM Bremen / Westerhold et al. 2020 | 0.0134 (0.0030) | 0.0160 (0.0035) | 0.0172 (0.0037) | 0.0172 (0.0037)
“Franz & Tiedemann2002” | 0.1167 (0.0055) | 0.1209 (0.0056) | 0.1391 (0.0063) | 0.1396 (0.0063)
“Tiedemann & Franz 1997” | 0.0608 (0.0064) | 0.0631 (0.0064) | 0.0744 (0.0071) | 0.0748 (0.0072)
“Billups et al. 1998” | 0.0641 (0.0042) | 0.0698 (0.0044) | 0.0711 (0.0044) | 0.0712 (0.0044)
“Bickert et al. 1997” | 0.1406 (0.0043) | 0.1625 (0.0049) | 0.2077 (0.0059) | 0.2073 (0.0058)
“deMenocal et al. 1997” | 0.1175 (0.0145) | 0.1364 (0.0163) | 0.2215 (0.0252) | 0.2201 (0.0249)
log-likelihood | 8445.04 | 6212.34 | 4982.08 | 4043.34
BIC | -16487.11 | -12021.72 | -9561.19 | -7683.71
Table 11: Maximum likelihood estimates of model (11). Standard errors in parentheses. LLH: log-likelihood.  $\sigma_{\eta,j}^{2}$ |  $\delta^{18}$O |
---|---|---
Period | 1 | 2 | 3 | 4 | 5 | 6 | LLH
$m=2$ | 0.0034 (0.0007) | 0.0107 (0.0019) | 0.0007 (0.0002) | 0.0425 (0.0049) | 0.0824 (0.0094) | 0.5175 (0.0493) | 9265.22
$m=4$ | 4.6e-8 (1.3e-8) | 3.2e-7 (1e-7) | 5e-9 (1.7e-9) | 3.4e-9 (8.2e-10) | 1.2e-12 (6.5e-13) | 4.7e-4 (7e-5) | 7234.54
$m=6$ | 4.1e-13 (1.8e-13) | 4.7e-14 (2.5e-14) | 7.4e-14 (2.7e-14) | 8.9e-15 (2.3e-15) | 2.9e-20 (1.5e-20) | 2.8e-7 (5.1e-8) | 6072.71
$m=8$ | 2.2e-18 (8.1e-19) | 2.2e-21 (8.3e-22) | 7.7e-19 (3.8e-19) | 2.9e-20 (1.0e-20) | 9.7e-30 (4.7e-30) | 1.3e-10 (4.1e-11) | 5033.15
|  $\delta^{13}$C |
Period | 1 | 2 | 3 | 4 | 5 | 6 | LLH
$m=2$ | 0.03024 (0.0036) | 0.0726 (0.0103) | 0.0041 (0.0011) | 0.0108 (0.0014) | 0.0423 (0.0059) | 0.1158 (0.0197) | 8445.04
$m=4$ | 1.8e-5 (4.6e-6) | 4.8e-6 (1.2e-6) | 2.0e-7 (6.2e-8) | 7.2e-9 (1.5e-9) | 8.1e-8 (2.4e-8) | 1.6e-6 (4.6e-7) | 6212.34
$m=6$ | 7.1e-12 (2.9e-12) | 1.3e-10 (4.7e-11) | 6.2e-13 (2.9e-13) | 3.0e-14 (6.0e-15) | 1.8e-15 (6.1e-16) | 6.2e-20 (5.7e-20) | 4982.08
$m=8$ | 1.9e-17 (5.4e-18) | 3.7e-15 (1.5e-15) | 2.3e-20 (1.0e-20) | 1.5e-19 (3.2e-20) | 5.6e-23 (4.4e-23) | 4.1e-27 (2.5e-27) | 4043.34

Figure 4 plots the estimated optimal gain functions for $\delta^{18}$O (left
panel) and $\delta^{13}$C (right panel) for different orders $m$. The figure
and Table 8 show that orders $m$ beyond 3 result in a cutoff very close to the
zero frequency. Nevertheless, the resulting estimated latent processes $\mu$
still show a good fit and a discernible difference in the degree of smoothing
for orders $m$ between 4 and 6, as we now demonstrate in the estimation that
distinguishes measurement equation variances according to studies and
transition equation variances according to climate states.

To differentiate measurement equation variances according to studies, we
employ the measurement equation specification (2). The transition equation
differentiates variances according to climate states as defined in Table 4. We
modify Equation (11) such that the last element in the state vector becomes

$\mu_{t+\Delta
t}^{(1)}=\mu_{t}^{(1)}+\eta_{t},\;\eta_{t}\stackrel{{\scriptstyle
i.i.d.}}{{\sim}}\mathsf{N}(0,\sigma_{\eta,m,j}^{2}),\;j=1,\ldots,6,\;m\in\\{2,4,6,8\\}.$
(13)

Tables 9 and 10 report the estimated measurement equation variances for
$\delta^{18}$O and $\delta^{13}$C, respectively. Table 11 reports the
estimated transition equation variances. Figures 5 and 6 show the estimated
latent state processes $\mu$ for $\delta^{18}$O and $\delta^{13}$C,
respectively, in the left column. The right column in both figures shows the
corresponding standardized residuals. Figures 5 and 6 show that there is no
visually discernible difference in the estimated latent processes between
orders $m=6$ and $m=8$. Tables 9 and 10 show that the differences in the
estimated measurement equation variances between $m=6$ and $m=8$ are small and
within the variation indicated by the estimated standard errors. Table 11
shows differences in transition equation variances between $m=6$ and $m=8$ of
several orders of magnitude, but in very high decimal places, $10^{-8}$ being
the biggest order that appears.

Between $m=4$ and $m=6$, in contrast, there is a visually discernible
difference in the estimated latent state processes $\mu$ for $\delta^{18}$O
and $\delta^{13}$C, and the estimated variances, both in measurement and
transition equation, differ substantially. We conclude that integration orders
up to $m=6$ are meaningful for $\delta^{18}$O and $\delta^{13}$C. We recommend
using specifications (13) for the transition equation and (2) for the
measurement equation, as Tables 9, 10, and 11 show that the differences in the
estimated variances are substantial. Comparing log-likelihoods across $m$’s is
not meaningful, since $m$ is a design choice. For given $m$, however,
comparing log-likelihoods in Tables 8 on the one hand and 9 and 10 on the
other, reveals that allowing for study- and climate states-dependent variances
substantially improves the model fit.

Figure 5: Left column: $\delta^{18}$O data (blue) and smoothed state
$\mathds{E}[\mu_{t+\Delta t}|\\{y_{s},\;s\leq\sum\Delta t\\}]$ (red) according
to model (11) for $m=2,4,6,8$. Right column: standardized residuals.

Figure 6: Left column: $\delta^{13}$C data (blue) and smoothed state
$\mathds{E}[\mu_{t+\Delta t}|\\{y_{s},\;s\leq\sum\Delta t\\}]$ (red) according
to model (11) for $m=2,4,6,8$. Right column: standardized residuals.

##  5 Imputing regularly time-stamped values

State space models handle missing observations seamlessly, since the Kalman
recursions compute the predicted values $\mathds{E}[\mu_{t+\Delta
t}|\\{y_{s},\;s\leq t\\}]$, the filtered values $\mathds{E}[\mu_{t+\Delta
t}|\\{y_{s},\;s\leq t+\Delta t\\}]$, and the smoothed values
$\mathds{E}[\mu_{t+\Delta t}|\\{y_{s},\;s\leq\sum\Delta t\\}]$ for all time
stamps, also those that have NAs in the $y_{t}$ vector in the measurement
equation (Durbin and Koopman, 2012). We obtain smoothed estimates
$\mathds{E}[\mu_{t+\Delta t}|\\{y_{s},\;s\leq\sum\Delta t\\}]$ of the
unobserved component $\mu$ from any of the models discussed in this paper at
an equidistantly spaced grid of time stamps in the following manner.

We define a constant time increment $\Delta$ (in years) and generate an
equidistant partition of the interval $[-67,0]$ (MY) at mesh $\Delta$. This
results in $N_{g}=[67000000/\Delta]$ time stamps for which there are no
observations (unless by coincidence), and so we fill the corresponding
elements in $y_{t}$ with “NA.” Here, $[\cdot]$ denotes the Gauss-bracket. This
partition is merged with the irregularly spaced data and ordered according to
time stamps, resulting in $N+N_{g}=23722+[67000000/\Delta]$ time stamps. The
maximum likelihood estimate of the parameter vector of the chosen model
estimated on the $N$ available data points is used as input to the Kalman
recursions, which are then run on the extended data matrix of $N+N_{g}$ rows.
The Kalman recursions compute the predicted, filtered, and smoothed values for
the $N+N_{g}$ observations, and we extract the $N_{g}$ regularly time stamped
smoothed values.

This state space method of obtaining regularly time stamped values can
therefore be seen as an alternative to the nonuniform discrete inverse Fourier
transform (Bagchi and Mitra, 1996). The state space method in the time domain
has the advantage of allowing for a higher resolution of other features of the
data, in our application a differentiation of variances in the measurement
equation according to studies or species and in the transition equation
according to climate states. It also provides estimates of the variances, and
thus signal-to-noise ratios and cutoff frequencies, that are optimal in the
maximum likelihood sense. This leaves only the order of integration (order of
the filter) to be chosen by the researcher.

Based on both the Bayes information criterion and the log-likelihood of the
fitted model, we arrive at a preference for the bivariate random-walk-plus-
noise model with measurement equation variances differentiated according to
studies of origin and transition equation variances differentiated according
to climate states. We therefore use this model to illustrate the imputation of
regularly time stamped values via the Kalman smoother, which may be used by
researchers in further analyses employing standard time series techniques.
Figure 7 shows equispaced time-stamped smoothed values for the unobserved
components for $\delta^{18}$O (top row) and $\delta^{13}$C (bottom row) for
three different mesh values $\Delta=(100000,\,10000,\,1000)$.

Figure 7: Regularly time-stamped smoothed values from Model (5) with
measurement variances according to studies and transition variances according
to climate states. Upper row: $\delta^{18}$O with mesh $\Delta=100000$ years
(left panel), $\Delta=10000$ years (middle panel), $\Delta=1000$ years (right
panel). Bottom row: $\delta^{13}$C with mesh $\Delta=100000$ years (left
panel), $\Delta=10000$ years (middle panel), $\Delta=1000$ years (right
panel).

##  6 Conclusion

We specified state space models for the time series of $\delta^{18}$O and
$\delta^{13}$C provided in Westerhold et al. (2020). The models addressed the
central challenges of the time series: non-stationarity with trends in varying
directions, irregular spacing of the time stamps of the observations, multiple
observations at certain time stamps, different studies of origin of the data,
different climate states covered by the data, and covariation in the two
series. The models were specified and estimated by maximum likelihood in the
time domain, allowing for a direct interpretation as Butterworth filters in
the frequency domain. We showed that accounting for the different studies of
origin of the data by way of differentiating measurement equation variances
substantially improves the model fit. Another improvement is achieved by
differentiating the variances in the state equation by climate states, as they
were identified in Westerhold et al. (2020). Bivariate specifications are
fruitful in that they use the covariation in the $\delta^{18}$O and
$\delta^{13}$C series to inform the estimation of the models. We showed that
integration orders of the unobserved component, which at the same time are the
orders of the corresponding Butterworth filters, up to six are meaningful for
the two time series. They yield visually distinct estimated latent components
with different degrees of smoothing. Therefore, the order is a design choice
made by the researcher reflecting their preference for a degree of smoothing.
Given this choice, and the data, the maximum likelihood estimation yields
optimal estimates of the measurement and state equation variances. These
correspond to optimal signal-to-noise ratios (and thus cutoff frequencies) for
the corresponding Butterworth filters. In summary, well developed statistical
theory, both in the time and in the frequency domains, can be employed to
specify time series models for $\delta^{18}$O and $\delta^{13}$C that can
resolve many of the rich features of these fascinating data. In future
research, we will address uncertainty in the time stamps of the observations
and relate the $\delta^{18}$O and $\delta^{13}$C measurements to paleo-
estimates of atmospheric CO2 concentrations.

##  Appendix A System matrices of the state space models

###  A.1 Univariate random walk plus noise specifications

The measurement equation relates the data points from the data file of
Westerhold et al. (2020) to the unobserved components $\mu_{t}$ that are the
central part of our modeling approach. A fact about the data set that has
direct relevance for the dimension chosen for the measurement equation is that
there are at most four data points at any given time stamp for each of the two
time series, and they can originate from different studies. Denote the data
points by $y$, the studies by $i$, the unobserved components by $\mu$, the
climate states by $j$, and time by $t$.

The set of unique time stamps is given by $t_{\nu}$, $\nu\in\\{1,\ldots,N\\}$,
$N=23722$. The time stamps $t_{\nu}$ fall into the interval
$[-67.10113,\,-0.000564]$, which is read as 67.10113 to 0.000564 million years
ago. At time $t_{\nu}$,

$\underbrace{\left[\begin{array}[]{c}y_{i_{1},t_{\nu}}\\\ y_{i_{2},t_{\nu}}\\\
y_{i_{3},t_{\nu}}\\\
y_{i_{4},t_{\nu}}\end{array}\right]}_{y_{t_{\nu}}\in\mathds{R}^{4\times
1}}=\underbrace{\left[\begin{array}[]{c}1\\\ 1\\\ 1\\\
1\end{array}\right]}_{Z\in\mathds{R}^{4\times
1}}\mu_{t_{\nu}}+\underbrace{\left[\begin{array}[]{c}\varepsilon_{i_{1},t_{\nu}}\\\
\varepsilon_{i_{2},t_{\nu}}\\\ \varepsilon_{i_{3},t_{\nu}}\\\
\varepsilon_{i_{4},t_{\nu}}\end{array}\right]}_{\varepsilon_{t_{\nu}}\in\mathds{R}^{4\times
1}},\;\mu_{t_{\nu}}\in\mathds{R}.$

The index $i_{\\{1,2,3,4\\}}\in\\{1,\,\ldots,\,34\\}$ denotes the 34 different
studies of origin coded in the data file. The zero-mean measurement equation
error vector is assumed to be identically and independently normally
distributed:

$\varepsilon_{t_{\nu}}\stackrel{{\scriptstyle
i.i.d.}}{{\sim}}\mathsf{N}(0,H_{t_{\nu}}),$

with covariance matrix

$H_{t_{\nu}}=\text{diag}\left(\sigma^{2}_{\varepsilon,i_{1}},\sigma^{2}_{\varepsilon,i_{2}},\sigma^{2}_{\varepsilon,i_{3}},\sigma^{2}_{\varepsilon,i_{4}}\right).$

The relation between $t_{\nu}$ and $i_{(\cdot)}$ is the mapping “study
$i_{(\cdot)}$ provided the data point $y_{i_{\cdot},t_{\nu}}$ for time stamp
$t_{\nu}$.”

The transition equation for the unobserved component is given as follows.

$\mu_{t_{\nu}+\Delta t_{\nu}}=\mu_{t_{\nu}}+\sqrt{\Delta
t_{\nu}}\,\eta_{j,t_{\nu}},$

where $j\in\\{1,\ldots,6\\}$ indicates the six different climate states
identified in Westerhold et al. (2020). The zero-mean transition equation
error is assumed to be identically and independently normally distributed:
$\eta_{j,t_{\nu}}\stackrel{{\scriptstyle
i.i.d.}}{{\sim}}\mathsf{N}(0,Q_{t_{\nu}})$, $Q_{t_{\nu}}=\sigma_{\eta,j}^{2}$.
The relation between $t_{\nu}$ and $j$ is the mapping “time $t_{\nu}$ falls
into period $j$ according to Table 4.”

###  A.2 Bivariate random walk plus noise specifications

For the bivariate specifications, let (1) denote $\delta^{18}$O and (2)
$\delta^{13}$C.

The measurement equation is given as follows.

$\underbrace{\left[\begin{array}[]{c}y_{i_{1},t_{\nu}}^{\delta^{18}O}\\\
y_{i_{2},t_{\nu}}^{\delta^{18}O}\\\ y_{i_{3},t_{\nu}}^{\delta^{18}O}\\\
y_{i_{4},t_{\nu}}^{\delta^{18}O}\\\ y_{i_{1},t_{\nu}}^{\delta^{13}C}\\\
y_{i_{2},t_{\nu}}^{\delta^{13}C}\\\ y_{i_{3},t_{\nu}}^{\delta^{13}C}\\\
y_{i_{4},t_{\nu}}^{\delta^{13}C}\end{array}\right]}_{y_{t_{\nu}}\in\mathds{R}^{8\times
1}}=\underbrace{\left[\begin{array}[]{cc}1&0\\\ 1&0\\\ 1&0\\\ 1&0\\\ 0&1\\\
0&1\\\ 0&1\\\ 0&1\end{array}\right]}_{Z\in\mathds{R}^{8\times
2}}\underbrace{\left[\begin{array}[]{c}\mu_{t_{\nu}}^{\delta^{18}O}\\\
\mu_{t_{\nu}}^{\delta^{13}C}\end{array}\right]}_{\mu_{t_{\nu}}\in\mathds{R}^{2\times
1}}+\underbrace{\left[\begin{array}[]{c}\varepsilon_{i_{1},t_{\nu}}^{\delta^{18}O}\\\
\varepsilon_{i_{2},t_{\nu}}^{\delta^{18}O}\\\
\varepsilon_{i_{3},t_{\nu}}^{\delta^{18}O}\\\
\varepsilon_{i_{4},t_{\nu}}^{\delta^{18}O}\\\
\varepsilon_{i_{1},t_{\nu}}^{\delta^{13}C}\\\
\varepsilon_{i_{2},t_{\nu}}^{\delta^{13}C}\\\
\varepsilon_{i_{3},t_{\nu}}^{\delta^{13}C}\\\
\varepsilon_{i_{4},t_{\nu}}^{\delta^{13}C}\end{array}\right]}_{\varepsilon_{t_{\nu}}\in\mathds{R}^{8\times
1}}.$

The distribution assumption is as above for the univariate model, simply
extended to the longer error vector:

$H_{t_{\nu}}=\text{diag}\left(\sigma^{2}_{\varepsilon,i_{1},\delta^{18}O},\sigma^{2}_{\varepsilon,i_{2},\delta^{18}O},\sigma^{2}_{\varepsilon,i_{3},\delta^{18}O},\sigma^{2}_{\varepsilon,i_{4},\delta^{18}O},\sigma^{2}_{\varepsilon,i_{1},\delta^{13}C},\sigma^{2}_{\varepsilon,i_{2},\delta^{13}C},\sigma^{2}_{\varepsilon,i_{3},\delta^{13}C},\sigma^{2}_{\varepsilon,i_{4},\delta^{13}C}\right),$

and the relation between $t_{\nu}$ and $i_{(\cdot)}$ is as before the mapping
“study $i_{(\cdot)}$ provided the data point $y_{i_{\cdot},t_{\nu}}$ for time
stamp $t_{\nu}$.”

The transition equation for the unobserved component $\mu$ is given as
follows.

$\underbrace{\left[\begin{array}[]{c}\mu_{t_{\nu}+\Delta
t_{\nu}}^{\delta^{18}O}\\\ \mu_{t_{\nu}+\Delta
t_{\nu}}^{\delta^{13}C}\end{array}\right]}_{\mu_{t_{\nu}+\Delta
t_{\nu}}\in\mathds{R}^{2\times 1}}=\underbrace{\left[\begin{array}[]{cc}1&0\\\
0&1\end{array}\right]}_{T\in\mathds{R}^{2\times
2}}\underbrace{\left[\begin{array}[]{c}\mu_{t_{\nu}}^{\delta^{18}O}\\\
\mu_{t_{\nu}}^{\delta^{13}C}\end{array}\right]}_{\mu_{t_{\nu}}\in\mathds{R}^{2\times
1}}+\underbrace{\left[\begin{array}[]{c}\eta_{j,t_{\nu}}^{\delta^{18}O}\\\
\eta_{j,t_{\nu}}^{\delta^{13}C}\end{array}\right]}_{\eta_{t_{\nu}}\in\mathds{R}^{2\times
1}}.$

In the transition equation, we allow for correlation of the transition error
from the unobserved component in $\delta^{18}$O with the one in
$\delta^{13}$C: $\eta_{t_{\nu}}\stackrel{{\scriptstyle
i.i.d.}}{{\sim}}\mathsf{N}(0,Q_{t_{\nu}})$,

$Q_{t_{\nu}}=\left[\begin{array}[]{cc}\sigma_{\eta,\delta^{18}O,j}^{2}\Delta
t_{\nu}^{\delta^{18}O}&\rho_{j}\sigma_{\eta,\delta^{18}O,j}\sigma_{\eta,\delta^{13}C,j}\min\\{\Delta
t_{\nu}^{\delta^{18}O},\Delta t_{\nu}^{\delta^{13}C}\\}\\\
\rho_{j}\sigma_{\eta,\delta^{18}O,j}\sigma_{\eta,\delta^{13}C,j}\min\\{\Delta
t_{\nu}^{\delta^{18}O},\Delta
t_{\nu}^{\delta^{13}C}\\}&\sigma_{\eta,\delta^{13}C,j}^{2}\Delta
t_{\nu}^{\delta^{13}C}\end{array}\right],$

where the $j$ indicate the climate state according to Table 4, and

$\Delta t_{\nu}^{(\delta^{18}O,\delta^{13}C)}=\left\\{\begin{array}[]{l}\Delta
t_{\nu},\qquad\qquad\qquad\rule{17.07182pt}{0.0pt}\text{if preceding obs at
$t_{\nu-1}$ is not NA},\\\ \Delta t_{\nu}+\Delta
t_{\nu-1}^{(\delta^{18}O,\delta^{13}C)},\quad\text{if preceding obs at
$t_{\nu-1}$ is NA}.\end{array}\right.$

##  Appendix B Measurement equation variances by benthic foraminifera species

Figure 8: Top left panel: $\delta^{18}$O data (blue) and smoothed state $\mathds{E}[\mu_{t+\Delta t}|\\{y_{s},\;s\leq\sum\Delta t\\}]$ (red) according to model (4) from Kalman smoothing recursions. Top right panel: same for $\delta^{13}$C. Note that the $y$-axis for $\delta^{18}$O is reversed, following common practice. Bottom left panel: standardized residuals for $\delta^{18}$O, bottom right panel: standardized residuals for $\delta^{13}$C. Table 12: Maximum likelihood estimates of model (4). Standard errors in parentheses. “Species” refers to the entry in column “benthic Species” in tab “Table S33” of the data file aba6853_tables_s8_s34.xlsx of Westerhold et al. (2020). BIC: Bayes information criterion.  Variance | Species | No. obs. |  $\delta^{18}$O | No. obs. |  $\delta^{13}$C
---|---|---|---|---|---
$\sigma_{\varepsilon,1}^{2}$ | “NTRUE” | 6949 | 0.0120 (0.0003) | 6949 | 0.0093 (0.0002)
$\sigma_{\varepsilon,2}^{2}$ | “OUMB” | 240 | 0.0510 (0.0052) | 56 | 0.0524 (0.0113)
$\sigma_{\varepsilon,3}^{2}$ | “CSUBS or CEOCEA, not specified in publication” | 301 | 0.0125 (0.0016) | 301 | 0.0162 (0.0019)
$\sigma_{\varepsilon,4}^{2}$ | “CPRAE” | 136 | 0.0165 (0.0030) | 136 | 0.0203 (0.0034)
$\sigma_{\varepsilon,5}^{2}$ | “CSPP” | 2310 | 0.0188 (0.0007) | 2310 | 0.0332 (0.0011)
$\sigma_{\varepsilon,6}^{2}$ |  CSPP $>$250 |  | 0.0177 (0.0017) |  | 0.0073 (0.0008)
|  “CSPP, $>$250” | 353 |  | 353 |
|  “CSPP, specimen $>$250 $\mu$m” | 10 |  | 10 |
$\sigma_{\varepsilon,7}^{2}$ | CSPP other |  | 0.0209 (0.0042) |  | 0.0231 (0.0045)
| “CSPP, whole specimen” | 25 |  | 25 |
| “CSPP, 150-250” | 2 |  | 2 |
|  “CSPP, $>$250, Reruns” | 58 |  | 58 |
$\sigma_{\varepsilon,8}^{2}$ | “CGRIM” | 318 | 0.0294 (0.0028) | 318 | 0.0203 (0.0019)
$\sigma_{\varepsilon,9}^{2}$ | “CHAVA” | 92 | 0.0246 (0.0046) | 92 | 0.0111 (0.0024)
$\sigma_{\varepsilon,10}^{2}$ | “CMUND” | 4380 | 0.0105 (0.0003) | 4380 | 0.0126 (0.0003)
$\sigma_{\varepsilon,11}^{2}$ | “PWUEL; CMUND” | 4249 | 0.0114 (0.0004) | 4249 | 0.0084 (0.0003)
$\sigma_{\varepsilon,12}^{2}$ | “PWUEL (eitherCKULL or PWUEL)” | 807 | 0.0125 (0.0007) | 807 | 0.0129 (0.0007)
$\sigma_{\varepsilon,13}^{2}$ | “PWUEL” | 3405 | 0.0592 (0.0017) | 3407 | 0.1458 (0.0039)
$\sigma_{\varepsilon,14}^{2}$ | “CKULL” | 255 | 0.0551 (0.0054) | 254 | 0.0678 (0.0074)
$\sigma_{\varepsilon,15}^{2}$ | “UVIG” | 149 | 0.0498 (0.0064) | 133 | 0.0591 (0.0010)
$\sigma_{\varepsilon,16}^{2}$ | “GORB” | 30 | 0.0604 (0.0167) | 0 | —
$\sigma_{\varepsilon,17}^{2}$ | “PMURR” | 91 | 0.0719 (0.0115) | 0 | —
$\sigma_{\varepsilon,18}^{2}$ | “NUMB” | 65 | 0.0741 (0.0139) | 65 | 0.4352 (0.0789)
$\sigma_{\varepsilon,19}^{2}$ | “CBRA” | 22 | 0.0408 (0.0139) | 22 | 0.0733 (0.0242)
$\sigma_{\varepsilon,20}^{2}$ | “CCIC” | 12 | 0.0657 (0.0284) | 12 | 0.1171 (0.0499)
$\sigma_{\eta}^{2}$ |  |  | 1.6051 (0.0485) |  | 1.1171 (0.0499)
log-likelihood |  | 24259 | 9169.53 | 23939 | 8855.73
BIC |  |  | -18127.50 |  | -17520.06

##  Appendix C Unobserved components and residuals plots for model (5)

Figure 9: Top left panel: $\delta^{18}$O data (blue) and smoothed state
$\mathds{E}[\mu_{t+\Delta t}|\\{y_{s},\;s\leq\sum\Delta t\\}]$ (red) according
to model (5) for the first climate state from Kalman smoothing recursions. Top
right panel: same for $\delta^{13}$C. Note that the $y$-axis for
$\delta^{18}$O is reversed, following common practice. Bottom left panel:
standardized residuals for $\delta^{18}$O, bottom right panel: standardized
residuals for $\delta^{13}$C.

Figure 10: Top left panel: $\delta^{18}$O data (blue) and smoothed state
$\mathds{E}[\mu_{t+\Delta t}|\\{y_{s},\;s\leq\sum\Delta t\\}]$ (red) according
to model (5) for the second climate state from Kalman smoothing recursions.
Top right panel: same for $\delta^{13}$C. Note that the $y$-axis for
$\delta^{18}$O is reversed, following common practice. Bottom left panel:
standardized residuals for $\delta^{18}$O, bottom right panel: standardized
residuals for $\delta^{13}$C.

Figure 11: Top left panel: $\delta^{18}$O data (blue) and smoothed state
$\mathds{E}[\mu_{t+\Delta t}|\\{y_{s},\;s\leq\sum\Delta t\\}]$ (red) according
to model (5) for the third climate state from Kalman smoothing recursions. Top
right panel: same for $\delta^{13}$C. Note that the $y$-axis for
$\delta^{18}$O is reversed, following common practice. Bottom left panel:
standardized residuals for $\delta^{18}$O, bottom right panel: standardized
residuals for $\delta^{13}$C.

Figure 12: Top left panel: $\delta^{18}$O data (blue) and smoothed state
$\mathds{E}[\mu_{t+\Delta t}|\\{y_{s},\;s\leq\sum\Delta t\\}]$ (red) according
to model (5) for the fourth climate state from Kalman smoothing recursions.
Top right panel: same for $\delta^{13}$C. Note that the $y$-axis for
$\delta^{18}$O is reversed, following common practice. Bottom left panel:
standardized residuals for $\delta^{18}$O, bottom right panel: standardized
residuals for $\delta^{13}$C.

Figure 13: Top left panel: $\delta^{18}$O data (blue) and smoothed state
$\mathds{E}[\mu_{t+\Delta t}|\\{y_{s},\;s\leq\sum\Delta t\\}]$ (red) according
to model (5) for the fifth climate state from Kalman smoothing recursions. Top
right panel: same for $\delta^{13}$C. Note that the $y$-axis for
$\delta^{18}$O is reversed, following common practice. Bottom left panel:
standardized residuals for $\delta^{18}$O, bottom right panel: standardized
residuals for $\delta^{13}$C.

Figure 14: Top left panel: $\delta^{18}$O data (blue) and smoothed state
$\mathds{E}[\mu_{t+\Delta t}|\\{y_{s},\;s\leq\sum\Delta t\\}]$ (red) according
to model (5) for the sixth climate state from Kalman smoothing recursions. Top
right panel: same for $\delta^{13}$C. Note that the $y$-axis for
$\delta^{18}$O is reversed, following common practice. Bottom left panel:
standardized residuals for $\delta^{18}$O, bottom right panel: standardized
residuals for $\delta^{13}$C.

## References

  * Bagchi and Mitra (1996) Bagchi, S. and S. K. Mitra (1996).  The nonuniform discrete Fourier transform and its applications in filter design. I. 1-D.  IEEE Transactions on Circuits and Systems II: Analog and Digital Signal Processing 43(6), 422–433.
  * Bell (1984) Bell, W. (1984).  Signal extraction for nonstationary time series.  The Annals of Statistics 12(2), 646–664.
  * Doornik (2021) Doornik, J. (2021).  An Introduction to OxMetrics 9.  Timberlake Consultants.
  * Durbin and Koopman (2012) Durbin, J. and S. J. Koopman (2012).  Time Series Analysis by State Space Methods (2 ed.).  Oxford University Press.
  * Gómez (1999) Gómez, V. (1999).  Three equivalent methods for filtering finite nonstationary time series.  Journal of Business & Economic Statistics 17(1), 109–116.
  * Gómez (2001) Gómez, V. (2001).  The use of butterworth filters for trend and cycle estimation in economic time series.  Journal of Business & Economic Statistics 19(3), 365–373.
  * Harvey and Koopman (2000) Harvey, A. and S. J. Koopman (2000).  Signal extraction and the formulation of unobserved components models.  The Econometrics Journal 3(1), 84–107.
  * Harvey and Koopman (2009) Harvey, A. and S. J. Koopman (2009).  Unobserved components models in economics and finance.  IEEE Control Systems Magazine 29(6), 71–81.
  * Harvey (1989) Harvey, A. C. (1989).  Forecasting, Structural Time Series Models and the Kalman Filter.  Cambridge: Cambridge University Press.
  * Harvey and Trimbur (2003) Harvey, A. C. and T. M. Trimbur (2003).  General model-based filters for extracting cycles and trends in economic time series.  Review of Economics and Statistics 85(2), 244–255.
  * Helske (2017) Helske, J. (2017).  KFAS: Exponential family state space models in R.  Journal of Statistical Software 78, 1–38.
  * Helske (2019) Helske, J. (2019).  Package “KFAS”: Kalman filter and smoother for exponential family state-space models.  Technical report, CRAN repository, version 1.3.4.
  * Karatzas and Shreve (1991) Karatzas, I. and S. E. Shreve (1991).  Brownian Motion and Stochastic Calculus (second ed.).  Springer.
  * Turner (2014) Turner, S. K. (2014).  Pliocene switch in orbital-scale carbon cycle/climate dynamics.  Paleoceanography 29(12), 1256–1266.
  * Westerhold et al. (2020) Westerhold, T., N. Marwan, A. J. Drury, D. Liebrand, C. Agnini, E. Anagnostou, J. S. Barnet, S. M. Bohaty, D. De Vleeschouwer, F. Florindo, et al. (2020).  An astronomically dated record of earth’s climate and its predictability over the last 66 million years.  Science 369(6509), 1383–1387.
  * Whittle (1983) Whittle, P. (1983).  Prediction and Regulation by Linear Least-Square Methods (second ed.).  University of Minnesota Press.
  * Zachos et al. (2001) Zachos, J., M. Pagani, L. Sloan, E. Thomas, and K. Billups (2001).  Trends, rhythms, and aberrations in global climate 65 ma to present.  Science 292(5517), 686–693.