In the last two decades the physicists have become more and more interested in complex systems with hierarchical structure extending over many scales costa11 . They are modeled by graphs including modules of different sizes which are grouped in a hierarchical manner baraba03 . A time series generated by such a complex system contains in a compressed form all the intricacies of its complex evolution. Hence one needs special algorithms to identify, separate, and analyze the evolution at different scales of a complex multiscale time series gao07 .

Usually the spectral analysis is realized in a nonadaptive manner, the time series being projected on a predetermined frequency grid, most often equidistant. This is the case even with the most popular multiresolution analysis supplied by the wavelet decomposition percival00 . A better approach would be a data-driven numerical algorithm which should separate several intrinsic spectral components arranged in an increasing hierarchy of frequencies. For a correct characterization of the dynamics over different time scales, an essential issue is the proper separation of time scales between the intrinsic components. Otherwise, we risk to assign to a specific time scale a mixture of unrelated dynamics. The only algorithm known to us which adapts to the characteristics of nonlinear and nonstationary time series is the empirical mode decomposition (EMD) method huang98 . But the EMD components are not constrained with respect to their spectral overlapping.

In this paper we present a numerical algorithm which is not only adaptive, but also decomposes a complex multiscale time series into several intrinsic components with minimum spectral superposition. Our decomposition algorithm is based on the condition that the spectral width of each intrinsic component should be minimum and disjoint from that of the neighboring components.

We test our decomposing algorithm by applying it to a financial time series which is one of the most known example of hierarchical multiscale time series dacoro2001 . In financial markets the time scales range continuously from minutes, where the market makers act, up to several years, intervals characteristic for investment funds. All these investors interact with each other. Asset prices are also influenced by many exterior economic and financial factors with even larger time scales (see for example Chap. 9 in tsay2010 ). Our decomposition algorithm allows us to show that the financial volatility contains a hierarchy of time scales with a heterogeneous scale structure and different dynamical properties. The small-scale intrinsic components have a random nature while the large-scale intrinsic components have a deterministic behavior rendering the financial time series nonstationary.

The statistical mechanics of the complex systems with a hierarchy of temporal or spatial scales is given by a superposition of several statistics, known as superstatistics beck03 ; beck11 . The simplest form of the superstatistical model has only two separate scales with independent random dynamics. For example, the small-scale dynamics may describe a system in local thermodynamical equilibrium with a Boltzmann-Gibbs factor $e^{-\beta E}$. If the intensive parameter $\beta$ slowly fluctuates at much larger scales, then the statistics of the whole system is given by the mixture of the microscopic statistics

where $f(\beta )$ is the probability distribution function (pdf) of the parameter $\beta$. The resulting statistics generally have non-Gaussian properties with heavy tails beck03 . The essential supposition in this approach is the independence of the two scales, not their number. If there are several independent time scales, the integral in Eq. (1) is replaced by a multiple integral sobyanin11 .

Many complex systems in various areas of science have been successfully modeled by this simple superstatistical model: turbulence, astrophysics, quantitative finance, climatology, random networks, medicine, etc. beck09 . Remarkable results are obtained when the local equilibrium hypothesis holds as, for example, in turbulence beck07 . But for biological or social systems this hypothesis is unlikely to stand all the time and then the global pdf is no more a simple integral as in Eq. (1) because we have to take into account the interactions between the different time scales. In this paper we show that in most of the time the random small-scale intrinsic components of the financial volatility have an independent evolution, but they become synchronized during financial crashes. Therefore a financial market seems to be a more complex system than a statistical system in local equilibrium in a slowly evolving environment.

The basic quantity in quantitative finance is the logreturn defined as $r_n=\log (S_{n+1}/S_n)$, where $\{S_n,1\le n\le N+1\}$ are the prices of a financial asset. The basic properties of daily financial returns, the so-called stylized facts cont01 , are captured by the heteroscedastic model

where $\{z_n\}$ is a Gaussian white noise with zero mean and unit variance and $\{v_n\}$ is a stochastic process called volatility with strictly positive values and independent of $\{z_n\}$ taylor07 . It is a particular case of the superstatistical model. If in Eq. (2) the volatility fluctuates much slowly than the Gaussian noise, then the pdf of the logreturns is a mixture of Gaussians

where $f(v)$ is the pdf of the volatility. The superstatistical model has been successfully applied to several financial time series straet09 ; gerig09 ; vamos2010 ; vamos2013 .

In Fig. 1(a) we plot the daily returns computed from the daily closing values of the S&P500 index from 1 January 1950 to 28 July 2012, a time series containing $N=15724$ values. The volatility is obtained as the moving average (MA) of the absolute values of the logreturns liu99

where $K$ is the semilength of the averaging window and . If $n⩽K$ ($n>N-K$), then the average is taken over the first $n+K$ (the last $N-n+K+1$) values of $|r_n|$. The properties of this MA with normalized padding with zero at boundaries are analyzed in vamos2012 .

Figures 1(b) and (c) show the volatility and the noise $z_n=r_n/v_n$ estimated by an MA with $K=3$ such that the small scale fluctuations of the volatility are preserved vamos2010 ; vamos2013 . The estimated noise is Gaussian [Fig. 2(a)] but it is uncorrelated only for lags larger than 5 days [Fig. 2(b)] due to the spurious correlations introduced by the MA (see Appendix B in vamos2012 ). Then the estimated volatility is accurate only for time scales larger than 5 days and its intrinsic components may be determined only for these time scales.

To measure the components spectral width we need a flexible local measure suited for a complex time series. Time scales are usually defined as the periods of the constant-amplitude trigonometric components of the Fourier decomposition. Instead of such a global averaged description, we use the characteristic times defined as the spacings between successive local extrema of the time series, also equal to the lengths of its monotonic segments huang98 . For an arbitrary time series we denote the characteristic times by ${\lambda }_j$, $1\le j\le J$. In the simplest cases, the set $\{{\lambda }_j\}$ supplies a full description of the time series spectrum. For instance, a harmonic oscillation has all ${\lambda }_j$ equal to half the oscillation period.

The situation is more complex in the case of a noisy time series. The characteristic times distribution is dominated by the small-scale fluctuations of the noise and it does not distinguish the large-scale variations possibly existing in the time series. Quantitatively, a time series cannot contain monotonic variations smaller than ${\lambda }_{\min }=\min \{{\lambda }_j\}$, but there may exist components with monotonic variations lasting longer than ${\lambda }_{\max }=\max \{{\lambda }_j\}$. For instance, the histogram of the characteristic times of the estimated volatility of the S&P500 index is dominated by fluctuations with characteristic times smaller than ${\lambda }_{\max }=14$ [Fig. 2(c)], which are much smaller than the characteristic times of the large-scale components of the volatility. Our goal is to decompose the volatility into a hierarchy of spectral components with characteristic times larger than 5 business days, such that the intervals $I_{\lambda }=[{\lambda }_{\min },{\lambda }_{\max }]$ associated to them should be disjoint.

We will decompose the volatility into components $v_n^{(c)}$

where $V$ is the mean volatility and $w_n^{(C)}$ are the volatility residuals of order $C$. This is a generalization of the superstatistical model with two independent scales described in the previous section. First, the number of the time scales is not limited to two. If they are independent, then the integral in Eq. (3) is replaced by a multiple integral sobyanin11 . Here we allow that the evolutions at different scales may be correlated and then the global pdf cannot be computed as an integral. Also, the residuals $w_n^{(C)}$ may be non-Gaussian. Finally, we make no assumption on the nature of the components (deterministic or random). In this section we describe the numerical decomposition algorithm and in the next section we exemplify how it works in the case of the S&P500 index.

Similar to the SiZer algorithm chaudhuri99 , we compute the components $v_n^{(c)}$ by successive smoothings, which gradually eliminate the variations with the largest characteristic times from the previous component. Let us assume that the first $c$ components have been computed. We denote by $I_{\lambda }^{(c)}$ the interval containing its characteristic times ${\lambda }_j^{(c)}$ ($1\le j\le J^{(c)}$) equal to the lengths of the monotonic segments. The next component is determined such that the interval $I_{\lambda }^{(c+1)}$ should have minimum width and should be as close as possible to the previous interval $I_{\lambda }^{(c)}$. In this way the time scales of the components are not predetermined, but they are controlled by data themselves and the resulting algorithm has an adaptive nature.

For smoothing we use the repeated central MA (RCMA) with normalized padding with zero at boundaries (Chap. 4 in vamos2012 ) obtained by repeating the smoothing (4) several times. Its parameters are the semilength $K$ of the averaging window and the number of repetitions $i$ of the MA. Roughly the same smoothing is realized by averaging many times with a small $K$ or several times with a large $K^{\prime }>K$ vamos2012 . Therefore we maintain $i$ between 15 and 20 repetitions and we vary mainly the semilength $K$ from one component to the next.

We denote by $\{u_n^{(c+1;1)}\}$ the time series obtained after one application of the RCMA to the residuals of order $c$, $w_n^{(c)}=v_n/(V{v}_n^{(1)}\mathrm{\dots }{v}_n^{(c)})$ and by $\{u_n^{(c+1;2)}\}$ that obtained after a second application to $w_n^{(c)}/u_n^{(c+1;1)}$. It is possible that a single smoothing with the RCMA cannot extract all the variations with the characteristic times in $I_{\lambda }^{(c+1)}$. If the most part of the monotonic segments with characteristic times in $I_{\lambda }^{(c+1)}$ has not been extracted, i.e., if $\Vert u_n^{(c+1;2)}\Vert /\Vert u_n^{(c+1;1)}\Vert >0.1$, where $||\cdot ||$ is the quadratic norm, then we apply once again the RCMA to the new partial residuals $w_n^{(c)}/(u_n^{(c+1;1)}{u}_n^{(c+1;2)})$. We stop the averagings after $k$ iterations when the quadratic norm of $u_n^{(c+1;k+1)}$ is one order of magnitude smaller than that of the product $u_n^{(c+1;1)}{u}_n^{(c+1;2)}\mathrm{\dots }{u}_n^{(c+1;k)}$.

After several applications of the RCMA, there remain some monotonic segments with reduced lengths and amplitudes which are not related to the characteristic times within of the $c+1$ component. A large-scale monotonic variation may be separated into two shorter segments with the same monotony by a small inflexion with opposite monotony. We replace these three segments by a single monotonic segment with the length equal to the sum of the three lengths using the ACD algorithm designed to determine the monotonic component of an arbitrary time series vamos2007 ; vamos2012 . In this way the three initial characteristic times are replaced by a single one equal to their sum so that the interval $I_{\lambda }^{(c+1)}$ is shifted towards the larger values of the characteristic times. We eliminate the small monotonic segments in ascending order of their lengths until $I_{\lambda }^{(c+1)}$ does not overlap $I_{\lambda }^{(c)}$. We choose the parameters of the RCMA and the number of the ACD iterations such that the resulting interval $I_{\lambda }^{(c+1)}$ has a minimum width.

Now we check if the random intrinsic components of the volatility are independent, as the superstatistical model (3) assumes. We can destroy the possible correlations between the three components of $v_n^{(\mathrm{rand})}$ by shuffling their values and then multiplying them. The majority of the values of the initial and the shuffled pdfs coincide, proving that most of the time the three random intrinsic components fluctuate independently [Fig. 9(a)]. Their tails differ from each other only for $v_n^{(\mathrm{rand})}>3$, for 120 values representing 0.76% of all the volatility values. Hence more than 99% of the time the random intrinsic components fluctuate independently and the central part of the resulting pdf can be computed as a multiple integral generalizing Eq. (3). We obtained the same result for two other financial indices, DJI and NIKKEI.

It is well-known that volatility has heavy tails liu99 . The question is if the tail of the volatility random part $v_n^{(\mathrm{rand})}$ remains heavy. A pdf with heavy tails has an asymptotic behavior following a power law and its moments of order higher than the absolute value of the exponent are infinite resnick2007 . Because the existing statistical methods to recognize and to measure power law tails have significant errors resnick2007 , we compare the volatility tails with the tails of the Laplace distribution $e^{-\mu |x|}/(2\mu )$, where $\mu$ is a positive real parameter. We use this distribution as a reference because, unlike the power law distributions, all its moments are finite although it has heavier tails than the Gaussian distribution.

We need several thousands of observations to establish if a sample pdf has heavy or exponential tails heyde04 , condition satisfied by our volatility time series. We have numerically generated 1000 uncorrelated time series with Laplace distribution of unit variance ($\mu =1/\sqrt{2}$) having the same length as the volatility time series. Then we have computed the cumulative distribution function (cdf) for each artificial time series and we have established the limits containing 95% of the complementary cdf (1 - cdf) [dashed lines in Fig. 9(b)]. We have also generated 1000 shuffled versions of $v_n^{(\mathrm{rand})}$, we have standardized them by extracting their mean and then dividing the result by their standard deviation, and we have established the 95% interval of variation of their complementary cdfs [continuous lines in Fig. 9(b)].

The right tail of the volatility random part $v_n^{(\mathrm{rand})}$ is significantly heavier than those of the Laplace distribution [Fig. 9(b)]. This tail is reduced by shuffling, but the tails of the standardized shuffled time series remain almost in all cases heavier than the exponential tails. Hence, even if the random intrinsic components of the volatility were independent, their product would have heavy tails. The correlation between them only enhances their heaviness. In conclusion, the extreme volatility is the result of both the superstatistical mechanism (mixture of different statistics) and the correlation between the random volatility components.

The random intrinsic components of the volatility become correlated during financial crashes when the financial markets structure is reorganized and all the categories of investors act in the same way, i.e., they become synchronized onnela2003b ; mcdonald2008 ; peron2011 . Synchronization of dynamically interacting units in a complex network was observed in a large variety of domains arenas2008 . Hence we have shown that the synchronization of the network structure is accompanied by the temporal synchronization of the dynamics at different time scales.

There are many possible applications of our adaptive decomposition algorithm due to its ability to separate nonlinear and nonstationary time series into disjoint intrinsic spectral components. They can be analyzed by usual statistical tests and the possible spectral heterogeneities may be discovered. For instance, the existence of the slowly varying trend in a noisy time series could be separated from a random noise with characteristic times much smaller than those of the trend. Also, one could analyze the intrinsic components of the concentration in an anomalous diffusion process generated by the superposition of a microscopic Brownian motion over a transport through a random porous medium suciu .

Here we have shown that the well-separated spectral intrinsic components of the financial volatility allow the search for correlations between the dynamics at different scales. Generalizing this result we conjecture that the temporal multiscale synchronization of the communities of a complex network may be a general mechanism to generate heavy tails different from those proposed till now as the combination of exponentials, the Yule process, critical phenomena, self-organized criticality newman2005 , the superstatistical model beck03 or the mixture of Gaussians gherghiu04 . The fact that the volatility heavy tails result from the simultaneous action of two distinct such mechanisms indicates that any of these mechanisms could collaborate to generate heavy tails when the system has a hierarchical complex structure.

In accordance with Eq. (2), the complex structure of the volatility is conveyed to logreturns. The heterogeneity of the intrinsic components and the synchronization of small-scale components might provide an explanation for the difficulty to model the logreturn distribution voit2005 . In particular, even if the financial returns have a universal tail index around 3 gopi99 , there exist significant variations between individual companies plerou99 . These variations are small for intraday returns, but they increase for daily returns. The variations are even larger for exchange rates between different currencies dacoro2001 . Such variations can be analyzed in relation with the complex evolution of the intrinsic components obtained with the adaptive decomposition algorithm presented here.