Bayesian estimation of early warning signals

Background

After my postdoc position at PIK finished I returned to UiT for a 13 month position as a Researcher, funded by the TiPES project. During this time I started working on early warning signals, focusing on applications in climate. When a dynamical system approaches a bifurcation point, which can induce tipping points, the system will exhibit changes to statistical indicators, including increased variance and autocorrelation. Typically, these early warning signals are estimated using a sliding window approach, and my aim was to develop an alternative approach based on Bayesian hierarchical modeling. Through this work I have collaborated Luc Hallali (UiT), Martin Rypdal (UiT), Sigrunn Holbek Sørbye (UiT), Håvard Rue (KAUST), Christian Franzke (IBS), Clara Hummel (UiT), Alessandro Cotronei (UiT) and Niklas Boers (PIK/TUM).

1. Introduction

Many systems can undergo abrupt and sometimes irreversible transitions when critical thresholds, known as tipping points, are crossed. In climate, anthropogenic forcing has put several important components of the Earth system at risk of crossing tipping points, with potentially severe consequences for life on the planet (McKay et al. 2022). Examples include a potential collapse of the Atlantic Meridional Overturning Circulation (AMOC; Boers 2021), destabilization of the Greenland Ice sheets (Boers and Rypdal YEAR) and dieback of the Amazon rainforest (CITE). In order to prevent possibly irreversible changes it is of crucial to know when a system is approaching a tipping point, before it is too late to course-correct.

Tipping points are naturally studied within the framework of dynamical systems, often formulated in terms of stochastic differential equations that capture the interaction between deterministic dynamics, external forcing and random fluctuations. Let, for example, $x(t)$ denote the time-dependent state variable, representing for example temperature or ice sheet volume, $dB(t)$ denote a white noise process describing the noise in the system and $F(x)$ denote the forcing of the system. The dynamical system can be expressed by

dx(t) = F(x(t)) dt + \sigma dB(t).

One could think of this equation as describing the motion of a particle in the presence of a potential $V(x)$ , with drift $F(x)=-V'(x)$ and diffusion term $\sigma dB(t)$ describing its movement. Points where $dx(t)/dt=0$ are called fixed-points. These are considered stable if a small perturbation (in any direction) near the fixed point will decay over time, and unstable otherwise. The following figure shows the potential $V(x)$ for an example of a system with two stable fixed points (valleys), $x_1$ and $x_3$ , separated by an unstable fixed point (peak), $x_2$ . If the particle is located in the valley of an equilibrium any perturbation will be met with an opposing force pushing the particle back towards the equilibium. If the perturbation is so large that the particle crosses the unstable fixed point separating the two domains of attraction, the particle will then be pushed towards the other stable fixed point instead. This phenomenon, where the tipping point is crossed solely from perturbations caused by the internal variation of the system, is called noise-induced tipping.

On the other hand, if the dynamics of the system depend on some slowly changing control parameter $\mu(t)$ , then fixed points may shift, vanish or spawn, depending on the value of $\mu(t)$ . This means that the stability of an equilibrium may diminish over time and eventually be lost, forcing the state variable to another equilibrium. Points where the stability or number of fixed points changes are called bifurcation points, and a critical transition caused by the control parameter approaching a bifurcation point is called a bifurcation-induced tipping process.

The potential over the state variables before (a), at (b) and after (c) the control parameter μ has reached the bifurcation point. The stable fixed points are shown in blue, while the unstable fixed points (tipping point) are shown in red. In panel (b) the unstable and first stable fixed points coincide as a single unstable saddle point.

Potential (a) and bifurcation diagram (c) of a noise-induced tipping process with constant control parameter μ=1, against the potential (b) and bifurcation diagram (d) of a bifurcation-induced tipping process where the control parameter μ increases linearly from -3 to 3.

The nature of an equilibrium can be studied by investigating the linearization of the dynamical system around some stable fixed point $x_s$

dx(t) = -\lambda (x(t)-x_s) dt + \sigma dB(t),

where $\lambda = -F'(x_s)$ is the restoring rate. This is known as the Langevin stochastic differential equation and has the solution

x(t) = x_0 +\int_{0}^t e^{-\lambda t}dB(t) ,

which forms an Ornstein-Uhlenbeck (OU) process, and discretization yields an AR(1) process

x_t = \phi x_{t-1} + \varepsilon_t, \qquad \varepsilon_t \sim \mathcal{N}\left( 0, \frac{1-\phi^2}{2\lambda}\sigma^2 \right).

When a system approaches a bifurcation that induces destabilization, its restoring rate will decrease, causing it to recover more slowly from perturbations. This phenomenon is called critical slowing down and causes perturbations to persist for longer, leading to increased variance and autocorrelation that can serve as early warning signals of an impending critical transition (CITECITECITE).

The most common approach to estimating early warning signals is to estimate each indicator in a sliding window (CITECITECITE), and then to use statistics like Kendall's tau to perform a hypothesis test on whether or not a statistically significant increase in the early warning signal indicators has been detected.

Sliding window estimates of the variance and autocorrelation of the observed data up to the tipping point (vertical line). Increased variance and autocorrelations suggests that the system may be approaching a bifurcation point.

While simple, sliding windows presents some limitations as methods for detecting early warning signals. First, it requires a choice on the window length. Large windows offer higher estimation accuracy, but will be less representative of the momentary state of the indicator as it essentially applies a high pass filter, while smaller windows better capture the changes of the indicators, but will provide less accurate estimates. A good trade-off can be difficult to achieve. Moreover, if the dataset contains trends sliding windows require a more careful approach to detrending.

2. Model-based early warning signal detection

In Myrvoll-Nilsen et al. (2025) we circumvent sliding windows by introducing a model-based approach to early warning signals. In short, we introduce a time-dependent AR(1) process where the lag-one autocorrelation parameter is a linear function of time

\phi(t) = a+bt , \qquad t\in (0,1) .

We assign appropriate priors on $a$ and $b$ and incorporate this into a hierarchical Bayesian model and perform Bayesian inference using the R-INLA framework. This ensures efficient inference and gives us uncertainty quantification expressed through the posterior distributions. We say that early warning signals is detected if the posterior marginal probability of a positive slope exceeds some threshold, $1-\alpha$ , i.e.

P(b>0\mid \bm y) \geq 1-\alpha .

The inferred evolution of the time-dependent lag-one autocorrelation parameter (a), the fitted 2nd order polynomial trend (b) and the posterior marginal distribution of the memory slope parameter (c) with area above 0 shaded.

However, Boers (2021) showed that these indicators can be biased if the system is driven by external noise that itself has increasing autocorrelation or variance,leading to false positive alarms. To account for such bias, Boettner and Boers (2022) and Morr and Boers (2024) suggest that the OU process should be driven by correlated noise rather than white noise. Based on this idea, we assume in Hallali et al. (2025) that the internal noise term is itself a time-dependent AR(1) process.

x_{t+1} = \phi(t) x_t + \varepsilon_{t+1}

\varepsilon_{t+1} = \rho(t) \varepsilon_t + \xi_{t+1}

where $x$ is the state variable, with time-dependent autocorrelation parameter $\phi(t)$ , and $\varepsilon$ is the internal noise, with time-dependent autocorrelation parameter $\rho(t)$ .

In Sørbye et al. (2025) we have also looked at extending model-based approaches to long-range dependent processes, utilizing the computationally efficient AR(1) mixture approximation of Sørbye et al. (2019).

3. Applications

4. Relevant publications