Critical slowing down in dynamical systems driven by nonstationary correlated noise

Precursor signals for bifurcation-induced critical transitions have recently gained interest across many research ﬁelds. Common indicators, including variance and autocorrelation increases, rely on the dynamical system being driven by white noise. Here, we show that these metrics raise false alarms for systems driven by time-correlated noise, if the autocorrelation of the noise process increases with time. We introduce an indicator for systems driven by nonstationary short-term memory noise, and show that this indicator performs well in situations where the classical methods fail.


I. INTRODUCTION
Many theoretical, experimental, and real-world systems undergo abrupt shifts from one dynamical regime to another. Anticipating or predicting such regime shifts is of vital importance in many instances, such as for major Earth system components or ecosystems, but also in physiology, neurology, or finance [1][2][3][4][5][6]. Such regime shifts can often be described in terms of critical transitions of low-dimensional nonlinear dynamical systems, induced by codimension-1 bifurcations.
The theory behind this type of transition goes back to the fluctuation-dissipation theorem [7,8] and is understood well [9]. Transitions triggered by codimension-1 bifurcations in stochastic dynamical systems with additive noise forcing are associated with a characteristic widening of the underlying basin of attraction on the way to the bifurcation [9,10]. This leads to a weakening of the restoring force to perturbations, which causes the phenomenon of critical slowing down (CSD) [11][12][13]. CSD changes the statistical properties of the system and, in particular, leads to an increase of the variance and lag-1 autocorrelation (AC1) before the critical threshold is reached. A data-driven approach to identify such statistical precursor signals for forthcoming critical transitions thus relies on measuring the variance and AC1 of a time series encoding the dynamics of the system in sliding windows, and searching for consistent increases in both indicators. This specific method has already found a number of applications to models and real-world data [14][15][16][17].
Precursor signals for critical transitions in terms of rising variance and AC1 have been identified in ecological and geological time series, controlled biological experiments, and a multitude of model simulations of different systems [10,12,22]. Besides these classical CSD indicators, which are most widely used in practical applications, a number of additional indicators have been proposed, including changes in skewness [9], diverging susceptibility [23], mean exit time [24], and increases in complexity of ARMA models [25]. In recent years there have been considerable efforts to improve the performance of these statistical methods by reducing false and raising true positive detection rates. Existing work has mainly focused on systems driven by uncorrelated (white) noise, which makes the methods prone to failure if this assumption is not fulfilled-which is likely in the climate system [26])-although some work has been done to construct methods for systems with correlated noise using detrended fluctuation analysis [27,28].
Here, we propose an analytical early-warning indicator for systems driven by noise with exponentially decaying autocorrelation structure. We show that the classical CSD indicators produce false alarms if the autocorrelation of the driving noise increases, while the indicator proposed here remains reliable in such cases. We thereby provide the analytical basis and confirmation for recent work that followed a purely data-driven approach to account for changing temporal correlation in the driving noise [29]. The method introduced here also suits itself to generalizations of the noise model and constitutes a flexible approach for situations when the classical CSD indicators are not applicable.
The structure of this paper is as follows: We first describe the occurrence of the classical CSD indicators for traditional one-dimensional systems driven by white noise. Thereafter, we show that the classical indicators cannot accurately predict transitions in cases of correlated noise with increasing autocorrelation, as described by a nonstationary AR(1) process. We then propose an early-warning indicator based on estimating the varying AR(1) parameters. This estimator reduces to the AC1 in the white-noise case. We then demonstrate the performance of our method, and compare it to existing alternatives, on simple mathematical models.

A. Classical critical slowing down
The occurrence of CSD can be formally explained on the basis of stochastic dynamics. Assume that the system is well described by a deterministic function f and driven by high-frequency forcing that can be represented by additive stochastic noise. The typical noise model assumed in the context of CSD is additive, stationary white noise with variance σ 2 , so that the dynamical equation reads for system states x, control parameter r, and a Wiener process W t . If the system is close to equilibrium, locally it suffices to consider the linearized dynamics since f (x * ) = 0 by definition. Of course, the equilibrium position x * can itself change in an open system, hence x * = x * (r(t ), t ). By performing a coordinate change x → (x − x * ) this part of the dynamics can be absorbed. Under the assumption that the control parameter changes sufficiently slowly, so that the system state can track the equilibrium position (quasistatic assumption), the linear approximation holds after the transformation. With λ:= − f (x * , r(t ), t ), Eq. (2) takes the formẋ where the negative sign causes λ to be a positive number for stable solutions. For constant (or slowly varying) λ, this equation defines the Ornstein-Uhlenbeck process; in its stationary state, the parameter λ is related to the variance Var(x t ) and the autocorrelation AC(x s , x t ) by [30] Var In this framework, the loss of (linear) stability when approaching the bifurcation manifests as a decrease in λ; by Eq. (4) this leads to the increases in variance and AC1 that are widely used to identify CSD in a given time series. Note that it has also previously been proposed to use the restoring rate λ as a parametric indicator for an approaching bifurcation [29,31,32], while some nonparametric approaches have been devised to analyze the full nonlinear system [33,34]. In practice, the statistical properties are usually calculated on discrete time series data. Sampling a realization of the Ornstein-Uhlenbeck process X t := x n t , n ∈ Z, in time steps of size t = const, leads to a discrete-time process which is a first-order autoregressive (AR1) model. In the stationary case, Eq. (4) implies that (see [35] and Supplemental Material [36]). Therefore, a decrease in λ, which leads to an increase of ϕ, in turn leads to an increase in variance and AC1 also when estimated from a discrete time series. The above derivation already reveals some limitations of the classical approach based on increasing variance and AC1. In particular, an increase of the variance may not only be caused by a decrease in the restoring rate λ, but also by an increase in σ . An increase in both variance and AC1 is therefore necessary for reliable detection of critical slowing down, even if the assumptions on white additive noise are justified; this has been prominently discussed in [13]. In cases where the assumption of white noise cannot be justified, however, further theoretical development is needed, as we will show in the following.

B. Critical slowing down in systems driven by correlated noise
In addition to assuming σ = const, the classical CSD indicators in terms of rising variance and AC1 rely on the assumption that the autocorrelation of the process only increases due to a weakening in the restoring force. A natural extension to the process is a model in which the autocorrelation may change independently of critical slowing down. In fact, for any real-world natural system, it would be very hard to justify that the driving noise is white, and correlated noise structures should be considered much more realistic. A straightforward generalization of Eq. (3) is thus an Ornstein-Uhlenbeck process driven by another Ornstein-Uhlenbeck process [37]. By increasing the autocorrelation parameter μ of the noise term, one directly manipulates both the autocorrelation and the variance of the process x t . Since the noise process is described in itself by an Ornstein-Uhlenbeck process, the covariance of the noise term decreases exponentially in time. This model can thus be considered as a linear model driven by a short-term memory process, a type of noise that is very often considered for (paleo-)climate records [38]. A discrete analog for this model is which has the following statistical properties (see Supplemental Material [36]): From these equations, it is clear that a loss of stability (and increase in ϕ) still leads to the increases in variance and AC1, just as for the case of white-noise forcing. However, in this model an increase in ρ also increases both quantities. This implies that classical CSD indicators, i.e., rising variance and AC1, are prone to false positives if there are no strong reasons 013230-2 to assume that the noise forcing is indeed white, since in that case the autocorrelation of the noise term may increase. An example is given in Fig. 1.

Estimating ϕ as an indicator of critical slowing down
Casting CSD into the above framework permits the definition of an early-warning indicator. The increase of the AR(1) parameter ϕ in Eq. (8) is a direct measure of the system's stability, and the motivation to use the variance and AC1 relies on the behavior of ϕ. We therefore propose to estimate ϕ directly, and to use it as a robust indicator for CSD for which the likelihood of false positives is reduced.
One can estimate ϕ from time series data using an ordinary least-squares (OLS) regression. Given the time series data y = (y 0 , y 1 , . . . , y n ), the OLS fit function is The best estimator for ϕ is given by the value ϕ minimizing S, given by The OLS estimator is unbiased if i.e., if the noise is uncorrelated.

White noise
In the classical situation with white-noise forcing, Eq. (13) holds and the estimator is unbiased. Since E(X t ) = 0 for an AR(1) process, Eq. (12) coincides with the sample estimation procedure for the linear (Pearson) correlation coefficient r(X t+1 , X t ) and is therefore identical to the AC1 estimation, consistent with Eq. (6). Thus, direct estimation of ϕ is equivalent with the classical CSD indicators if the time series itself can be described well by an AR(1) process driven by white noise.

AR(1) noise
For the model described by Eq. (8), direct estimation of ϕ using Eq. (12) yields a biased result. It is, however, possible to construct an unbiased estimator [39][40][41]. One begins with the biased least-squares estimator for a time series y: For large n this quantity converges to a biased estimator [39]: which only coincides with ϕ for the white-noise case ρ = 0.
One can see that an increase in ρ alone leads to an increase in ϕ b (for ϕ, ρ ∈ (0, 1)), which makes ϕ b unsuited as an indicator of CSD. However, the convergence properties of the least-squares estimator for ρ, are also known to be [39] Using both of these properties, one can construct an adjusted, unbiased estimator [40] for ϕ > ρ, given by or alternatively for ρ > ϕ by This estimation procedure works as long as the processes are stationary (−1 < ϕ, ρ < 1). For a detailed discussion on the convergence properties of these estimators, see [40,42].

III. RESULTS
To demonstrate the benefit of our method, the performance of the classical CSD indicators, variance and AC1, are compared to the simple parameter ϕ and the adjusted parameter ϕ a for two models. The first model is a linear model without the potential for abrupt transitions, and the second is a nonlinear model with a double-well potential, which thus allows for a critical transition from one fixed point to the other. Both are driven by white noise and nonstationary AR(1) noise for comparison. The equations are discretized using the Euler-Maruyama method and the different CSD indicators are calculated in rolling windows (see Supplemental Material for details on the employed parameters [36]). For the AR(1) noise driving the systems, the autocorrelation parameter ρ is increased linearly from 0 to 0.6 over the course of the simulation. The linear trend of the quantities is calculated using a linear least-squares regression (for the double-well model only data prior to the transition are considered).
An increase in one of the CSD indicators might of course also happen due to random fluctuations. To evaluate the significance of these trends we use a surrogate test based on Fourier surrogates with randomized phases. For a given time series, 1000 surrogates are created. These surrogates coincide with the original time series in their relevant statistical properties; in particular, the total variance and total lag-1 autocorrelation are the same. However, any true trend will be removed in the surrogates, as desired for such a surrogate null model. The likelihood that a trend in the original time series is real and not caused by random fluctuations can then be estimated by comparing the magnitude of the trend in the original time series to the trends in the surrogates. This is done by calculating the proportion of trends that have a larger magnitude in the surrogates than the original time series. The test can be turned into a binary classifier by choosing a significance threshold for this test statistic-the proportion of surrogates that have a larger trend than the original time series. For example, a threshold of 0.05 would mean a maximum of 5% of the surrogate trend magnitudes are allowed to be larger than the trend in the original time series for the trend to be considered statistically significant. For more information on this procedure see [43] as well as Appendix, and for application to CSD estimation see [15,44,45].
The performance of the different CSD indicators is tested for varying values of this threshold. To quantify the performance of the indicators, the number of significant positives trends is calculated for ensembles of 500 realizations. The proportions of these ensembles that yield significant trendsaccording to the test prescribed above-as a function of the threshold are shown in Figs. 2 and 3. the white-noise case and Eq. (7) for the nonstationary AR(1) noise case. This linear model is used as a baseline for the performance of the methods, for which no indication of CSD should occur. The trends in the estimator are calculated for 10 4 and 10 5 total data points in the time series (see Fig. 2 and Supplemental Material [36]). For the white-noise case, the classical indicators as well as ϕ and ϕ a method all detect a number of significant (false) positive trends that matches the significance threshold (Fig. 2). This is the expected behavior for this null model in case of a correctly working CSD indicator. In contrast, based on variance, AC1, or ϕ, a significant excess of false positive detections occurs in the case of nonstationary AR(1) noise. As predicted, the changes in the noise term are misinterpreted by these indicators as signs of CSD; increases in the data availability only increase the false positive rate since these effects are more clearly resolved. In contrast, the adjusted ϕ a does not respond to the changes in the noise forcing and hence performs exactly as expected for the AR(1) noise, and similarly as the classical methods in the white-noise case.

B. Double-well model
Next, we compare the CSD indicator proposed here, i.e., ϕ a , to the previously proposed indicators for time series obtained with the nonlinear double-well model, given in the continuous case bẏ (r evolves linearly from −1 to 1; see Supplemental Material [36] for further details). For a varying control parameter r, this model undergoes a saddle-node bifurcation preceded by a widening of the basin of attraction. To calculate the CSD indicators in this case, the resulting time series has to be detrended beforehand in order to avoid biases induced by any underlying nonlinear trends. In the case at hand, we detrend the individual stochastic realizations by subtracting the corresponding deterministic realization (σ = 0). For this model, significantly positive trends constitute a true positive detection of CSD and hence a justified alarm of a forthcoming critical transition. The classical indicators perform well in the case of white-noise forcing and exhibit a high number of true positives prior to the transition. Their performance increases with the size of the time series (see Supplemental Material [36]). For the AR(1) noise model, the number of significant trends increases faster than for the white-noise model, indicating changes in the noise term are misattributed as signs of CSD. In contrast, ϕ a performs similarly in the white-noise case and in the AR(1) noise case, indicating that ϕ a is able to disentangle effects in the noise term from actual CSD on the way to a bifurcation-induced transition. Its performance in the white-noise case is however slightly lower than for the classical methods by sometimes misattributing some of the effects of CSD to noise effects. Therefore, using just this method slightly increases the chance for false negative detections.

IV. DISCUSSION
Investigating the performance of the classical CSD indicators on simple models reconfirms that CSD provides powerful tools for the prediction of bifurcation-induced transitions. Under the right circumstances, these tools yield near 100% effective precursor signals.
However, we also demonstrated that the classical CSD indicators fail in cases where correlated noise forcing cannot be excluded, which is likely the case in many real-world situations, such as when analyzing climate data; blindly applying the classical CSD indicators, i.e., variance and AC1, likely leads to incorrect conclusions. We demonstrated this by the performance of variance and AC1 on a linear model driven by nonstationary AR(1) noise; in this case the two classical CSD indicators lead to a near 100% false positive rate. These results call attention to the need for further research into indicators for CSD in the presence of correlated noise.
Directly estimating the restoring rate of the system, or equivalently the closely related AR(1) parameter ϕ, presents strong candidates for improvement. We have demonstrated that one can readily adjust the estimation to obtain unbiased values ϕ a for different and especially for nonstationary noise models. Using this unbiased estimate of the AR(1) parameter, i.e., ϕ a , directly as a CSD indicator works robustly in cases of short-memory noise with varying variance and autocorrelation. We emphasize that using the unbiased estimate ϕ a as introduced above efficiently decreases the probability of false positive detections, which is very high if the classical indicators are used in the-in practical cases very likely-situation of noise forcing with changing temporal correlation strength. At the same time, our results indicate that using this method might slightly increase the false negative detection rate for actual transitions, by misattributing trends in the autocorrelation to changes in the noise properties. For real-world time series data, the choice of method should therefore be informed by additional information about the most probable noise model.
Since the estimation procedure for ϕ a is more involved than for the classical methods, it is likely that more data points will be needed in practical applications of this method compared to the classical methods. This limits the usability of this method in situations where only sparse data are available, e.g., when analyzing (paleo-)climate time series. We emphasize that the generalization made here is also applicable to other noise models. For example, Brouste and Iacus [46] propose a method to estimate the restoring rate for time series data derived from an Ornstein-Uhlenbeck process driven by a fractional Gaussian noise process. This quantity can be used as a CSD indicator in cases where the noise of the system is better modeled by a long-term memory noise process, and can be complementary to the method proposed here.

APPENDIX: RANDOM PHASE FOURIER SURROGATE TEST
To test the significance of a (linear) trend within a time series, a surrogate data approach is taken. For a given time series y, a number of pseudodatasets is used. These surrogate sets are created using the method of random phase Fourier surrogates. For a dataset y = (y 0 , y 1 , . . . , y n ), the algorithm works as follows.
(1) Calculate the discrete Fourier transform, FT (y) k = Y k = n−1 j=0 y n e − i2π n k j . (2) Create a set of uniformly distributed random phases, e iφ k .
This method is used since the surrogate dataset has the same variance and (linear) correlation structure as the original time series, but since the components are randomly shuffled any real (linear) trend is removed. The linear trend's significance is estimated by fitting a linear function f (x) = ax + b to the original and surrogate data using an ordinary least-squares regression; the trend is indicated by the magnitude and sign of the parameter a. If there is a significant (positive) trend in the original time series this parameter will be much larger for the real dataset than for most of the surrogate set, since any trend in the surrogate data is a purely random effect stemming from the variance and autocorrelation of the time series. Consequently, if the trends in the surrogate data are of similar or larger magnitude than for the original time series, it is an indicator that the trend is insignificant and can be explained purely by stochastic effects.
Following this logic, a hypothesis test can be established: Create a large number of surrogates (1000 in practice) and calculate the trend parameter a for the real (a ref ) and all surrogate datasets (a s ). The probability that the trend in the original time series can be explained by stochastic effects alone is then p = 1 − number of a s for which a s a ref total number of a s .
This can be turned into a binary classifier by choosing an appropriate significance threshold value for p (which is varied for Figs. 2 and 3 in the main text).