# Changepoint Detection: An Analysis of the Central England Temperature Series

Xueheng Shi
Department of Statistics,
University of California, Santa Cruz
&Claudie Beaulieu
Department of Ocean Sciences,
University of California, Santa Cruz
<EMAIL_ADDRESS>
&Rebecca Killick
Department of Statistics,
Lancaster University
&Robert Lund
Department of Statistics,
University of California, Santa Cruz

###### Abstract

This paper presents a statistical analysis of structural changes in the
Central England temperature series. This series contains one of the longest
surface temperature records available and a changepoint analysis of it reveals
several interesting aspects. Regression functions with structural breaks,
including mean and trend shifts, are fitted to the series and compared via two
commonly used multiple changepoint penalized likelihood criteria. In the end,
the optimal model is judged to be one containing three location and trend
shifts, with a transition to a rapidly warming regime circa 1989. The
variability of the series is not found to be significantly changing, and shift
features are judged to be more plausible than short- or long-memory
autocorrelations. The analysis serves as a walk-through tutorial of different
changepoint techniques, illustrating what can statistically be inferred from
different models.

##  1 Introduction

Climate time series often contain structural changes and other nonlinear
features. Sudden structural changes are called changepoints and may shift
series characteristics, including means, variances, trends, and
autocorrelations. It is crucial to know changepoint times in climate series,
especially when assessing long-term trends [1, 2, 3, 4]. On average,
approximately six changepoints occur over a century of operation in a randomly
selected US climate station [5, 6].

Structural changes may be attributed to changes in climate (e.g., an increase
in the magnitude of a trend), related to decadal variability, or may be
induced by artificial changes in measurement procedures such as station
relocations or instrumentation changes. Climate time series are often
correlated, inducing memory at time scales longer than the measurement
frequency [7]. This memory is often modeled as a first-order autoregressive
(AR(1)) process in climate studies [1, 8, 9], but alternative descriptions are
possible [10].

Statistical methods to detect changepoints have rapidly evolved over the last
few decades. These include methods to detect a single shift in the series’
mean [11], in its variance [12], or in its general mean structure [13, 14].
Much of the changepoint detection literature assumes independent and
identically distributed model errors (termed white noise here). Pre-whitening
techniques to mitigate the effects of autocorrelation have been previously
used [8, 15]. It is important to allow autocorrelation and mean shifts
features in tandem in models since both features can inject similar run
patterns into a time series. [4, 16, 17] show how misleading changepoint
inferences can be should one ignore autocorrelation in a series. This point
will be reinforced as we proceed.

Multiple changepoint methods often arise with climate series. Methods designed
to detect at most one changepoint have been applied iteratively to estimate
multiple changepoint configurations through a process known as binary
segmentation [18]. Binary segmentation is now known to poorly perform at
detecting multiple changepoints [16]. See [19] for an interesting attempt to
fix binary segmentation. Penalized likelihood methods, the approach taken
here, were developed in [20, 21, 22, 23] and tend to perform better [16]. Many
of these methods allow for autocorrelation. The references above are by no
means exhaustive; indeed, the changepoint literature is currently vastly
expanding.

As most methodological papers are not written with user comprehension in mind,
the technical changepoint literature may seem impenetrable to non-
statisticians, making it challenging for the practitioner to select an
appropriate approach. Compounding difficulties, [24, 4] show that a
changepoint setup can produce spurious inferences should prominent data
features be ignored — the choice of model and method is critical. Indeed,
changepoint techniques can produce vastly different results when models and
assumptions are slightly changed.

The aim of this paper is to present, through an example, the steps needed to
conduct a comprehensive changepoint analysis of a climate series. To this end,
we analyze the central England temperature (CET) series by fitting a hierarchy
of changepoint models, from simpler models designed to detect mean shifts
only, to models capable of detecting shifts in means, linear trends, and
variances. We focus on penalized likelihood multiple changepoint techniques,
which enable us to walk the reader through a layered analysis. As our analysis
increases in structural complexity with successive layers, the underlying
assumptions of each layer are justified. Emphasis is placed on implementation
and interpretation rather than the theoretical foundations of penalized
likelihoods. Nonetheless, references to the formal statistical literature are
provided.

The rest of this paper proceeds as follows. The CET series used here is
introduced in the next section. Section 3 provides some rudimentary background
on changepoint models, describing the penalized likelihood methods used here.
The next three sections present fits of various multiple changepoint models.
Results for each fit motivate subsequent fits. Remarks about the optimal model
are made in the final section along with concluding comments.

##  2 The CET Series

The CET time-series is perhaps the longest instrumental record of surface
temperatures in the world, commencing in 1659 and now spanning 362 years
(through 2020). This record has been previously analyzed for long term changes
[25, 26, 27]; however, to our knowledge, no detailed changepoint analysis of
it has been previously conducted. Changepoints are physically plausible in the
CET record for a variety of reasons. First, artificial shifts near the
record’s onset may exist when data quality was lower [28]. Second, cooling
associated with the Maunder minimum in the early 1700s may be present [25].
Finally, an increase in the pace of climate change arising globally circa the
1960s [4] may be locally detectable. While we will not delve into attribution
of found changepoint times in detail here, the length of the CET record
affords us the opportunity to explore a variety of features.

The CET series, available at https://www.metoffice.gov.uk/hadobs/hadcet, was
provided by the UK Met Office. Measurements began in 1659 and were mostly
compiled by [29, 30] until 1973, then continued and updated to 1991 in [28].
The series is now kept by the Hadley Centre, Met Office. The annual CET time
series is a composite of 15 stations in the UK, located over a roughly
triangular area bounded by Lancashire, London, and Bristol. The series is thus
representative of the climate of the English Midlands. The locations of
stations used to form the composite time series is depicted in the top graphic
in Figure 1. The CET temperatures, presented in the bottom graphic in Figure
1, have been corrected for inhomogeneities due to changes in measurement
practices through time [29, 30, 28], and for urban warming since 1960 [31]. As
such, the CET time series can be analyzed in the context of long-term climate
change.

Figure 1: Station locations and annual average temperatures of Central
England.

##  3 Structural Change Models

To explore structural changes in the CET series, a layered changepoint
analysis, gradually building on past findings, will be conducted. Let $X_{t}$
denote the annual temperature observed at time $t$ and suppose that data from
the years $1,\ldots,N$ are available. In general, a changepoint analysis
partitions the series into $m+1$ distinct regimes, each regime having
homogeneous characteristics. The number of changepoints $m$ is unknown and
needs to be estimated from the series. Let $\tau_{i}$ denote the $i$th
changepoint time. The changepoint times partition the series into $m+1$
distinct regimes (segments); boundary conditions take $\tau_{0}=0$ and
$\tau_{m+1}=N$.

All regression models in this paper have the time series form

$X_{t}=f(t)+\epsilon_{t},\qquad t=1,2,\ldots,N,$ (1)

where $f(t)=E[X_{t}]$ is the mean of the series at time $t$. The structural
form of $f$ will vary as we proceed, generally containing location and/or
trend parameters and their shifts. The model errors
$\\{\epsilon_{t}\\}_{t=1}^{N}$ have a zero mean and may be correlated in time.
We work with AR(1) errors for simplicity, but more complex time series models
are possible. While it is important to allow for some correlation in annual
data, the form of the error correlation structure is typically not as crucial.

The AR(1) difference equation governing the errors $\\{\epsilon_{t}\\}$ is

$\epsilon_{t}=\phi\epsilon_{t-1}+Z_{t},$

where $\phi\in(-1,1)$ and $\\{Z_{t}\\}$ is zero mean white noise (WN) with
unknown variance $\sigma^{2}$. Solutions to the AR(1) equation have
geometrically decaying correlations:
$\mbox{Corr}(\epsilon_{t},\epsilon_{t+h})=\phi^{h}$ for $h\geq 0$. Because the
data are annually averaged, a Gaussian errors $\\{\epsilon_{t}\\}$ is
statistically realistic. An implication of this is that future model
likelihoods will be Gaussian based.

To estimate the changepoint structure and model parameters from the data,
penalized likelihood methods will be used. The idea behind a penalized
likelihood is to choose the changepoint configuration that minimizes the
penalized likelihood

$-2\log(L^{*}(m;\tau_{1},\ldots,\tau_{m}))+P(m;\tau_{1},\ldots,\tau_{m}).$ (2)

The notation here is as follows: $L^{*}(m;\tau_{1},\ldots,\tau_{m})$ is the
optimal time series likelihood that can be achieved from a model with $m$
changepoints that occur at the times $\tau_{1},\ldots,\tau_{m}$ for the fixed
data sample $X_{1},X_{2},\ldots,X_{N}$. To determine
$L^{*}(m;\tau_{1},\ldots,\tau_{m})$, one must estimate all parameters in the
mean function $f$ and the AR(1) model errors assuming that $m$ changepoints
occur at the times $\tau_{1},\ldots,\tau_{m}$. This procedure will be
discussed further below. The quantity $P(m;\tau_{1},\ldots,\tau_{m})$ is the
penalty for having a model with $m$ changepoints at the times
$\tau_{1},\ldots,\tau_{m}$. As more and more changepoints are added to the
model, the overall fit gets better ($L^{*}$ gets larger); the penalty term
prevents the model fit from being overparametrized (having too many
changepoints).

Many penalty structures have been proposed in the statistics and climate
literature. These include Akaike information criterion (AIC), Bayesian
information criterion (BIC), modified Bayesian information criterion (mBIC),
and Minimum description lengths (MDL). The two that we will use here are BIC
and MDL. These two penalties were judged as "winners" in a changepoint
detection context in the recent statistical comparison in [16]. AIC penalties
are not considered here because they often estimate an excessive number of
changepoints [16]. The BIC penalty for having $m$ changepoints at the times
$\tau_{1},\ldots,\tau_{m}$ is simply $m\log(N)$; additional parameters are
penalized at the rate of $\log(N)$ per model parameter.

The BIC penalty is proportional to the number model parameters. In contrast,
the MDL penalty is more complex in form, also accounting for the changepoint
location times $\tau_{1},\ldots,\tau_{m}$. The MDL penalty depends on the form
of $f$ and is rooted in information theory, which quantifies the computer
memory needed to store model specifications (the best model uses the minimal
space). MDL penalties have proven useful in changepoint detection [20, 23]).
Other climate penalty literature for changepoint tasks includes [32]. Our
penalized likelihood objective functions for structural changes are summarized
in Table 1. When $m=0$, the MDL penalty is zero in all cases since no shift
parameters enter the model.

Table 1: Penalized likelihoods for our models

Criteria | Objective function
---|---
BIC | $N\log(\hat{\sigma}^{2})+N+N\log(2\pi)+(2m+3)\log(N)$
MDL | $N\log(\hat{\sigma}^{2})+N+N\log(2\pi)+2\log(N)+2\log(m)+\sum_{i=1}^{m+1}\log(\tau_{i}-\tau_{i-1})+2\sum_{i=2}^{m}\log(\tau_{i})$
(a) Penalized likelihood functions for $E[X_{t}]=\mu_{r(t)}$ assuming AR$(1)$
errors

Criteria | Objective Function
---|---
BIC | $N\log(\hat{\sigma}^{2})+N+N\log(2\pi)+(3m+4)\log(N)$
MDL | $N\log(\hat{\sigma}^{2})+N+N\log(2\pi)+2\log(N)+2\log(m)+2\sum_{i=1}^{m+1}\log(\tau_{i}-\tau_{i-1})+2\sum_{i=2}^{m}\log(\tau_{i})$
(b) Penalized likelihood functions for $E[X_{t}]=\mu_{r(t)}+\beta_{r(t)}t$
assuming AR$(1)$ errors

Criteria | Objective Function
---|---
BIC | $N\log(\hat{\sigma}^{2})+N+N\log(2\pi)+(3m+3)\log(N)$
MDL | $N\log(\hat{\sigma}^{2})+N+N\log(2\pi)+\log(N)+2\log(m)+2\sum_{i=1}^{m+1}\log(\tau_{i}-\tau_{i-1})+2\sum_{i=2}^{m}\log(\tau_{i})$
(c) Penalized likelihoods for $E[X_{t}]=\mu_{r(t)}+\beta_{r(t)}t$ assuming
white noise errors

Criteria | Objective function
---|---
BIC | $N\log(\hat{\sigma}^{2})+N+N\log(2\pi)+(2m+4)\log(N)$
MDL | $N\log(\hat{\sigma}^{2})+N+N\log(2\pi)+3\log(N)+2\log(m)+\sum_{i=1}^{m+1}\log(\tau_{i}-\tau_{i-1})+2\sum_{i=2}^{m}\log(\tau_{i})$
(d) Penalized likelihood functions for $E[X_{t}]=\mu_{r(t)}+\beta t$ assuming
AR$(1)$ errors

Methods for handling multiple changepoint analyses without penalized
likelihoods exist. One popular technique is termed binary segmentation [33].
Binary segmentation works with any single changepoint technique, termed an at
most one change (AMOC) method. Many AMOC tests have been developed, including
cumulative sums (CUSUM) [34], likelihood ratios [35], Chow tests [36], and sum
of squared CUSUM tests [16]. Binary segmentation first analyzes the entire
series for a changepoint. If a changepoint is found, the series is split into
subsegments at the changepoint time and the two subsegments are further
scrutinized for additional changepoints. The procedure is repeated iteratively
until no subsegments are deemed to have additional changepoints. While simple
and computationally convenient, binary segmentation is one of the poorer
performing multiple changepoint techniques [16], often being fooled by
changepoints that occur close to one another or shifts that move the series in
antipodal directions. Attempts to fix binary segmentation exist — see the wild
binary segmentation and related methods in [19] and [37]. Unfortunately, these
techniques often assume independent model errors or are restricted to single
parameter changes per regime (for example, mean shifts only).

A drawback of penalized likelihood methods involves computation time. There
are $\binom{N-1}{m}$ distinct changepoint configurations having $m$
changepoints. Summing this over all $m$ shows that there are $2^{N-1}$
distinct changepoint configurations that need to be searched in an exhaustive
optimization of a penalized likelihood, a daunting task for $N=200$ on even
the world’s fastest computers. As a solution, genetic algorithms (GA) will be
used to optimize our penalized likelihoods. GAs are random walk searches that
are unlikely to visit suboptimal changepoint configurations. While [23] show
how to devise a GA in climate changepoint applications, generally available
GAs have now become savvy enough to confidently handle our needs. The GA
optimizations performed here use the R package GA [38]. We mention that [22]
and [39], two rapid dynamic programming based multiple changepoint
configuration optimizers, currently cannot handle our needs: [39] assumes
independence of model errors and [22] assumes all parameters change at each
changepoint time (including the correlation and error variances). GAs are the
only optimization method that reasonably handle all models considered here.

##  4 Mean Shifts

The simplest regression form in a changepoint analysis is that of mean shifts.
This model posits $f(\cdot)$ to have form

$\displaystyle f(t)=\begin{cases}\mu_{0},\quad&1\leq t\leq\tau_{1},\\\
\mu_{1},\quad&\tau_{1}+1\leq t\leq\tau_{2},\\\ \;\quad\vdots\\\
\mu_{m},\quad&\tau_{m}+1\leq t\leq N.\end{cases}$ (3)

The model’s mean structure is compactly written as $f(t)=E[X_{t}]=\mu_{r(t)}$,
where $r(t)\in\\{0,\ldots,m\\}$ denotes the regime index in force at time $t$;
for example, $r(t)=0$ for $1\leq t\leq\tau_{1}$.

Given $m$ and the changepoint times $\tau_{1},\ldots,\tau_{m}$, the means are
first estimated via segment averages:

$\hat{\mu}_{k}=\frac{1}{\tau_{k}-\tau_{k-1}}\sum_{t=\tau_{k-1}+1}^{\tau_{k}}X_{t},\quad
k=0,1,\ldots,m.$

While sample means are not the exact maximum likelihood estimators of the mean
parameters for correlated series, they are typically very close and are easy
to compute (unlike maximum likelihood estimators). A detailed discussion of
least squares versus maximum likelihood estimator differences for time series
is contained in [40]. Next, a regime-wise mean is subtracted from the series
via $D_{t}=X_{t}-\hat{f}(t)=X_{t}-\hat{\mu}_{r(t)}$. The AR$(1)$ parameter is
then estimated via

$\hat{\phi}=\frac{\sum_{t=1}^{N-1}D_{t}D_{t+1}}{\sum_{t=1}^{N}D_{t}^{2}}.$ (4)

One-step-ahead predictions of the time series are now computed by

$\hat{D}_{t}=\hat{\phi}\hat{D}_{t-1},\quad t\geq 2,$ (5)

with the start-up condition $\hat{D}_{1}=0$. The white noise variance in the
AR(1) model is estimated as

$\hat{\sigma}^{2}=\frac{1}{N}\sum_{t=1}^{N}\hat{D}_{t}^{2}.$ (6)

Calculations will reveal a negative Gaussian log-likelihood of [23]

$-2\log(L^{*}(m;\tau_{1},\ldots,\tau_{m}))=N\log(\hat{\sigma}^{2})+\underbrace{N+N\log(2\pi).}_{\text{Constant}}$
(7)

The constant term does not change over distinct changepoint configurations —
it can be neglected in the changepoint configuration comparisons. The above
equations show how to evaluate model likelihoods given the changepoint
configuration; the optimal changepoint configuration is found by a GA search.

The fitted mean shift estimated for the CET series found by our GA is
summarized in Table 2. The two penalties yield different changepoint
configurations, with MDL estimating seven changepoints and BIC five. It is
common to get different segmentations with different penalties; changepoint
analyses crucially depend on the choice of penalty. These differences aside,
both penalties pinpoint 1989 as a changepoint time. MDL segmentations often
place changepoints close together, signifying an outlier or outlying run in
the data. Here, both MDL and BIC deemed the “cold year" in 1740 an outlier,
bracketing this time by two changepoints. Because MDL methods are based on
information theory and not large sample statistical asymptotics, they often
flag outliers. Estimates of the AR(1) autocorrelation parameter are
$\hat{\phi}=0.13$ for BIC and $\hat{\phi}=0.06$ for MDL. While neither value
represents heavy autocorrelation, the differences in the two estimates bring
us to a salient point: as the number of changepoints in a model fit increases,
the estimated level of autocorrelation tends to decrease. In an AR(1) series,
$\hat{\phi}$ is the estimated correlation between two successive years of
observations. Indeed, positive autocorrelation and mean shifts can inject
similar run patterns into a time series — it is important to put both
components into a model and let the likelihood decide which is the
statistically preferable structure. Should autocorrelation be neglected, one
risks flagging spurious changepoints.

Figures 2 and 3 plot the BIC and MDL estimated mean shift configurations
against the series; both fits are visually appealing. The MDL fit places more
shifts at the beginning of the record, perhaps suggesting that the data during
these times is less reliable. Evident in this fit is that the last two (BIC)
or three (MDL) regimes are successively moving the series higher in a
"staircase", which may be due to the presence of a long-term warming trend. To
investigate this, we move to a changepoint analysis where trends are allowed.

Figure 2: Estimated CET mean shift structure. BIC flags five mean shifts in
1706, 1740, 1741, 1893, and 1989.
Figure 3: Estimated CET mean shift structure. MDL flags seven mean shifts in
1691, 1699, 1727, 1740, 1741, 1911, and 1989.

##  5 Trend Shifts

A model that allows for a linear trend in $f(\cdot)$ is the piece-wise linear
form

$f(t)=\begin{cases}\mu_{0}+\beta_{0}t,\qquad&1\leq t\leq\tau_{1},\\\
\mu_{1}+\beta_{1}t,\qquad&\tau_{1}+1\leq t\leq\tau_{2},\\\ \qquad\vdots\\\
\mu_{m}+\beta_{m}t,\qquad&\tau_{m}+1\leq t\leq N,\end{cases}.$ (8)

One can compactly write this as $E[X_{t}]=f(t)=\mu_{r(t)}+\beta_{r(t)}t$ using
our definition of $r(t)$.

While statistical methods for trend estimation are well developed, most
changepoint theory for models with trends only allows for shifts in the
intercept parameter [41, 42]. Even in the single changepoint case, it is not
clear whether the trend slopes in the two regimes should be constrained as
identical; see the debate in [24] and [43]. Compared to the mean shift case,
research on multiple trend shift models is sparse. However, [44] present a
dynamic programming approach to identify trend shifts using a penalty based on
absolute distances that is neither the MDL nor BIC. Their $\\{\epsilon_{t}\\}$
must be white noise (uncorrelated) with a zero mean and constant variance. See
[45] and [46] and their related R package strucchange by [47] for more
details.

The least squares estimators for the $i$th regime’s parameters are computed
from data in this regime only:

$\hat{\beta}_{i}=\frac{\sum_{t=\tau_{i}+1}^{\tau_{i+1}}(X_{t}-\bar{X}_{i})(t-\bar{t}_{i})}{\sum_{t=\tau_{i}+1}^{\tau_{i+1}}(t-\bar{t}_{i})^{2}},\quad\hat{\mu}_{i}=\bar{X}_{i}-\hat{\beta}_{i}\bar{t}_{i},$

where
$\bar{X}_{i}=(\sum_{t=\tau_{i}+1}^{\tau_{i+1}}X_{t})/(\tau_{i+1}-\tau_{i})$
and $\bar{t}_{i}=(\tau_{i+1}+\tau_{i}+1)/2$. While again these are not the
exact maximum likelihood estimators in correlated settings, they again are
typically very close to them [40].

One next computes the detrended series via

$D_{t}=X_{t}-\hat{f}(t)=X_{t}-(\hat{\mu}_{r(t)}+\hat{\beta}_{r(t)}t).$ (9)

The equations in (4)-(7) apply without modification. One caveat is that since
regression lines contain two parameters, we require all regimes to be at least
three units long (so that fits in any single regime are not perfect).

GA optimizations of the BIC and MDL penalized likelihoods yield identical
configurations, both flagging three breaks at the times $1700,1739$, and
$1988$. This methodological agreement is convenient, but is not typical in
changepoint analyses. Table 2 lists the other parameters in the model fit.
Cooling occurs during the first 39 years, followed by an increasing-trend
second regime, with subsequent shifts to two warming trend regimes. The last
regime, entered in 1989, is warming with a trend of $1.1^{\circ}$C per
century.

Visually, the piecewise linear function, depicted in Figure 4, looks to better
describe the series than the simple mean shift model. The AR(1) correlation
estimate moves down to $\hat{\phi}=0.06$. This parameter tests as
insignificantly different from zero with standard time series methods [10].
When $\phi=0$, an AR(1) model reduces to white noise. This point is
emphasized: our model fits prefer the trend shift mean structure over
structures involving correlated errors. The estimated white noise variance is
$\hat{\sigma}^{2}=0.29$.

Figure 4: Estimated CET trend shift structure. Both BIC and MDL flag three
shifts in 1700, 1739, and 1988.

We now examine residuals from the fitted model in Figure 5. For AR(1) errors
with $\phi=0$, these are simply the $\\{D_{t}\\}$ in (9). The conclusion that
$\phi=0$ implies that $\\{D_{t}\\}$ is uncorrelated in time $t$. The residuals
also appear to have an approximate zero mean (in fact, the sample mean of
residuals fitted from least squares methods will be exactly zero). What is not
clear is whether or not $D_{t}$ has a constant variance in time $t$.

Figure 5: Residuals of the trend shift model.

[4] explore issues related to changes in the second moment of a climate
series. To further examine the variance structure of $\\{\hat{D}_{t}\\}$, we
develop a model that allows for variance shifts in a zero mean independent
Gaussian series. This can be constructed by positing that $\hat{D}_{t}\sim
N(0,\sigma^{2}(t))$, where $\sigma^{2}(t)$ has the shift structure

$\sigma^{2}(t)=\begin{cases}v_{0}^{2},\qquad&1\leq t\leq\tau_{1},\\\
v_{1}^{2},\qquad&\tau_{1}+1\leq t\leq\tau_{2},\\\ \qquad\vdots\\\
v_{m}^{2},\qquad&\tau_{m}+1\leq t\leq N,\end{cases}.$ (10)

We do not require $m$ or $\tau_{1},\ldots,\tau_{m}$ to be the same as those in
the previous analyses.

The maximum likelihood estimators of the regime variances are

$\hat{v}_{\ell}^{2}=\frac{1}{\tau_{\ell+1}-\tau_{\ell}}\sum_{t={\tau_{\ell}}}^{\tau_{\ell+1}}\hat{D}_{t}^{2}.$

Plugging this back into the Gaussian likelihood gives

$-2\log(L(m;\tau_{1},\ldots,\tau_{m}))=\sum_{k=0}^{m}(\tau_{k+1}-\tau_{k})\hat{v}_{k}^{2}+N\log(2\pi)+N.$

Moving to penalties, the BIC and MDL penalties are

$\displaystyle\mbox{BIC}(m;\tau_{1},\tau_{2},\cdots,\tau_{m})$
$\displaystyle=\sum_{k=0}^{m}(\tau_{k+1}-\tau_{k})\hat{v}_{k}^{2}+N\log(2\pi)+N+(2m+1)\log(N),$
$\displaystyle\mbox{MDL}(m;\tau_{1},\tau_{2},\cdots,\tau_{m})$
$\displaystyle=\sum_{k=0}^{m}(\tau_{k+1}-\tau_{k})\hat{v}_{k}^{2}+N\log(2\pi)+N+$
$\displaystyle\quad+2\log(m)+\sum_{i=1}^{m+1}\log(\tau_{i}-\tau_{i-1})+2\sum_{i=2}^{m}\log(\tau_{i}).$

A GA was again used to optimize the penalized likelihoods. Here, we restricted
each segment to have length two or more so that no regime achieves a perfect
fit ($\hat{v}_{k}>0$ for each segment $k$). In this case, both BIC and MDL
prefer models with no changepoints. Hence, we do not find evidence of changing
variabilities in the Central England annual temperatures.

For completeness, we consider a model with a common trend slope in all regimes
as in [41]:

$X_{t}=\mu_{r(t)}+\beta t+\epsilon_{t},$

where $\beta$ is the trend slope applicable to the whole series, $\mu_{r(t)}$
is as in (3), and the errors $\\{\epsilon_{t}\\}$ are an AR$(1)$ process.

The ordinary least square estimators of $\beta$ and
$\mu_{0},\mu_{1},\ldots,\mu_{m}$ have the explicit form

$\hat{\beta}=\frac{\sum_{i=0}^{m}\sum_{\tau_{i}+1}^{\tau_{i+1}}(X_{t}-\bar{X}_{i})(t-\bar{t}_{i})}{\sum_{i=0}^{m}\sum_{\tau_{i}+1}^{\tau_{i+1}}(t-\bar{t}_{i})^{2}},\quad\hat{\mu}_{i}=\bar{X}_{i}-\hat{\beta}\bar{t}_{i},\quad
i=0,1,\ldots,m,$

where $\bar{X}_{i}$ and $\bar{t}_{i}$ are as before. These are again very
close to the maximum likelihood estimators [40].

A GA was used to find the best fitting model, which is plotted against the
data in Figure 6. Both BIC and MDL flag a single mean shift in 1988. Since the
BIC and MDL penalized likelihoods in Table 2 are larger for the constant slope
model than for regime varying trend slopes, the inference is that regime
varying slopes are needed.

Figure 6: Estimated CET trend shift structure imposing a constant regime trend
slope. Both BIC and MDL flag a single changepoint in 1988.

##  6 Additional Models

This section briefly considers some other models used to describe temperature
trends. While our purpose is to be comparative and complete, we do not provide
full statistical detail.

###  6.1 Joinpin models

In the climate literature, there is debate over whether trend models should
impose continuity in $E[X_{t}]$ at the changepoint times [48]. These so-called
joinpin models require that $E[X_{t}]=f(t)$ is continuous in time $t$. Here,
we compare a joinpin model to the mean and trend shift models fitted in
previous sections. Unfortunately, it is not clear what an appropriate MDL
penalty is for this case, nor does this seem to be an easy matter to rectify.
A BIC penalty is available and is simply $3\log(N)$.

To fit a joinpin model, the package in [44] was used. The estimated error
variance of $0.29$ from the discontinuous model and a BIC penalty were used to
select a changepoint configuration. This model fit assumes IID errors, which
does not seem implausible from our past work. The fitted model flags a single
changepoint in 1972, as depicted in Table 2 and Figure 7. The resulting fit is
visually appealing at the end of the series; however, the fit at the start of
the time series is poor, with long periods of overestimated and underestimated
temperatures. The fit is stable against changes in the estimator variance.
When compared to other fits, the joinpin model has a much higher BIC than the
mean and trend shift models (see Table 2). As such, we do not comment further
about joinpin structures.

Figure 7: Estimated CET joinpin shift structure. BIC flags one shift in 1972.

###  6.2 Long-memory models

A body of climate literature argues that climate time series exhibit long-
memory, where the autocorrelation function decays slowly, often via a power
law [49, 50]. Long-memory and changepoint features can inject similar run
properties into a climate series, which is appreciated in the statistical and
econometric literatures [51, 52, 53, 54]. [55] suggest that the daily CET time
series exhibits long-memory.

To compare our results to a long-memory model, we fit an autoregressive
fractionally integrated moving-average (ARFIMA) model to the CET series. In
particular, ARFIMA models with no moving-average component, an integration
parameter $d$ with $0<d<0.5$, and an autoregressive component of orders zero
and one, are considered. These are legitimate long-memory models.

To fit the ARFIMA models, the R package fracdiff [56] was used. A BIC penalty
was calculated and is listed in Table 2. An MDL penalty is not informative
since this model does not have any changepoints. The fitted long-memory model
has the largest BIC amongst all models compared here, suggesting that
changepoints, rather than long-memory, is more plausible in the CET series.
For additional evidence that changepoints are preferred to long-memory
features, we applied the time varying wavelet spectrum methods in [57] to the
CET series. These methods were used on surface temperatures in [58] and shown
to well discriminate changepoint versus long-memory models in long time
series. The results confirm that a changepoint model is more appropriate than
a long-memory model. The fitted model of order zero was also preferred to the
fitted model of order one, reinforcing that correlation aspects in the CET
series are minimal.

##  7 Comments, Conclusions, and Discussion

This study compared several common changepoint model fits (mean shifts, trend
shifts, and joinpin), as well as a long-memory model, to the CET time series.
To our knowledge, this is the first time a detailed changepoint analysis has
been conducted on this long record. Starting with a simple mean shift model,
several different changepoint model fits, illustrating the techniques and
salient points of changepoint analyses, were presented.

Among all fitted models, the winner (best fitting model), which has the
smallest BIC and MDL scores, has trend and mean shifts in 1700, 1739, and
1988. Table 2 presents the log-likelihood, BIC, and MDL scores of all model
fits along with some additional information. Table 3 provides estimates of our
best fitting model’s intercept and slope parameters. Having both BIC and MDL
penalties agree on model type and changepoint configuration adds robustness to
our conclusions. Models with shifts were deemed preferable to models having
correlation.

While a formal attribution of these three changes is not our emphasis, some
brief words about the detected shifts is provided. First, shifts flagged
before 1772 may be due to inferior data quality over this period [27]. As
such, the first two shifts may represent artificial changes in the record.
This said, the shift in 1700 also coincides with the rapid cooling of the
Maunder minimum [25]. This feature of the CET data is known and should be
flagged by most changepoint approaches. A third changepoint flagged in 1988 is
not surprising given that the rate of warming has increased at the global
level in the 1960/70s in a range of surface temperature records, as discussed
in studies using both trend shift and joinpin models [3, 4, 48, 59].

A joinpin model has been argued as more physically plausible for representing
temperatures [48] and was included in this comparison. Our fitted joinpin
configuration flags a single changepoint in 1972. While any flagged
changepoint seems plausible towards the end of the record due to increased
warming, the joinpin fit to the earliest data is very poor. This can be
visually seen in Figure 7, but is also reflected by the higher MDL and BIC
scores for that model in Table 2. A joinpin model should be used when a
discontinuous mean function is unlikely or physically implausible. With the
CET series, it is not evident whether the estimated mean function should be
continuous or discontinuous. Elaborating, for series containing “only a single
station", mean discontinuities are physically expected. However, when more and
more station records are averaged, mean function discontinuities are reduced,
becoming more implausible with an increasing number of stations. Should a
discontinuous mean function be deemed possible, a trend shift model provides
greater flexibility since it can simultaneously approximate a joinpin mean
structure and detect discontinuous shifts [4]. Also, a model containing only
mean shifts will flag a sequence of shifts in an attempt to follow a long term
trend should the data have a trend and it not be included in the model. Should
there be any suspicion about a trend in the record, we recommend using
techniques that allow for trends.

As an aside, this paper did not detour into target minus reference series
aspects. Indeed, a more detailed changepoint analysis could be done by
subtracting all $\binom{15}{2}=105$ pairs of series from one another and
examining these differences for changepoints. Once all individual series are
homogenized (changepoint adjusted), the 15 series could be reaveraged for a
better composite CET series. See [6] for more on this procedure.
Unfortunately, we did not have the 15 original series at our disposal;
moreover, the technical methods presented here would not change with this
data.

A residual analysis was conducted to ensure that the underlying assumptions of
the model were met. With the CET series, residuals of the trend shifts were
judged to be uncorrelated (white noise). However, climate time series often
exhibit autocorrelation that should be taken into account. We stress the
importance of verifying the underlying assumptions in any changepoint model.
Indeed, neglecting positive autocorrelation raises the risk of detecting
spurious shifts. Also, the series’ autocorrelation may be more complex than
AR(1) and may itself contain shifts [2, 4]. Some climate series may also
contain long-memory autocorrelations [60]. An additional challenge lies with
the ambiguity between long-memory and changepoint models as both features can
yield “similar looking" series. Because of this, a long-memory model was
included as part of our comparison. We found that the CET time series is best
represented by a multiple trend shift changepoint structure and not a long-
memory model. Such a comparison is not possible for all climate series since
lengthy records are required to analyze any long-memory series [58]. The CET
time series, which is the longest publicly available surface temperature
series, enables this comparison.

Table 2: Model fitting results

Model | Penalty | Flagged Changepoints | Log-likelihood | Penalized Score
---|---|---|---|---
Mean shifts+AR(1) | BIC | 1706, 1740, 1741, 1893, 1989 | -287.91 | 652.41
MDL | 1691,1699,1727,1740,1741,1911,1989 | -276.26 | 649.00
Trend shifts+AR(1) | BIC | 1700,1739,1988 | -289.42 | 655.43
MDL | 1700,1739,1988 | -289.42 | 645.98
Trend shifts+WN | BIC | 1700,1739, 1988 | -290.52 | 651.74
MDL | 1700,1739, 1988 | -290.52 | 642.30
Mean shifts+fixed slope+AR(1) | BIC | 1988 | -310.22 | 655.79
MDL | 1988 | -310.22 | 653.30
Joinpin | BIC | 1972 | -321.19 | 654.17
Long memory+AR(1) | BIC | NA | -316.59 | 656.75
Long memory | BIC | NA | -319.31 | 655.93

Table 3: Parameters estimates of the best fitting model: trend shifts and white noise errors Segment |  Intercept (∘C) |  Slope (∘C/yr)
---|---|---
1659-1699 |  9.25 |  -0.027
1700-1738 |  7.78 |  0.026
1739-1987 |  8.78 |  0.002
1988-2020 |  6.50 |  0.011

## Acknowledgement

Rebecca Killick gratefully acknowledges funding from EP/R01860X/1 and
NE/T006102/1. Robert Lund and Xueheng Shi thank funding from NSF DMS-2113592.

## Data statement

The Central England data used in this study is available at
https://www.metoffice.gov.uk/hadobs/hadcet/. We used the annual records from
1659-2020.

## References

  * [1] Robert Lund, Xiaolan L Wang, Qi Qi Lu, Jaxk Reeves, Colin M. Gallagher, and Yang Feng.  Changepoint detection in periodic and autocorrelated time series.  Journal of Climate, 20(20):5178–5190, 2007.
  * [2] Claudie Beaulieu, Jie Chen, and Jorge L. Sarmiento.  Change-point analysis as a tool to detect abrupt climate variations.  Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 370(1962):1228–1249, 2012.
  * [3] Niamh Cahill, Stefan Rahmstorf, and Andrew C Parnell.  Change points of global temperature.  Environmental Research Letters, 10(8):084002, 2015.
  * [4] Claudie Beaulieu and Rebecca Killick.  Distinguishing trends and shifts from memory in climate data.  Journal of Climate, 31(23):9519 – 9543, 2018.
  * [5] J. Murray Mitchell Jr.  On the causes of instrumentally observed secular temperature trends.  Journal of Atmospheric Sciences, 10(4):244–261, 1953.
  * [6] Matthew J Menne and Claude N Williams Jr.  Homogenization of temperature series via pairwise comparisons.  Journal of Climate, 22(7):1700–1717, 2009.
  * [7] Klaus Hasselmann.  Stochastic climate models part I. Theory.  Tellus, 28(6):473–485, 1976.
  * [8] Michael Robbins, Colin Gallagher, Robert Lund, and Alexander Aue.  Mean shift testing in correlated data.  Journal of Time Series Analysis, 32(5):498–511, 2011.
  * [9] D.L. Hartmann, A.M.G. Klein Tank, M. Rusticucci, L.V. Alexander, S. Bronnimann, Y. Charabi, F.J. Dentener, E.J. Dlugokencky, D.R. Easterling, A. Kaplan, B.J. Soden, P.W. Thorne, M. Wild, and P.M. Zhai.  Observations: Atmosphere and Surface, book section 2, page 159–254.  Cambridge University Press, Cambridge, United Kingdom and New York, NY, USA, 2013.
  * [10] P. J. Brockwell and R. A. Davis.  Time Series: Theory and Methods.  Springer-Verlag, 2nd edition edition, 1991.
  * [11] Herman Chernoff and Shelemyahu Zacks.  Estimating the current mean of a normal distribution which is subjected to changes in time.  The Annals of Mathematical Statistics, 35(3):999–1018, 1964.
  * [12] Der-Ann Hsu.  Tests for variance shift at an unknown time point.  Journal of the Royal Statistical Society: Series C, 26(3):279–284, 1977.
  * [13] Richard E Quandt.  The estimation of the parameters of a linear regression system obeying two separate regimes.  Journal of the American Statistical Association, 53(284):873–880, 1958.
  * [14] Michael W Robbins, Colin M. Gallagher, and Robert Lund.  A general regression changepoint test for time series data.  Journal of the American Statistical Association, 111(514):670–683, 2016.
  * [15] Francesco Serinaldi and Chris G Kilsby.  The importance of prewhitening in change point analysis under persistence.  Stochastic Environmental Research and Risk Assessment, 30(2):763–777, 2016.
  * [16] Xueheng Shi, Colin Gallagher, Robert Lund, and Rebecca Killick.  A comparison of single and multiple changepoint techniques for time series data.  Computational Statistics and Data Analysis, Minor revision in review:arXiv preprint arXiv:2101.01960, 2021.
  * [17] Colin M. Gallagher, Rebecca Killick, Robert Lund, and Xueheng Shi.  Autocovariance estimation in the presence of changepoints.  Submitted to Journal of Time Series Analysis, arXiv preprint arXiv:2102.10669, 2021.
  * [18] Sergei N Rodionov.  A sequential algorithm for testing climate regime shifts.  Geophysical Research Letters, 31(9), 2004.
  * [19] Piotr Fryzlewicz.  Wild binary segmentation for multiple change-point detection.  Annals of Statistics, 42(6):2243–2281, 2014.
  * [20] Richard A. Davis, Thomas C. M. Lee, and Gabriel A. Rodriguez-Yam.  Structural break estimation for nonstationary time series models.  Journal of the American Statistical Association, 101(473):223–239, 2006.
  * [21] QiQi Lu, Robert Lund, and Thomas C. M. Lee.  An MDL approach to the climate segmentation problem.  Annals of Applied Statistics, 4:299–319, 2010.
  * [22] Rebecca Killick, Paul Fearnhead, and Idris A Eckley.  Optimal detection of changepoints with a linear computational cost.  Journal of the American Statistical Association, 107(500):1590–1598, 2012.
  * [23] Shanghong Li and Robert Lund.  Multiple changepoint detection via genetic algorithms.  Journal of Climate, 25(2):674–686, 2012.
  * [24] Robert Lund and Jaxk Reeves.  Detection of undocumented changepoints: A revision of the two-phase regression model.  Journal of Climate, 15(17):2547–2554, 2002.
  * [25] G. Plaut, M. Ghil, and R. Vautard.  Interannual and interdecadal variability in 335 years of central England temperatures.  Science, 268(5211):710–713, 1995.
  * [26] David I. Harvey and Terence C. Mills.  Modelling trends in Central England temperatures.  Journal of Forecasting, 22(1):35–47, 2003.
  * [27] Eric Hillebrand and Tommaso Proietti.  Phase changes and seasonal warming in early instrumental temperature records.  Journal of Climate, 30(17):6795–6821, 2017.
  * [28] David E Parker, Tim P Legg, and Chris K Folland.  A new daily Central England temperature series, 1772–1991.  International Journal of Climatology, 12(4):317–342, 1992.
  * [29] Gordon Manley.  The mean temperature of Central England, 1698–1952.  Quarterly Journal of the Royal Meteorological Society, 79(340):242–261, 1953.
  * [30] Gordon Manley.  Central England temperatures: monthly means 1659 to 1973.  Quarterly Journal of the Royal Meteorological Society, 100(425):389–405, 1974.
  * [31] David Parker and Briony Horton.  Uncertainties in Central England temperature 1878–2003 and some improvements to the maximum and minimum series.  International Journal of Climatology, 25(9):1173–1188, 2005.
  * [32] Henri Caussinus and Olivier Mestre.  Detection and correction of artificial shifts in climate series.  Journal of the Royal Statistical Society: Series C (Applied Statistics), 53(3):405–425, 2004.
  * [33] Andrew John Scott and M Knott.  A cluster analysis method for grouping means in the analysis of variance.  Biometrics, 30:507–512, 1974.
  * [34] Ewan S Page.  Continuous inspection schemes.  Biometrika, 41(1):100–115, 1954.
  * [35] Venkata Jandhyala, Stergios Fotopoulos, Ian MacNeill, and Pengyu Liu.  Inference for single and multiple change-points in time series.  Journal of Time Series Analysis, 34(4):423–446, 2013.
  * [36] Gregory C Chow.  Tests of equality between sets of coefficients in two linear regressions.  Econometrica: Journal of the Econometric Society, 28:591–605, 1960\.
  * [37] Birte Eichinger and Claudia Kirch.  A MOSUM procedure for the estimation of multiple random change points.  Bernoulli, 24(1):526–564, 2018.
  * [38] Luca Scrucca.  GA: a package for genetic algorithms in R.  Journal of Statistical Software, 53(4):1–37, 2013.
  * [39] Robert Maidstone, Toby Hocking, Guillem Rigaill, and Paul Fearnhead.  On optimal multiple changepoint algorithms for large data.  Statistics and Computing, 27(2):519–533, 2017.
  * [40] Jaechoul Lee and Robert Lund.  A refined efficiency rate for ordinary least squares and generalized least squares estimators for a linear trend with autoregressive errors.  Journal of Time Series Analysis, 33(2):312–324, 2012.
  * [41] QiQi Lu and Robert Lund.  Simple linear regression with multiple level shifts.  Canadian Journal of Statistics, 35(3):447–458, 2007.
  * [42] Colin M. Gallagher, Robert Lund, and Michael Robbins.  Changepoint detection in climate time series with long-term trends.  Journal of Climate, 26(14):4994–5006, 2013.
  * [43] Xiaolan L. Wang.  Comments on “Detection of undocumented changepoints: A revision of the two-phase regression model”.  Journal of Climate, 16(20):3383–3385, 2003.
  * [44] Robert Maidstone, Paul Fearnhead, and Adam Letchford.  Detecting changes in slope with an ${L}_{0}$ penalty.  ArXiv preprint arXiv:1701.01672, 2017.
  * [45] Jushan Bai and Pierre Perron.  Estimating and testing linear models with multiple structural changes.  Econometrica, 66:47–78, 1998.
  * [46] Jushan Bai and Pierre Perron.  Computation and analysis of multiple structural change models.  Journal of Applied Econometrics, 18(1):1–22, 2003.
  * [47] Achim Zeileis, Friedrich Leisch, Kurt Hornik, Christian Kleiber, Bruce Hansen, Edgar C Merkle, and Maintainer Achim Zeileis.  Package ‘strucchange’.  R package version, pages 1–5, 2015.
  * [48] Stefan Rahmstorf, Grant Foster, and Niamh Cahill.  Global temperature evolution: recent trends and some pitfalls.  Environmental Research Letters, 12(5):054001, 2017.
  * [49] Naiming Yuan, Minghu Ding, Yan Hiuang, Zuntao Fu, Elena Xoplaki, and Juerg Luterbacher.  On the long-term climate memory in the surface air temperature records over antarctica: A nonnegligible factor for trend evaluation.  Journal of Climate, 28(15):5922 – 5934, 2015.
  * [50] Richard Blender and Klaus Fraedrich.  Long time memory in global warming simulations.  Geophysical Research Letters, 30(14), 2003.
  * [51] Francis X. Diebold and Atsushi Inoue.  Long memory and regime switching.  Journal of Econometrics, 105(1):131–159, 2001.
  * [52] Clive W. J. Granger and Namwon Hyung.  Occasional structural breaks and long memory with an application to the S&P 500 absolute stock returns.  Journal of Empirical Finance, 11(3):399–421, 2004.
  * [53] Terence C. Mills.  Time series modelling of two millennia of Northern Hemisphere temperatures: long memory or shifting trends?  Journal of the Royal Statistical Society: Series A (Statistics in Society), 170(1):83–94, 2007.
  * [54] Chun Yip Yau and Richard A. Davis.  Likelihood inference for discriminating between long-memory and change-point models.  Journal of Time Series Analysis, 33(4):649–664, 2012.
  * [55] Joanna Syroka and Ralf Toumi.  Scaling of Central England temperature fluctuations?  Atmospheric Science Letters, 2(1):143–154, 2001.
  * [56] Martin Maechler.  fracdiff: Fractionally Differenced ARIMA aka ARFIMA(p,d,q) Models, 2020.  R package version 1.5-1.
  * [57] Ben Norwood and Rebecca Killick.  Long memory and changepoint models: a spectral classification procedure.  Statistics and Computing, 28(2):291–302, 2018.
  * [58] Claudie Beaulieu, Rebecca Killick, David Ireland, and Ben Norwood.  Considering long-memory when testing for changepoints in surface temperature: A classification approach based on the time-varying spectrum.  Environmetrics, 31(1):e2568, 2020.
  * [59] Eric Ruggieri.  A Bayesian approach to detecting change points in climatic records.  International Journal of Climatology, 33(2):520–528, 2013.
  * [60] D. I. Vyushin, P. J. Kushner, and Francis Zwiers.  Modeling and understanding persistence of climate variability.  Journal of Geophysical Research: Atmospheres, 117(D21), 2012.