Spatiotemporal dynamics analysis and systemic risk measurement of energy price system based on complex network

Physica A 526 (2019) 120702

Contents lists available at ScienceDirect

Physica A journal homepage: www.elsevier.com/locate/physa

Spatiotemporal dynamics analysis and systemic risk measurement of energy price system based on complex network ∗

Hua Xu a,c , Minggang Wang b,c , , Weiguo Yang a a

Faculty of Science, Jiangsu University, Zhenjiang, 212013, Jiangsu, China School of Mathematical Science, Nanjing Normal University, Nanjing, 210042, Jiangsu, China c Taizhou College, Nanjing Normal University, Taizhou 225300, Jiangsu, China b

highlights • • • •

We propose a method to characterize the spatiotemporal dynamics of energy price system. We propose a new method to measure the systemic risk based on complex network. The ‘‘risk hopping" phenomenon of energy price system is revealed. The topological structure of energy price correlated modal transition network is revealed.

article

info

Article history: Received 11 October 2018 Received in revised form 16 January 2019 Available online 12 April 2019 Keywords: Energy price Correlation mode Systemic risk Complex network

a b s t r a c t Multivariate time series contain rich information of corresponding complex systems. This paper proposes a new method to characterize the evolution and transition characteristics of the correlation modes for multivariate time series(ETCCMMTS). Based on ETCCMMTS, the characteristics of energy price systems are analyzed. The results show that energy price correlated modes exhibit co-movement, group-occurrence and the small-world property. The results prove that futures prices have the leading role, the crude oil futures price affects the stability of an energy price system, and the crude oil spot price has the highest clustering coefficient, whereas the heating oil futures price has the largest betweenness coefficient. The topological absorptive is introduced to measure the systemic risk, and the "risk hopping" phenomenon is revealed. The recursive graph indicates that the fluctuations of energy price display not only short-range correlation for 1–2 months but also long-range correlation for 7–8 years. Additionally, the long-range correlation will lead to high systematic risk. In the course of energy price fluctuations, there is a total of 41 different types of correlation modes. The cumulative time of different correlation modes may result in a ‘‘cumulative time jump" and shows a linear growth in general. In addition, the topological structure of energy price correlated modal transition network (EPCMTN) is revealed. © 2019 Elsevier B.V. All rights reserved.

1. Introduction Time series can be regarded as the output of a system. Its rich information can reflect the microstructure of the system. Through mining and analyzing data from the time series, we can obtain much valuable information on the system in ∗ Correspondence to: School of Mathematical Science, Nanjing Normal University, Nanjing, 210042, Jiangsu, China. E-mail address: [email protected] (M. Wang). https://doi.org/10.1016/j.physa.2019.03.067 0378-4371/© 2019 Elsevier B.V. All rights reserved.

2

H. Xu, M. Wang and W. Yang / Physica A 526 (2019) 120702

question. Complex networks have been a research hotspot in recent years, and their main idea is to take the relation between parts of a real system as a complex network and then characterize the relationship between those parts in the form of network to better understand the essential characteristics of real systems. With the advent of BA scale-free networks [1], WS small-world networks [2], NW networks [3], ER random networks [4] and other complex networks, complex networks have been applied in more fields, which provides us with a new perspective and a new method for studying complexity problems. In recent years, the application of complex network theories to the analysis of nonlinear time series has attracted much attention. The main idea of such a method is to transform nonlinear time series into the corresponding complex networks with certain algorithm(s) and then use the topological structure of complex networks to analyze the properties of nonlinear time series. At present, there are generally four ways to convert univariate time series into complex networks. (1) Pseudo-periodic time series transitions [5], which are mainly targeted at analyzing pseudo-periodic time series. Their main idea is to map each period of a time series into one node. If the phase space distance or correlation coefficient of two periodic fragments satisfies certain conditions, the corresponding nodes of the two sequence segments are connected, and a network is constructed. (2) The visibility graph (VG) method, which was first put forward by Lucas Lacasa et al. [6]. The main idea of this method is to map all of the data in the time series onto the nodes. If all nodes between two nodes fall below the line that connects the two nodes, that is, the two nodes are ‘‘visible’’ then the connection between the two nodes is set up, thus forming a complex network. The advantages of this network are that it can be used to analyze all types of time series and maintain most of the properties of the original time series. Through this method, periodic sequences, random sequences and fractal sequences are transformed into regular networks, random networks and scale-free networks, respectively. Recently, many different types’ visibility algorithms are proposed and have been successfully applied to many fields, such as human behavior analysis, stock index analysis, foreign exchange rate, and gold price time series [7–14]. (3) The phase space reconstruction method [15,16]. It starts from the phase space reconstruction of time series analysis, maps fixed-length time series segments into nodes of a network, and then uses the correlation coefficients (or distances) between these nodes to determine whether these nodes are connected or not. If the absolute value (or distance) of one correlation coefficient is larger than a certain threshold value, it is considered that the two nodes are connected. This method is also widely used at present [17–19]. (4) The coarse graining method [20–23]. In this method, the fluctuations of time series are transformed into specific signal sequences. A fixed-length signal sequence is regarded as a network node to connect nodes of time series in chronological order and to construct a weighted complex network with direction. The method is widely used in international oil trade series, energy price series and consumer price index series [24–28]. In recent years, some scholars have studied high-dimensional time series with complex network theories [29,30]. This approach has received wide application [31–41]. In summary, research studies on time series with complex network theories have achieved good results, resulting in a number of new algorithms that can transform time series into complex networks. However, the application of complex network theories to the analysis of time series is still an open problem. Further study is still needed to develop more efficient methods to transform time series into complex networks. Current research studies mainly focus on analyzing univariate time series with complex networks [5–28]. Although there are some researches study on the analysis of multivariate time series with complex networks [29–38], those research studies are generally still at an early stage and defective with regard to solving practical problems. In [29], although the global dynamics of high-dimensional systems are analyzed, the study of the local evolution of high-dimensional systems is not sufficient. In [30], multiple complex networks are used to characterize multivariate time series. The complexity of its algorithm is higher, but the reliability of the characterization of the actual system decreases. To analyze multivariate time series with complex networks, this paper presents a method to convert multiple time series into complex networks. The basic idea is as follows: Firstly, we divide multiple time series by constructing sliding windows, obtain a series of data matrices that change continuously, and construct sequences of correlation modes based on a complex network by calculating the correlation of data matrices. Then, we define a series of dynamic indexes of correlation modes and analyze the local and global dynamics of the highdimensional system by using the evolutionary relation between correlation modes. Finally, we construct the transition network of correlation modes to find core correlation modes and analyze the rules of correlation modes transition. Compared with previous studies, our method has several advantages. (1) The analysis method proposed in this paper is simpler and more practical. (2) We can use stochastic matrix theory by constructing sliding windows and divide multivariate time series. These guarantee that we can obtain better continuity in segmented time series. (3) We define a dynamic topological index to characterize the evolution of correlation modes. The index can be used to not only characterize the local evolutionary features of multivariate time series but also measure system risks. (4) The method proposed in this paper can be used to study the correlation mode transition characteristics of multivariate time series at each small-scale time window, which makes up for the shortcomings of previous studies. This method can be used to analyze the correlation mode evolution of multiple energy price time series. Additionally, the effectiveness of our method has been verified by empirical analysis. 2. Methods In this paper, our method of analyzing the evolution and transition characteristics of the correlation modes for multivariate time series includes four steps:

H. Xu, M. Wang and W. Yang / Physica A 526 (2019) 120702

3

Fig. 1. Internal connections within correlation modes.

Step 1: Construct Correlation Modes. A multivariate time series is defined as X = xi t ′ , i = 1, 2, . . . , N, t ′ = 1, 2, . . . , M. In the formula, N represents the number of time series, which is greater than or equal to three in this paper. M represents the length of time series. The time series on the entire large scale can be divided into small-scale segments by using the sliding time windows. Assume that the length of the sliding window is L, whose value is usually greater than that of N according to random matrix theory. Let l be the step length between sliding time windows. To ensure that segments of the time ] [ Lsmall-scale + 1 small-scale time windows, series are continuous, we require that l < L. In this manner, we can obtain T = M − l where [] refers to the rounding function. For every small-scale time window, we use Eq. (1) to calculate the correlation.

{ ( )}

(t )

cij

[

⟨ ⟩] [ ⟨ ⟩] t t t (k) − Xi( ) · Xj( ) (k) − Xj( ) = √ ⟨ ⟩]2 √∑ [ ⟨ ⟩]2 ∑L [ (t ) t L (t ) (t ) · (k) − Xi( ) k=1 Xi k=1 Xj (k) − Xj ∑L

k=1

(t)

(t)

(t )

Xi

(t)

(t)

(t)

(t)

(1)

∑L

(t)

(t)

∑L

(t)

where cij ∈ C(t) , Xi (k) ∈ Xi , Xj (k) ∈ Xj , k = 1, 2, . . . , L, ⟨Xi ⟩ = k=1 Xi (k)/L, and ⟨Xj ⟩ = k=1 Xj (k)/L. As we know, Eq. (1) is Pearson correlation coefficient which can only measure the linear correlation among variables. In order to overcome the defect of Pearson correlation coefficient, many methods have been developed to measure the nonlinear correlation, for example detrended cross-correlation analysis (DCCA) [42], multifractal detrended cross-correlation analysis (MF-DCCA) [43] and so on. In this paper, we mainly use the network topology to analyze the characteristics of the system, for the sake of simplicity, we still use Pearson correlation coefficient. In future work, we will test the impact of using different correlation measures on the results. Thus, within each time window, we can obtain one correlation coefficient matrix C(t ) . Assigning appropriate threshold values to rc and using Eq. (2), we can convert C(t ) into the adjacency matrix A(t ) : (t )

aij

⎧ (t ) ⎪ ⎨ 1, cij ≥ rc = 0, −rc < cij(t ) < rc ⎪ ⎩ −1, c (t ) ≤ −r c ij

(2)

(t)

where aij ∈ A(t) . We can use the adjacency matrix A(t) in every time window t to obtain the correlation mode NC (t) . The internal connections within each correlation mode NC (t) are as follows: at time t, if time series i and j are positively (t) correlated, i.e., aij = 1, then there is a connected edge between nodes i andj in correlation modes, which are shown in (t)

Fig. 1(a). If time series i and j are negatively correlated, i.e., aij = −1, there is also a connected edge between nodes i (t) aij

andj in correlation modes, which are shown in Fig. 1(b). If = 0, there is no connected edge between nodes i and j in correlation modes. For example, if there are four nodes, A, B, C and D, in the correlation modes, where A and B, as well as B and C, are positively correlated, A and C are negatively correlated, and D has no connection with A, B and C, then, the internal connections within correlation modes are shown in Fig. 1(c). Step 2: Topological Structure Evolution and Systematic Risk Measurement. According to the topological structure of a complex network [38], we define the dynamic topological index of { index } (t)

the correlation mode sequence Nc = Nc

, where t = 1, 2, . . . , T . Because there are positive and negative connections (t)+

in the structures of correlation modes, we take the degree of the positive connection of node i at time t as ki . The (t)+ larger ki is, the greater the degree of synchronous change between node i and its neighboring nodes is. We take the (t)− (t)− degree of the negative connection of node i at time t as ki . The larger ki is, the greater degree of asynchronous (t) (t)+ − change between node i and its neighboring nodes is. Therefore, at moment t, node i’s degree ki = ki + k(t) , which i (t) means that the greater ki is, the more highly connected nodei is with its neighboring nodes, and this node is under greater influence of its neighboring nodes. The average positive connection degree, average negative connection degree − + and average degree of every node are set as ⟨k+ i ⟩, ⟨ki ⟩ and ⟨ki ⟩, respectively. Take Fig. 1(c) as an example: ⟨kA ⟩ = 1, − + − + − + − ⟨kA ⟩ = 1, ⟨kA ⟩ = 2, ⟨kB ⟩ = 2, ⟨kB ⟩ = 0, ⟨kB ⟩ = 2, ⟨kC ⟩ = 1, ⟨kC ⟩ = 1, ⟨kC ⟩ = 2, and ⟨kD ⟩ = ⟨kD ⟩ = ⟨kD ⟩ = 0. At time t, the average positive connection degree, average negative connection degree and average degree of correlation modes

4

H. Xu, M. Wang and W. Yang / Physica A 526 (2019) 120702

are ⟨k(t)+ ⟩, ⟨k(t)− ⟩ and ⟨k(t) ⟩, respectively. Similar to the calculation method used in [38], we can define the clustering coefficient, the betweenness centrality coefficient, the average path length of correlation modes of time t, respectively. We can use the dynamic topological index of correlation modes to measure the changes of various elements in the complex system. The structural characteristics of correlation modes reflect the interaction between different nodes from different aspects, which varies with time. Thus, we can use the change of structural characteristics of correlation modes to measure the risk evolution of the entire system. In Ref. [39], the absorption rate (AR) is defined to measure the evolution of systemic risks: En =

n ∑

λi /N

(3)

i=1

where λi is the eigenvalue of the correlation matrix, i = 1, 2, . . . , n, and N denotes the dimension of the correlation matrix. It can be seen that the absorption rate simply measures the evolution of the risk evolution of a complex system with the change of the number of n largest eigenvalues that contain market information but ignores the influence of other factors in the system on systematic risks. To characterize the risks of a complex system more comprehensively, we take into consideration the node characteristics, community structure and global characteristics of relevant modes and use topological indexes of correlation modes as defined above to design the dynamic typological absorption rate (DTAR) to measure systematic risks: (t )

Me

( (t )+ ) ( )−1 ( ) β1 k + β2 k(t )− + α2 β1 ℓ(t )+ + β2 ℓ(t )− + α3 β1 L(t )+ + β2 L(t )− = α1 N −1

(4)

where βi represents the weight of positive and negative correlation topological indexes between the elements of correlation modes in a complex system, and αi refers to the weight of different topological indexes in correlation modes. (t ) According to the definition of topological indexes in correlation modes and Eq. (4), it can be worked out that 0 ≤ Me ≤ 1. (t ) (t ) In order words, the closer Me is to 1, the higher the risk of the system is at timet t. Conversely, the closer Me is to 0, the lower the risk of the system is at timet t. Step 3: The Global Evolution Characteristics of Correlation Modes. The topological structure of correlation modes changes with time. To describe this process from the global perspective, we use Euclidean distance to measure the relationship between correlation modes. The Euclidean distance between correlation modes Nc (tm ) and Nc (tn ) is defined as follows:

  N N [ ]2 ∑ ∑ (t ) (t ) d [Nc (tm ) , Nc (tn )] = √ aij m − aij n ,

(t )

aij k ∈ A(tk )

(5)

i=1 j=1

Assign an appropriate threshold value to d˜ and use the following formula to construct a recursive graph:

{

} R (tm , tn ) = Θ d [Nc (tm ) , Nc (tn )] − d˜ ,

Θ (x) =

{

1, 0,

x≥0 x<0

(6)

where Θ (x) is the Heaviside function. We can use Eq. (6) to plot the relationship between correlation modes in twodimensional coordinates (both the abscissa and the ordinate are at time t). In the recursive graph, if and only if the Euclidean distance between correlation modes Nc (tm ) and Nc (tn ) is close enough, i.e., when R(tm , tn ) = 0, do we plot the sign ‘‘■’’ at (tm , tn ) and (tn , tm ). It is obvious that at (tm , tn ) and (tn , tm ), i.e., the main diagonal, the two ‘‘■’’ signs are always there. We can use it to characterize the global dynamics characteristics of changes in the correlation between the various elements of a complex system. Step 4: Build the Transitions Network of Correlation Modes. To characterize the rules of correlation modes transitions, we construct a correlation modes transitions network by using the relations between the correlation modes with time obtained in step { 2(and )} step 3. For a clear description of the construction process, we take six-dimensional random time series, i.e., X = xi t ′ , where i = 1, 2, . . . , 6, for example. Assume that the length of the random time series M = 800 and that the length of the sliding window L = 300. According to step 1, we can get a total of six time windows, as shown in Fig. 2(a) . We construct correlation modes in each time window and obtain six correlation modes, as shown in Fig. 2(b). It can be seen from Fig. 2(b) that correlation mode 2 and 5 have the same structure. Therefore, according to the chronological order of the correlation modes, the associated modal-switched directional network can be constructed, the results of which are shown in Fig. 2(d). Of course, we can also choose an appropriate threshold value according to the Euclidean distance between correlation modes to construct a directional network for correlation modes. The entire process of constructing the directional network for correlation modes is shown in Fig. 2 (a, b, c, d). Using the topological structure of the directional network, such as the node strength and clustering coefficient, we can analyze the dynamics characteristics of the correlation transitions between the modes.

H. Xu, M. Wang and W. Yang / Physica A 526 (2019) 120702

5

Fig. 2. The entire process of constructing the directional network. (a) Divide the time windows. (b) Construct correlation modes. (c) Generate a recursive graph (d) Construct a directional network of correlation modes.

3. Data and results 3.1. Data sources Daily prices of crude oil futures, WTI spot price, gasoline futures, gasoline spot price, heating oil futures, and heating oil spot price during the period between October 3, 2005 and February 8, 2016 were selected as the sample data, and they were labeled FO, FS, FG, SG, FH, and SH, respectively. The daily prices of crude oil futures, gasoline futures and heating oil futures were download from NYMEX Futures Prices database (https://www.eia.gov/dnav/pet/PET_PRI_FUT_S1_D.htm). The daily prices of WTI spot price, gasoline spot price and heating oil spot price were download from Spot Prices database (https://www.eia.gov/dnav/pet/pet_pri_spt_s1_d.htm). These six energy prices play the important role in the energy market. By analyzing the correlation of the selected six energy prices, we can reveal the relationship between the energy futures market and the spot market, and then reveal the characteristics of the energy price market. Therefore, this paper is focused on the correlations among different data, but the data selected here have different units, so they are first standardized. The results are as follows: It can be seen from Fig. 3 that there is a relatively high correlation between every two energy prices (the smallest Pearson correlation coefficient is 0.8674), which indicates that fluctuations of different energy prices follow the same pattern. However, it needs to be noticed that the correlation among different energy prices in different periods is not a constant, it changes with time, which is nonlinear and unstable. For instance, the correlation coefficient of the futures price of crude oil and the spot price of gasoline is 0.6410 during the period between October 3, 2005 and February 6, 2008, which shows that the two prices are weakly correlated. However, during the period between December 14, 2006 and July 21, 2008, the correlation coefficient is 0.9099, which indicates that they are strongly correlated. It can be seen that, during the process of fluctuations of energy prices, the correlation among them continues changing. This change reflects the complexity of the energy price system. Study of the change can effectively reveal the nature of the energy price system. In the following subsection, the evolution features of the correlations among energy prices will be studied based on the method of ETCCMMTS. 3.2. Construction of correlation modes for the energy price system. Since we want to analyze the evolution characteristics of correlations among energy prices, the length of the sliding window L and the step l are setted neither too big nor too small [31–37]. Because the data selected in this paper are daily energy prices and there are five price data each week, 2 months (40 days) was selected as the length of a time window, and 3 weeks (15 days) was selected as the time step, i.e., L = 40l = 15, a small-scale time window can be constructed.

[ T =

M −L l

] [ ] 2602 − 40 +1 = + 1 = 171 15

(7)

6

H. Xu, M. Wang and W. Yang / Physica A 526 (2019) 120702

Fig. 3. Fluctuation charts of standardized energy prices.

In each time window, the correlation coefficient of energy prices was calculated according to Eq. (1). The results show that different energy prices have positive correlations in each time window. Therefore, the adjacent matrices of correlation modes constructed in each time window consist only of the elements 0 and 1. That is, if there are connected edges among different energy prices, the connected edges are in the form of Fig. 1(a). In the following parts, to visualize the drawing of correlation modes for energy prices, a line without arrows is used to represent a connected edge. The threshold value of rc in step 1 is defined as follows: rc =

T ∑ ⟨ t =1

C (t )

⟩/

T , C (t ) =

⟨

⟩

N N ⏐ ⏐/ ∑ ∑ ⏐ (t ) ⏐ ⏐Cij −Iij ⏐ (N × N )

(8)

j=1 i=1

where Iij ∈ I and I represents the N-dimensional }unit matrix. The correlation modes of energy prices during the fluctuation { process in each time window are NC = NC (t ) , t = 1, 2, . . . , 171. Correlation modes in some of the time windows are ’’ in the correlation modes mean connection ‘‘ ’’ shown in Fig. 4, where ‘‘ From Fig. 4, it can be inferred that the structure of correlation modes during the fluctuation of energy prices in different small-scale time windows changes with time. The change indicates changes in correlation among different energy prices. During the period between October 24, 2014 and December 10, 2014, the correlation mode of energy prices has 8 connected edges (as shown in the 7th chart in Fig. 4). In this mode, the gasoline features share edges with the crude oil futures price, the crude oil spot price and the gasoline spot price, which means that the gasoline features price is highly correlated to the crude oil futures price, the crude oil spot price and the gasoline spot price. In other words, the fluctuation of the gasoline features price is more affected by fluctuations of the crude oil futures price, the crude oil spot price and the gasoline spot price and less by the fluctuations of the heating oil futures price and the heating spot price. However, during the period between December 30, 2015 and January 21, 2016, the number of connected edges in the correlation mode was reduced to 5 (as shown in the 9th chart in Fig. 4), which shows that the overall correlation among energy prices is low. The mode is divided into two independent parts in this period. In one part, the heating oil futures price and the heating oil spot price share no connected edges with other energy prices, which indicate that fluctuations of the two prices are slightly affected by those of other energy prices. In another part, the futures price of gasoline shares 3 connected edges with other energy prices, which indicates that, although the overall correlation among energy prices in the energy price system becomes weak, some correlations remain stable or even become stronger. This also shows that evolution features of correlation modes of energy prices can be used to characterize some evolution features of the energy price system. The whole time series set was divided into small fragments by the sliding windows with the sliding step length l. The advantage of utilizing the method is that the fragments have the feature of memory and transitivity [31–37]. With different L and l, we can obtain the different correlation modes. The results of the sensitivity analysis onL and l are shown in Fig. 5. We can set different L and l depending on different needs of analysis. 3.3. Evolution characteristics of correlation modes in the energy price system To get a clearer analysis of evolution features of correlation modes in the energy price system, an analysis of the stochastic system was first made based on dynamic typological indexes of correlation modes defined in step 2. To make a

H. Xu, M. Wang and W. Yang / Physica A 526 (2019) 120702

7

Fig. 4. Correlation modes in some of the small-scale time windows.

Fig. 5. Sensitivity analysis on L and l. (a) Percentage of positively correlated. (b) Percentage of uncorrelated. (c) Percentage of negatively correlated for different values of scale L and l. (d) Percentage of negatively correlated for different values of scale L.

8

H. Xu, M. Wang and W. Yang / Physica A 526 (2019) 120702

Fig. 6. Evolution graphs of the topological indexes in the stochastic system. (a) The average positive degree. (b) The average positive clustering coefficient. (c) The average positive path length. (d) The average negative degree. (e) The average negative clustering coefficient. and (f) The average negative path length.

Table 1 Average values of typological indexes of energy prices and stochastic system. Energy price system

⟨k⟩ ⟨ℓ⟩ ⟨B⟩ ⟨L⟩

Stochastic system

FO

SO

FG

SG

FH

SH

Average

4.1345 0.9160 0.8881 –

4.1170 0.9228 0.8878 –

4.0000 0.8559 1.0536 –

3.3977 0.8298 0.1725 –

4.0292 0.7554 1.9237 –

3.6667 0.7345 1.0978 –

3.8908 0.8357 1.0039 1.5952

1.0117 0.4234 1.6667 1.8538

comparison between the stochastic system and the energy price system, six random time series with the same lengths as the energy price series were generated using MATLAB. To make the time series more general, a computer simulation was run 1000 times to get six average random time series. Then, the same parameters were taken out of the six time series by using the method of constructing correlation modes for energy prices to construct correlations modes in the stochastic system. Unlike the energy price system, there are both positive and negative correlations in the stochastic system. After calculation, the average positive topological indexes are shown in Fig. 6(a, b, c) and the negative topological indexes are shown in Fig. 6(d, e, f). Green solid lines in Fig. 6 represent the 95% confidence interval, and blue lines refer to the average values of corresponding indexes. It can be seen from Fig. 6 that correlation modes in the stochastic system feature small correlation values, low clustering coefficients and long average paths. These features are typical of a stochastic system. The above calculation shows that there are only positive correlation modes in the energy system when we set L = 40, l = 15. To make a comparison among correlation modes in the energy system and the stochastic system, we constructed the correlation modes in the stochastic system: (t )

aij

⏐ ⏐ ⎧ ⎨ 1, ⏐⏐cij(t ) ⏐⏐ ≥ rc ⏐ ⏐ = ⎩ 0, ⏐⏐c (t ) ⏐⏐ < r c ij

(9)

After calculation, the evolution graphs of the topological indexes in the energy price system at time t are shown in Fig. 7(a, b, c), and those in the stochastic system are shown in Fig. 7(d, e, f). Fig. 7 shows the trends of evolution concerning the average correlations, the average clustering coefficients and the average path lengths of both the energy price system and the stochastic system. It can be inferred that the evolution of dynamic indexes in the energy price system is more regular than that of dynamic indexes in the stochastic system. According to Fig. 7, the average values of typological indexes in the energy price system and in the stochastic system are calculated and shown in Table 1. First, we made a comparison of the dynamic evolution indexes between energy price system and stochastic system. It was found that the average correlation for energy prices is 3.8908, which shows that there is an average of four types of energy prices that are highly correlated with the fluctuations at any time. Meanwhile, the average correlation for random time series is 1.0117, which shows that the fluctuations of random time series are independent of each other. The average

H. Xu, M. Wang and W. Yang / Physica A 526 (2019) 120702

9

Fig. 7. Evolution graphs of the topological indexes in the energy price system. (a) The average degree. (b) The average clustering coefficient. (c) The average path length, and in the stochastic system. (d) The average degree. (e) The average clustering coefficient. (f) The average path length.

clustering coefficient in energy price system is 0.8357, much larger than 0.4234, which is the average clustering coefficient in the stochastic system. This indicates that there is a clustering effect in the correlations among energy prices during their fluctuations. The average path length in the energy price system is 1.5952, slightly shorter than that in the stochastic system, which indicates that the correlations among energy prices take on the small-world property. What are the roles of different energy prices in the analysis of their correlations during fluctuation? This question can be answered through the analysis of dynamic typological indexes. It can be seen from Table 1 that almost dynamic typological indexes of energy futures prices are larger than those of corresponding energy spot prices, indicating that the energy futures price is in a leading position in the energy price system. Previous research studies on the relationship between the energy futures market and energy spot market have drawn the same conclusion. However, they have seldom explored the roles of each energy price in the whole entire energy price system. By analyzing the results in Table 1, the roles of different energy prices in the energy price system can be differentiated. It is found that the average correlation of the crude oil futures price (FO) is the largest, showing that it is most closely correlated to other energy prices in the entire energy price system and that it plays an important role in maintaining the stability of the system. The average clustering coefficient of crude oil spot prices (SO) is the largest, which indicates that it has high polymerization in the energy price system. The average betweenness coefficient of the heating oil futures price (FH) is the largest, which shows that it plays an important intermediary role in the energy price system. The sensitivity analysis on L and l for the topological indexes of correlation modes in the energy price system are shown in Fig. 8(a, b, c). With an increase in scaleL and l, the average degree and the average clustering coefficient increase (see Fig. 8(a) and (b)), the average path length decrease (see Fig. 8(c) and (d)). Generally, the correlation modes in the energy price system feature bigger degree values, higher clustering coefficients and shorter average paths compared with the stochastic system. To summarize, based on ETCCMMTS, we not only get the evolution features of correlation modes in the energy price system, but also we can obtain the roles of different energy prices in energy price system. 3.4. Measurement of risks in the energy price system Measurement of the risk in the energy price system has been a research hotspot in the field of energy economics. The dynamic typological absorption rate RE (t ) (Eq. (4)) defined in step 2 will be used to measure system risks in the energy price system. Because there are only positive correlations in the correlation modes of the energy price system, Eq. (4) can be reduced to: Me(t ) = α1 k(t ) /(N − 1) + α2 ℓ(t ) + α3 L(t )

⟨

⟩

⟨

⟩

⟨

⟩−1

(10)

In the equation, system risks of the stochastic system can be measured by calculating the dynamic typological absorption rate and absorption rate defined in Ref. [39] separately to get the evolution of the stochastic system risk, which is shown in Fig. 9(a, b). The comparison between Figs. 9(a) and 9(b) shows that risks measured by the two methods follow a similar trend of evolution with both small values. This indicates that risks in the stochastic system are kept at a low level. In addition, compared with the absorption rate defined in Ref. [39], the dynamic typological absorption rate is more disperse, differential, and precise in characterizing the energy price system. The results of system risks of the energy price system measured by the above two methods are shown in Fig. 9 (c, d).

10

H. Xu, M. Wang and W. Yang / Physica A 526 (2019) 120702

Fig. 8. Sensitivity analysis on L and l. (a) The average degree. (b) The average clustering coefficient. (c) The average path length for different values of scale L and l. (d) The average path length for different values of scale L.

Fig. 9. The results of system risk. (a) The dynamic typological absorption rate of the stochastic system. (b) The absorption rate of the Stochastic System. (c) The dynamic typological absorption rate of the energy price system. (d) The absorption rate of the energy price system.

The comparison between Fig. 9(c) and (d) shows that the dynamic typological index can better reflect the risks evolution futures of the energy price system than the absorption rate. Fig. 9(c) shows that there was a ‘‘risk jump’’ in the energy price system on November 21, 2008 caused by the financial crisis in early 2008, and then the system risk kept increasing. It did not decrease until another ‘‘risk jump’’ appeared on March 14, 2012. Thanks to the development

H. Xu, M. Wang and W. Yang / Physica A 526 (2019) 120702

11

of various clean energies in recent years, the system risk of energy price was maintained at a low level. Fig. 9(d) and (c) follow similar trends of evolution, but the former is not differential enough and fails to reveal the ‘‘risk jump’’ in the energy price system. The above empirical analysis shows that the dynamic typological absorption rate (DTAR) we defined in this paper can measure system risk effectively. 3.5. Global characteristics of the change of energy price fluctuation correlation modes With step 3, the distance matrices between stochastic systems and energy price system are first calculated, and the results are shown in Fig. 10(a) and (b). In Fig. 10(a), apart from the diagonal, the data is uniform distributed, which shows that the distances between the correlation modes for stochastic systems are almost the same. In contrast, in Fig. 10(b), the data is distributed in blocks, which indicates that the distance distribution between the correlation modes for the energy price system is more complicated. To further extract this complexity, we took the threshold d˜ and calculated and drew the long-range global recursive graphs of the stochastic system and the energy price system, respectively. The results are shown in Fig. 10(c) and (d). The principle of selecting the threshold of d˜ is that the points on the long-range global recursive graph of the stochastic system fall roughly near the main diagonal (as shown in Fig. 10(c)), which means that the distances between the correlation modes for stochastic systems are very large and that there is hardly any correlation between the correlation modes. In contrast, according to the mapping method of the long-range recursive graph of step 3, the long-range global recursive graph of the energy price system is shown in Fig. 10(d), which indicates that the correlation modes of the energy price system have both short-range correlation and long-range correlation. To obtain the cycle length of the correlation, we randomly extracted the correlation length of the correlation modes of the stochastic systems and the energy price systems in Fig. 10(c) and (d). The results are shown in Fig. 10(e), in which the blue bar represents the correlation length distribution of the correlation modes for the energy price systems and the red bar shows the correlation length distribution of the correlation modes for the stochastic systems. It can be seen that, for the stochastic systems, there is almost no correlation between the correlation modes, and there is only a short-range correlation between several modes. However, for the energy price systems, the correlation modes have both 1–2 months’ short-range correlation and 7–8 years’ long-range correlation, indicating that the correlation between the fluctuations of various energy prices in the energy price systems show complex global change characteristics. According to the change of the time interval between the correlation modes for the energy price system, the energy price fluctuation system can be divided into three different periods. The results are shown in Fig. 10(f), where, because the time interval represented by the blue ‘‘□’’ is too short (from December 16, 2009 to January 12, 2010), it cannot be listed as a separate period. The correlation patterns in the energy price systems show short-range correlations before November 30, 2008 and long-range correlations between November 30, 2008 and December 16, 2009. Later, it is followed by short-range correlations in short-time and long-range correlations. Since March 13, 2012, it showed short-range correlations again. Combined with the energy price system risk analysis above, it can be concluded that, when it shows long-range correlations between the various energy price fluctuations in the energy price system, the energy price system will face a higher systemic risk; meanwhile, when it shows short-range correlations, the energy price system will face a lower systemic risk. Therefore, it is possible to effectively identify the risks in the energy price system by making use of the changes in the correlations between the correlation modes of the energy price system. 3.6. Energy price correlation modes transmission network (EPCMTN) In the above analysis, 171 correlation modes of energy prices are obtained, but in fact, the structures of energy price correlation modes at different periods might be the same. For instance, the correlation modes of energy prices between October 21, 2008 and December 26, 2008 (as shown in the fourth graph of Fig. 4) and those between October 29, 2010 and December 27, 2010 (as shown in the sixth graph of Fig. 4) are the same. Only 41 correlation modes have different structures, and their structures are shown in Fig. 11. Are there temporal patterns concerning the appearance of correlation modes with different structures? To answer this question, the cumulative time intervals of correlation modes with different structures in the energy price fluctuation process have been examined, as shown in Fig. 12. It can be seen that the time intervals of correlation modes with different structures in the process of energy price fluctuation are not evenly spaced, but the time intervals between new modes gradually increase with time. When the modes accumulate to a certain number (the 24th mode in the figure), the next new mode will appear via a ‘‘cumulative time jump’’ (as shown by the red dotted circle in Fig. 12), which indicates that the appearance time of correlation modes for the energy price system is heterogeneous. Fig. 12 also indicates that the appearance time of the different modes before and after the ‘‘jump’’ is distributed in a straight line. With the least squares regression, we carried out regression on the appearance time of correlation modes with different structures before and after the ‘‘jump.’’ The corresponding regression equations are y = 2.9304x − 9.2971 and y = 3.8382x + 106.8676, and the corresponding R2 are 0.9315 and 0.9805, respectively. These results mean that the reliability of the results is high, which indicates that the cumulative time of correlation modes with different structures in the process of energy price fluctuations before and after the ‘‘jump’’ shows a linear growth trend. However, in application, for an unknown system, it is difficult to determine when the ‘‘jump’’ will appear due to the absence of relevant data. If the difference before and after the ‘‘jump’’ is put aside and the cumulative time of correlation modes with different structures in the process of energy

12

H. Xu, M. Wang and W. Yang / Physica A 526 (2019) 120702

Fig. 10. The results of global characteristics. (a) Distance distribution of the correlation modes for a stochastic system. (b) Distance distribution of the correlation modes for an energy price system. (c) Long-range global recursive graph of a stochastic system. (d) Long-range global recursive graph of an energy price system. (e) Correlation length distribution of a stochastic system and energy price system. (f) Division of energy price fluctuation periods.

price system is taken as a period of time, with the least squares, the regression equation will be y = 4.9143x−28.5659, and the corresponding R2 is 0.9479. The accumulation time of correlation modes with different structures in the energy price system over the entire period also obeys the linear growth trend. In summary, the accumulation time of correlation modes with different structures in the energy price system is not chaotic but subject to a linear growth trend, and a ‘‘jump’’ phenomenon appears in the growth process. Forty-one correlation modes with different structures are transformed into each other over time, and the energy price correlation modes transmission network (EPCMTN) is obtained, which is shown in Fig. 13. The numbers of nodes in Fig. 13 correspond to those in Fig. 11. It can be seen clearly from Fig. 13 that, in EPCMTN, few nodes have high strength, while most nodes have low strength, which indicates that EPCMTN is a scale-free network. 3.6.1. Node strength and strength distribution of EPCMTN Because EPCMTN is a directed weighted network, the node strength and the strength distribution will be discussed in this section [38]. For an unweighted network, the degree of a node reflects the connection of the node with other nodes. In a weighted network, the strength of a node has a similar nature. The strength of a node reflects not only the connection between the node and other nodes but also the weight of the connected edges. According to the construction of the EPCMTN network, the nodes are connected in a time sequence. Therefore, apart from the first and last nodes, the strength towards and outwards the nodes is the same. In this session, only the intensities towards the nodes are calculated

H. Xu, M. Wang and W. Yang / Physica A 526 (2019) 120702

13

Fig. 11. Energy price fluctuation correlation modes with different structures.

Fig. 12. Cumulative time distribution of energy price correlation modes with different structures.

and analyzed, and for the sake of narration, the strength towards a node will be directly referred to as the strength. The nodes strength of the EPCMTN are shown in Fig. 14(a). It is clear that few nodes are of high strength, while most nodes are of low strength, and the strength of 21.95% nodes accounts for 64.12% of the total strength. The double logarithm graph of the EPCMTN node strength distribution is shown in Fig. 14(b). It can be seen from Fig. 14(b) that the EPCMTN node strength follows a power-law distribution with a power index of γ = 0.7031.

14

H. Xu, M. Wang and W. Yang / Physica A 526 (2019) 120702

Fig. 13. Energy price correlation modes transmission network (EPCMTN).

Fig. 14. Node strength and strength distribution of EPCMTN. (a) Node strength. (b) Double logarithmic distribution of node strength.

The strength of a node reflects the importance of the node in the entire network. We can obtain the top six nodes in strength in EPCMTN are 19, 18, 22, 34, 8 and 14. The fully connected mode of ‘‘crude oil-gasoline-heating oil’’ (Mode 19) has the highest node strength in EPCMTN, indicating that the fully connected mode enjoys a core status in EPCMTN and a high degree of linkage between energy prices during energy price fluctuations. 3.6.2. Other topologies of correlation modes in EPCMTN Clustering coefficients are usually used to describe the local structural properties of a network, and the clustering coefficients of nodes reflect the degree of connection between adjacent nodes [38]. The nodes with the highest clustering coefficient in EPCMTN are 22, 1, 10, 20, 21 and 27. The relationship between clustering coefficients and node strength is shown in Fig. 15(a). It can be seen that small nodal strength nodes in EPCMTN tend to have high clustering coefficients. Mode 22 has both the highest clustering coefficient and the highest node strength, which indicates that Mode 22 plays an important role in EPCMTN. How long does it take for one correlation mode in EPCMTN to transform into another? To answer this question, the shortest path length of EPCMTN should be calculated. With the Floyd algorithm, the distance between any two nodes

H. Xu, M. Wang and W. Yang / Physica A 526 (2019) 120702

15

Fig. 15. The results of other topologies. (a) Relationship between clustering coefficient and node strength. (b) Distance distribution among correlation modes in EPCMTN. (c) Path length distribution in EPCMTN. (d) Relationship between betweenness and node strength.

in EPCMTN has been calculated, and the results are shown in Fig. 15(b). A different path length distribution is shown in Fig. 15(c). Fig. 15(b) shows that the EPCMTN has a diameter of 16, indicating that the path length between any two modes in the EPCMTN will not exceed 16. As shown in Fig. 15(c), the distance distribution between any two modes in the EPCMTN nearly obeys the normal distribution, and the distance ranging from five to ten between two modes accounts for 60% of the total modes in EPCMTN. The average distance is 6.9527, indicating that the energy price correlation shows short-range correlation. On average, a correlation mode needs to go through six to seven modes before being transformed into any other correlation mode, which indicates that there is relatively frequent conversion between modes in EPCMTN. In addition, it has been found that there are a few core modes in EPCMTN. Even if the current mode is not a core mode, it can go back to the core mode after six to seven transition modes, which has important practical significance for predicting the change of the correlation between energy price systems. To analyze the transition modes in EPCMTN, the top six nodes in node betweenness in EPCMTN are 19, 16, 18, 14, 11 and 1, and the relationship between the node betweenness and the node strength in EPCMTN is shown in Fig. 15(d). There are three common modes, namely the modes of 19, 16 and 14, among which Mode 19 and Mode 16 have both high node strength and betweenness. Fig. 15(d) shows that node with the highest strength in EPCMTN has the highest node betweenness. To summarize, the core correlation modes in the energy price system are the fully connected modes of ‘‘crude oil gasoline - heating oil’’ (Mode 19), the fully connected mode of ‘‘crude oil - gasoline’’ (Mode 16) and the fully connected mode of ‘‘gasoline - heating oil’’ (Mode 22). The occurrence probability of these three core modes is 19.9406 times, 3.0099 times and 3.0099 times, respectively, the appearance of other non-core modes, which indicates that the correlation modes in the energy price fluctuation process transform into each other with the three modes as the focus. The ‘‘gasoline-heating oil full linkage’’ mode (Mode 22) has the largest polymerization effect, and the ‘‘crude oil - gasoline - heating oil full linkage’’ mode (Mode 19) and the ‘‘crude oil - gasoline full linkage’’ mode (Mode 16) play an important mediating role. These three core modes transform into each other in the energy price fluctuation process. If the current correlation mode is not a core mode, it can go back to the core mode after six to seven transition modes. 4. Discussion The main idea of this paper is to reveal the inherent characteristics of complex systems from the evolution of multiple time series correlation. Therefore, the main study object of this paper is the change of multiple time series correlation

16

H. Xu, M. Wang and W. Yang / Physica A 526 (2019) 120702

rather than time series. Based on the complex network theory, the authors propose the new four-step method for analyzing the evolution and transformation characteristics of correlation modes of multiple time series (ETCCMMTS). The first step is to construct a sliding window to divide the multiple time series into a series of small modules and convert the series of small modules into a series of correlation modes. Second, introduce the dynamic topological structure index and the dynamic topological entropy to describe the evolution of the time series correlation and the evolution of the system risk. Third, introduce a global recursive graph to show the evolution of multiple time series correlation. Finally, construct a directed weighted modes transformation network to reveal the transformation characteristics of the correlation of multiple time series. The numerical results verify the practicability of the method proposed in this paper. With the method of ETCCMMTS, the characteristics of the energy price system from the perspective of the correlation between energy price fluctuations are as follows. (1) The average degree and the average clustering coefficient of energy price system correlated modes are 3.8908 and 0.8357, respectively, and they are larger than the average degree (1.0117) and average clustering coefficient (0.4234) of the correlation modes for stochastic systems, respectively, which indicates that there is linkage and concentration in the process of energy price fluctuation. However, the average betweenness and path length of the correlation modes are 1.0039 and 1.5952, respectively, which are smaller than the average betweenness (1.6667) and average path length (1.8358) of the correlation modes for the stochastic systems. This shows that the energy price fluctuation has small-world characteristics. (2) Through empirical calculation, it is found that the topological indexes of the energy futures prices are all larger than the energy spot prices, and the leading role of the energy futures price in the price system is explained from the network topological structure. The heating oil futures price has the largest average degree (3.1965), the crude oil futures price enjoys the largest average clustering coefficient (0.7334), and the gasoline futures price enjoys the largest average betweenness (1.3266), indicating that the heating oil futures price affects the energy price system stability, the crude oil futures price has a polymerization effect in the energy price system, and the gasoline futures price plays the intermediary function in the energy price system. (3) The absorptivity of the dynamic topology is introduced to measure the energy price risk. It is found that there are two ‘‘risk jumps’’ in the evolution of the energy price system risk. One occurred on November 21, 2008, and the other occurred on March 14, 2012, and since then, the risk of the energy price system began to decline. Especially in recent years, with the development of various clean energies, the risk level in the energy price system has remained relatively low. (4) The global change of correlation modes for stochastic systems is evenly distributed. In contrast, the global change of the energy price correlation modes presents an even block distribution, which shows that the energy price fluctuation has both long-range and short-range correlation. It is also found that long-range correlation of source price fluctuation can lead to high risks of the energy price system. (5) The 41 different correlation modes for energy price fluctuation are obtained. The results show that the accumulation time of different correlation modes is not random but instead shows a linear growth trend. With the construction of EPCMTN, it is found that the node strength of EPCMTN is subject to a power-law distribution. Although the energy price fluctuation is highly complex, the transformation of the energy price fluctuation correlation focuses mainly on three core correlation modes. It takes six to seven non-core transmission correlation modes, on average, to return to a core mode. These results are of great practical significance to reveal the essential characteristics of the energy price system and to construct a new energy price forecasting system. Acknowledgments The Research was supported by the following foundations: The Philosophy and Social Science Research Project in Colleges and Universities of Jiangsu Province, China (2018SJA2123), The National Natural Science Foundation of China (71503132, 71811520710, 11571142), Qing Lan Project of Jiangsu Province, China (2017), Six talent peaks project in Jiangsu Province, China (JY-055). References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13]

A.L. Barabási, R. Albert, Emergence of scaling in random networks, Science 286 (5439) (1999) 509–512. D.J. Watts, S.H. Strogatz, Collective dynamics of ‘small-world’networks, Nature 393 (6684) (1998) 440–442. M.E.J. Newman, D.J. Watts, Renormalization group analysis of the small-world network model, Phys. Lett. A 263 (1999) 341–346. P. Erdős, A. Rényi, On the existence of a factor of degree one of a connected random graph, Acta Math. Hungar. 17 (3–4) (1966) 359–368. J. Zhang, M. Small, Complex network from pseudoperiodic time series: Topology versus dynamics, Phys. Rev. Lett. 96 (23) (2006) 238701. L. Lacasa, B. Luque, F. Ballesteros, et al., From time series to complex networks: The visibility graph, Proc. Natl. Acad. Sci. 105 (13) (2008) 4972–4975. B. Luque, L. Lacasa, F. Ballesteros, et al., Horizontal visibility graphs: Exact results for random time series, Phys. Rev. E 80 (4) (2009) 046103. I.V. Bezsudnov, A.A. Snarskii, From the time series to the complex networks: The parametric natural visibility graph, Physica A 414 (2014) 53–60. Z.K. Gao, Q. Cai, Y.X. Yang, et al., Multiscale limited penetrable horizontal visibility graph for analyzing nonlinear time series, Sci. Rep. 6 (2016) 35622. X. Li, M. Sun, C. Gao, et al., The parametric modified limited penetrable visibility graph for constructing complex networks from time series, Physica A 492 (2018) 1097–1106. M. Wang, A.L.M. Vilela, R. Du, et al., Exact results of the limited penetrable horizontal visibility graph associated to random time series and its application, Sci. Rep. 8 (1) (2018) 5130. Y. Yang, J. Wang, H. Yang, et al., Visibility graph approach to exchange rate series, Physica A 388 (20) (2009) 4431–4437. M. Stephen, C. Gu, H. Yang, Visibility graph based time series analysis, PLoS One 10 (11) (2015) e0143015.

H. Xu, M. Wang and W. Yang / Physica A 526 (2019) 120702

17

[14] M.C. Qian, Z.Q. Jiang, W.X. Zhou, Universal and nonuniversal allometric scaling behaviors in the visibility graphs of world stock market indices, J. Phys. A 43 (33) (2010) 335002. [15] Y. Yang, H. Yang, Complex network-based time series analysis, Physica A 387 (5–6) (2008) 1381–1386. [16] X. Xu, J. Zhang, M. Small, Superfamily phenomena and motifs of networks induced from time series, Proc. Natl. Acad. Sci. 105 (50) (2008) 19601–19605. [17] R.V. Donner, Y. Zou, J.F. Donges, et al., Recurrence networks—a novel paradigm for nonlinear time series analysis, New J. Phys. 12 (3) (2010) 033025. [18] N. Marwan, J.F. Donges, Y. Zou, et al., Complex network approach for recurrence analysis of time series, Phys. Lett. A 373 (46) (2009) 4246–4254. [19] Q. Tang, J. Liu, H. Liu, Comparison of different daily streamflow series in US and China, under a viewpoint of complex networks, Modern Phys. Lett. B 24 (14) (2010) 1541–1547. [20] W.D. Chen, H. Xu, Q. Guo, Dynamic analysis on the topological properties of the complex network of international oil prices, Acta Phys. Sin. 59 (7) (2010) 4514–4523. [21] M. Wang, L. Tian, From time series to complex networks: The phase space coarse graining, Physica A 461 (2016) 456–468. [22] X. Gao, H. An, W. Fang, et al., Characteristics of the transmission of autoregressive sub-patterns in financial time series, Sci. Rep. (2014) 4. [23] X. Gao, H. An, W. Fang, et al., Transmission of linear regression patterns between time series: From relationship in time series to complex networks, Phys. Rev. E 90 (1) (2014) 012818. [24] M. Wang, Y. Chen, L. Tian, et al., Fluctuation behavior analysis of international crude oil and gasoline price based on complex network perspective, Appl. Energy 175 (2016) 109–127. [25] H. An, X. Gao, W. Fang, et al., Research on patterns in the fluctuation of the co-movement between crude oil futures and spot prices: A complex network approach, Appl. Energy 136 (2014) 1067–1075. [26] X. Gao, W. Fang, F. An, et al., Detecting method for crude oil price fluctuation mechanism under different periodic time series, Appl. Energy 192 (2017) 201–212. [27] M. Wang, L. Tian, Regulating effect of the energy market—Theoretical and empirical analysis based on a novel energy prices–energy supply–economic growth dynamic system, Appl. Energy 155 (2015) 526–546. [28] H. Chen, L. Tian, M. Wang, et al., Analysis of the dynamic evolutionary behavior of American heating oil spot and futures price fluctuation networks, Sustainability 9 (4) (2017) 574. [29] K. Iwayama, Y. Hirata, K. Takahashi, et al., Characterizing global evolutions of complex systems via intermediate network representations, Sci. Rep. (2012) 2. [30] L. Lacasa, V. Nicosia, V. Latora, Network structure of multivariate time series, Sci. Rep. (2015) 5. [31] Z.K. Gao, X.W. Zhang, N.D. Jin, et al., Recurrence networks from multivariate signals for uncovering dynamic transitions of horizontal oil-water stratified flows, Europhys. Lett. 103 (5) (2013) 50004. [32] Z.K. Gao, S.S. Zhang, W.D. Dang, et al., Multilayer network from multivariate time series for characterizing nonlinear flow behavior, Int. J. Bifurcation Chaos 27 (04) (2017) 1750059. [33] R. Du, Y. Wang, G. Dong, et al., A complex network perspective on interrelations and evolution features of international oil trade, 2002–2013, Appl. Energy 196 (2017) 142–151. [34] M. Wang, L. Tian, H. Xu, et al., Systemic risk and spatiotemporal dynamics of the consumer market of China, Physica A 473 (2017) 188–204. [35] Z.K. Gao, Y.X. Yang, P.C. Fang, et al., Multiscale complex network for analyzing experimental multivariate time series, Europhys. Lett. 109 (3) (2015) 30005. [36] Z.K. Gao, M. Small, J. Kurths, Complex network analysis of time series, Europhys. Lett. 116 (5) (2017) 50001. [37] M. Wang, L. Tian, R. Du, Research on the interaction patterns among the global crude oil import dependency countries: A complex network approach, Appl. Energy 180 (2016) 779–791. [38] X.F. Wang, X. Li, G.R. Chen, The Theories and Application of Complex Network, Tsinghua Unversity Press, Beijing, 2016 (in Chinese). [39] H. Meng, W.J. Xie, Z.Q. Jiang, et al., Systemic risk and spatiotemporal dynamics of the US housing market, Sci. Rep. 4 (2014) 3655. [40] X. Zhang, L. Feng, Y. Berman, et al., Exacerbated vulnerability of coupled socio-economic risk in complex networks, Europhys. Lett. 116 (1) (2016) 18001. [41] X. Zhang, B. Podobnik, D.Y. Kenett, et al., Systemic risk and causality dynamics of the world international shipping market, Physica A 415 (2014) 43–53. [42] B. Podobnik, H.E. Stanley, Detrended cross-correlation analysis: A new method for analyzing two nonstationary time series, Phys. Rev. Lett. 100 (8) (2008) 084102. [43] W.X. Zhou, Multifractal detrended cross-correlation analysis for two nonstationary signals, Phys. Rev. E 77 (6) (2008) 066211.