Network structure of multivariate time series

Our understanding of a variety of phenomena in physics, biology and economics crucially depends on the analysis of multivariate time series. While a wide range tools and techniques for time series analysis already exist, the increasing availability of massive data structures calls for new approaches for multidimensional signal processing. We present here a non-parametric method to analyse multivariate time series, based on the mapping of a multidimensional time series into a multilayer network, which allows to extract information on a high dimensional dynamical system through the analysis of the structure of the associated multiplex network. The method is simple to implement, general, scalable, does not require ad hoc phase space partitioning, and is thus suitable for the analysis of large, heterogeneous and non-stationary time series. We show that simple structural descriptors of the associated multiplex networks allow to extract and quantify nontrivial properties of coupled chaotic maps, including the transition between different dynamical phases and the onset of various types of synchronization. As a concrete example we then study financial time series, showing that a multiplex network analysis can efficiently discriminate crises from periods of financial stability, where standard methods based on time-series symbolization often fail.

where x ∈ R d and S ∈ {s 1 , s 2 , ..., s p } d . Once the original (continuous state) series {x(t)} is symbolized into a (discrete state) sequence of symbols {S(t)}, several measures such as the ones involving frequency histograms of finite size samples (as unbiased estimators of probability densities) such as complexity or information [1,5] measures can be computed and used to analyse the system. A large majority of the methods that analyse empirical time series require such preprocessing, although it is fair to say that this is not always clarified or explicitly stated. However, such preprocessing is far from trivial, in the sense that results usually depend on such ad hoc procedure. First, how many symbols p should we define? Note that standard option, widely used in the field of symbolic dynamics, is to makes use of two symbols (p = 2), where s 1 ≡ L and s 2 ≡ R (alternatively and without loss of generality, s 1 ≡ 1 and s 2 ≡ 2). In general, each symbol indeed corresponds to the label of a different, non-overlapping cell c i , such that the set C = ∪ p i=1 c i is a partition or tile of the phase space under study. For instance, for p = 2 the function f : [a, b] → [a, b] is usually symbolized according to a homogeneous partition: c 1 = [a, b/2), c 2 = [b /2, b]. Now, again, is this is best choice for the definition of cells? The response is also not unique, and as a matter of fact, each dynamical system will typically require a different symbolisation and partition. For instance, suppose that the distribution of phase space visits of a given map is not uniform but it is peaked around some neighbourhood. Intuitively, one would need to refine the partition in that neighbourhood, to capture fine-grained details which would otherwise be lost. These kind of situations inevitably require a full exploration of the map's phase space prior to any symbolization.
More dramatically, note that the dynamics induced by the symbolization process can in general be very different under different phase space partitions and symbolisations. Concretely, only for so called generating partitions the dynamics after symbolization remain equivalent, although which are the generating partitions is a nontrivial question that lacks a general solution. As a toy example, let us consider the dynamics of a logistic map f (x) = rx(1 − x) in [0, 1], where standard symbolic dynamics uses p = 2. For r = 3.6 slightly below the first Misiurewicz point, the map is chaotic, but the chaotic attractor is splitted into two disjoint sets. Orbits in this case densely fill the attractor, although they make 'jumps' between each chaotic band in a periodic manner. Therefore, if we the induced dynamics would become totally periodic, and every trace of the chaotic dynamics would be lost (see figure 1 for an illustration if this toy example). An autocorrelation function or a frequency analysis (power spectrum) -options that do not require a series symbolisation-would also misleadingly conclude that the system is in a periodic state. Interestingly enough, in such a (pathological) toy model, the visibility approach is not misleaded [2].
Note also that finite size effects can play an important role in the symbolisation preprocessing. If the series under study is too short, an excessive number of symbols would inevitably introduce spurious finite size effects, and accordingly, every measure based on the statistics of the symbolised series would be artificially biased. The optimal number of symbols therefore usually depends as well on the size of the series under study.
Finally, it is worth recalling that value of measures computed from symbolised series, such as the well known Shannon entropy, depend on the particular symbolisation, and hence are not invariant under smooth coordinate changes [1] (or under simple changes of experimental units!). All these issues are usually seen as intrinsic and inevitable drawbacks of the symbolisation method, as this is at least dependent on both p and the specific partition in an ad hoc way. Note at this point that the visibility method also induces a symbolisation in the original (multivariate) time series, as a mapping can be straightforwardly defined between the series and the multiplex network's vector integer degree sequence, where k(t) ≡ (k [1] (t), k [2] (t), ..., k [d] (t)) ∈ Z d and by construction k [α] (t) > 1. Within this 'network symbolization', note that the number of symbols p is not a free parameter that needs to be tuned anymore, much on the contrary, the number of symbols emerge naturally by construction and varies from map to map (from series to series). Also, the distribution of symbols is not directly coarse-grained from the distribution of visits to the different regions of the attractor. As a matter of fact, by construction the visibility algorithm makes use of the whole time series of (continuous) data to generate the degree sequence, and as there is no symbolization prior to this mapping, in principle fine-grained fluctuations are taken into account at all scales. Furthermore, no ad hoc phase space partition needs to be defined within the visibility method. This in principle further removes the second source of ambiguity found for standard symbolization. Moreover, note also that visibility algorithms are invariant under several changes of scale in the original time series [3,4], hence removing yet another source of ambiguity present in the symbolisation procedure. Finally, note that whereas the degree sequence can be understood as a time series (global) symbolisation, other structures involving more sophisticated properties of the visibility graphs (degree-degree correlations, spectral properties, etc) can be directly retrieved and hence, at least in principle this method can yield more (or other) information of the dynamical process not described by the symbolised series.
In what follows we compare the performance of the (multiplex) visibility methodology with the results obtained via symbolisation, in the context of diffusively coupled chaotic maps.

II. ORDER PARAMETER AND PHASE DIAGRAM IN DIFFUSIVELY COUPLED CHAOTIC MAPS: A COMPARISON BETWEEN THE MULTIPLEX AND THE TIME SERIES SYMBOLISATION APPROACHES.
We focus on the first system of five diffusively coupled chaotic logistic maps considered in the main text, whose evolution is described by We generate multivariate orbits {x(t)} N t=1 , x ∈ R 5 of size N and proceed to compute the interlayer averaged mutual information I HVG via the associated visibility multiplex, as defined in the main text. We then compare this proposed order parameter with an analogous measure I SYMB directly performed in the symbolized series. Without loss of generality, for a symbolized series with p integer symbols, let us first define the pairwise mutual information such that In [5], p = 2 and a homogeneous phase space partition was used, but other choices can be made as well. In what follows we compare the performance of I HVG and I SYMB as order parameters that encapsulate the rich dynamical transitions that the system of CMLs undergo as we increase the coupling constant. Concretely, we investigate I SYMB for (i) different number of symbols p, and (ii) different degrees of heterogeneity of the partition. Finally we investigate the effect of shortening the size of the multivariate series under study, N . Such assessment is based in three different criteria, namely (i) capacity of capturing the monotonic increase of synchronisation for weak coupling ( < 0.1, inside the Fully Developped Turbulent (FDT) phase), (ii) accuracy to detect major dynamical transitions (such as the transition from FDT to a randomly selected pattern (PS), or the FDT to Spatio-temporal Intermittency (STI)) through sharp changes in the order parameter, and (iii) accuracy to detect secondary dynamical transitions, such as the onset of multiband chaotic attractors inside the STI phase.
In what follows, we show that the visibility approach is at least as accurate as the optimal symbolisation (i.e., the optimal selection of symbols and partition). Since this optimal selection is not known a priori, we conclude that the visibility approach is, at least in the case under study, more efficient.
A. Effect of the number of symbols for a homogeneous phase space partition Let us first fix an homogeneous partition of the phase space [0, 1] d for which cells are hypercubes of dimension d with size 1/p, where p is the number of symbols considered, and let us explore the effect of varying p on the performance of I SYMB . Note at this point that, if we only consider the 'symbolization' aspect of the visibility algorithm, remind that noisy series have an associated HVG with a mean degree k = 4. Therefore we might consider that, in a first approximation, the HVG method should indeed be comparable to a symbolization with p ≈ 4. We summarise For ≈ 0.34, the chaotic attractor splits into different bands, and such secondary transition is captured by a sharp change in the visibility multiplex mutual information (blue). Note that in the case of symbolized series, both for small number of symbols (p = 2) as well as for p > 10, this transition is not captured. (Second row ) In the left panel, we plot the behavior of I SYMB with p = 2 in the weak coupling regime ( < 0.1), where the system displays spatio-temporal chaos (FDT phase) and the mutual information among maps increase monotonically (but weakly) with . We find that I SYMB is not monotonic in this region. In the right panel we plot the same results for different symbols p, recovering the correct monotonic increasing found with I HVG when the number of symbols is increased. (Third row ) In the left panel, we plot I SYMB with p = 100 symbols again the weakly coupling regime. The monotonic increasing shape is lost again if the number of symbols is too large. This is confirmed in the right panel, concluding that for large number of symbols I SYMB doesn't capture the subtle mutual information increase with and misleadingly predicts constant mutual information.
our results in figure 2. As expected, qualitatively speaking, the visibility approach yields similar results as the symbolisation method for p = 4, 5. For the standard symbolisation (the one used in [5]) with p = 2, symbolisation fails to accurately describe several properties such as the increase of synchronisation for weak coupling and the onset of multiband attractors inside STI. Also, for large values p > 10, these properties are not captured anymore by the symbolised series.    results dramatically depend on the selection of c and are summarised in table I and in figure 3. We haven't found, for the case p = 2, a proper partition that gathers better or equal results than the visibility method. Interestingly, the standard (homogeneous) phase space partition (c = 0.5), used in [5], is shown to be suboptimal when compared to other selections such as c = 0.7 − 0.8. Notice that the HVG method is quite robust against shortening, as for sizes of 2 8 we already find the correct profile, including a monotonically increasing function for weak coupling (FDT phase), the location of primary transitions, and the location of onset of multiband chaotic attractor. Similar behaviour is found for p = 5, whereas the effect of shortening the size is more acute for larger number of symbols.

C. Effect of time series size for different number of symbols
Finally, we study the effect of shortening the size of the time series. Usually, the shorter the time series, the smaller the upper bound for the number of symbols that can be used before the lack of statistics takes a predominant role and ruins the accuracy of the measurements. Hence the question, which is the 'optimal' number of symbols conditioned to the size of the series under study? Our results are summarised in the panels of figure 4. Roughly speaking, the robustness of both the visibility method and the symbolisation with p ≈ 5 (the optimal symbolisation according to our previous analysis) against series shortening is similar (the same qualitative results are found in both cases up to short series of size N = 2 8 data).

III. EDGE OVERLAP AS AN ALTERNATIVE MEASURE TO DESCRIBE PHASES
We can extend the analysis performed using the pairwise layer averaged mutual information to another multiplex measure. The so called average edge overlap o is defined as: where o ij is the overlap of the edges between node i and node j at the different layers, and K is the total number of pairs of nodes connected on at least one of the M layers [6]. Notice that o ij = 0 if a In figure 5 we plot the average edge overlap for the system of CMLs studied in the main body of the paper, as a function of the coupling constant . The shape of this purely multiplex measure is qualitatively similar to the pairwise layer averaged mutual information, although subtle dynamical transitions, such as the onset of multiband chaotic attractors, is not clearly seen via this particular metric.

IV. DETECTING THE ONSET OF SYNCHRONISATION IN GLOBALLY COUPLED CHAOTIC MAPS.
To round off the validation part presented in the main body of this paper, we consider here a system of Globally Coupled Maps (GCMs) [7]: ∀α = 1, . . . , M , where the dynamics of each unit is governed by the logistic map f (x) = µx(1 − x), µ ∈ [0, 4]. This system can be indeed considered as a mean-field version of CMLs. In particular we consider values of µ for which each map is chaotic, i.e. µ > µ ∞ = 3.56995..., excluding periodic windows. In those cases, complete synchronization of the GCM is only reached when the system is in the simplest chaotic attractor, where x [α] = x [β] ∀α, β, and the dynamics effectively reduces to that of a single logistic map. It can be proved [7] that this is indeed the stable regime for those values of for which the Lyapunov exponent λ 0 of one isolated logistic map satisfies the inequality λ 0 + log(1 − ) < 0. The onset of complete synchronization is therefore reached at a critical value of the coupling strength c = 1 − exp(λ 0 ) −1 . For example at µ = 4, by making use of the the analytic expression λ 0 = log 2, we get that the onset of complete synchronization occurs at c = 1/2, whereas for other values of µ the value of c can be derived from the numerical evaluation of λ 0 . The multiplex visibility graph approach is able to predict the position of the onset of complete synchronization. We conjecture that the average mutual information I attains its maximum at the onset of complete synchronization, and accordingly we propose the quantity multiplex c ≡ min (argmax(I( ))) as a measure of c . In Fig. 6 we plot multiplex c versus c for a system of 5 globally coupled chaotic maps and for different values of µ ∈ [µ ∞ , 4], finding a remarkable agreement in every case.

V. ANALYSIS OF MULTIVARIATE FINANCIAL SERIES
We have analyzed a large dataset of financial stocks comprising stock evolution of the 35 major american companies from the New York Stock Exchange (NYSE) and Nasdaq in the period 1998-2012, the majority of which belong to the Dow Jones Industrial Average (see table II for a detailed list of all the companies). The NYSE is the largest and most liquid cash equities exchange in the world by market capitalization. The series have very high resolution (one data per minute), yielding O(2 · 10 6 ) data per company.