Exact Multivariate Amplitude Distributions for Non-Stationary Gaussian or Algebraic Fluctuations of Covariances or Correlations
EExact Multivariate Amplitude Distributions forNon–Stationary Gaussian or Algebraic Fluctuationsof Covariances or Correlations
Thomas Guhr and Andreas Schell
Fakult¨at f¨ur Physik, Universit¨at Duisburg–Essen, Duisburg, GermanyE-mail: [email protected]
Abstract.
Complex systems are often non–stationary, typical indicators arecontinuously changing statistical properties of time series. In particular, thecorrelations between different time series fluctuate. Models that describe themultivariate amplitude distributions of such systems are of considerable interest.Extending previous work, we view a set of measured, non–stationary correlationmatrices as an ensemble for which we set up a random matrix model. We use thisensemble to average the stationary multivariate amplitude distributions measuredon short time scales and thus obtain for large time scales multivariate amplitudedistributions which feature heavy tails. We explicitly work out four cases, combiningGaussian and algebraic distributions. The results are either of closed forms or singleintegrals. We thus provide, first, explicit multivariate distributions for such non–stationary systems and, second, a tool that quantitatively captures the degree of non–stationarity in the correlations.
Keywords : non–stationarity, multivariate statistics, algebraic heavy tails, randommatrix theory
Submitted to:
J. Phys. A: Math. Gen.
1. Introduction
A wealth of complex systems show non–stationarity as characteristic feature, i.e. theylack any kind of equlibrium [1, 2, 3, 4]. Non–stationarity shows up in many differentways [5, 6], for example, the increasing share of wind energy fed into the power gridis greatly appreciated from an environmental viewpoint, but causes problems for thestability of the grid due to non–stationarity and the lack of predictability for wind speedand wind direction. In electroencephalography (EEG) electrical currents are recorded atdifferent positions on the scalp to measure the brain activity. The correlations stronglydepend on the overall state of the brain [7, 8]. Wave packets traveling through disordered a r X i v : . [ c ond - m a t . d i s - nn ] N ov xact Multivariate Amplitude Distributions i.e. modelsbased on the assumption of randomness for some parts or aspects of the systems? —Sometimes other approaches, e.g. from many–body physics carry over or provide usefulinspiration. This is especially so for the interplay between coherent, collective motionand incoherent motion of the individual particles [23, 24, 25].Here, we have two goals. First, we construct analytical results for the multivariatedistributions of amplitudes, measured as time series in correlated, non–stationarysystems. Second, we thereby provide quantitative measures for the degree of non–stationarity in the correlations. We use the word “amplitude” to refer to any observable,directly measured or inferred from a measurement, at a given position and a given time.The word “position” is meant in a general sense, a geographical point, a location,a specific stock and so on. The amplitudes at all times of measurement at a givenposition form the corresponding time series. The challenge we master is to capture thenon–stationarity of the correlations between these time series in terms of a statisticalensemble. We considerably extend a random matrix approach to tackle these issues thatwe presented some years ago [26], see also some modification in Ref. [27]. It also turnedout useful in the context of credit risk and its impact on systemic stability [28, 29].Relations to other systems were discussed in Ref. [30]. In Ref. [26], we showed that thefluctuations of the correlations lift the tails of the multivariate amplitude distributions,making them heavy–tailed. Our new results extend that.Our model may be viewed as a justification of compounding or mixture techniquesin mathematics and superstatistics in physics [31, 32, 33, 34, 34, 35, 36]. By tracingheavy–tailed distributions back to fluctuations of the correlations, we explain such ad–hoc approaches. However, our previous studies were based on a Gaussian assumptionfor the statistical ensemble. To avoid misunderstandings, we underline that thisGaussian assumption does not at all imply Gaussian form of the multivariate amplitudedistribution, rather, its tails are exponential [26]. Here, we extend this by alsoconsidering ensembles of algebraic covariance or correlation matrices in the form ofcertain determinants. This is relevant, as many complex systems, not only financialmarkets [37], are known to have algebraically heavy–tailed distributions. Our earlierresults [26] show a discrepancy from the exponential behavior way out in the tailsbetween data and model which is due to the above mentioned Gaussian assumption. xact Multivariate Amplitude Distributions
2. Posing the Problem and Developing the Strategy
To set the tone and to introduce our notation, we briefly review basics on covariancesand correlations in Sec. 2.1. After a discussion of non–stationarity in Sec. 2.2, we presentthe key ideas of our random matrix model in Sec. 2.3 and choose in Sec. 2.4 explicitforms of the distributions which are the ingredients of our model. xact Multivariate Amplitude Distributions Suppose we measured amplitudes A k ( t ) at K positions labeled k = 1 , . . . , K and atequidistant times t = 1 , . . . , T . These amplitudes can be any type of data, temperatures,water levels, electrical voltages, stock prices and so on. The term “position” is used ina very general way, it can be a geographical point, when e.g. water levels of riversare measured or a location on the scalp of a person in EEG measurements, or in anabstract sense a company when stock prices are considered. We always assume thattime is normalized such that t is an integer and T is the total number of points in time.Keeping with a common notation, we write the postition as an index and the time as theargument of the data A k ( t ), other notations such as A ( k, t ) or A kt would be equivalent.Our data may be ordered in the rectangular K × T matrix A = A (1) A (2) · · · A ( T ) A (1) A (2) · · · A ( T )... ... . . . ... A K (1) A K (2) · · · A K ( T ) . (1)The rows of this matrix are the well known time series A k ( t ) , t = 1 , . . . , T which containthe amplitudes at the same position for all times. Importantly, the columns may beinterpreted in a similar, but dual, way: they are the position series which collect theamplitudes at the same time for all positions. These two types of series provide differentinformation, particularly on covariances and correlations. To obtain the covariances ofthe time series, we first have to substract their mean values in time, (cid:104) A k (cid:105) T = 1 T T (cid:88) t =1 A k ( t ) , (2)which yields the normalized time series B k ( t ) = A k ( t ) − (cid:104) A k (cid:105) T , t = 1 , . . . , T (3)that are the rows of the K × T matrix B . The covariances of the time series are thenthe normalized scalar productsΣ kl = (cid:104) B k B l (cid:105) T = 1 T T (cid:88) t =1 B k ( t ) B l ( t ) . (4)In the statistics literature, one often uses the prefactor 1 / ( T −
1) and refers to thecovariances (4) as biased. In physics, the present choice seems more common. The Σ kl are the elements of the K × K sample covariance matrix of the time seriesΣ = 1 T BB † . (5) xact Multivariate Amplitude Distributions σ k = Σ kk with σ = diag ( σ , . . . , σ K ) (6)being a diagonal matrix. All variance are positive definite and C = σ − Σ σ − (7)is the correlation matrix of the time series.However, when studying the covariances of the position series, we must applyanother normalization employing the mean values of the columns, (cid:104) A ( t ) (cid:105) K = 1 K K (cid:88) k =1 A k ( t ) , (8)leading to the normalized position series (cid:101) B k ( t ) = A k ( t ) − (cid:104) A ( t ) (cid:105) K , k = 1 , . . . , K . (9)These are the colums of the K × T matrix (cid:101) B . The covariances of the position series arethe normalized scalar productsΞ( t, s ) = (cid:104) (cid:101) B ( t ) (cid:101) B ( s ) (cid:105) K = 1 K K (cid:88) k =1 (cid:101) B k ( t ) (cid:101) B k ( s ) (10)which form the T × T covariance matrix of the position seriesΞ = 1 K (cid:101) B † (cid:101) B . (11)Introducing the variances of the position series ξ ( t ) = Ξ( t, t ) , and ξ = diag ( ξ (1) , . . . , ξ ( T )) , (12)we may define D = ξ − Ξ ξ − (13)as the correlation matrix of the position series. The covariances Σ kl quantify thedifferences or similarity of two different time series at equal times, while the covariancesΞ( t, s ) give information on what happens at the same positions for different times andthus on possible non–Markovian behavior of the system under consideration. By theirvery definition, Σ and Ξ and thus C and D are positive semi–definite, they are positivedefinite for K ≤ T or T ≤ K , respectively. For uncorrelated time series, one has C = K , i.e. the K × K unit matrix, Markovian data are characterized by D = T . It is worthmentioning that the covariance matrices Σ and Ξ have physical dimensions, namely the xact Multivariate Amplitude Distributions C and D are dimensionless. As this will become relevant in the later discussion, we continueto work for the time being with both, the covariances and the correlations.In the statistics literature, the sample covariance and correlation matrices areviewed as estimators of the population covariance or correlation matrices. The ideabehind is that every sample is a subset of the full data set, e.g. the income distributionfor a whole country is not obtained by asking everyone who receives a salary, it isestimated from a sample, i.e. from a subset of the population that is consideredsufficiently large. This issue, however, is only of limited relevance in the present context.We have systems in mind for which all time series are at our disposal, i.e. the samplematrices yield very reliable estimators. Rather, it is the issue of non–stationarity thatwe will tackle here. In a non–stationary complex system, crucial parameters or distributions of observableschange in an erratic, unpredictable way over time. Among the examples mentionedin Sec. 1, correlations in financial markets are particulary illustrative. Figure 1 showstwo large correlation matrices of stock prices for companies in the Standard and Poor’s
Figure 1.
Correlation matrices of K = 306 companies in the Standard and Poor’s500 index for the fourth quarter of 2005 and the first quarter of 2006, the darker, thestronger the correlation. The companies are sorted according to industrial sectors.Taken from Ref. [26].
500 index ordered according to the industrial sectors. The time series were measuredin subsequent quarters. Two important features are visible. First, the matrices lookclearly different, because the business relations between the companies and the marketexpectations of the traders change in time. Second, the coarse structure remainssimilar, indicating some stability of the industrial sectors. Subsequently measuredcorrelation matrices for other complex systems would have similar properties. Thequestion is how the non–stationarity, i.e. the fluctuations of the correlations on shortertime scales influence the statistics of a complex system studied over long time scales. xact Multivariate Amplitude Distributions r k , k = 1 , . . . , K ordered in the K component vector r = ( r , . . . , r K ). In our notation,we distinguish the amplitudes A k ( t ), measured at a particular time t from the amplitudes r k , sampled over some time interval. As in our previous study [26], we assume that wemay divide the data measured over a large time scale T tot into epochs of length T ep within which the system is to a good approximation stationary, while it changes fromepoch to epoch as in the example of the financial market in Fig. 1. This time scaleshould be chosen in such a way that the number N ep = T tot T ep (14)of epochs is integer. While a proper determination of the time scale T ep is essentialfor data analysis, we emphasize that this time scale does not directly enter our modelto be set up in the sequel. Figure 2 illustrates our construction which interprets the Figure 2.
Large time interval of length T tot , divided into N ep epochs of length T ep . non–stationary system studied over the long time T tot as assembled of N ep subsequentsystems which are all approximately stationary on their time scale T ep . Random Matrix Theory usually relies on a concept which is sometimes referred to assecond ergodicity, these are the models in class (I), as defined in Sec. 1. Ergodicity is anindispensable feature of statistical mechanics, as it states the equivalence of time andensemble averages. It is worth underlining that the ensembles in statistical mechanicsare fictitious, a powerful mathematical construction to facilitate, by means of ergodicity,the physically relevant time average. When using random matrices to model spectralstatistics of individual systems such as a chaotic billiard or the celebrated Hydrogenatom in a strong magnetic field, a very similar line of reasoning is applied. Empirically,the spectral statistics of the individual system is obtained by spectral averages of onesingle spectrum which, to allow for a meaningful statement, has to contain a largenumber of levels or resonances. Second ergodicity, first proven in Ref. [53], statesthat such a spectral average is equivalent to an average over an ensemble of randommatrices, provided their dimension is very large. Again, this ensemble is fictitious as xact Multivariate Amplitude Distributions N ep different matrices Σ ep or C ep . We may view either set of N ep matrices as a matrixensemble, and our goal is to model it by an ensemble of random matrices. Thisobservation will lead us to a random matrix model in class (II). We model the trulyexisting ensemble of positive definite matrices Σ ep or C ep by random ones,Σ ep −→ N XX † or C ep −→ N XX † , (15)where the form of the right hand sides ensures positive semi–definiteness. The randommatrices X , the model data matrices, have dimension K × N where the number of rows K is determined by the fact that the matrices Σ ep and C ep have dimension K × K . Thenumber of columns N , however, is not fixed by the relation (15), it may be interpretedas length of the random model time series which form the rows of X . It is a freeparameter in our construction and is, at the end of our calculations, not even restrictedto integer values. To provide already now an intuitive and plausibel interpretation andanticipating the later discussion, we mention that the amount of information on a systemgrows with T , the longer the time series, the lower are the fluctuations when measuringaverages. Accordingly, N in our model characterizes the strength of the fluctuations inthe ensemble of the matrices Σ ep or C ep around the matrix mean values Σ tot or C tot ,respectively, for the long time interval. We notice that the truly measured matrices arestrongly correlated, the stronger, the closer in time their epochs are. Our model is notmeant to capture those correlations directly, rather, it is set up to model the fluctuationsof the matrices Σ ep or C ep in different epochs around their mean values, which implicitlytakes care of the mentioned correlations across the epochs. This is fully consistent withthe basic idea of statistical mechanics and in the present context best understood whencomparing with chiral random matrices.Motivated by our data analysis to be published elsewhere, we assume that theanalytical forms of the multivariate amplitude distributions p ( r | Σ ep ) or p ( r | C ep ) for agiven epoch have the same functional form for all epochs. Nevertheless, they differ fromepoch to epoch. We assume that this variation may be fully captured by the positivedefinite matrices Σ ep , which have the dimension of the amplitudes squared, or C ep , whichare dimensionless, respectively, and differ from epoch to epoch. Importantly, these xact Multivariate Amplitude Distributions p ( r | Σ ep ) or p ( r | C ep ) have Gaussian form. For other functional forms, Σ ep or C ep cannotdirectly be obtained in this way, but it is essential to construct our model in such a waythat they are related to the sample covariance or correlation matrices and it must bepossible to determine them from those.In the sequel, we focus on the correlation matrices, but later on we will showhow everything easily carries over to covariance matrices. We draw the random datamatrices X and thus the model matrices XX † /N from a distribution w ( X | C, D ). Itparametrically depends on a positive definite matrix C which, similary to the abovediscussion, only coincides with the K × K sample correlation matrix of time seriesmeasured over the long time T tot for a Gaussian choice of w ( X | C, D ), but has to berelated to the sample correlation matrix in other cases. It is an intrinsic property ofour model that the second input matrix, the matrix D , modeling the correlations of theposition series, has dimension N × N . The crucial concept now consists of carrying outthe random matrix ensemble average (cid:104) p (cid:105) ( r | C, D ) = (cid:90) p (cid:18) r (cid:12)(cid:12)(cid:12) N XX † (cid:19) w ( X | C, D )d[ X ] , (16)in which the replacement (15) is used in the multivariate distributions p ( r | C ep ). Theinvariant measure or volume elementd[ X ] = K (cid:89) k =1 N (cid:89) n =1 d X k ( n ) (17)is simply the product of the differentials of all independent matrix elements. Theensemble average (16) yields the distribution for all amplitude data measured over thelong time T tot subjected to the non–stationarity of the covariances. Motivated by data analyses, we choose two forms of the multivariate amplitudedistribution in each of the epochs. As in Ref. [26], we employ the multivariate Gaussian p G ( r | C ep ) = 1 (cid:112) det 2 πC ep exp (cid:18) − r † C − r (cid:19) (18)and, as a new choice to model heavy tails, the algebraic distribution p A ( r | C ep ) = α K lm (cid:18) m r † C − r (cid:19) l α K lm = (cid:114) m K Γ( l )Γ( l − K/
2) 1 (cid:112) det 2 πC ep (19) xact Multivariate Amplitude Distributions l and m . It is normalizable if l > K/
2, thenormalization constant α K lm depends on the number of the positions K , i.e. on thedimension of the amplitude vector r and on l and m , see Sec. 3.1. Not all moments of thisdistribution exist, and the convergence has to be guaranteed when doing integrals. Thisimposes conditions on the parameters, we come back to this point. We notice that bothdistributions depend on the quadratic form r † C − r which is known as the Mahalanobisdistance [56]. Importantly, the distribution (19) converges to the Gaussian (18), whenboth parameters l and m are taken to infinity under the conditionlim l,m →∞ ml = 2 . (20)In the model, the expectation value (cid:104) rr † (cid:105) serves as an estimator for the samplecovariances or correlations. We find (cid:104) rr † (cid:105) Y = β Y C ep (21)with β Y = , if Y = Gm l − K − , if Y = A . (22)Due to its very definition, we have β G = 1 for the multivariate Gaussian, but a differentvalue β A for the algebraic distribution. The relation in the latter case then suggests theuseful fixing m = 2 l − K − l and m . Only with this choice, C ep can be estimated by thesample correlation matrix, otherwise only up to some factor. The distributions (18) and(19) are not sensitive to non–Markovian effects, even if they are in the data. This isso, because, in a data analysis, they are obtained by first sampling all amplitudes at all T ep times, and then lumping together all these T ep distributions. Similarly, one couldalso consider a distribution (cid:101) p ( (cid:101) r | D ep ) of the position amplitudes (cid:101) r = ( (cid:101) r (1) , . . . , (cid:101) r ( T ep )),which is tested by sampling for all K positions, and then lumping together. Such adistribution is only sensitive to non–Markovian effects, but not to the correlations ofthe time series in the usual sense.Both types of correlations, of time and position series, are accounted for in thedistributions that we now choose to model the truly existing ensemble of correlationmatrices, fluctuating around the positive definite matrices C and D . They generalizethe amplitude distributions (18) and (19) from N = 1 to arbitrary N . A natural choicefor the distribution of the random data matrix X is the multimultivariate Gaussian w G ( X | C, D ) = 1 √ det 2 πD ⊗ C exp (cid:18) −
12 tr D − X † C − X (cid:19) , (24) xact Multivariate Amplitude Distributions D ⊗ C as the N × N matrix with the K × K matrixentries D nm C . The distribution (24) is known in the literature as doubly correlatedWishart distribution [57, 58, 59, 60, 61]. In Refs. [26, 27] we only considered theMarkovian case, i.e. , D = N , here we go beyond that by allowing arbitrary positivenon–Markovian correlation matrices D . To take care of heavy tails in the distributionof the truly existing ensemble, we also choose the algebraic, determinant distribution w A ( X | C, D ) = α KNLM det L (cid:18) N + 1 M D − X † C − X (cid:19) α KNLM = (cid:114) M KN N (cid:89) n =1 Γ( L − ( n − / L − ( K + n − /
2) 1 √ det 2 πD ⊗ C (25)depending on two shape parameters L and M . This distribution is a generalization ofthe ones introduced by Forrester and Krishnapur [62] and by Wirtz, Waltner, Kieburgand Kumar [63]. Existence of the normalization integral is guaranteed if L > K + N − . (26)Formally, the amplitude distribution (19) is a special case for N = 1 and D = 1. Thenormalization constant α KNLM includes the normalization constant α K lm for N = 1,its derivation is sketched in Sec. 3.1. Accordingly, the distribution (25) converges to theGaussian (24) if L and M are taken to infinity under the conditionlim L,M →∞ ML = 2 . (27)The first matrix moment is given by (cid:28) N XX † (cid:29) Y = B Y C , (28)where B Y = , if Y = GM L − K − N − , if Y = A . (29)In the Gaussian case, the confirmation of the value B G = 1 is natural, see Ref. [60].However, in the algebraic case, the first marix moment only exists if L > K + N + 12 . (30)As to be expected, the lower threshold is larger by one as in the condition (30) forthe existence of the normalization. Of course, the value B A is different from B G = 1.Because the matrix moments in the algebraic case are also of interest in a general xact Multivariate Amplitude Distributions M = 2 L − K − N − . (31)Only by such a fixing, which puts B A to one, we guarantee that the first matrix momentof the distribution (25) is equal to C , facilitating the comparison of tails with theGaussian case. If instead of the correlation matrices C and D the covariance matricesΣ and Ξ are employed, we have to replace C in the relation (28) with Σ and to multiplyits right hand side with tr Ξ /N [60, 64].Two observations are worth mentioning. First, the positive definiteness of C and D allows us to write the matrix in the trace and in the determinant of Eqs. (24) and (25)in real–symmetric, positive semidefinite form, D − X † C − X −→ ( D − / X † C − / )( C − / XD − / ) (32)= ( C − / XD − / ) † ( C − / XD − / ) . (33)We will use this in some of the calculations to follow, we then write, for Y = A, G , w Y ( (cid:101) X † (cid:101) X ) with (cid:101) X = C − / XD − / instead of w Y ( X | C, D ). As the distributions ofthe random data matrices only involve matrix invariants, we also have the remarkableidentities w Y ( (cid:101) X † (cid:101) X ) = w Y ( (cid:101) X (cid:101) X † ) , Y = G, A , (34)where the determinant in the algebraic case Y = A is N × N on the left–hand and K × K on the right–hand side.Second, as the correlation matrices C ep are dimensionless, the amplitudes r in thedistributions (18) and (19) must be dimensionless as well. If necessary, one can quicklygo over to the distributions involving the covariance matrix Σ ep = σ ep C ep σ ep with thevariances σ ep as in Eq. (7), we have p Y ( r (cid:48) | Σ ep )d[ r (cid:48) ] = p Y ( r | C ep )d[ r ] (35)for Y = G, A with the new amplitudes r (cid:48) = σ ep r carrying a physical dimension.Similarly, we have w Y ( X (cid:48) | Σ , Ξ)d[ X (cid:48) ] = w Y ( X | C, D )d[ X ] (36)with the rescaled, dimension carrying random matrices X (cid:48) = σXξ .
3. Methods
Although the calculations to be carried out in Sec. 4 are based on a straightforwardstrategy, the details become rather involved and complex. To still render the derivations xact Multivariate Amplitude Distributions
Facilitating the computations to follow, we notice that the algebraic distribution (19) isthe integral transform p A ( r | C ep ) = ∞ (cid:90) χ l − K/ ( z ) p G (cid:16) r (cid:12)(cid:12)(cid:12) mz C ep (cid:17) d z (37)of the multivariate Gaussian (18), involving the χ distribution χ q ( z ) = 12 q/ Γ( q/ z q/ − exp( − z/ z ) (38)of q degrees of freedom. Formula (37) has a generalization to the space of N × N real–symmetric, positive definite matrices Z , indicated by the notation Z >
0. We have theintegral transform w A ( (cid:101) X † (cid:101) X ) = (cid:114) M KN N (cid:89) n =1 Γ( L − ( n − / L − ( K + n − / (cid:90) i NL ( Z ) w G (cid:18) M Z (cid:101) X † (cid:101) X (cid:19) d[ Z ] (39)with (cid:101) X = C − / XD − / . The volume element d[ Z ] is the product of the differentials ofall independent variables. We introduce the Ingham–Siegel distribution i Nq ( Z ) = 1 π N ( N − / (cid:81) Nn =1 Γ( q − ( n − /
2) det q − ( N +1) / Z exp( − tr Z )Θ( Z ) (40)which is, apart from a scaling factor in the exponent, a matrix generalization of the χ distribution (38). The matrix Heaviside distribution Θ is one, whenever all eigenvalues ofthe real–symmetric matrix in its argument are positive and zero otherwise. Formula (39)is an application of the Ingham–Siegel integral [65, 66] (cid:90) Z> exp( − tr ZR )det q − ( N +1) / Z d[ Z ] = π N ( N − / (cid:81) Nn =1 Γ( q − ( n − / q R , (41)where the matrix R is real–symmetric as well. Convergence is guaranteed if q ≥ ( N + 1) /
2. Using Eq. (41) and the normalization of the distribution (24), one easilyderives the normalization constant α KNLM in the distribution (25). Owing to Eq. (34), xact Multivariate Amplitude Distributions w A ( (cid:101) X (cid:101) X † ) = (cid:114) M KN K (cid:89) k =1 Γ( L − ( k − / L − ( N + k − / (cid:90) i KL ( Z ) w G (cid:18) M Z (cid:101) X (cid:101) X † (cid:19) d[ Z ] , (42)where Z is now a positive K × K integration matrix.Integrals over the multimultivariate Gaussian (24) and related expressions areconveniently done as integrals over multivariate Gaussian–like functions by observingthat the trace in the exponent may, by virtue oftr F X † GX = x † ( F † ⊗ G ) x , (43)be written as quadratic form of the KN component vector x , constructed from thecolumns X ( n ) , n = 1 , . . . , N , of the K × N matrix X according to x = X (1)... X ( N ) . (44)Formula (43) holds for real matrices F and G of dimensions N × N and K × K ,respectively.To carry out the ensemble averages, we always work with the Fourier transform, i.e. with the charateristic function of the amplitude distribution, depending on the K component vector ω , (cid:104) ϕ (cid:105) ( ω | C, D ) = (cid:90) (cid:104) p (cid:105) ( r | C, D ) exp(i ω · r )d[ r ]= (cid:90) ϕ (cid:18) ω (cid:12)(cid:12)(cid:12) N XX † (cid:19) w ( X | C, D )d[ X ] , (45)where we used Eq. (16). Here, ϕ is the characteristic function of p . The correspondingFourier backtransform is (cid:104) p (cid:105) ( r | C, D ) = 1(2 π ) K (cid:90) (cid:104) ϕ (cid:105) ( ω | C, D ) exp( − i ω · r )d[ ω ] . (46)In view of relation (37), it suffices to employ the Gaussian amplitude distributions, p G (cid:18) r (cid:12)(cid:12)(cid:12) N XX † (cid:19) = 1(2 π ) K (cid:90) ϕ G (cid:18) ω (cid:12)(cid:12)(cid:12) N XX † (cid:19) exp( − i ω · r )d[ ω ] ϕ G (cid:18) ω (cid:12)(cid:12)(cid:12) N XX † (cid:19) = exp (cid:18) − N ω † XX † ω (cid:19) . (47)While the inverse of XX † appears in the exponent of the distribution, it is XX † itselfin the charateristic function, implying that the ensemble average in Eq. (45) can, byusing all above results, be reduced to Gaussian integrals. xact Multivariate Amplitude Distributions All results on special functions to be summarized here as well as notations andconventions stem from Ref. [67]. Due to the invariances of the random matrix ensembles,we often have to calculate Fourier back transforms (46) of characteristic functions whichdepend, due to the structure of our model, on ω only via a quadratic form involvingthe correlation matrix C . To indicate that, we write occasionally (cid:104) ϕ (cid:105) ( ω † Cω ) instead of (cid:104) ϕ (cid:105) ( ω | C, D ). Changing variables according to ω → C / ω , we have1(2 π ) K (cid:90) (cid:104) ϕ (cid:105) ( ω † Cω ) exp( − i ω · r )d[ ω ] = 1(2 π ) K √ det C (cid:90) (cid:104) ϕ (cid:105) ( ω ) exp( − i ω · C − / r )d[ ω ] , (48)introduce hyperspherical coordinates with radius ρ = | ω | and recognize the angularintegrals as spherical Bessel function of zeroth order in K dimensions. More explicitly,we choose ϑ as the angle between ω and C − / r , find a factor of sin K − ϑ in the Jacobianand use the integral π (cid:90) exp( − i z cos ϑ ) sin µ ϑ d ϑ = √ π (cid:18) z (cid:19) µ/ Γ (cid:18) µ + 12 (cid:19) J µ/ ( z ) , (49)where J µ/ ( z ) is the Bessel function of the first kind and of order µ/
2. The remainingangular integrals give a group volume. Collecting everything, we arrive at1(2 π ) K (cid:90) (cid:104) ϕ (cid:105) ( ω † Cω ) exp( − i ω · r )d[ ω ]= 1 √ det 2 πC √ r † C − r ( K − / ∞ (cid:90) (cid:104) ϕ (cid:105) ( ρ ) J ( K − / ( ρ √ r † C − r ) ρ K/ d ρ , (50)which reduces the K dimensional Fourier backtransform to a certain one–dimensionalHankel transform. Furthermore, we need the confluent hypergeometric or Kummerfunction F ( µ, ν, − z ) = 2Γ( ν )Γ( µ ) √ z ν − ∞ (cid:90) exp( − t ) t µ − ν J ν − (2 √ zt )d t (51)as well as the Tricomi function U ( µ, ν, z ) = 1Γ( µ ) ∞ (cid:90) exp( − zt ) t µ − (1 + t ) ν − µ − d t . (52)We will also employ the integral ∞ (cid:90) J ν ( bt ) t ν +1 ( a + t ) µ +1 d t = a ν − µ b µ µ Γ( µ + 1) K ν − µ ( ab ) (53) xact Multivariate Amplitude Distributions a and b which relates the ordinary and the modified Bessel function ofthe second kind or Macdonald function K µ ( w ) = w µ ∞ (cid:90) exp (cid:18) − t − w t (cid:19) t − µ − d t (54)of order µ . The latter appears in the integral U ( µ, ν, z ) = 2 z ( ν − / Γ( µ )Γ( µ − ν + 1) ∞ (cid:90) exp( − t ) t µ − ( ν +1) / K ν − (2 √ zt )d t (55)which connects Tricomi and Macdonald function. The Macdonald function may also beused to simplify the following double integral over two χ distributions and an arbitraryfunction h ( zz (cid:48) ), depending on the arguments of the former, ∞ (cid:90) d zχ κ ( z ) ∞ (cid:90) d z (cid:48) χ λ ( z (cid:48) ) h ( zz (cid:48) ) = ∞ (cid:90) d uh ( u ) ∞ (cid:90) d zχ κ ( z ) ∞ (cid:90) d z (cid:48) χ λ ( z (cid:48) ) δ ( u − zz (cid:48) )= ∞ (cid:90) d uh ( u ) ∞ (cid:90) d zz χ κ ( z ) χ λ (cid:16) uz (cid:17) = 2 − κ − λ Γ( κ )Γ( λ ) ∞ (cid:90) h ( u ) K λ − κ ( √ u ) u ( λ + κ ) / − d u (56)owing to the integral representation (54). Finally, we need the Gaussian hypergeometricfunction F ( λ, µ, ν, z ) = Γ( ν )Γ( µ )Γ( ν − µ ) (cid:90) t µ − (1 − t ) ν − µ − (1 − tz ) − λ d t (57)which is related to the ordinary and modified Bessel functions by the integrals ∞ (cid:90) t − λ K µ ( at ) J ν ( bt )d t = b ν Γ(( ν − λ + µ + 1) / ν − λ − µ + 1) / λ +1 a ν − λ +1 Γ( ν + 1) F (cid:18) ν − λ + µ + 12 , ν − λ − µ + 12 , ν + 1 , − b a (cid:19) (58)and ∞ (cid:90) t − λ K µ ( at ) K ν ( bt )d t = b ν Γ(( ν − λ + µ + 1) / ν − λ − µ + 1) / λ +2 a ν − λ +1 Γ((1 − λ + µ − ν ) / − λ − µ − ν ) / − λ ) F (cid:18) ν − λ + µ + 12 , ν − λ − µ + 12 , − λ, − b a (cid:19) (59)with constants a and b . All formulas given here are valid within certain parameterranges which can be found in Ref. [67]. xact Multivariate Amplitude Distributions
4. Derivations and Results for the Multivariate Ensemble AveragedAmplitude Distributions
After some preparatory remarks in Sec. 4.1, we study the cases Gaussian–Gaussian,Gaussian–Algebraic, Algebraic–Gaussian and Algebraic–Algebraic in Secs. 4.2 to 4.5,respectively.
The amplitude distributions within the epochs depend on the amplitudes r and acorrelation matrix C ep considered to be fixed within the epoch. Its fluctuations aremodeled by the distributions of the random data matrix X . We calculate the ensembleaverages (cid:104) p (cid:105) Y Y (cid:48) ( r | C, D ) = (cid:90) p Y (cid:18) r (cid:12)(cid:12)(cid:12) N XX † (cid:19) w Y (cid:48) ( X | C, D )d[ X ] (60)as amplitude distribution for the long time interval for all combinations Y, Y (cid:48) = G, A .We also compute the estimators for the sample correlations or covariances,respectively. They are given by (cid:104) rr † (cid:105) Y Y (cid:48) = (cid:90) rr † (cid:104) p (cid:105) Y Y (cid:48) ( r | C, D )d[ r ] . (61)Inserting formula (60) and interchanging the order of integration, we find with Eqs. (21)and (28) (cid:104) rr † (cid:105) Y Y (cid:48) = (cid:90) d[ X ] w Y (cid:48) ( X | C, D ) (cid:90) d[ r ] rr † p Y (cid:18) r (cid:12)(cid:12)(cid:12) N XX † (cid:19) = (cid:90) d[ X ] w Y (cid:48) ( X | C, D ) β Y N XX † = β Y B Y (cid:48) C (62)with β Y and B Y (cid:48) given in Eqs. (22) and (29). As in Sec. 2.4, if instead of the correlationmatrices C and D the covariance matrices Σ and Ξ are used, C in relation (62) has tobe replaced with Σ and its right hand side must be multplied with tr Ξ /N . Generalizing the results of Ref. [26], we include non–Markovian effects bei considering anon–trivial correlation matrix D (cid:54) = N . We find for the ensemble averaged characteristic xact Multivariate Amplitude Distributions (cid:104) ϕ (cid:105) GG ( ω | C, D ) = (cid:90) ϕ G (cid:18) ω (cid:12)(cid:12)(cid:12) N XX † (cid:19) w G ( X | C, D )d[ X ]= 1 √ det 2 πD ⊗ C (cid:90) exp (cid:18) − x † (cid:18) N ⊗ ωω † N + D − ⊗ C − (cid:19) x (cid:19) d[ X ]= 1 (cid:112) det( N ⊗ K + D ⊗ Cωω † /N )= 1 (cid:112) det( N + Dω † Cω/N ) , (63)where we drastically reduce the dimension of the determinant in the last step with thehelp of Sylvester’s theorem. To calculate the ensemble averaged distribution (cid:104) p (cid:105) GG ( r | C, D ) = 1(2 π ) K (cid:90) exp( − i ω · r ) (cid:112) det( N + Dω † Cω/N ) d[ ω ] , (64) i.e. , the Fourier backtransform (46), we may use formula (50) and arrive at (cid:104) p (cid:105) GG ( r | C, D ) = 1 √ det 2 πC √ r † C − r ( K − / ∞ (cid:90) J ( K − / ( ρ √ r † C − r ) (cid:112) det( N + Dρ /N ) ρ K/ d ρ , (65)which depends, due to the invariances of the ensemble, on the amplitudes only viathe quadratic form r † C − r . Furthermore, the matrix D enters formula (65) only byits eigenvalues Υ = diag (Υ(1) , . . . , Υ( N )) such that det( N + Dρ /N ) = (cid:81) Nn =1 (1 +Υ( n ) ρ /N ). In the non–Markovian case D = N , the integral over ρ in Eq. (65) can bedone with the help of Eq. (53) (cid:104) p (cid:105) GG ( r | C, N ) = 12 N/ − Γ( N/ (cid:112) det 2 πC/N K ( K − N ) / ( √ N r † C − r ) √ N r † C − r ( K − N ) / , (66)confirming the result of Ref. [26].For data analysis, the asymptotic behavior for large arguments √ N r † C − r isimportant. The Macdonald functions has a well–known exponential decay, whichdetermines the determines the asymptotics of (cid:104) p (cid:105) GG ( r | C, N ). There is no reason tobelieve that this exponential behavior turns to an algebraic one for D (cid:54) = N , we presentsome arguments in Appendix A. To calculate the ensemble average of the characteristic function, we employ the integralrepresentation (39) of the algebraic distribution, interchange the order of the matrix xact Multivariate Amplitude Distributions (cid:104) ϕ (cid:105) GA ( ω | C, D ) = (cid:90) ϕ G (cid:18) ω (cid:12)(cid:12)(cid:12) N XX † (cid:19) w A ( X | C, D )d[ X ]= (cid:114) M KN N (cid:89) n =1 Γ( L − ( n − / L − ( K + n − / (cid:90) d[ Z ] i NL ( Z ) (cid:90) ϕ G (cid:18) ω (cid:12)(cid:12)(cid:12) N XX † (cid:19) w G (cid:18) M Z (cid:101) X † (cid:101) X (cid:19) d[ X ] . (67)Using Eq. (43), we write the X dependent terms in the exponent as − N ω † XX † ω − M tr Z (cid:101) X † (cid:101) X = − N tr X † ωω † X − M tr D − / ZD − / X † CX = − x † (cid:18) N ⊗ ωω † N + 2 M D − / ZD − / ⊗ C − (cid:19) x (68)and do the Gaussian integral over X which yields √ π KN (cid:112) det ( N ⊗ ωω † /N + 2 D − / ZD − / /M ⊗ C − )= √ π KN det N/ C det ( K − / (2 D − / ZD − / /M ) (cid:112) det( ω † Cω N /N + 2 D − / ZD − / /M )= (cid:114) M N ( K − √ π KN √ det D ⊗ C det ( K − / Z (cid:112) det( ω † CωD/N + 2
Z/M )= √ πM N ( K − √ det D ⊗ C det ( K − / Z (cid:90) exp (cid:18) − ζ † (cid:18) ω † Cω N D + 1
M Z (cid:19) ζ (cid:19) d[ ζ ] . (69)Again, it was possible to considerably reduce the dimension of the determinant, theresulting determinant is written as a Gaussian integral over an N component vector ζ .Inserting this into Eq. (67), we have (cid:104) ϕ (cid:105) GA ( ω | C, D ) = 1 √ πM N N (cid:89) n =1 Γ( L − ( n − / L − ( K + n − / (cid:90) d[ Z ] i NL ( Z )det ( K − / Z (cid:90) exp (cid:18) − ζ † (cid:18) ω † Cω N D + 1
M Z (cid:19) ζ (cid:19) d[ ζ ] . (70)Observing ζ † Zζ = tr ζζ † Z , the matrix integral is conveniently recognized as an Ingham–Siegel integral (41), (cid:90) i NL ( Z )det ( K − / Z exp (cid:18) − ζ † M Zζ (cid:19) d[ Z ]= N (cid:89) n =1 Γ( L − ( K − / − ( n − / L − ( n − /
2) 1det L − ( K − / ( N + ζζ † /M ) . (71) xact Multivariate Amplitude Distributions L − ( K − / ( N + ζζ † /M ) = 1(1 + ζ † ζ/M ) L − ( K − / = ∞ (cid:90) χ L − ( K − / ( z ) exp (cid:18) − ζ † ζ M z (cid:19) d z , (72)where we used Eqs. (37) and (38). Many of the remaining Γ functions cancel. Collectingeverything we observe that the ζ integral can be done and (cid:104) ϕ (cid:105) GA ( ω | C, D ) = (cid:114) M N Γ( L − ( K − / L − ( K + N − / ∞ (cid:90) χ L − ( K − / ( z )d z (cid:112) det( z N /M + ω † CωD/N ) (73)is the final expression for the ensemble averaged characteristic function. By virtue ofEq. (50), we find (cid:104) p (cid:105) GA ( r | C, D ) = (cid:114) M N Γ( L − ( K − / L − ( K + N − /
2) 1 √ det 2 πC √ r † C − r ( K − / ∞ (cid:90) ∞ (cid:90) J ( K − / ( ρ √ r † C − r ) χ L − ( K − / ( z ) (cid:112) det( z N /M + ρ D/N ) d zρ K/ d ρ , (74)which expresses the ensemble averaged amplitude distribution as twofold integral. Aspointed out in Sec. 2.4, the algebraic distribution w A ( X | C, D ) converges to the Gaussian w G ( X | C, D ) in the limit
L, M → ∞ under the condition (27). Consequently, this limitof Eq. (74) also reproduces the result (65). This is easily shown by rescaling z → zM and then carrying out a saddlepoint approximation for the z integral, the saddlepointlies at z = 1.One of the integrals in Eq. (74) can be carried out by making the change of variables z → ρ z which removes the ρ dependence from the determinant. The ρ integral isthen of the form (51), yielding the Kummer function. Changing variables according to z = M/N u we eventually end up with (cid:104) p (cid:105) GA ( r | C, D ) = Γ( L − ( N − / L − ( K + N − / K/
2) 1 (cid:112) det 2 πCM/N ∞ (cid:90) F (cid:18) L − N − , K , − uN M r † C − r (cid:19) u ( K − / d u (cid:112) det ( N + uD ) . (75)In the Markovian case D = N , the determinant becomes a power and formula (52)allows us to do the integral, (cid:104) p (cid:105) GA ( r | C, N ) = Γ( L − ( N − / L − ( K − / L − ( K + N − / N/
2) 1 (cid:112) det 2 πCM/NU (cid:18) L − N − , K − N , N M r † C − r (cid:19) , (76) xact Multivariate Amplitude Distributions (cid:104) p (cid:105) GA ( r | C, D ) is algebraic, not exponentialas in the case of a Gaussian distribution of the random data matrices. More precisely,we have (cid:104) p (cid:105) GA ( r | C, D ) ∼ r † C − r ) L − ( N − / (77)for large values of √ r † C − r . Of course, this behavior comes from the algebraicdistribution of the random data matrices, the derivation is given in Appendix A. Major steps in the calculation can be traced back to the Gaussian–Gaussian case. UsingEq. (37) in Eq. (60) and interchanging integrals, we find (cid:104) p (cid:105) AG ( r | C, D ) = ∞ (cid:90) d zχ l − K/ ( z ) (cid:90) p G (cid:18) r (cid:12)(cid:12)(cid:12) mz N XX † (cid:19) w G ( X | C, D )d[ X ]= ∞ (cid:90) d zχ l − K/ ( z ) (cid:114) zm K (cid:90) p G (cid:18)(cid:114) zm r (cid:12)(cid:12)(cid:12) N XX † (cid:19) w G ( X | C, D )d[ X ]= ∞ (cid:90) χ l − K/ ( z ) (cid:114) zm K (cid:104) p (cid:105) GG (cid:18)(cid:114) zm r (cid:12)(cid:12)(cid:12) C, D (cid:19) d z (78)for the ensemble averaged amplitude distribution and (cid:104) ϕ (cid:105) AG ( ω | C, D ) = ∞ (cid:90) χ l − K/ ( z ) (cid:104) ϕ (cid:105) GG (cid:18)(cid:114) mz ω (cid:12)(cid:12)(cid:12) C, D (cid:19) d z (79)for the characteristic function, which both are one–dimensional χ transforms of thecorresponding function in the Gaussian–Gaussian case. Inserting formula (63) we arriveat (cid:104) ϕ (cid:105) AG ( ω | C, D ) = (cid:114) m N Γ( l − ( K − N ) / l − K/ ∞ (cid:90) χ l − ( K − N ) / ( z )d z (cid:112) det( z N /m + ω † CωD/N ) (80)for the characteristic function, which coincides with the result (73) in the Gaussian–Algebraic case if we replace L , M and N with proper combinations of l , m and N ,which are, however, not easy to guess as the different dependencies on these parametershave different origin. With formula (50), we find (cid:104) p (cid:105) AG ( r | C, D ) = (cid:114) m N Γ( l − ( K − N ) / l − K/
2) 1 √ det 2 πC √ r † C − r ( K − / ∞ (cid:90) ∞ (cid:90) J ( K − / ( ρ √ r † C − r ) χ l − ( K − N ) / ( z ) (cid:112) det( z N /m + ρ D/N ) d zρ K/ d ρ (81) xact Multivariate Amplitude Distributions (cid:104) p (cid:105) AG ( r | C, D ) = Γ( l )Γ( l − K/ K/
2) 1 (cid:112) det 2 πCm/N ∞ (cid:90) F (cid:18) l, K , − uN m r † C − r (cid:19) u ( K − / d u (cid:112) det ( N + uD ) (82)for arbitrary D as well as (cid:104) p (cid:105) AG ( r | C, N ) = Γ( l )Γ( l − ( K − N ) / l − K/ N/
2) 1 (cid:112) det 2 πCm/NU (cid:18) l, K − N , N m r † C − r (cid:19) (83)in the Markovian case D = N . The similarity of the distributions to those given inSec. 4.3 does not come as a surprise, but it is a purely mathematical one. From theviewpoint of physics and data analysis, the cases Gaussian–Algebraic and Algebraic–Gaussian are very different and only identical for the irrelevant parameter value N = 1,corresponding to model time series of length one.For arbitrary D , the asymptotic behavior of (cid:104) p (cid:105) AG ( r | C, D ) can be infered fromSec. 4.3 and Appendix A. We find (cid:104) p (cid:105) AG ( r | C, D ) ∼ r † C − r ) l (84)for large values of √ r † C − r , which results from the asymptotic relation (77) by replacing L − ( N − / l . This is plausibel, because the algebraic behavior stems onlyfrom the amplitude distributions within the epochs which formally coincides with thedistribution of the random data matrices for l = L and N = 1. As in the Algebraic–Gaussian case we use a shortcut, here by observing that the desiredfunctions can be written as integrals over the corresponding ones in the Gaussian–Algebraic case. The steps in Eq. (78) carry over to an algebraic distribution of therandom data matrix, yielding for the ensemble averaged amplitude distribution theintegral (cid:104) p (cid:105) AA ( r | C, D ) = ∞ (cid:90) χ l − K/ ( z ) (cid:114) zm K (cid:104) p (cid:105) GA (cid:18)(cid:114) zm r (cid:12)(cid:12)(cid:12) C, D (cid:19) d z (85) xact Multivariate Amplitude Distributions (cid:104) p (cid:105) GA calculated in Sec. 4.3. Similarly, we have (cid:104) ϕ (cid:105) AA ( ω | C, D ) = ∞ (cid:90) χ l − K/ ( z ) (cid:104) ϕ (cid:105) GA (cid:18)(cid:114) mz ω (cid:12)(cid:12)(cid:12) C, D (cid:19) d z (86)for the characteristic function. Pluging in the result (73) and rearranging terms, we findfor the latter (cid:104) ϕ (cid:105) AA ( ω | C, D ) = 2 N √ M m N Γ( l − ( K − N ) / L − ( K − / l − K/ L − ( K + N − / ∞ (cid:90) ∞ (cid:90) χ L − ( K − / ( z ) χ l − ( K − N ) / ( z (cid:48) )d z d z (cid:48) (cid:112) det( zz (cid:48) N / ( M m ) + ω † CωD/N ) , (87)which may, by virtue of Eq. (56), be cast into the form of a certain one–dimensionalBessel transform, (cid:104) ϕ (cid:105) AA ( ω | C, D ) = 2 K +( N +1) / − L − l √ M m N Γ( l − K/ L − ( K + N − / ∞ (cid:90) K l − L +( N − / ( √ u ) √ u L + l − K +( N +1) / − (cid:112) det( u N / ( M m ) + ω † CωD/N ) d u . (88)Applying formula (50) yields after some algebra the ensemble averaged amplitudedistribution (cid:104) p (cid:105) AA ( r | C, D ) = 2 K +( N +1) / − L − l √ det 2 πC √ r † C − r ( K − / √ M m N Γ( l − K/ L − ( K + N − / ∞ (cid:90) ∞ (cid:90) J ( K − / ( ρ √ r † C − r ) K l − L +( N − / ( √ u ) √ u L + l − K +( N +1) / − ρ K/ (cid:112) det( u N / ( M m ) + ρ D/N ) d u d ρ . (89)Analogously to Sec. 4.3, the ρ integral can be performed by making the change ofvariables u → uρ , which removes ρ from the determinant. The ρ integral is then of theform (58) and gives the Gaussian hypergeometric function. With the change of variables u = M m/N v , we arrive at (cid:104) p (cid:105) AA ( r | C, D ) = Γ( l )Γ( L − ( N − / (cid:112) det πCM m/N Γ( l − K/ L − ( K + N − / K/ ∞ (cid:90) F (cid:18) l, L − N − , K , − N r † C − rM m v (cid:19) v ( K − / d v (cid:112) det( N + vD ) . (90)Even in the Algebraic–Algebraic case, a reduction to a one–dimensional integral can becarried out for arbitrary D . In the Markovian case D = N , we start from Eq. (89), use xact Multivariate Amplitude Distributions ρ integral which allows us to apply formula (59) to the remaining u integral, (cid:104) p (cid:105) AA ( r | C, N ) = Γ( l )Γ( l − ( K − N ) / (cid:112) det πCM m/N Γ( l − K/ L − ( N − / L − ( K − / L − ( K + N − / L + l − ( K − / N/ F (cid:18) l, L − N − , L + l − K − , − N r † C − rM m (cid:19) , (91)which is a complicated, but still a closed form expression as in the previous cases.As both, the amplitude distributions within the epochs with parameter l as well asthe distribution of the random data matrices with parameters L and N , are algebraic,the asymptotic behavior may be governed by either of them, (cid:104) p (cid:105) AA ( r | C, D ) ∼ r † C − r ) l , if l < L − N − ln r † C − r ( r † C − r ) l , if l = L − N − r † C − r ) L − ( N − / , if l > L − N − . (92)The derivation is given in Appendix A.
5. A Technique and Formulae for the Analysis of Data
We begin with general considerations in Sec. 5.1, before we collect the results for thefour cases in Secs. 5.2 to 5.5. In Sec. 5.6, some figures illustrate our results.
To compare our K variate distributions with data, we further extend the methodintroduced in Ref. [26] which we developed in the spirit of aggregation. The crucialidea is to construct K univariate distributions out of the K variate one which arethen overlaid, i.e. , lumped together, or, if meaningful, analyzed individually. We takeadvantage of the fact that all K variate distributions depend on the amplitudes r viathe Mahalanobis distance r † C − r only. As anticipated in Sec. 3 and explicitly shown inSec. 4, the Fourier backtransform (46) which yields the ensemble averaged amplitudedistribution (cid:104) p (cid:105) Y Y (cid:48) ( r | C, D ) has always the form (48). To decouple the amplitudes, werotate the vektor r into the eigenbasis of the correlation matrix C . More precisely, weuse the diagonalization C = U Λ U † such that C − / = U Λ − / U † , (93)where U is an orthogonal K × K matrix and Λ is the diagonal matrix of the eigenvaluesΛ k . As they are positive definite, the square roots Λ / k are real, we choose them positive. xact Multivariate Amplitude Distributions r = U † r (94)as new arguments of the ensemble averaged amplitude distribution. The correspondingJacobi determinant is one and the functional form of the distribution is thus not altered.In the exponential function on the right hand side of Eq. (48), we have ω · C − / r = ω · U Λ − / ¯ r = ( U † ω ) · Λ − / ¯ r . (95)As the characteristic function (cid:104) ϕ (cid:105) Y Y (cid:48) depends on the vector ω via its length only, thechange ω → U † ω of integration variables fully removes U from the integrand. We find (cid:104) p (cid:105) Y Y (cid:48) (¯ r | C, D ) = 1(2 π ) K √ det C (cid:90) (cid:104) ϕ (cid:105) Y Y (cid:48) ( ω ) exp( − i ω · Λ − / ¯ r )d[ ω ] . (96)We are now ready to define K univariate distributions by integrating out all rotatedamplitudes but the k –th one, (cid:104) p (cid:105) ( k ) Y Y (cid:48) (¯ r k | C, D ) = (cid:90) (cid:104) p (cid:105) Y Y (cid:48) (¯ r | C, D )d[¯ r ] (cid:54) = k k = 1 , . . . , K . (97)Inserting Eq. (96), we can do all these K − (cid:81) l (cid:54) = k π Λ / l δ ( ω l ).Hence, all ω integrals except the k –th one can be carried out and we arrive at a one–dimensional Fourier transform, (cid:104) p (cid:105) ( k ) Y Y (cid:48) (¯ r k | Λ k , D ) = 1 π √ Λ k ∞ (cid:90) (cid:104) ϕ (cid:105) Y Y (cid:48) ( ω k ) cos ω k ¯ r k √ Λ k d ω k , (98)which reduces to a Fourier cosine transform because (cid:104) ϕ (cid:105) Y Y (cid:48) is an even function. Wethus obtain K univariate distributions of the same form, but scaled with √ Λ k . It ishelpful to rewrite Eq. (98) by usingcos z = (cid:114) πz J − / ( z ) (99)which yields (cid:104) p (cid:105) ( k ) Y Y (cid:48) (¯ r k | Λ k , D ) = 1 √ π Λ k (cid:18) ¯ r k Λ k (cid:19) / ∞ (cid:90) (cid:104) ϕ (cid:105) Y Y (cid:48) ( ρ ) J − / (cid:115) ¯ r k Λ k ρ √ ρ d ρ . (100)This formula greatly simplifies the ensuing derivations that involve algebraicdistributions, as we can proceed exactly as in Sec. 4, we only have to set K = 1 andreplace r † C − r with ¯ r k / Λ k in the term J ( K − / ( ρ √ r † C − r ) ρ K/ / √ r † C − r ( K − / relatedto the spherical Bessel function. We notice that (cid:104) ϕ (cid:105) Y Y (cid:48) ( ρ ) is unchanged and, in general,depends on K directly as well as implicitly. Furthermore, the arguments leading to theasymptotic behavior carry over to the present discussion. xact Multivariate Amplitude Distributions (cid:104) p (cid:105) ( k ) Y Y (cid:48) (¯ r k | Λ k , D ) follow directly fromthe result (62). Due to the symmetry of the distributions, they coincide with the secondmoments (cid:104) ¯ r k (cid:105) ( k ) Y Y (cid:48) . Consider the left hand side of Eq. (62) in the rotated coordinates (94), (cid:104) rr † (cid:105) Y Y (cid:48) = U (cid:104) ¯ r ¯ r † (cid:105) Y Y (cid:48) U † (101)and thus (cid:104) ¯ r ¯ r † (cid:105) Y Y (cid:48) = U † (cid:104) rr † (cid:105) Y Y (cid:48) U = U † β Y B Y (cid:48) CU = β Y B Y (cid:48) Λ . (102)Naturally, the correlation matrix in the eigenbasis of C is diagonal and proportional toΛ. This implies (cid:104) ¯ r k (cid:105) ( k ) Y Y (cid:48) = β Y B Y (cid:48) Λ k (103)with β Y and B Y (cid:48) given in Eqs. (22) and (29).When analyizing data, K is given, we obtain the matrices C and D or Σ and Ξ byusing the the originally measured amplitudes for sampling over the long time interval T tot . The data are normalized in different ways for the computation of the time andthe position series, see Sec. 2.1, implying that C and D have different eigenvalues.The result (103) is important as (cid:104) ¯ r k (cid:105) Y Y (cid:48) serves as sample variance for the univariatedistributions (cid:104) p (cid:105) ( k ) Y Y (cid:48) (¯ r k | Λ k , D ) in the eigenbasis of the correlation matrix, i.e. for therotated amplitudes. This variance allows one to determine an unknown parameter in theunderlying algebraic distributions, if the fixings (23) and (31) have not been made. In theGaussian–Algebraic and Algebraic–Gaussian cases, the result (103) provides relationsbetween the unknown parameters M , L , N and m , l , respectively. In the Algebraic–Algebraic case, there is only one relation for all of these parameters, hence, only one,not both, can be inferred from the sample variance. The relation between l and m canbe obtained within the epochs by sampling variances or by fitting. In all cases, theparameter N is a fit parameter, measuring the strength of the fluctuations. It might behelpful to use higher moments to obtain further relations. However, experience tells, that N sensitively determines the shape and is best obtained by fitting the whole distribution (cid:104) p (cid:105) ( k ) Y Y (cid:48) (¯ r k | Λ k , D ) to the data. Once more, if instead of the correlation matrices C and D the covariance matrices Σ and Ξ are used, Λ k in Eq. (103) has to be interpreted asthe k –th eigenvalue of Σ and the right hand side of Eq. (103) must be multiplied withtr Ξ /N . In a general non–Markovian case with D (cid:54) = N , we obtain the one–dimensional integral (cid:104) p (cid:105) ( k ) GG (¯ r k | Λ k , D ) = 1 π √ Λ k ∞ (cid:90) N + Dρ /N ) cos ρ ¯ r k √ Λ k d ρ (104) xact Multivariate Amplitude Distributions D = N to (cid:104) p (cid:105) ( k ) GG (¯ r k | Λ k , N ) = 12 ( N − / Γ( N/ (cid:112) π Λ k /N (cid:115) N ¯ r k Λ k ( N − / K (1 − N ) / (cid:115) N ¯ r k Λ k (105)in agreement with Ref. [26]. This case is exceptional, because its Gaussian nature leadsto a very simple relation bewteen the results (105) and (66). The former follows fromthe latter by just putting K = 1. This is not so in the other three cases. Carrying out steps analogous to the ones in Sec. 4.3, we obtain the formulae for a generalMarkovian case with D (cid:54) = N , (cid:104) p (cid:105) ( k ) GA (¯ r k | Λ k , D ) = Γ( L − ( K + N ) / L − ( K + N − / π (cid:112) k M/N ∞ (cid:90) F (cid:18) L − K + N , , − uN M ¯ r k Λ k (cid:19) d u (cid:112) u det ( N + uD ) , (106)as well as for a Markovian situation with D = N , (cid:104) p (cid:105) ( k ) GA (¯ r k | Λ k , N ) = Γ( L − ( K + N ) / L − ( K − / L − ( K + N − / N/ (cid:112) π Λ k M/NU (cid:18) L − K + N , − N , N M ¯ r k Λ k (cid:19) . (107)We also find (cid:104) p (cid:105) ( k ) GA (¯ r k | Λ k , D ) ∼ (cid:18) ¯ r k Λ k (cid:19) − L +( K + N ) / − (108)as asymptotic result. This case is, as already mentioned in Sec. 4, to some extent similar to the previous oneand we arrive for arbitrary D at (cid:104) p (cid:105) ( k ) AG (¯ r k | Λ k , D ) = Γ( l − ( K − / l − K/ π (cid:112) k m/N ∞ (cid:90) F (cid:18) l − K − , , − uN m ¯ r k Λ k (cid:19) d u (cid:112) u det ( N + uD ) (109) xact Multivariate Amplitude Distributions D = N at (cid:104) p (cid:105) ( k ) AG (¯ r k | Λ k , N ) = Γ( l − ( K − / l − ( K − N ) / l − K/ N/ (cid:112) π Λ k m/NU (cid:18) l − K − , − N , N m ¯ r k Λ k (cid:19) . (110)Moreover, (cid:104) p (cid:105) ( k ) AG (¯ r k | Λ k , D ) ∼ (cid:18) ¯ r k Λ k (cid:19) − l +( K − / (111)is the asymptotic behavior. Finally, we provide the results for the Algebraic–Algebraic case. In a non–Markoviansituation with D (cid:54) = N , we find after steps analogous to the ones in Sec. 4.5 (cid:104) p (cid:105) ( k ) AA (¯ r k | Λ k , D ) = Γ( l − ( K − / L − ( K + N ) / π (cid:112) Λ k M m/N Γ( l − K/ L − ( K + N − / ∞ (cid:90) F (cid:18) l − K − , L − K + N , , − vNM m ¯ r k Λ k (cid:19) d v (cid:112) v det( N + vD ) , (112)while (cid:104) p (cid:105) ( k ) AA (¯ r k | Λ k , N ) = Γ( l − ( K − / l − ( K − N ) / (cid:112) π Λ k M m/N Γ( l − K/ L + l − ( K − L − ( K − / L − ( K + N ) / L − ( K + N − / N/ F (cid:18) l − K − , L − K + N , L + l − ( K − , − NM m ¯ r k Λ k (cid:19) , (113)results in the Markovian case D = N . We also obtain (cid:104) p (cid:105) ( k ) AA (¯ r k | Λ k , D ) ∼ (cid:18) ¯ r k Λ k (cid:19) − l − ( K − / , if l < L − N − (cid:18) ¯ r k Λ k (cid:19) − l − ( K − / ln ¯ r k Λ k , if l = L − N − (cid:18) ¯ r k Λ k (cid:19) − L +( K + N ) / − , if l > L − N − (114)for the asymptotic behavior. xact Multivariate Amplitude Distributions To render possible a comparison of the results involving algebraic distributions withthose in the Gaussian–Gaussian case, we choose values of L and l which ensure theexistence of the first matrix moment, see Sec. 2.4. We notice that the conditions onthe existence of the algebraic distributions, i.e. , of their normalizations, are slightlyweaker. We also use the fixings (23) and (31), implying that β A and B A become one andmaking sure that the ensemble averaged correlation matrix, corresponding to the samplecorrelation matrix, coincides in all cases with the input matrix C . Thus, according toEq. (103) the variances (cid:104) ¯ r k (cid:105) ( k ) Y Y (cid:48) are simply given by Λ k . The functional form of alldistributions (cid:104) p (cid:105) ( k ) Y Y (cid:48) (¯ r k | D ) then allows us to normalize the rotated amplitude ¯ r k by thestandard deviation, ˜ r = ¯ r k √ Λ k , (115)such that all K distributions in this variable coincide and the corresponding variances areall given by one. For the graphical representation, it is useful to look at the Markoviancase. We consider K = 100 positions. In Figs. 3 and 4, we use the shape parameters L = 55 and l = 55, as well as N = 5 which is a typical value from an empirical viewpoint. AAGA AGGG (cid:45) (cid:45) r (cid:142) pd f Figure 3.
Probability densities (cid:104) p (cid:105) ( k ) Y Y (cid:48) in the Markovian case versus the rotatedamplitudes ˜ r , normalized to unit standard deviation. The four cases Gaussian–Gaussian, Gaussian–Algebraic, Algebraic–Gaussian, Algebraic–Algebraic are labeledYY’ = GG, GA, AG, AA, respectively. Number of positions K = 100, shapeparameters L = 55 and l = 55, strength parameter for fluctuations of correlations N = 5. Linear scale, abscissa between -2 and +2. As seen in Fig. 3, the more algebraic, the stronger peaked is the distribution, and asFig. 4 shows, the heavier are the tails, corroborating the construction of our model. InFigs. 5 and 6, we look at the Gaussian–Algebraic case for K = 100 more closely by xact Multivariate Amplitude Distributions AA GA AGGG (cid:45) (cid:45) (cid:45) (cid:45) r (cid:142) pd f Figure 4.
As Fig. 3. Logarithmic scale, abscissa between -8 and +8. varying the dependence on L and N . As L is increased for fixed N = 5, the convergenceof the ensemble distribution for the random data matrix to the Gaussian competeswith the algebraic tail. Put differently, the exponential tail known from the Gaussian– GA L (cid:61)
55, 60, 65, 70, 75 (cid:45) (cid:45) (cid:45) (cid:45) r (cid:142) pd f Figure 5.
Probability densities (cid:104) p (cid:105) ( k ) GA in the Markovian case versus the rotatedamplitudes ˜ r , normalized to unit standard deviation, K = 100, N = 5, L varyingas indicated, the smaller L , the heavier the tail. Logarithmic scale. Gaussian case seems to win for larger L . But this is not so, it depends on the rangein ˜ r considered. On a sufficiently large ˜ r scale, the algebraic tail wins. Only in thehard limit L → ∞ for fixed ˜ r , there is an exponential tail everywhere. Furthermore,as seen in Fig. 6, the smaller the parameter N for fixed L = 65, the stronger are the xact Multivariate Amplitude Distributions GA N (cid:61)
3, 6, 9, 12 (cid:45) (cid:45) (cid:45) (cid:45) r (cid:142) pd f Figure 6.
Probability densities (cid:104) p (cid:105) ( k ) GA in the Markovian case versus the rotatedamplitudes ˜ r , normalized to unit standard deviation, K = 100, L = 65, N varyingas indicated, the smaller N , the heavier the tail. Logarithmic scale. non–stationary fluctuations of the random data matrices. Of course, this must have alarge impact on the distribution, but it is interesting to observe that the effect saturatesquickly, the curves for N = 9 and N = 12 are almost on top of each other.
6. Conclusions
We put forward an approach to model multivariate amplitude distributions in thepresence of non–stationarities that imply fluctuations of correlations. This is a frequentlyencountered situation in complex systems. Our model is based on a division of timescales, the idea behind and the necessary assumption is the existence of a shorter timescale on which the system is in good approximation stationary. Longer time scales arethen viewed as assembled from several or many intervals of the shorter time scale. Weargued that the correlation or covariance matrices measured in all these shorter timeintervals, referred to as epochs, form a truly existing ensemble that we modeled by anensemble of random matrices. Averaging the multivariate amplitude distributions ofthe epochs over this ensemble yields heavy tailed amplitude distributions for the largertime scales on which the non–stationarity is crucial. In short, the fluctations of thecorrelations lift the tails of the multivariate amplitude distributions.We evaluated our model explicitly for four cases, combining Gaussian andalgebraic multivariate amplitude distributions in the epochs with Gaussian and algebraicdistributions for the random correlation or covariance matrices. It is a welcome featurethat we arrive in all cases at closed form results if the system is Markovian, but even ifit is not, the resulting formulae are single integrals. A highly appreciated effect, as inevery ensemble approach, is the reduction in the number of relevant parameters. In all xact Multivariate Amplitude Distributions
Acknowledgement
We thank Mario Kieburg and Shanshan Wang for helpful discussions.
Appendix A. Asymptotics of the Ensemble Averaged AmplitudeDistributions
We always write w = √ r † C − r , we begin with the Gaussian–Algebraic case of Sec. 4.3.To study the integral in Eq. (75), to which we refer as Q GA , we write the determinantas Gaussian integral using an N component vector ζ and employ the power series of theKummer function Q GA = ∞ (cid:90) d uu K/ − (cid:112) det( N + uD ) F (cid:18) L − N − , K , − N w M u (cid:19) = 1 π N/ ∞ (cid:88) i =0 Γ( L − ( N − / i )Γ( K/ L − ( N − / K/ i ) i ! (cid:18) − N w M (cid:19) i (cid:90) d[ ζ ] exp( − ζ † ζ ) ∞ (cid:90) d uu K/ i − exp( − uζ † Dζ ) . (A.1) xact Multivariate Amplitude Distributions u integral can now be calculated, and we resum the resulting series Q GA = Γ( K/ π N/ (cid:90) d[ ζ ] exp( − ζ † ζ )( ζ † Dζ ) − K/ (cid:18) N w M ζ † Dζ (cid:19) − L + N − = Γ( K/ √ det πD (cid:18) N w M (cid:19) ( N − K ) / (cid:90) d[ η ] exp (cid:18) − N w M η † D − η (cid:19) ( η † η ) L − ( N + K − / (1 + η † η ) − L +( N − / , (A.2)where we changed variables according to ζ = (cid:112) N w / M D − / η in the second equation.Introducing hyperspherical coordinates η = ρe η with d[ η ] = ρ N − d ρ dΩ N where dΩ N isthe infinitesimal solid angle in N dimensions, we identify the ρ integral as the integralrepresentation (52) of the Tricomi function Q GA = Γ( K/ L − ( K − / √ det πD (cid:18) N w M (cid:19) ( N − K ) / (cid:90) U (cid:18) L − K − , N − K , N w M e † η D − e η (cid:19) dΩ N . (A.3)The Tricomi function has an asymptotic expansion in its last argument which yields asleading term the inverse of this last argument raised to the power of its first argument.Thus we find (cid:104) p (cid:105) GA ( r | C, D ) −→ Γ( L − ( N − / L − ( K − / L − ( K + N − / (cid:112) det 2 πCM/N √ det πD (cid:18) N w M (cid:19) − L + N − (cid:90) ( e † η D − e η ) − L +( K − / dΩ N . (A.4)Most conveniently, the angular average appears only as a factor independent of w and we arrive, as claimed, at the asymptotic behavior (76) of the ensemble averagedamplitude distribution for large w = √ r † C − r . Indirectly, this derivation shows thatthe ensemble avearged amplitude distribution (27) in the case of Gaussian distributedrandom data matrices cannot be algebraic. As discussed in Sec. 4.3, the limit L, M → ∞ of (cid:104) p (cid:105) GA ( r | C, D ) under the condition (65) yields (cid:104) p (cid:105) GG ( r | C, D ). This limit, however, isin conflict with the above used asymptotic expansion of the Tricomi function. If suchan expansion cannot be applied, the angular average will never separate off as a simpleproduct, seriously hampering a direct derivation of the exponential asymptotics in theGaussian–Gaussian case for D (cid:54) = N .The algebraic tail in the Algebraic–Gaussian case of Sec. 4.4 is derived accordingly.We notice that in this case, too, the tail in the Markovian situation D = N immediatelyfollows from the results (76) and (83) due to the above mentioned asymptotic expansionof the Tricomi function which features its last argument to the power of its firstargument. With this in mind, we now turn to the Algebraic–Algebraic case of Sec. 4.5. EFERENCES D = N as we may expect that it yields thesame asymptotics for the non–Markovian one, too. Thus we consider in the result (91)the hypergeometric function to which we refer as Q AA . In the regime w (cid:29) M m/N , wemay use the approximation Q AA (cid:39) F (cid:18) l, L − N − , L + l − K − , − N r † C − rM m (cid:19) . (A.5)An asymptotic expansion can then be infered from the transformation formula [67]sin π ( b − a ) π F ( a, b, c, z )Γ( c ) = ( − z ) − a Γ( b )Γ( c − a )Γ( a − b + 1) F (cid:18) a, a − c + 1 , a − b + 1 , z (cid:19) − ( − z ) − b Γ( a )Γ( c − b )Γ( b − a + 1) F (cid:18) b, b − c + 1 , b − a + 1 , z (cid:19) . (A.6)The first term in the defining Gauss series of the hypergeometric function is one,implying that the two hypergeometric functions on the right hand side become onefor large z . The asymptotics is hence determined by ( − z ) − a or ( − z ) − b , depending onwhether a < b or b > a . The sine function on the left hand side produces a sign whichensures the positivity of the distribution. To apply this particular line of reasoningto Q AA , we must assume that N is even, for odd N cancelations have to prevent thefunctions Γ( a − b + 1) and Γ( a − b + 1) from diverging. If b = a , we bring the sinefunction on the right hand side of Eq. (A.6) and use l’Hospital’s rule which yields F ( a, b, c, z ) ∼ ( − z ) − a Γ( a )Γ( c ) (cid:18) − z ) − ψ (1) − ψ ( a ) − ψ ( c − a )Γ( c − a ) z + . . . (cid:19) , (A.7)resulting in a logarithmic correction to the asymptotic behavior. Here ψ ( z ) = Γ (cid:48) ( z ) / Γ( z )is the digamma function. References [1] J. B. Gao. Recurrence Time Statistics for Chaotic Systems and Their Applications.
Physical Review Letters , 83:3178–3181, Oct 1999.[2] Rainer Hegger, Holger Kantz, Lorenzo Matassini, and Thomas Schreiber. Copingwith Nonstationarity by Overembedding.
Physical Review Letters , 84:4092–4095,May 2000.[3] Pedro Bernaola-Galv´an, Plamen Ch. Ivanov, Lu´ıs A. Nunes Amaral, and H. EugeneStanley. Scale Invariance in the Nonstationarity of Human Heart Rate.
PhysicalReview Letters , 87:168105, Oct 2001.[4] Christoph Rieke, Karsten Sternickel, Ralph G. Andrzejak, Christian E. Elger, PeterDavid, and Klaus Lehnertz. Measuring Nonstationarity by Analyzing the Loss ofRecurrence in Dynamical Systems.
Physical Review Letters , 88:244102, May 2002.
EFERENCES
Journal of Physics A , 37, 02 2004.[6] R K P Zia and B Schmittmann. A possible classification of nonequilibrium steadystates.
Journal of Physics A , 39(24):L407–L413, may 2006.[7] Jan Pieter Pijn, Jan Van Neerven, Andr´e Noest, and Fernando H Lopesda Silva. Chaos or noise in EEG signals; dependence on state and brain site.
Electroencephalography and Clinical Neurophysiology , 79(5):371–381, 1991.[8] Markus M¨uller, Gerold Baier, Andreas Galka, Ulrich Stephani, and HiltrudMuhle. Detection and characterization of changes of the correlation structure inmultivariate time series.
Physical Review E , 71:046116, Apr 2005.[9] R. H¨ohmann, U. Kuhl, H.-J. St¨ockmann, L. Kaplan, and E. J. Heller. Freak Wavesin the Linear Regime: A Microwave Study.
Physical Review Letters , 104:093901,Mar 2010.[10] Jakob J. Metzger, Ragnar Fleischmann, and Theo Geisel. Statistics of ExtremeWaves in Random Media.
Physical Review Letters , 112:203903, May 2014.[11] Henri Degueldre, Jakob Metzger, Theo Geisel, and Ragnar Fleischmann. RandomFocusing of Tsunami Waves.
Nature Physics , 12, 11 2015.[12] Geert Bekaert and Campbell R. Harvey. Time-Varying World Market Integration.
Journal of Finance , 50(2):403–444, 1995.[13] Francois Longin and Bruno Solnik. Is the correlation in international equity returnsconstant: 1960-1990?
Journal of International Money and Finance , 14(1):3–26,1995.[14] J.-P. Onnela, A. Chakraborti, K. Kaski, J. Kert´esz, and A. Kanto. Dynamicsof market correlations: Taxonomy and portfolio analysis.
Physical Review E ,68:056110, Nov 2003.[15] Yiting Zhang, Gladys Hui Ting Lee, Jian Cheng Wong, Jun Liang Kok,Manamohan Prusty, and Siew Ann Cheong. Will the US economy recover in 2010?A minimal spanning tree study.
Physica A , 390(11):2020–2050, 2011.[16] Dong-Ming Song, Michele Tumminello, Wei-Xing Zhou, and Rosario N. Mantegna.Evolution of worldwide stock markets, correlation structure, and correlation-basedgraphs.
Physical Review E , 84:026108, Aug 2011.[17] Leonidas Sandoval and Italo De Paula Franca. Correlation of financial markets intimes of crisis.
Physica A , 391(1):187–208, 2012.[18] Michael C. M¨unnix, Takashi Shimada, Rudi Sch¨afer, Francois Leyvraz, Thomas H.Seligman, Thomas Guhr, and H. Eugene Stanley. Identifying States of a FinancialMarket.
Scientific Reports , 2:644, Sep 2012.[19] F. Ghasemi, Muhammad Sahimi, Joachim Peinke, and M Tabar. Analysis ofNon-stationary Data for Heart-Rate Fluctuations in Terms of Drift and DiffusionCoefficients.
Journal of Biological Physics , 32:117–28, 11 2006.
EFERENCES
Physical Review E , 87:062139, Jun 2013.[21] Fatemeh Ghasemi, Muhammad Sahimi, J. Peinke, R. Friedrich, G. Reza Jafari,and M. Reza Rahimi Tabar. Markov analysis and Kramers-Moyal expansion ofnonstationary stochastic processes with application to the fluctuations in the oilprice.
Physical Review E , 75:060102, Jun 2007.[22] Rudi Sch¨afer and Thomas Guhr. Local normalization: Uncovering correlations innon-stationary financial time series.
Physica A , 389(18):3856–3865, 2010.[23] A. Bohr and B.R. Mottelson.
Nuclear Structure: Volume 2, Nuclear deformations .Nuclear Structure. Benjamin, 1969.[24] Vladimir Zelevinsky. Quantum Chaos and Complexity in Nuclei.
Annual Reviewof Nuclear and Particle Science , 46(1):237–279, 1996.[25] Thomas Guhr, Axel M¨uller-Groeling, and Hans A. Weidenm¨uller. Random MatrixTheories in Quantum Physics: Common Concepts.
Physics Reports , 299:189–425,1998.[26] Thilo A. Schmitt, Desislava Chetalova, Rudi Sch¨afer, and Thomas Guhr. Non-stationarity in financial time series: Generic features and tail behavior.
EurophysicsLetters , 103(5):58003, 2013.[27] Frederik Meudt, Martin Theissen, Rudi Sch¨afer, and Thomas Guhr. ConstructingAnalytically Tractable Ensembles of Non-Stationary Covariances with anApplication to Financial Data.
Journal of Statistical Mechanics , 2015(11):P11025,03 2015.[28] Thilo Schmitt, Desislava Chetalova, Rudi Sch¨afer, and Thomas Guhr. Credit Riskand the Instability of the Financial System: an Ensemble Approach.
EurophysicsLetters , 105, 09 2013.[29] Thilo Schmitt, Rudi Sch¨afer, and Thomas Guhr. Credit Risk: Taking FluctuatingAsset Correlations into Account.
Journal of Credit Risk , 11:73–94, 09 2015.[30] Rudi Sch¨afer, Sonja Barkhofen, Thomas Guhr, Hans-J¨urgen St¨ockmann, and UlrichKuhl. Compounding Approach for Univariate Time Series with NonstationaryVariances.
Physical Review E , 92:062901, Dec 2015.[31] Satya Dubey. Compound gamma, beta and F distributions.
Metrika: InternationalJournal for Theoretical and Applied Statistics , 16(1):27–31, 1970.[32] O. Barndorff-Nielsen, J. Kent, and M. Sørensen. Normal Variance-Mean Mixturesand z Distributions.
International Statistical Review / Revue Internationale deStatistique , 50(2):145–159, 1982.[33] C. Beck and E.G.D. Cohen. Superstatistics.
Physica A , 322(C):267–275, 2003.
EFERENCES
Journal of Physics A ,42(17):175207, 2009.[35] Anthony Paul Doulgeris and Torbjørn Eltoft. Scale Mixture of Gaussian Modellingof Polarimetric SAR Data.
EURASIP Journal on Advances in Signal Processing ,2010, 2010.[36] Florence Forbes and Darren Wraith. A new family of multivariate heavy-taileddistributions with variable marginal amounts of tailweight: Application to robustclustering.
Statistics and Computing , 24, 11 2014.[37] JP Bouchaud and M Potters.
Theory of Financial Risks, From Statistical Physicsto Risk Management . Cambridge University Press, New York, 2000.[38] M.L. Mehta.
Random Matrices . ISSN. Elsevier Science, 2004.[39] J. Wishart. The Generalised Product Moment Distribution in Samples From aNormal Multivariate Population.
Biometrika , 20A(1-2):32–52, 1928.[40] Laurent Laloux, Pierre Cizeau, Jean-Philippe Bouchaud, and Marc Potters. NoiseDressing of Financial Correlation Matrices.
Physical Review Letters , 83:1467–1470,Aug 1999.[41] Laurent Laloux, Pierre Cizeau, Marc Potters, and Jean-Philippe Bouchaud.Random Matrix Theory and Financial Correlations.
International Journal ofTheoretical and Applied Finance , 03(03):391–397, 2000.[42] Vasiliki Plerou, Parameswaran Gopikrishnan, Bernd Rosenow, Lu´ıs A. Nunes Ama-ral, and H. Eugene Stanley. Universal and Nonuniversal Properties of Cross Cor-relations in Financial Time Series.
Physical Review Letters , 83:1471–1474, Aug1999.[43] Vasiliki Plerou, Parameswaran Gopikrishnan, Bernd Rosenow, Luis A. NunesAmaral, Thomas Guhr, and H. Eugene Stanley. Random Matrix Approach toCross Correlations in Financial Data.
Physical Review E , 65:066126, Jun 2002.[44] Szil´ard Pafka and Imre Kondor. Estimated correlation matrices and portfoliooptimization.
Physica A , 343(C):623–634, 2004.[45] Marc Potters, Jean-Philippe Bouchaud, and Laurent Laloux. FinancialApplications of Random Matrix Theory: Old Laces and New Pieces.
Acta PhysicaPolonica B , 36:2767–2784, 09 2005.[46] S Drozdz, J Kwapien, and P Oswiecimka. Empirics Versus RMT in FinancialCross-Correlations.
Acta Physica Polonica B , 39(1):4027–4039, 2008.[47] Jaroslaw Kwapien, Stanislaw Drozdz, and P. Oswiecimka. The bulk of the stockmarket correlation matrix is not pure noise.
Physica A , 359:589–606, 05 2006.[48] Giulio Biroli, Jean-Philippe Bouchaud, and Marc Potters. The Student ensembleof correlation matrices: eigenvalue spectrum and Kullback-Leibler entropy.
ActaPhysica Polonica B , 38:4009–4026, 10 2007.
EFERENCES
Physica A , 343:694–700, 04 2001.[50] Zdzis(cid:32)law Burda, Romuald A. Janik, Jerzy Jurkiewicz, Maciej A. Nowak, GaborPapp, and Ismail Zahed. Free random L´evy matrices.
Physical Review E , 65:021106,Jan 2002.[51] Gernot Akemann and Pierpaolo Vivo. Power law deformation of Wishart–Laguerreensembles of random matrices.
Journal of Statistical Mechanics , 2008(09):P09002,2008.[52] Zdzislaw Burda, Andrzej Jarosz, Maciej A. Nowak, Jerzy Jurkiewicz, G´abor Papp,and Ismail Zahed. Applying free random variables to random matrix analysis offinancial data. Part I: The Gaussian case.
Quantitative Finance , 11(7):1103–1124,Jul 2011.[53] J.B. French, P.A. Mello, and A. Pandey. Ergodic behavior in the statistical theoryof nuclear reactions.
Physics Letters B , 80:17–19, 1978.[54] R. U. Haq, A. Pandey, and O. Bohigas. Fluctuation Properties of Nuclear EnergyLevels: Do Theory and Experiment Agree?
Physical Review Letters , 48:1086–1089,Apr 1982.[55] J.J.M. Verbaarschot and T. Wettig. Random Matrix Theory and Chiral Symmetryin QCD.
Annual Review of Nuclear and Particle Science , 50(1):343–410, 2000.[56] P.C. Mahalanobis. On the generalised distance in statistics.
Proceedings of theNational Institute of Science of India , 2:49–55, 1936.[57] Steven H. Simon and Aris L. Moustakas. Eigenvalue density of correlated complexrandom Wishart matrices.
Physical Review E , 69:065101, Jun 2004.[58] Zdzis(cid:32)law Burda, Jerzy Jurkiewicz, and Bart(cid:32)lomiej Wac(cid:32)law. Spectral moments ofcorrelated Wishart matrices.
Physical Review E , 71:026111, Feb 2005.[59] M. R. McKay, A. J. Grant, and I. B. Collings. Performance Analysis ofMIMO-MRC in Double-Correlated Rayleigh Environments.
IEEE Transactionson Communications , 55(3):497–507, 2007.[60] Daniel Waltner, Tim Wirtz, and Thomas Guhr. Eigenvalue Density of the DoublyCorrelated Wishart Model: Exact Results.
Journal of Physics A , 48:175204, 122014.[61] Zdzislaw Burda and Jerzy Jurkiewicz. Chapter 13, Heavy-tailed Random Matrices.In Gernot Akemann, Jinho Baik, and Philippe Di Francesco, editors,
OxfordHandbook of Random Matrix Theory , page 270. Oxford University Press, 09 2015.[62] Peter Forrester and Manjunath Krishnapur. Derivation of an eigenvalue probabilitydensity function relating to the Poincare disk.
Journal of Physics A , 42, 2009.[63] Tim Wirtz, Daniel Waltner, Mario Kieburg, and Santosh Kumar. The CorrelatedJacobi and the Correlated Cauchy-Lorentz ensembles.
Journal of StatisticalPhysics , 162:495, 01 2016.
EFERENCES submitted for publication , 2020.[65] Carl Ludwig Siegel. ¨Uber Die Analytische Theorie Der Quadratischen Formen.
Annals of Mathematics , 36(3):527–606, 1935.[66] Yan V. Fyodorov. Negative moments of characteristic polynomials of randommatrices: Ingham–Siegel integral as an alternative to Hubbard–Stratonovichtransformation.
Nuclear Physics B , 621(3):643–674, Jan 2002.[67] I.S. Gradsteyn and I.M.Ryzhik.