An early warning system for multivariate time series with sparse and non-uniform sampling

TL;DR

This paper introduces the ratio of deviations (RoD), a new early warning statistic for multivariate systems, effective even with sparse, non-uniform, and high-frequency data, predicting bifurcations before they occur.

Contribution

The paper proposes RoD, a novel early warning test statistic, and demonstrates its effectiveness in predicting Hopf bifurcations in complex, irregularly sampled multivariate time series.

Findings

RoD is asymptotically related to autocorrelation.

RoD effectively predicts bifurcations in synthetic sparse data.

RoD performs well as a classifier on high-frequency time series.

Abstract

In this paper we propose a new early warning test statistic, the ratio of deviations (RoD), which is defined to be the root mean squared of successive differences divided by the standard deviation. We show that RoD and autocorrelation are asymptotically related, and this relationship motivates the use of RoD to predict Hopf bifurcations in multivariate systems before they occur. We validate the use of RoD on synthetic data in the novel situation where the data is sparse and non-uniformly sampled. Additionally, we adapt the method to be used on high-frequency time series by sampling, and demonstrate the proficiency of RoD as a classifier.

Tables

Table 1. Table 1 : Values used for experiments generating observations at random times.

Experiment	$\alpha$	$\beta$	Mean Sampling Time = $\left(\alpha +\beta \right)/2$
1	20	40	30
2	25	50	37.5
3	25	75	50
4	50	100	75

Table 2. Table 2 : Values used for experiments generating observations at random times.

a	AUC	$\alpha$	$\beta$	Window	a	AUC
1	0.867	25	50	500	10	0.830
1	0.980	25	75	500	10	0.937
1	0.876	50	100	750	10	0.911
1	0.840	50	100	1000	10	0.800

Equations

$${\nu }^2(t_n)$$

$${\sigma }^2(t_n)$$

$$\text{RoD}(t_n)=\frac{\nu (t_n)}{\sigma (t_n)}.$$

$$\underset{n\to \infty }{\lim }{\text{RoD}}^2(t_n)=2(1-{\rho }_X(1)),$$

$${\rho }_X(1)=\frac{\text{E}[X_{t_i}{X}_{t_{i-1}}]}{{\sigma }^2}.$$

$$\underset{n\to \infty }{\lim }{\text{RoD}}^2(t_n)$$

$$=\underset{n\to \infty }{\lim }\frac{\frac{1}{n-1}\sum _{i=2}^n{X}_{t_i}^2+\frac{1}{n-1}\sum _{i=1}^{n-1}{X}_{t_i}^2-2\frac{1}{n-1}\sum _{i=2}^n{X}_{t_{i-1}}{X}_{t_i}}{\frac{1}{n}\sum _{i=1}^n{X}_{t_i}^2}$$

$$=\underset{n\to \infty }{\lim }\left(\frac{\frac{1}{n-1}\sum _{i=2}^n{X}_{t_i}^2}{\frac{1}{n}\sum _{i=1}^n{X}_{t_i}^2}+\frac{\frac{1}{n-1}\sum _{i=1}^{n-1}{X}_{t_i}^2}{\frac{1}{n}\sum _{i=1}^n{X}_{t_i}^2}-2\frac{\frac{1}{n-1}\sum _{i=2}^n{X}_{t_{i-1}}{X}_{t_i}}{\frac{1}{n}\sum _{i=1}^n{X}_{t_i}^2}\right)$$

$$=2\frac{{\sigma }^2}{{\sigma }^2}-2{\rho }_X(1)$$

$$=2(1-{\rho }_X(1)).$$

$$X_{t_k}=a_{t_k}{X}_{t_{k-1}}+{\xi }_{t_k}.$$

$$\vec{X}(t_k)=A(t_k)\vec{X}(t_{k-1})+\vec{\xi }(t_k),$$

$$\begin{array}[]{rl}dx&=[\lambda(t)x-y+2\eta x(x^{2}+y^{2})-x(x^{2}+y^{2})^{2}]dt+\eta dW_{1}\\ dy&=[x+\lambda(t)y+2\eta y(x^{2}+y^{2})-y(x^{2}+y^{2})^{2}]dt+\eta dW_{2},\end{array}$$

$$\begin{array}[]{rl}dx&=\frac{1}{a}(3x-x^{3}-y)dt+\sigma dW_{1}\\ dy&=[x-\lambda(t)]dt+\sigma dW_{2}\\ dz&=(x-z)dt+\sigma dW_{3},\end{array}$$

Full text

1 Introduction

The prospect of identifying a tipping point in an observed system prior to its occurrence has recently excited the scientific community. Deterioration of a stable state may indicate the proximity of a tipping point, although neither necessarily implies the other [2]. The recent emphasis on early detection of tipping points has its roots in applications that are sensitive to anthropogenic activity, such as climate and ecology [4, 5]. Applications where that human intervention could prevent an undesirable regime shift makes the idea of early warning particularly enticing, and these concepts are spreading to other applications where intervention is critical, e.g. medicine [3, 12].

Popular early warning signals include observing an increasing trend in a test statistic such as autocorrelation, variance, or skewness [2, 8, 10, 12], with the theory based in dynamical systems. As a dynamical system approaches a tipping point, the equilibrium state loses stability, and the destabilization is reflected in the eigenvalues of the Jacobian. The time it takes for small perturbations to return to the equilibrium state increases as the leading eigenvalue approaches the boundary between stable and unstable in the complex plane—either the imaginary axis for systems modeled in continuous time or the unit circle for map dynamics. This phenomenon is known as critical slowing down, referring to the decreased rate of return. The requisite details on the theory of dynamical systems can be found in most introductory textbooks on the subject such as those by Strogatz [11] or Glendinning [7].

Developing early warning signals to detect catastrophic changes, with the most common being a fold or saddle-node bifurcation, has been the primary focus of the field to date. Catastrophic changes are necessarily irreversible, and thus are accompanied by a ’point of no return.’ Other bifurcations (e.g., supercritical Hopf bifurcations) are reversible, and thus are not always considered tipping points. However, in many applications it is still desirable to identify preventable or predictable changes in stability, even if they are reversible.

In this paper we present a new test statistic for use as an early warning signal, called the ratio of deviations (RoD), relating the root mean squared of successive deviations to the standard deviation of a time series. We show that for a time series generated from a weakly stationary lag-1 autoregressive process, RoD is asymptotically related to lag-1 autocorrelation. According to this relationship, as autocorrelation decreases the RoD increases. Inspired by the asymptotic relationship, we examined what an increase in the RoD indicates about the future values of a time series. We propose that an increase in RoD is indicative of a change in the nature of oscillations. Thus, the RoD should be able to detect passage through a Hopf bifurcation.

Supercritical Hopf bifurcations are reversible, and as such have not been a primary focus in the early warning literature. Subcritical Hopf bifurcations may be accompanied by a catastrophic saddle-node bifurcation of periodic orbits, with a region of bistability where the two attractors are a stable equilibrium point and a stable periodic orbit. However, even in the supercritical case, there is a hysteresis effect with a dynamic Hopf bifurcation in that a trajectory can remain near the newly unstable equilibrium for quite some time [1]. Figure 1 depicts the phenomenon. Thus, even if a Hopf bifurcation is not detected until after it has occurred, there may still be enough time left to prevent any negative effects of a regime change. Conversely, ‘reversing’ a Hopf bifurcation after it is visually apparent (i.e., after a trajectory is near the periodic orbit), could require a significant shift in the parameter.

RoD is designed to be used on short time series with irregularly sampled observations, and we extend its use to longer, high-frequency time series through sampling. To demonstrate the predictive skill of RoD, we simulated a system known to exhibit a supercritical Hopf bifurcation as a parameter passes through a critical value, plus control runs where the bifurcation parameter remains constant. We ran two types of tests: (a) to evaluate use on short time series and (b) to evaluate the extension to longer time series. RoD (paired with another test statistic) performed adequately as a predictor on the short time series. Its performance on longer time series is exceptional.

The paper will proceed as follows: Section 2 describes the theory behind the ratio of deviations as an early warning signal. Section 3 discusses the performance on synthetic data. We conclude with a discussion in Section 4.

2 Ratio of Deviations

Let ${X_{t_{k}}}$ be a univariate time series where $k=1,\ldots,n,\ldots$. The RoD relates the root mean squared of successive deviations (RMSSD), $\nu$, and standard deviation, $\sigma$. By definition, we have

[TABLE]

where $\mu=\text{E}[X]$. RoD is defined as

[TABLE]

Proposition 2.1**.**

If ${X_{t_{i}}}$ is a weakly stationary process, then

[TABLE]

where $\rho_{X}(1)$ is the lag-1 autocorrelation of $X_{t_{i}}$.

Proof.

Both standard deviation and RMSSD are independent of the mean, so we assume that $\text{E}[X]=0$ zero without loss of generality. Since $X_{t_{i}}$ is weakly stationary, we have that

[TABLE]

Thus,

[TABLE]

∎

Remark 2.2**.**

We can compute an AR(1) coefficient and RoD for any finite time series, however the preceding result relies on $\sigma^{2}$ (eventually) being time-independent.

2.1 Detecting Changes in an Underlying System

By examining a ratio of two distinct measures of dispersion, RoD is designed to detect a change in the nature of deviations in a time series. We say RoD detects a change in the variable $X_{t}$ at time $t_{k}$ both $\text{RoD}(t_{k})>\text{RoD}(t_{k-1})$. Note that we only require a single increase in RoD rather than an increasing trend, in contrast with the standard practice for other early warning signals such as autocorrelation.

It is possible that RoD increases even if both standard deviation and RMSSD decrease. To remedy this, we recommend pairing RoD with a tandem value, similar to using autocorrelation in conjunction with standard deviation [6]. We have explored requiring an increase in standard deviation or an increase in RMSSD. Requiring an increase in standard deviation reduces the number of positives more than requiring RMSSD to increase does. A third option is to set restrictions on the range of values a given variable can take, and require a value not meeting those restrictions. In the third case, the reduction in positives will clearly be dependent upon the range. The tandem metric that makes the most sense may depend on the application, a priori knowledge of the underlying system, or properties of the observations (especially on short time series).

For a univariate time series $X_{t}$, we assume the underlying model is of the form

[TABLE]

A change in $a_{t}$ will affect the stability of the system, possibly causing an increase in RoD depending on the magnitude and direction of the change. Alternatively, RoD could detect a change due to a rare event where $\xi_{t_{k}}$ takes a value far from its mean.

Because detection can be triggered due to a random event rather than deterioration of a stable state, RoD is best applied to multivariate systems. Consider the multivariate linear system

[TABLE]

where $\vec{X}(t)=(X_{1}(t),\ldots,X_{n}(t))\in\mathbb{R}^{n}$, $\vec{\xi}(t)=(\xi_{1}(t),\ldots,\xi_{n}(t))\in\mathbb{R}^{n}$, and $A(t)=\left(a_{ij}(t)\right).$ By calculating the RoD for each variable $X_{i}(t)$ individually, we are able to detect changes in particular subsystems. If we only detect a change in one variable, we conclude that it is likely due to the noise term. However, if we detect a change in all variables (or all variables of a subsystem), it likely indicates a structural change in the system. Furthermore, because we will use RoD on short time series with long times between observations, we would expect the effects of noise to subside before the next observation if the system were stable.

Note that we are detecting changes by observing an increase in the RoD of each univariate time series. Proposition 2.1 suggests that an increase in the RoD roughly corresponds with a decrease in autocorrelation. If the underlying system loses stability due to decreasing autocorrelation, we should expect the trajectory to develop oscillations of increasing amplitude. Thus we anticipate that RoD is a mechanism for detecting Hopf bifurcations.

3 Validation on Synthetic Data

We evaluated RoD as a test statistic on a suite of systems that are known to exhibit Hopf bifurcations as a parameter $\lambda$ passes through a critical value $\lambda_{c}$, including an additive white noise term in each equation. All systems were simulated according to the Euler-Maruyama method using the ‘deSolve’ package in R.

3.1 A Simple Example: Hopf Normal Form

First we demonstrate the method on a simple example: the Hopf normal form with additive white noise. The equations of the Hopf normal form are

[TABLE]

and we take $\eta=0.25.$ The system was simulated from $t=0$ to $t=100$ using time steps of 0.05. We ramped $\lambda(t)$ linearly from -1 at $t=0$ to 1 at $t=100$, so the system undergoes a dynamic Hopf bifurcation at $(t,\lambda)=(50,0)$. The time series $x(t)$ is plotted in Figure 2. We then sampled the time series for both $x$ and $y$ at times $t_{i}$, where $t_{0}=0$ and $t_{i+1}=t_{i}+\text{unif}(4,8)$ to generate sample time series $(x_{i},y_{i})$. The $x_{i}$ are depicted by triangles in Figure 2. Let $X_{n}=\{x_{0},x_{1},\ldots,x_{n}\}$ and define $Y_{n}$ similarly. We computed the RoD for $X_{n}$ and $Y_{n}$ for each $n$. We say we detect a change in the variable $x$ if both $\text{RoD}(X_{n+1})>\text{RoD}(X_{n})$ and $\sigma(X_{n+1})>\sigma(X_{n})$ (similarly for $y$). The red triangle in Figure 2 indicates the observation when a change was detected in both variables.

3.2 Experiments on a 3D Variation of the Van der Pol System

A Hopf bifurcation can be observed in a 2-dimensional subsystem, at least locally (possibly after a change of coordinates). However, in practice it is not necessarily clear which 2-dimensional subsystem one should consider nor what coordinate transformation should be used. We tested the performance of RoD as a warning signal on two different parametrizations of 3D variation of the Van der Pol system.

In certain parameter regimes, the Van der Pol system is an excitable system, meaning that a small perturbation (in a particular direction) can lead to a big oscillation [7, 11]. Excitability is related to a separation of time scales and phenomenon known as canard explosion. For our purposes, it is sufficient to think of a canard explosion as Hopf bifurcation in which the amplitude of the periodic orbits grows exponentially in terms of distance in parameter space from the bifurcation. A reader interested in delving deeper into the rich literature of canard-related phenomena can begin with works by Krupa and Szmolyan (e.g., [9]) or those of Wechselberger (e.g., [13]).

The variant of the Van der Pol system we examine is

[TABLE]

where $a$ is a time-scale parameter that affects the growth rate of the amplitude of periodic orbits. We will refer to system (8) with $a=10$ as ‘normal’ because the amplitudes grow as expected according to a normal Hopf bifurcation. We will refer to the system with $a=1$ as ‘excitable’ or ‘with canard explosion’ because the amplitudes grow much faster than normal. Figure 3 shows how the amplitudes of the periodic orbits grow after bifurcations in each case.

We ran 200 simulations of each parametrization for each of five prescribed noise intensities ($\sigma=0,0.01,0.05,0.1,0.25$.) In half of the runs $\lambda$ was ramped to pass through a Hopf bifurcation, and the other half were control runs where $\lambda$ was held constant (i.e., $\lambda(t)\equiv\lambda_{0}$). The initial value $\lambda(0)=\lambda_{0}=1.2$ was chosen so that the system had an attracting equilibrium point, and we used the equilibrium point as the initial condition for the model. Each simulation ran for 2000 dimensionless time units with steps of 0.05 units. For the simulations when $\lambda$ was ramped, we chose $d\lambda/dt$ so that $\lambda(1000)=\lambda_{c}=1$, i.e., $\lambda$ passed through the critical value exactly half-way through the simulation.

As stated earlier, we do not want to compute the RoD on the whole trajectory with high-frequency observations. Instead, we observe the trajectories iteratively at random times, so that $t_{i+1}=t_{i}+\Delta T$ where $t_{0}=0$, $\Delta T=unif(\alpha,\beta),$ and $\beta>\alpha>0.$ For each system we ran experiments varying $(\alpha,\beta)$, generating 100 time series from each trajectory for each of the values in Table 1. Additionally, we experimented with the window over which we used observations to calculate the RoD, using windows of 250, 500, 750, and 1000 time units.

Since we want to detect bifurcations before they occur at $t=1000$, the window of 1000 time units represents an experiment using the entire trajectory. Figures 4-7 break down the results of these experiments by parametrization, window length, and observation rate when using RoD with RMSSD.

3.3 Experiments with High-frequency Data

Next we turn our attention to experiments with high-frequency data. We use the same trajectories as those simulated in the previous subsection. Even though the goal is now to detect a bifurcation using high-frequency data, we still do not want to use all of the observations, or even all of the observations in a given window. Instead, we sample the data and compute the RoD on the random samples.

For this experiment, we use the same random samples of observations as in the previous subsection, we just utilize them differently. Previously, we treated each sample of observations as if it were the only data we had to generate a binary prediction as to whether or not the system would undergo a Hopf bifurcation. Now, we know that we are actually generating 100 predictions for each trajectory, and we use them to determine a probability that the system will undergo a bifurcation.

Since we are now determining a ‘probability’ for each trajectory, we can measure performance using area under the curve (AUC) where ‘the curve’ refers to the receiver operating characteristic (ROC) curve. The ROC curve depicts the performance of a binary classifier as the threshold used to separate positive predictions from negative is varied. We are not required to use a probability of 0.5 to determine which trajectories we expect undergo bifurcation. This observation means we no longer have to worry about false positives from individual samples, but rather on trajectories as a whole, allowing us to use RoD alone as a test statistic (e.g., without RMSSD or SD). We found that RoD was best as a classifier when using a window length of 500 (half of the full trajectory before bifurcation) with $(\alpha,\beta)=(25,75)$, giving us an average of 10 observations in each RoD calculation. The AUCs with these parameters were 0.937 for the ‘normal’ system and 0.98 for the ‘excitable’ system. Figure 8 depicts the ROC curves for the optimal case.

In general, RoD performed well as a classifier in each of the experiments with between 8 and 14 observations expected to lie in each window. We already discussed that a window length of 500 with an average of 10 observations per window gave the best results, but the second-best performance also came with an average of 10 observations per window (and window length of 750). Table 2 lists the AUC for each of the experiments where the expected number of observations per window was in the interval $[8,14]$.

4 Discussion

In this paper we presented a new early warning test statistic, RoD, that can be used with low-frequency data and non-uniform sampling rates. RoD adequately detects Hopf bifurcations in situations where there is one relatively short time series with sparse observations. It is interesting to note the role noise plays in the performance. No false positives are detected in the absence of noise. In the presence of noise, however, the magnitude of the noise does not affect the false positive rate. On the other hand, stronger noise improves the true positive rate before the bifurcation occurs. Recalling Figure 1, we emphasize that in the absence of noise (or when noise has sufficiently small influence), there will be a significant delay between when the bifurcation occurs and when it is observed. As such, we should think of the true positive rates in Figures 4-7 as lower bounds on the effectiveness of RoD as a detection metric. In some applications, we may even wish to detect Hopf bifurcations that have already occurred.

In situations where one would like to detect an upcoming Hopf bifurcation from a longer, high-frequency time series, we recommend producing a number of short, sparse time series by sampling the high-frequency observations at irregular intervals. When aggregating the results of RoD on the set of shorter time series, the performance is much stronger than on a single short time series.

Many of the difficulties in tipping point prediction for climatological and ecological applications arise from having only one time series. It would be interesting to see how the sampling technique affects performance of other early warning signals designed to detect catastrophic bifurcations.