Information dynamic spectrum characterizes system instability toward critical transitions

Abstract

This paper addresses the need of characterizing system instability toward critical transitions in complex systems. We propose a novel information dynamic spectrum framework and a probabilistic light cone method to automate the analysis. Our framework uniquely investigates heterogeneously networked dynamical systems with transient directional influences, which subsumes unidirectional diffusion dynamics. When the observed instability of a system deviates from the prediction, the method automatically indicates the approach of an upcoming critical transition. We provide several demonstrations in engineering, economics, and social systems. The results suggest that early detecting critical transitions of synchronizations, sudden collapse, and exponential growth is possible.

1 Introduction

Detection and prediction of emerging tipping points are major challenges in complex systems [1]–[5], because while self-organized emerging interactions can facilitate information exchange, they increase the risk of attack or failure. When a system operates in a high-risk, unstable region, a small perturbation can induce a critical transition that leads to catastrophic failures. Although studies of data-driven computational models, with prior knowledge of individual systems, have greatly advanced the understanding of such emerging phenomena, we remain unequipped to accurately detect and predict tipping points prior to critical transitions [2].

Current works in emerging tipping point detection focus on homogeneous structured dynamical systems. In [6], Scheffer et al. give a seminal review on early warning signals for detecting critical transitions in the ecological domain. Signals such as increased temporal correlation, skewness, and spatial correlation of population dynamics are used to quantify the phenomena of critical slowing down as early warning indicators. These methods do not address heterogeneously networked systems. Through the remainder of this paper, networked dynamical system will be abbreviated as network. Scheffer et al. in [7] highlight the need to anticipate critical transitions in heterogeneous networks with empirical indicators. Such heterogeneity typically leads to connectivity patterns that are statistically characterized by heavy tails, large fluctuations, scale-free properties, and non-trivial correlations, which cannot be described by traditional modeling techniques that rely on the homogeneity assumption of a network [8]. Consequently, network science has evolved from studying complex systems by modeling single, non-interacting networks to modeling interdependent networks [9]. In particular, Buldrev et al.[4] show that a network of networks is extremely vulnerable with respect to random failures: a random removal of a small fraction of the nodes from a network can trigger a catastrophic cascade of failures. The interconnectedness and interdependency of complex networks further pose hyper risks [1]. Understanding complex interactions and coupling dynamics in large-scale complex networks is undoubtedly a crucial step toward preventing dangerous system behaviors.

The spectral early warning signals (EWS) theory [10], [11], inspired by [6], is one of the few attempts to detect critical transitions in heterogeneous networks. The spectral EWS theory states that the covariance spectrum can quantify the phenomenon of critical slowing down in heterogeneous networks by mathematically proving the link between complex network structures and observed time series. In particular, the leading covariance eigenvalue provides a structure-invariant indicator of upcoming critical transitions. As the system reaches the tipping point, the self-organized structure reveals itself via the components of covariance spectrum. Although spectral EWS quantifies how much the elements of a system change together, the symmetric nature of covariance spectrum does not permit the analysis of directional influences among elements.

Transfer entropy (TE) [12] and symbolic transfer entropy (STE) [13] have been proposed to identify directional influence in complex systems. For instance, STE is used to analyze brain activity data for the detection and identification of asymmetric dependences of brain regions in epileptic seizure activity [13]. In [14], the transfer entropy matrix is used on financial market data to analyze the asymmetrical influence of mature markets on emerging markets. Although the transfer entropy quantified in [13], [14] showed the promising results in analyzing financial and neurophysiological data, the changing structures and dynamics of networks were not addressed. Therefore, there is a need to quantify local transfer entropy that is changing in time [15].

Our proposed information dynamic spectrum framework detects and indicates critical transitions by quantifying both global and local directional influences in heterogeneous networks. Our framework is based on a novel associative transfer entropy (ATE) measure which decomposes directional influence of transfer entropy into associative states of the influences. We transform multivariate time series of a complex system into the spectrum of the transfer entropy matrix (TEM) and the spectrum of the associative transfer entropy matrix (ATEM) in order to capture information dynamics of the system. We develop novel spectral radius measures of TEM and ATEM to detect early warning signs of source-driven instability and to reveal the sources and dynamics of directional influences. Motivated by the TE and ATE behaviors observed in different critical transition examples, we then develop a method to automate the generation of early warning indicators. In particular, the nature of convex growth of spectral radius of TEM and ATEM (TE and ATE trajectories), prior to critical transitions, enables us to generate probabilistic light cones of instability trajectories using natural logarithmic curve modeling. This method enables us to analyze system instability trajectories, which in turn provides early indications of critical transitions. We demonstrate our methods on four types of transitions: (1) oscillatory synchronizations of networked non-Foster circuits; (2) pitchfork bifurcations of chaotic systems with various canonical network structures; (3) abrupt falls of directional influences in Latin America stock indices and Dow Jones during the 2008 financial crisis; and (4) exponential growth in Wikipedia editing behaviors. Our results suggest that one can analyze the trajectories of system instability to indicate the nature of upcoming critical transitions prior to their tipping points (see Table 1).

Table 1 Examples of critical transitions demonstrated in this paper

2 Information dynamic spectrum

We propose an information dynamic spectrum framework for monitoring complex systems. The assumptions in this work are (1) the outputs of the elements in the complex system can be observed as time progresses, i.e. time series are available, and (2) the outputs are the results of internal dynamics within individual elements and external interactions with other elements. An example is the time series of stock indices over a period of time. The fluctuations of a stock index depend on its own dynamics, as well as the hidden interactions with other indices. One example application of the proposed framework is to monitor the dynamics of this complex network of stock indices.

We will here establish the notations. Given the time series of a system of m elements: $X(t)={[x_1(t),x_2(t),\dots ,x_m(t)]}^T$, where t is time and $t=1,2,\dots ,N$, $x_i(t)$ is the time series of the i th element, and T represents the transpose of a vector, an interesting problem is to analyze the information dynamics of the system. In the following subsections, we will first describe the definition of transfer entropy, and then introduce a novel measure, associative transfer entropy, followed by the proposed information dynamic spectrum framework.

2.1 Transfer entropy

Transfer entropy (TE) is a directional measure of information flow between a pair of time series. A time series is a sequence of time-ordered observations of a node in a network, or an element of a system. TE quantifies how much information is transferred from the current state into the future from one time series to another. The concept of TE extends from the concept of mutual information, which is a measure of mutual dependence between two time series and more mutual information means less uncertainty if knowing one time series. However, mutual information is symmetric and does not inform the direction the information flows. To address this, TE takes into account the dynamics of information transport and provides how much information is transferred from the current state into the future.

To give a formal definition of TE, let $x=\{x(1),x(2),x(3),\dots ,x(N)\}$ and $y=\{y(1),y(2),y(3),\dots ,y(N)\}$ be two time series of interest and denote $x^{(\ell )}(t)=(x(t-\ell +1),\dots ,x(t))$ the vector of the time series x at time t with length-ℓ history. The amount of TE from source x to destination y with a time lag τ in the future is defined as [12],

$${TE}_{x\to y}=\sum _{y(t+\tau ),y^{(k)}(t),x^{(\ell )}(t)}p(y(t+\tau ),y^{(k)}(t),x^{(\ell )}(t))log\frac{p(y(t+\tau )|y^{(k)}(t),x^{(\ell )}(t))}{p(y(t+\tau )|y^{(k)}(t))},$$

(1)

where $p(\cdot ,\cdot ,\cdot )$ represents the joint probability and $p(\cdot |\cdot ,\cdot )$ and $p(\cdot |\cdot )$ represent the conditional probabilities. This definition is derived from the Shannon entropy: $H=-{\sum }_{z}p(z)logp(z)$, a measure of uncertainty. The definition of ${TE}_{x\to y}$ in Eq. (1) can be expressed as $H(y(t+\tau )|y^{(k)}(t))-H(y(t+\tau )|y^{(k)}(t),x^{(\ell )}(t))$. This is the uncertainty in y at time $t+\tau$ given its history $y^{(k)}(t)$, minus the uncertainty in y at time $t+\tau$ given its history $y^{(k)}(t)$ and x’s history $x^{(\ell )}(t)$. Therefore, TE measures the net amount of information transferred from x to y. An important feature of TE is that it is asymmetric, in contrast to mutual information, which is a symmetric measure. More details about TE, including some interesting examples, can be found in the original paper [12].

2.2 Associative transfer entropy

Since TE only quantifies the net amount of information going from a source to a destination, it does not distinguish the types of effects the information transferred. The idea of the proposed associative transfer entropy (ATE) is to decompose TE by constraining the associated states of processes. It is often important to identify the types of the information flow, rather than the gross amount of information flow. For instance, the amount of influence for an increase in x that leads to a likely increase in y in the future can be the same as the amount of influence for an increase in x that leads to a likely decrease in y. There is no distinction between these two very different outcomes using TE. Our proposed ATE addresses such a difference and is defined by:

$${ATE}_{x\to y}^S=\sum _{(y(t+\tau ),y^{(k)}(t),x^{(\ell )}(t))\in S}p(y(t+\tau ),y^{(k)}(t),x^{(\ell )}(t))log\frac{p(y(t+\tau )|y^{(k)}(t),x^{(\ell )}(t))}{p(y(t+\tau )|y^{(k)}(t))},$$

(2)

where S is the associate state, a subset of the set of all possible states of $(y(t+\tau ),y^{(k)}(t),x^{(\ell )}(t))$. The set S represents a certain associated state between x and y. The purpose of ATE is to capture information transfer between two variables for a particular association of their states. For example, for binary time series, TE can be decomposed into a positive association ATE+ and a negative association ATE−, in order to distinguish two situations where the amount of information transferred may be the same but with opposite effects: in one case the value of the source x drives the value of the destination y in the same direction, and in the other case the source drives the destination in the opposite direction.

To illustrate this case of positive and negative influences, we simulate two time series $x(t)$ and $y(t)$ with a binary difference between the current value and the previous value, where 1 represents an increment and 0 represents a decrement. The probability of increment or decrement is conditioned by the previous increment or decrement. Figure 1 illustrates a few simulated examples with the following setup. Let $\dot{x}(t)=x(t)-x(t-1)$ and $\dot{y}(t)=y(t)-y(t-1)$ be the current time-step value minus the previous time-step value. We fix the conditional probabilities

$$p(\dot{x}(t+1)=a|\dot{x}(t)=b)=p(\dot{y}(t+1)=a|\dot{y}(t)=b)=p(\dot{x}(t+1)=a|\dot{y}(t)=b)=0.5,$$

where $a,b=0$ or 1. Given that there is an increment in x at t, the probability that there will be an increment at $t+1$ is equally likely. Keeping these conditional probabilities fixed at 0.5, the only bias will come from $p(\dot{y}(t+1)=0|\dot{x}(t)=0)=p(\dot{y}(t+1)=1|\dot{x}(t)=1)=P$. When $P<0.5$, an increment in x is more likely to cause a decrement in y, and a decrement in x is more likely to cause an increment in y. Figure 1(a)-(f) shows a few such simulated time series of length 1000 from arbitrary initializations with $P=0,0.2,0.4,0.6,0.8,1$, respectively.

Figure 1

Example simulations of positive and negative influence. The conditional probabilities of $p(\dot{x}(t+1)|\dot{x}(t))$, $p(\dot{y}(t+1)|\dot{y}(t))$, and $p(\dot{x}(t+1)|\dot{y}(t))$ are fixed and unbiased. The only bias is from $\dot{x}(t)\to \dot{y}(t+1)$: $p(\dot{y}(t+1)=0|\dot{x}(t)=0)=P$ and $p(\dot{y}(t+1)=1|\dot{x}(t)=1)=P$. The influence from x to y is negative when $P<0.5$, as seen in (a)-(c). Similarly, the influence from x to y is positive when $P>0.5$, as seen in (d)-(f).

Figure 2 plots TE, ATE+, and ATE− as functions of P, for $0\le P\le 1$, averaged over 100 trials for each P. In this binary case, we decompose TE into positive association ATE+ and negative association ATE−, summing over the sets $S_{+}=\{(\dot{y}(t+1),\dot{x}(t),\dot{y}(t))=(0,0,0),(0,0,1),(1,1,0),\text{ or }(1,1,1)\}$ and $S_{-}=\{(\dot{y}(t+1),\dot{x}(t),\dot{y}(t))=(0,1,0),(0,1,1),(1,0,0),\text{or }(1,0,1)\}$, respectively. Note that ATE+ sums over positive associations, which means that $\dot{y}(t+1)$ and $\dot{x}(t)$ are either both 0 or both 1. ATE− sums over negative associations, which means that one of $\dot{y}(t+1)$ and $\dot{x}(t)$ is 0 and the other is 1. In Figure 2(a), as x is the only influence for y in the future, we see that when $P>0.5$ and increases, ATE+ increases; similarly, when $P<0.5$ and decreases, ATE− increases. One the other hand, TE does not distinguish the types of influence between the two time series. In this case, the expectation values for ATE+ and ATE− can be derived as functions of P: $\text{ATE}+=P{log}_2(P/0.5)$ and $\text{ATE}-=(1-P){log}_2[(1-P)/0.5]$. These are consistent with the numerical results shown in the figure. To understand more about ATE+ and ATE−, we show two examples in which the influence on y from itself is positive, instead of a neutral influence in the case of Figure 2(a). In Figure 2(b), $p(\dot{y}(t+1)=a|\dot{y}(t)=a)=0.6$, while in Figure 2(c), $p(\dot{y}(t+1)=a|\dot{y}(t)=a)=0.9$, which has a stronger positive influence. Near $P=1$, one can observe that ATE+ in Figure 2(c) is much smaller than ATE+ in Figure 2(b), even though the positive influence from $\dot{x}(t)\to \dot{y}(t+1)$ is high. The explanation of this is that TE and ATE± negate its own influence, $\dot{y}(t)\to \dot{y}(t+1)$, and the self-influence in Figure 2(c) is 0.9, much larger than that in Figure 2(b), which is 0.6. In other words, TE and ATE± measures the net influence, which is the amount of external influence minus the internal influence.

Figure 2

ATE distinguishes positive and negative influence in information transfer. The TE, ATE+ and ATE− curves are functions of $P=p(\dot{y}(t+1)=0|\dot{x}(t)=0)=p(\dot{y}(t+1)=1|\dot{x}(t)=1)$.

2.3 ATE symbolization technique

The definition of TE and ATE in Eq. (1) and Eq. (2) respectively are based on continuous random variables. To estimate probability distributions from a finite number of observations, one can discretize the observed values. We adopt the symbolic transfer entropy method [13] that is based on the symbolization technique for permutation entropy [16]. We adopt this method to compute ATE because ATE sums over a subset of all the possible discrete states that TE sums over. Therefore, this is a natural extension of the computation. The key step is to transform the continuous-valued time series ${\{x(t)\}}_{t=1}^N$ into a symbolic time series with a suitable length n, where $n\ge 2$ is an integer. First, for each t, the n consecutive values $\{x(t),x(t+1),\dots ,x(t+n-1)\}$ are ordered in the ascending order. The sequence of the indices that correspond to such permutation is recorded and denoted by $\hat{x}(t)$, the associated permutation of order-n symbol of $x(t)$ at t. Then we estimate the ATE of ${\{x(t)\}}_{t=1}^N$ by calculating the ATE of ${\{\hat{x}(t)\}}_{t=1}^N$. This now can be done efficiently because the number of states to sum over in the ATE equation (2) is finite and the probabilities and conditional probabilities can be directly calculated based on the discrete-valued symbolic time series ${\{\hat{x}(t)\}}_{t=1}^N$.

2.4 Local TE and ATE

In complex systems, the amount of information exchange from one node to another is typically not constant. To properly handle dynamic data, we calculate the TE and ATE of time series in a local time window, so that TE and ATE become functions of time. For illustration, we consider a two-node network, which produces two time series $x(t)$ and $y(t)$, for $t=1,\dots ,1000$, with symbolic states $n=2$ for $\dot{x}(t)$ and $\dot{y}(t)$, where $\dot{x}(t)=1$ represents increment and $\dot{x}(t)=0$ represents decrement. We then simulate the data according to the conditional probabilities $p(\dot{y}(t+1)=0|\dot{x}(t)=0)=p(\dot{y}(t+1)=1|\dot{x}(t)=1)=P$ and therefore $p(\dot{y}(t+1)=1|\dot{x}(t)=0)=p(\dot{y}(t+1)=0|\dot{x}(t)=1)=1-P$. In the simulated data, the probabilities of associative influence change at $t=300$ and 600. For $1\le t<300$, $P=[0.9,1]$; for $300\le t<600$, $P=[0,0.2]$; and for $600\le t\le 1000$, $P=[0.9,1]$. Therefore, initially x has a strong positive influence, then a strong negative influence during the middle period, and finally a strong positive influence at the end. For simplicity, we fix $p(\dot{x}(t+1)=0|\dot{y}(t)=0)=p(\dot{x}(t+1)=1|\dot{y}(t)=1)=p(\dot{x}(t+1)=0|\dot{y}(t)=1)=p(\dot{x}(t+1)=1|\dot{y}(t)=0)=0.5$. The simulated data can be seen in the top graph of Figure 3. The local TE and ATE± curves are plotted in the bottom graph, with sliding window size $W=100$. One can see that initially ATE+ has large values, then ATE− has large values, and in the final stage ATE+ has large values. This is consistent with our setup of the simulated data. Near the switching of information flow at $t=300$ and $t=600$, the ATE± curves are not step functions, but instead have linear transitions. The reason is because the calculation of local ATE relies on a sliding window, which moves continuously and can include different dynamics within the window. Nonetheless, we observe the effectiveness of local ATE, compared to local TE, which does not distinguish between positive and negative influence.

Figure 3

TE and ATE as functions of time in the dynamic case. Top: Simulated time series with $p(\dot{x}(t+1)|\dot{x}(t),\dot{y}(t))$ fixed and unbiased, while $p(\dot{y}(t+1)|\dot{y}(t),\dot{x}(t))$ is switched at time $t=300$ and 600, indicated by the blue vertical lines. Bottom: The ATE curves are able to distinguish positive and negative influence, while the TE curve cannot.

2.5 Spectral radius of TEM and ATEM

Since TE and ATE are pairwise measurements, they do not give a single measurement of the total amount of entropy transferred within a system. To extend pairwise TE and ATE, we first define the ATE Matrix (ATEM) of a system of m elements as an $m\times m$ matrix with ij th entry ${(M^S)}_{ij}={ATE}_{x_i\to x_j}^S$, where the superscript S is the associative state that is defined for Eq. (2). Similarly, the ij th entry of the $m\times m$ TE Matrix (TEM) is ${(M)}_{ij}={TE}_{x_i\to x_j}$. TEM has been used in [14] to reveal the asymmetric influences from mature markets to emerging markets, but was not used to give the total amount of TE of the entire system, nor in a sliding window fashion to measure the transient influence dynamics. Note that since TE and ATE are directional, TEM and ATEM are nonsymmetric, and thus their eigenvalues are complex-valued.

We use the spectral radius of the TEM (resp. ATEM) to measure the total amount of TE (resp. ATE) in the entire network. The spectral radius of a matrix is the supremum of the absolute values of its eigenvalues. Furthermore, to use the spectral radius of a TEM or ATEM to monitor the system, we use a time sliding window to calculate local TEM or ATEM and therefore the spectral radius of a TEM or ATEM becomes a function that depends on time.

To demonstrate the usefulness of our local TE method, we show that for a variety of pitchfork bifurcations, the spectral radius function of the local TEM is able to detect the transition in advance. Figure 4 shows four different graph structures: an undirected 9-node chain, a directed 9-node chain, a downward 7-node binary tree, and an upward 7-node binary tree. We simulate supercritical pitchfork bifurcations on these network graphs and examine the simulated data with the spectral radius of each graph’s TEM. The data $x(t)$ is simulated according to the following ordinary differential equation [11] on a fixed network:

$$\frac{\partial x}{\partial t}=c(t)x-x^3+\alpha △x+\sigma d\omega ,$$

(3)

where △ represents the graph Laplacian of the network, $\alpha >0$ is a real number that scales the amount of diffusion across the nodes, and $\sigma d\omega$ is white noise with scale σ. The scaling function in our data simulation is chosen to be $c(t)=tanh(3t-10)$. We chose this pitchfork bifurcation equation because it generates data that will transition from one state to another rather quickly for individual nodes. With the graph Laplacian term in the equation, as an external force for each node, one can get different interesting bifurcation results. Our purpose of using this equation is to show that for a variety of pitchfork bifurcations, our local TE method is able to detect the transition in advance. For instance, the adjacency matrix A and Laplacian matrix L of the upward binary tree are:

$$A=\left[\begin{array}{ccccccc}0& 1& 1& 0& 0& 0& 0\\ 0& 0& 0& 1& 1& 0& 0\\ 0& 0& 0& 0& 0& 1& 1\\ 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0\end{array}\right]\text{and}L=\left[\begin{array}{ccccccc}2& -1& -1& 0& 0& 0& 0\\ 0& 2& 0& -1& -1& 0& 0\\ 0& 0& 2& 0& 0& -1& -1\\ 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0\end{array}\right],$$

(4)

respectively. In the adjacency matrix, $a_{ij}=1$ means there is a connection from node j to node i, and $a_{ij}=0$ means there is not a connection from node j to node i. The Laplacian matrix $L=D_{\mathrm{in}}-A$, where $D_{\mathrm{in}}$ is the in-degree matrix whose diagonal entry $d_{ii}$ is the sum of i th row of A, the number of connections going into node i. Figure 4 plots the simulated time series and their spectral radii of the local TEM. We observe that before transitioning or bifurcation, the spectral radius decreases rapidly, which provides an early indication of system transitioning. The spectral radius of TEM in these cases drops to the lowest point during transitioning due to the strong internal dynamics of individual nodes, and TE measures the net influence dynamics from one node to another.

Figure 4

Spectral radius of TEM as early indicator of phase transition. Simulated pitchfork bifurcations with different graph structures. The spectral radius of TEM decreases before bifurcation, which provides an early indication for phase transition.

2.6 Detect critical transitions with ATE and TE spectral radius curves

In this subsection, we explain how to use the TE and ATE curves defined in the previous subsection to characterize critical transitions prior to their occurrences (see Table 1). We give our observations on behaviors of TE and ATE in four different critical transition scenarios. The first critical transition scenario is from a stable state to oscillations in a network of non-Foster circuits. In this case, as a network starts to synchronize in phase, the pairwise TEs will likely start to increase because of the growing amount of interacting dynamics between circuits. The TE and ATE+ spectral radius curves may therefore start to increase. If the synchronization is in the negative mode, in which the phases differ in 180^∘, the ATE− spectral radius curve will start to increase.

The second scenario is the case where the time series of a system bifurcate, switching from one state to alternative states. The examples in Figure 4 demonstrate the nature of such transitions. Since the phenomenon of critical slowing down is observed prior to such transitions, the TE and ATE curves may decrease rapidly prior to the tipping point.

The third scenario is when there is a sudden drop in the time series of the system. It has been suggested that there is a critical slowing phenomenon before the tipping point of a crash in many real world systems, such as stock markets, housing market, etc. In this scenario, right before the crash, the interacting dynamics between nodes may start to decrease, which can be explained by the uncertainty of the information carried by other nodes. Therefore, the TE and ATE curves are likely to decrease rapidly before the crash.

The fourth scenario is the exponential growth in system activities. In this case, the TE and ATE+ curves are likely to increase, as a result of an increase in the interacting dynamics between nodes. Given one of the above scenarios, we can look out for critical transitions using TE and ATE curves.

3 Numerical results of information dynamic spectrum

3.1 ATE indicates instability of non-Foster circuit

3.1.1 Data

The circuit data used for TE and ATE evaluation is simulated from a non-Foster network described in [17]. A non-Foster network is an active network containing a power source, therefore it may not obey Foster’s reactance theorem. The non-Foster circuit data analyzed here contains three coupled circuits with internal feedback loops. The interface points include two antennas (Ant1 and Ant2) and two ports (i1 and i2) to record time series of the voltage. The circuit is operated in the stable (no oscillation) region, but near the unstable region. A perturbation is then added at timestep = 250 with 1 milli-Amp. When a small perturbation is added, the circuits may become oscillatory. The circuits will either become synchronized in oscillation, which is considered unstable, or return to the original stable state.

3.1.2 Analysis

Figure 5 shows our ATE analysis of the non-Foster network. At the top is a plot of the circuits in voltage over time at each interface point. The circuits are initially operated in the stable region, where there are no oscillations. The circuits become unstable, after a small perturbation is added. We perform the ATE analysis to determine whether the circuits will become unstable. The bottom plot shows the TE and ATE± curves over time. The curves are obtained from the absolute sum of spectrum of TEM and ATEM±, respectively. We found that the spectral radii of TEM and ATEM± also show similar trajectories, but since TEM is the sum of ATEM+ and ATEM−, the ATE+ spectral radius curve will never cross the TE spectral radius curve. Thus, the absolute sum of the spectrum in this example is more informative as an indicator of early detection of synchronizations. As shown in the plot, the ATE+ spectral absolute sum curve crosses over that of TE around 800, indicating that a significant positive influence is being transferred. Therefore, the system is likely to synchronize in the positive mode. On the other hand, if the negative influence had crossed the overall information transfer, the system is likely to sync in the negative mode. An example of this is illustrated in Additional file 4. As an early indication of critical transitions, the TE and ATE+ curves in Figure 5 increase rapidly before the full synchronization. Interestingly, TE and ATE curves reach their peaks at the start of the synchronization, while the ATE− curve flattens because there is no negative association. In addition, to show how the measures behave when there are no critical transitions, an example is included in Additional file 4 when there is no perturbation applied to the circuits. The parameter choices for this analysis are symbolic length $n=3$ and sliding window size $W=30$.

Figure 5

Spectral radii of TEM and ATEM as early indicators of synchronization. Top: Signals from non-Foster circuits in voltage. Bottom: The corresponding TE and ATE curves.

3.2 TEMs infer directional influences in Latin America stock indices

3.2.1 Data

Table 2 shows nine major Latin America stock market indices. Time series of these stock market indices were collected during a year span, centered at the October 2008 Crash, with the resolution of a day, except when the market is closed during the weekend. During the month of October 2008, the stock indices dropped dramatically. Figure 6 (top) shows the nine Latin America stock indices. Since the volumes of the nine stock indices are very diverse, we normalized the indices in order to show the drastic drop of the market in October 2008.

Figure 6

Spectral radii of ATEM as early indicators of the October 2008 Crash. Top: Latin America stock indices around the October 2008 Crash. Bottom: TE, ATE+, and ATE− curves decrease rapidly before the October 2008 Crash as early indicators.

Table 2 Latin America stock indices

3.2.2 Analysis

The bottom plot of Figure 6 shows the TE, ATE+, and ATE− curves. One can observe the rapid decrease of the curves before the October 2008 Crash, that is similar to the pitchfork bifurcation examples in Figure 4. Interestingly, two other decreases of the curves precede the October 2008 crash and two other bumps after the crash are not directly associated to any significant transition. They might be due to the general instability of the markets after these two major events. Additionally, there is another rapid decrease in the TE and ATE curves in February 2009, which corresponds to the market fall in March 2009. This market fall is due to that in March 2009, the stock in USA fell to its lowest level since 1997 for Dow and 1996 for S&P 500 [18]. We use TEMs of different periods to analyze the dynamics of the stock indices around a critical event. We identify the directional structures in different periods (see Figure 7), where each node represents a Latin America stock market index. The top row of Figure 7 shows the TEMs before, during, and after the October 2008 Crash from left to right. The red color corresponds to large values, while blue corresponds to small values. During the crash, the total amount of information transfer decreases. After the crash, the values of TEM return to the values before the crash. Visualization of the network structure shows that Panama is strongly influenced by Colombia and Brazil before and after the crash, yet during the crash Panama is primarily driven by Mexico and Venezuela. It is believed that [19] Panama’s economy has strong ties to shipping needs of Brazil (iron ore and soy from the north and Amazonia) and Colombia (particularly coal) that use the Panama Canal [20]. As for Mexico and Venezuela, there is thought that the general economic collapse of 2008 is attributed to the downturn of Panama’s real estate and offshore banking markets that were dominated by Mexican and Venezuelan investors associated with illicit trafficking and greatly affected by the drying up of funds from the US and European drug markets [21]. The parameter choices for this analysis are symbolic length $n=5$ and sliding window size $W=30$.

Figure 7

TEM analysis of the October 2008 Crash. TEMs infer directional influences in Latin America stock indices before, during, and after the October 2008 Crash.

3.3 ATE reveals Wikipedia motifs that drive the changes

3.3.1 Data

The dynamics of editing behaviors of Wikipedia’s content is explored in [22], in which temporal motifs, temporal bipartite graphs with multiple node and edge types for users and revisions, are proposed. The first two rows of Figure 8 show the twelve most frequent motifs from the years 2001 to 2011. For example, in the second row, the first motif at the left shows a minor edit of a Wikipedia page by an anonymous author, followed by a revert from a registered user, followed by a minor edit from an anonymous author. The time series of each motif is the number of times that motif occurs in each month. The top of Figure 9 plots the counts of the Wiki motifs for each month from January 2001 to April 2011. This example shows a different type of transition: an exponential growth in Wikipedia editing.

Figure 8

ATE reveals Wikipedia motifs that drive the growth of Wikipedia pages. Most frequent motifs of Wikipedia from left to right and first row to second row. (a) ATEM+ shows motifs 1, 8 and 12 have the most positive influence on other motifs. (b) ATEM− shows motif 9 has the most negative influence on other motifs. (c) TEM shows asymmetric influence among the motifs but the contrast is not as strong as ATEM+ and ATEM−.

Figure 9

Information dynamics of Wikipedia motifs. Top: Wikipedia motif occurrence in time. Bottom: TE, ATE+, and ATE− curves in blue, red, and green, respectively. There is a significant increase in the ATE+ curve in 2005, as an early indication of rapid growth in Wikipedia’s contents.

3.3.2 Analysis

The bottom row of Figure 8 from left to right shows global ATEM+, ATEM−, and TEM from 2001 to 2011 to identify important motifs that drive the changes in other motifs. In particular, ATEM+ shows that motifs 1, 8 and 12 have the most positive influence on other motifs, observing that consecutive ‘minor add’, ‘revert’, and ‘major add’ by registered users encourage Wikipedia’s content growth. ATEM− shows that motif 9 has the most negative influence on other motifs. TEM shows the asymmetric influence among the motifs but the contrast is not as strong as ATEM+ and ATEM−. The bottom of Figure 9 shows the TE, ATE+, and ATE−, i.e. spectral radii of local TEM and ATEM, in blue, red, and green, respectively. There is a significant increase in the ATE+ and TE curves near midyear 2005 and ATE+ crosses over ATE−. These suggest an early indication of the growth in Wikipedia’s contents. Interestingly, the peaks of TE and ATE+ in late 2008 occur at the beginning of the stabilization of the editing behaviors. The parameter choices for this analysis are symbolic length $n=5$ and sliding window size $W=30$.

3.4 TE and ATE curves as early indicators of stock market crash

3.4.1 Data

We collected the share prices of Dow Jones Industrial Average (DJIA) from June 1 to December 31, 2008, using 25 shares that persisted throughout the years near the October 2008 Crash. The stock indices are: MMM, AA, AXP, T, BA, CAT, KO, DD, XOM, GE, HPQ, HD, INTC, IBM, JNJ, JPM, MCD, MRK, MSFT, PFE, PG, UTX, VZ, WMT, and DIS. The time series resolution is one day, except when the market is closed during the weekend.

3.4.2 Analysis

Figure 10 shows the Dow Jones indices (top) and the TE and ATE curves (bottom) from June 2008 to December 2008. The indices are normalized and translated in the index value for visualization. During the October Crash, all the indices decrease rapidly. We see that the TE and ATE spectral radius curves behave similarly to the ones in Figure 4 before bifurcations and Figure 6 before the Latin America market crash. The rapid decrease of TE suggests a critical transition, where the dynamics within individual nodes become more significant than the dynamics between nodes. An additional observation is that the TE and ATE curves decrease quickly a second time after October 6, right before the second drop across all the indices. The parameter choices for this analysis are symbolic length $n=5$ and sliding window size $W=20$.

Figure 10

TE and ATE curves as early indicators of stock market crash. Top: Dow Jones indices from June 2008 to December 2008. Bottom: TE, ATE+, and ATE− curves in blue, red, and green, respectively. There are significant decreases in all three curves before the October 2008 Crash, as an early indicator of market crash.

3.5 Discussions

Our empirical data analysis suggests that early detection of critical transitions using local ATE curves is possible. The local ATE curve, spectral radius of the local ATEM over a sliding time window, is an indicator of the amount of net entropy transferred in the entire system at each time. On the other hand, by looking at the entries of ATEM, one can identify the most positive or negative influential nodes of the system, as well as the most positively or negatively influenced nodes. From the observations of the TE, ATE+, and ATE− curves in different examples that corresponds to different types of critical transitions, we found that (1) ATE+ increases rapidly before synchronization in the positive mode. ATE− increases rapidly before synchronization in the negative mode. (2) TE, ATE+, and ATE− drop rapidly before abrupt falls, because TE, ATE+, and ATE− are the net information flow and external information flow begins to decrease before the sudden collapse. (3) ATE+ increases rapidly before an exponential growth in a system.

4 Probabilistic light cones for trajectory prediction

We propose a model-based probabilistic light cone method for predicting the trajectories of the spectral radius of TEM or ATEM. As observed in the previous section, the behaviors of TE and ATE curves can be interpreted manually as early indications of critical transitions. Here, we provide an automated method to generate early indicators. Using a natural logarithmic model, our method first computes probabilistic light cones for TE and ATE curves. This gives a prediction value with a confidence interval. If the actual TE and ATE trajectories are outside the confidence interval, the method declares instability is detected, which serves as an early indicator of a critical transition.

Our computation of probabilistic light cone is inspired by [23], where an MCMC method is used to predict citation growth based on the preferential attachment. However, critical transition is not discussed in [23]. Our method, on the other hand, aims to generate early indicators for an upcoming critical transition as time progresses. We choose natural logarithmic curves to model the growth rate of the TE trajectory as information transfer approaches its maximum prior to a critical transition as observed in the non-Foster circuit example in the previous section. We apply a moving time window over the observed spectral radius (TE or ATE±) time series to derive the unknown coefficients and constants for natural logarithm curves. For a given prediction time point, we generate the probabilistic light cone based on 95% confidence intervals of predicted trajectories with fitted natural logarithm curves.

4.1 Model-based forecasting

We propose a statistical forecasting model to estimate the TE and ATE trajectories. An advantage of analyzing the TE and ATE curves is that TE and ATE remove the spikes of the raw data. We have observed that the TE curve can be well approximated by the natural logarithm with an unknown coefficient a and a constant c:

$$g(t)=aln(t)+c.$$

(5)

An example of this is shown in Figure 11. We see that this logarithmic function is able to approximate the TE curve of the non-Foster circuits data very well. The curve fitting is done on the interval from $t=700$ to $t=1100$ with the exponential function and logarithmic function in solid blue and solid red lines, respectively. The fitted curve is projected to a future time interval from $t=1100$ to $t=1700$ in dashed blue and dashed red lines, respectively. Our observations of the TE and ATE curves in other types of data also show that the curve approaches its maximum gradually, instead of a spiking or sudden approach. The idea here is to transform the original raw signals to better-behaved features (signals) for detection and prediction. Concave increases of curves are very common, such as autocorrelation, variances, skewness, flickering, state-transition-based local network entropy (SNE) [24], and information dissipation length (IDL) [25]. However, TE is a more general measure because it is an inherently non-symmetric measure that quantifies the net entropy transferred with appealing concave increase characteristics.

Figure 11

TE curves fitted with logarithm functions. The ATE curve (black) can be fitted quite well with the logarithm function (blue), in contrast to fitting with the exponential function.

To estimate the coefficient a of the logarithmic function in (5), instead of fitting the TE and ATE curves deterministically, we estimate the rate of change with various discrete time steps. This generates a probabilistic light cone, which consists of a collection of fitted curves using different discrete time steps. Taking the derivative of Eq. (5), we have $g^{\prime }(t)=\frac{a}{t}$. We obtain the following discretized version for the unknown a:

$$\frac{△g}{△t}=a\frac{1}{t}.$$

(6)

For a fixed timestep △t in a window $[T_{\mathrm{start}},T_{\mathrm{end}}]$, we use the least squares method to solve the unknown a. First, we write the following matrix equation:

$$\left[\begin{array}{c}\frac{1}{t_1+△t/2}\\ \frac{1}{t_2+△t/2}\\ ⋮\\ \frac{1}{t_k+△t/2}\end{array}\right]a=\left[\begin{array}{c}\frac{g(t_1+△t)-g(t_1)}{△t}\\ \frac{g(t_2+△t)-g(t_2)}{△t}\\ ⋮\\ \frac{g(t_k+△t)-g(t_k)}{△t}\end{array}\right],$$

(7)

where $t_1,\dots ,t_k+△t\in [T_{\mathrm{start}},T_{\mathrm{end}}]$. The approximation of $a_{△t}$ for the fixed timestep △t is then obtained by

$$a_{△t}={argmin}_a{\parallel \left[\begin{array}{c}\frac{1}{t_1+△t/2}\\ \frac{1}{t_2+△t/2}\\ ⋮\\ \frac{1}{t_k+△t/2}\end{array}\right]a-\left[\begin{array}{c}\frac{g(t_1-△t)-g(t_1)}{△t}\\ \frac{g(t_2-△t)-g(t_2)}{△t}\\ ⋮\\ \frac{g(t_k-△t)-g(t_k)}{△t}\end{array}\right]\parallel }_2.$$

(8)

The constant corresponding to this timestep is then $c_{△t}=g(T_{\mathrm{end}})-a_{△t}ln(T_{\mathrm{end}})$. Therefore, the predicted value for the future time $T_{\mathrm{end}}+t_d$ is

$$F_{△t}(T_{\mathrm{end}}+t_d)=a_{△t}ln(T_{\mathrm{end}}+t_d)+c_{△t}.$$

(9)

4.2 Probabilistic light cone and error estimation

Given the method to estimate the constants c and a in previous section, we generate a probabilistic light cone at each time as we vary △t to obtain multiple estimates of c and a. We define a probabilistic light cone at a given time consists of a collection of natural logarithm curves starting from that point.

To estimate the error of this prediction for the immediate next timestep, we calculate the following:

$${error}_{△t}(i)=g(t_i+△t)-[a_{△t}ln(t_i+△t)+c(i)],$$

(10)

where $c(i)=g(t_{i-1})-a_{△t}ln(t_{i-1})$.

Let E be the collection of all ${error}_{△t}$ over all timesteps △t and let $\sigma =\text{ standard deviation}\text{of }E$. Let G be the collection of all predicted value $G_{△t}(t_d)$ over all timesteps △t and $\mu =mean(G)$. Therefore, the 95% confidence interval for the future time $T_{\mathrm{end}}+t_d$ is

$$CI=[\mu -1.96\frac{\sigma }{\sqrt{L}},\mu +1.96\frac{\sigma }{\sqrt{L}}],$$

(11)

where L is the size of E.

4.3 Numerical results

Figure 12 shows the probabilistic light cone analysis on the non-Foster circuit data. The top plot shows prediction of the TE curve (blue) in order to find instability. The predicted values and their 95% confidence intervals are shown in green. In this probabilistic light cone analysis, the prediction leap time is $t_d=20$, the discrete timestep collection is $△t=1,\dots ,50$, and the curve fitting window size is $T_{\mathrm{start}}-T_{\mathrm{end}}=200$. Two (indicated in magenta diamonds) of the 36 points of the actual trajectory are outside the 95% confidence interval, just before the non-Foster circuits are fully synchronized around time $t=1500$. These two points indicate the instability of the TE curve and detect the full synchronization in advance. Figure 13 shows two snapshots of the probabilistic light cone video (see Additional file 1) produced according to the method described above and the predicted trajectory.

Figure 12

Trajectory prediction and error estimation of non-Foster circuit TE curve. Top: Actual TE trajectory in dark blue and the predicted trajectory and its 95% confidence interval in green. The predicted trajectory stays close to the actual trajectory. The magenta diamonds indicate that the actual trajectory points appear outside the 95% confidence interval of the prediction, which is used as early detection of synchronization, or circuits entering the unstable region. Bottom: non-Foster circuits in voltage.

Figure 13

Probabilistic light cone. Two snapshots of the probabilistic light cone and the predicted trajectory. The top panel of each snapshot shows the probabilistic light cone in red, the actual TE trajectory in blue and the predicted trajectory in green. The bottom panel of each snapshot shows the non-Foster circuit signals.

To see how far ahead in time $t_d$ one can predict the value in relation to error, the distribution of errors for each time $t_d$ is plotted in Figure 14. As the prediction time leap $t_d$ increases, the distribution of errors starts to flatten out, because the error increases. This is shown in the bottom right plot of Figure 14, where the curve is the mean squared error of the prediction values versus time leap $t_d$. The prediction method is designed to work for different values of $t_d$ within a reasonable range, since it is based on confidence values.

Figure 14

Error distribution for each prediction leap time step. As expected, a smaller prediction leap time step $t_d$ gives more accurate prediction, and as $t_d$ increases, the error distribution flattens out. We see that the mean squared error of prediction values increases as prediction timestep increases.

Figure 15 shows the probabilistic light cone analysis of the pitchfork bifurcations in Figure 4, where the leap time $t_d=10$ and the discrete timestep collection is $△t=1,\dots ,40$, and the curve fitting window size is $T_{\mathrm{start}}-T_{\mathrm{end}}=45$. The detected abrupt falls of the TE curve (red points) occur before the start time of the bifurcations. However, a couple drops that occur far before the transitions were also detected and gave false alarms.

Figure 15

TE trajectory prediction of pitchfork bifurcations. TE curves and the predicted trajectories of pitchfork bifurcations in Figure 4.

Figure 16 shows the probabilistic light cone analysis of the Latin America stock indices that is shown in Figure 6. The abrupt falls of the TE curve is detected and indicated in red, which occurs before the October 2008 Crash. The leap time is $t_d=5$ days and the discrete timestep collection is $△t=1,\dots ,20$, and the curve fitting window size is $T_{\mathrm{start}}-T_{\mathrm{end}}=25$.

Figure 16

Trajectory prediction and error estimation of Latin America stock index TE curve. Top: Actual TE trajectory in dark blue and the predicted trajectory and its 95% confidence interval in green. The magenta diamonds indicate that the actual trajectory points are outside the 95% confidence interval, as an early detection of the crash of the October 2008 crash. Bottom: Latin America stock indices.

Figure 17 shows the probabilistic light cone analysis of the Wikipedia motifs. The probabilistic light cone video is included in Additional file 2. At the top is the same TE curve (blue) of the Wikipedia motifs in Figure 8. The predicted values and their 95% confidence intervals are shown in green. In this probabilistic light cone analysis, the prediction leap time is $t_d=4$ months, the discrete timestep collection is $△t=1,\dots ,15$, and the curve fitting window size is $T_{\mathrm{start}}-T_{\mathrm{end}}=18$. The two magenta diamonds indicate that the actual TE values at those times are outside the prediction interval. These two points are above the respective prediction intervals, indicating rapid increase in TE or the instability of the TE curve. These can be used as early indications before the Wikipedia contents take off.

Figure 17

Trajectory prediction and error estimation of Wiki motif TE curve. Top: Actual TE trajectory in dark blue and the predicted trajectory and its 95% confidence interval in green. The magenta diamonds indicate that the actual trajectory points are outside the 95% confidence interval, as an early detection of the growth of the Wikipedia content. Bottom: Occurrence of the most frequent Wiki motifs.

Figure 18 shows the probabilistic light cone analysis of the Dow Jones data. A video is included in Additional file 3. At the top is the same TE curve (blue) of the Dow Jones stock indices in Figure 10. The predicted values and their 95% confidence intervals are shown in green. In this probabilistic light cone analysis, the prediction leap time is $t_d=5$ days, the discrete timestep collection is $△t=1,\dots ,15$, and the curve fitting window size is $T_{\mathrm{start}}-T_{\mathrm{end}}=18$. The two magenta diamonds show that the actual TE at those times is outside the prediction interval, indicating the instability of the TE curve. The magenta diamond on the left occurs at the bottom of the sudden drop of the TE curve, before the October crash, which suggests the probabilistic light cone method is able to early detect the market crash. The second magenta indicates a rapid increase in the TE. Since the objective here is to detect crash, we only look for a sudden drop in the TE curve. The phenomenon of a sudden increase in the TE curve in stock market will be left for future work. On the other hand, the second drop in TE that is observed in Figure 10 is not as steep as the first drop and was not caught by the probabilistic light cone method. We speculate that large-scale fluctuations after the October 2008 Crash may indicate that the system were on its paths to another region. Consequently, one may need to change parameters in our probabilistic light cone method (e.g. sliding time windows) to adapt toward a new system region.

Figure 18

Trajectory prediction and error estimation of Dow Jones stock market TE curve. Top: Actual TE trajectory in dark blue and the predicted trajectory and its 95% confidence interval in green. The magenta diamonds indicate that the actual trajectory points are outside the 95% confidence interval. The magenta diamond on the left indicates the detection of a critical transition. Bottom: Dow Jones stock market data.

5 Conclusions

We propose a novel information dynamic spectrum framework for automated detection of critical transitions and identification of directional influences. We have shown that the framework is able to: (1) provide an effective measure for quantifying associative, asymmetric directional influence rather than symmetric influence, (2) provide an effective formulation that captures system-wise directional influence, rather than pair-wise influence and (3) provide an effective measure for detecting instability in systems with directional influence dynamics. Within this framework, we further propose the probabilistic light cone method that predicts the ATE trajectories and indicates instabilities toward critical transitions. This model-based prediction method captures information that cannot be discerned by looking at individual time series signals alone and is novel compared to conventional deterministic methods or complicated probability state space predictions. The analyses performed in this work are based on historical data, and the rapid decrease or increase in TE is not a sufficient condition for all critical transitions. Rather, as shown in Table 1, the nature of upcoming critical transitions is revealed by the behaviors of ATE curves, in conjunction with examining the trends in the time series of a system. Future work requires a complete analysis of the proposed framework, as well as an adaptive parameter selection method for the probabilistic light cone analysis.