Towards a General Theory of Extremes for Observables of Chaotic Dynamical Systems.

In this paper we provide a connection between the geometrical properties of the attractor of a chaotic dynamical system and the distribution of extreme values. We show that the extremes of so-called physical observables are distributed according to the classical generalised Pareto distribution and derive explicit expressions for the scaling and the shape parameter. In particular, we derive that the shape parameter does not depend on the chosen observables, but only on the partial dimensions of the invariant measure on the stable, unstable, and neutral manifolds. The shape parameter is negative and is close to zero when high-dimensional systems are considered. This result agrees with what was derived recently using the generalized extreme value approach. Combining the results obtained using such physical observables and the properties of the extremes of distance observables, it is possible to derive estimates of the partial dimensions of the attractor along the stable and the unstable directions of the flow. Moreover, by writing the shape parameter in terms of moments of the extremes of the considered observable and by using linear response theory, we relate the sensitivity to perturbations of the shape parameter to the sensitivity of the moments, of the partial dimensions, and of the Kaplan-Yorke dimension of the attractor. Preliminary numerical investigations provide encouraging results on the applicability of the theory presented here. The results presented here do not apply for all combinations of Axiom A systems and observables, but the breakdown seems to be related to very special geometrical configurations.

Introduction

Extreme value theory (EVT) is gaining more and more prominence in a vast range of scientific fields because of its theoretical relevance in mathematical and physical sciences, and because it addresses the problem of understanding, modeling, evaluating risk factors such as those related to instabilities in the financial markets and to natural hazards related to seismic, climatic and hydrological extreme events. Even if the probability of extreme events is very low and decreases quickly with their magnitude, the associated risks can dominate over those coming from events belonging to the bulk of the statistics. An extensive account of recent results and relevant applications is given in [1].

The first comprehensive discussion of EVT dates back to the fundamental work by Gnedenko [2], who investigated the distribution of the maxima of a sample of independent identically distributed (i.i.d) stochastic variables. He showed that under very general conditions such maxima are distributed according to the so-called Generalised Extreme Value (GEV) distribution. The classic way of dealing with the statistical inference of extremes actually follows quite precisely the steps of the Gnedenko’s theorem. One partitions the experimental time series into bins of fixed length, extracts the maximum of each bin, and fits these data to the GEV distribution family using, typically, methods such as maximum likelihood estimation (MLE) or L-moments. See [3] for a detailed account of this methodology. It is possible to deal with extremes by taking a different point of view, i.e., by defining extremes as the events exceeding a given threshold. In the limit of very high threshold, we expect that the extremes are distributed according to the Generalized Pareto Distribution (GPD) introduced by Pickands III [4] and Balkema and De Haan [5]. In the case of i.i.d. variables, it is well known that a strong connections exists between the two methodologies. As shown in [6], we have that if block maxima obey the GEV distribution, then exceedances over some high threshold obey an associated GPD. Nonetheless, it is apparent that, whereas the two approaches are equivalent in the asymptotic limit, the GPD approach provides more reliable and more robust results when realistic, finite time series are considered (see, e.g., [7]).

A Brief Recapitulation of Extreme Value Theory

Gnedenko [2] studied the convergence of maxima of i.i.d. variableswith cumulative distribution function (cdf) where and are normalizing sequences and . Under general hypothesis on the nature of the parent distribution of data, Gnedenko [2] showed that the asymptotic distribution of maxima, up to an affine change of variable, belongs to a single family of generalized distribution called GEV distribution whose cdf can be written as:whereThis expression holds for , using (location parameter) and (scale parameter) as scaling constants, and is the shape parameter (also called the tail index). When , the distribution corresponds to a Gumbel type (Type 1 distribution). When the index is positive, it corresponds to a Fréchet (Type 2 distribution); when the index is negative, it corresponds to a Weibull (Type 3 distribution).

We briefly mention the Pareto approach to EVT. We define an exceedance as , which measures by how much the variable exceeds a given threshold . As discussed above, under the same conditions under which the block maxima of the i.i.d. stochastic variables obey the GEV statistics, the exceedances are asymptotically distributed according to the Generalised Pareto Distribution. Defining , we have that , where is the right endpoint of the distribution function for with [4, 6]:where the range of is if and if .

The relation between GEV and GPD parameters has been already discussed in literature in case of i.i.d variables [7-10]. It is first interesting to note thatfor , where the latter condition implies [11]. If we consider the upper range , we have that is only slightly smaller than 1, so that Eq. 4 implies thatso that the two distributions are asymptotically equivalent. This simple result is actually equivalent to the rather cumbersome formulas given in Coles [9] and Katz et al. [8] (and reported also by us in [12] for defining the correspondence between the parameters of the GPD and GEV distributions describing the statistics of extreme events extracted from the same data series.

Extreme Value Theory for Dynamical Systems

Recently, a great deal of attention has focused on understanding to what extent EVT can be applied to study the extreme of observables of deterministic dynamical systems. The main applications-driven motivation for this renewed interest probably comes from the spectacular development of numerical modeling in a geophysical fluid dynamical context and from the need to assess the ability of climate model to reproduce the observed statistical properties of extremes in present climate conditions and understand how they will change in an altered climate [13]. Other related applications include the numerical simulation of hydrological risk and of the production of electric energy from wind. It is clear that the matter is far from being trivial: numerical experiments on climate models of various degrees of complexity have shown that the speed of convergence (if any) of the statistical properties of the extremes definitely depends strongly on the chosen climatic variable of interest [3, 14–16].

Apart from these specific albeit very relevant applications, this problem has been addressed by the mathematical and statistical physical community. A first important result is that when a dynamical system has a regular (periodic of quasi-periodic) behaviour, we do not expect, in general, to find convergence to GEV distributions for the extremes of any observable, see Nicolis et al. [17], and Haiman [18]. Instead, if one chooses specific observables and considers dynamical systems obeying suitable mixing conditions, which guarantee the independence of subsequent maxima, it is possible to prove that the distribution of the block maxima of the observables converge to a member of the GEV family. The observables are expressed as , where is a function of the distance of the orbit from a point in the attractor , usually taken as the initial condition, such that has a global maximum for . The specific member of the GEV family (which is determined by the sign of the shape parameter) the maxima distribution converges to depends on the specific choice of . The paper by Collet [19] can be considered the cornerstone for the subsequent results obtained in the last few years [20-22]. The resulting parameters of the GEV distributions can be expressed as simple functions of the local dimension of the attractor. These results have been shown to be accurately detectable in numerical experiments when considering finite time series [23-25]. If, instead, the maxima are clustered, so that they feature a relative strong short-time correlation, the results have to be modified by introducing the extremal index [26].

Recently, it has been shown how to obtain results which are independent on whether the underlying dynamics of the system is mixing or, instead, regular . The key ingredient relies on using the Pareto rather than the Gnedenko approach. Such a shift in the point of view on extremes allows to derive results that do not dependent on whether extremes feature strong time-correlations or not. Assuming only that the local measures scales with the local dimension [27], it is possible to obtain by direct integration a GPD for the threshold exceedances of the observables introduced in [20–22, 26] when considering a generic orbit of a dynamical systems. In fact, the Pareto approach entails sampling all the points of the orbit that are very close to , thus sampling the local scaling of the invariant measure. With a suitable choice of , the resulting of the GPD is proportional to the inverse of local dimension [12]. The results obtained using the Pareto approach agree exactly with what was derived using the Gnedenko approach under the assumption of mixing dynamics in [20-25]. When the underlying system is mixing enough, the dynamical (Gnedenko approach) and the geometrical (Pareto approach) points of view on extremes give the same results, whereas differences emerge if the dynamics is such that strong time correlations exist between block-maxima. The selection of the extremes of the observables discussed above—using either the Gnedenko or the Pareto approaches—acts as magnifying lens near the initial condition, and that’s why one can extract information on the local dimension. Therefore, this provides a potentially viable alternative to e.g. the Grassberger–Procaccia algorithm [28] for the investigation of the scaling properties of the invariant measure of a chaotic attractor.

The results discussed above feature a major drawback when considering their relevance in many applications. Extreme events correspond to close returns of the orbit to to its initial condition. While relevant problems in natural sciences can be set in the framework of this class of observables (e.g. the classic problem of weather analogues in meteorology, already discussed by Lorenz in connection to the problem of predictability [29]), this is not the typical structure of the observables encountered in many applications, such as the case of total energy or enstrophy of a fluid flow. Recently, this problem has been addressed in Holland et al. [30], who have studied, using the GEV approach, whether EVT applies for physical observables of maps obeying the mixing conditions proposed in Freitas and Freitas [20], Freitas et al. [21, 26], Gupta et al. [22]. They consider a general observable reaching its maximum value restricted to the support of the invariant measure in (assuming for simplicity that such a point is unique), and assume that . Note that is indeed not of the form discussed above. They find that the block maxima of are asymptotically distributed according to a member of the GEV distributions, where the shape parameter is negative and can be written as a simple function of the partial dimensions along the stable and unstable manifolds at the point. This seems indeed a very relevant result, as it provides a universal property of extremes, regardless of the specific functional form of the observable .

Goals of the Paper

In this paper, we consider a GPD approach to EVT and try to complement and improve the results presented in [30] regarding the physical observables and those presented in [12] regarding the distance observables. We focus our attention on Axiom A systems [31], which are a special class of dynamical systems possessing a Sinai–Ruelle–Bowen (SRB) invariant measure [32] and featuring hyperbolicity in the attracting set. Such invariant measure coincides with the Kolmogorov’s physical measure, i.e. it is robust against infinitesimal stochastic perturbations. Another important property of Axiom A systems is that it is possible to develop a response theory for computing the change in the statistical properties of any observable due to small perturbations to the flow [33, 34]. Such response theory has recently been the subject of intense theoretical [35, 36], algorithmic [37] and numerical investigations [38-41] and is gaining prominence especially for geophysical fluid dynamical applications. Moreover, the response theory seems to provide powerful tools for studying multiscale systems and deriving parametrizations of the impact of the fast variables on dynamics of the slow [42, 43]. Finally, an important property of Axiom A systems is that, while the dynamics of natural or artificial systems is definitely not Axiom A in general, Axiom A systems can be considered as good ’effective’ models of actual systems with many degrees of freedom thanks to the so-called chaotic hypothesis, which is the extension in the non-equilibrium framework of the classic ergodic hypothesis [44]. Moreover, as discussed in [41], when we perform numerical simulations we implicitly assume that the system under investigation is Axiom A or Axiom A-equivalent. Therefore, considering Axiom A systems seems a good mathematical framework in view of providing results useful for a large spectrum of applications. The choice of considering Axiom A systems is instrumental in the derivation of various results on the relationship of EVT parameters to the dynamical and geometrical properties of the system, and will allow addressing the problem of the sensitivity of extremes to small perturbations of the system. The dependence of the properties of extremes of parametric modulations of the underlying dynamics is an issue of relevant theoretical as well as applicative interest. The practical interest stems from the fact that it is relevant to be able to control or predict variations in extreme events due to small perturbations to the dynamics. The theoretical interest comes from the fact that when considering extremes, universal parametric probability distributions can be defined, as opposed to the case of the bulk statistical properties. Because of this, we may hope to reconstruct the parameters descriptive of the EVT from simple moments of the distributions, express these in terms of observables of the system, and use the Ruelle response theory for expressing rigorously the sensitivity of extremes to small perturbations to the dynamics.

In Sect. 2 we show that by direct integration it is possible to derive the value of the two GPD parameters and , and, in particular, that the value of agrees with the GEV shape parameter obtained in [30]. We also show that, combining the results obtained using such physical observable and the distance observables considered in [12], it is possible to derive the estimates of the partial dimensions of the attractor along the stable and the unstable directions of the flow. In Sect. 3, we develop a linear response theory describing the impact of small time-independent -perturbations to the flow on the statistical properties of the extremes of the observable . We will first investigate the sensitivity of suitable defined observables describing above-threshold occurrences. We will focus on computing the changes of the shape parameter . We will find two equivalent expression for the sensitivity of with respect to . First, we will provide an expression for the sensitivity of in terms of the first two moments of the probability distribution of above-threshold events. Such expression entails a combination of observables of the Axiom A system, so that one can use Ruelle’s theory to compute the response to -perturbations to the dynamics. Nonetheless, we show that mathematical problems emerge when considering some limits. We follow the same approach for defining an expression for the sensitivity of the extremes of the distance observables considered in [12]. Then, we will relate the sensitivity of to the sensitivity of the Kaplan–Yorke dimension of the attractor. We will link our results to the well-established fact that both the extremes of observables and quantities like the Lyapunov exponents feature a good degree of regularity with respect to perturbations when one considers intermediate complexity to high-dimensional chaotic dynamical systems. In Sect. 4 we present the results of some numerical experiments performed using simple Hénon maps [45], aimed at providing support to our results. In Sect. 5 we briefly discuss the problems one faces when multiple time scales are present in the system. In Sect.  6 we comment our findings and present our conclusions.

Extreme Value Theory for Physical Observables of Axiom A Systems

Geometry of the Problem

Let us consider a continuous-time mixing Axiom A dynamical system on a compact manifold , where , with initial condition and evolution operator, is defined for all . Let us define as the attracting invariant set of the dynamical system , so that is the associated SRB measure with support . Let us now consider a smooth observable whose maximum restricted to the support of is unique, so that , , and is, moreover, not a critical point, so that , where the gradient is taken in . Therefore, we have that the the neutral manifold and the unstable manifold are tangent to the manifold in . We have that the intersection between the manifolds and is the empty set if . We define as the subset of included between the manifolds and . Furthermore, we define as the subset of included between the hyperplane tangent to the manifold in and the hyperplane , which is obtained by applying the translation given by the vector to the hyperplane , where . As , which is the limit of our interest, we have that and become undistinguishable up to leading order. See Fig. 1 for a a depiction of this geometrical construction.

Fig. 1

A low-dimensional cartoon of the geometrical construction used for deriving the EVL for exceedances above the threshold for the observable such that is realized for . a The manifolds and are depicted, together with the the attracting invariant set and the two hyperplanes and . is tangent to in and is obtained from via translation along . b The hyperplanes and delimit the region . Its intersection with is . c The manifolds and delimit the region . Its intersection with is . As , we have that

More in general, we denote as , with , the intersection between and the subset of included between the hyperplane and , where , is obtained from by applying to it the translation given by the vector .

It is now clear that we observe an exceedance of the observable above each time the systems visits a point belonging to .

We define the exceedances for the points as . An exceedance corresponds geometrically to a distance from , and the maximum exceedance corresponds to a distance between and .

Therefore, . The probability of observing an exceedance of at least given that an exceedance occurs is given by:where we have used the ergodicity of the system. Obviously, the value of the previous expression is 1 if . The expression contained in Eq. (5) monotonically decreases with (as if ) and vanishes when .

Derivation of the Generalised Pareto Distribution Parameters for the Extremes of a Physical Observable

We now wish to understand how to express the numerator and denominator of Eq. 5 as a function of , , and . We follow some of the ideas presented in Holland et al. [30] and use the fact that we are considering Axiom A systems. We define as the local dimension around , such thatIn order to proceed with the derivation of an extreme value law, such an asymptotics is not sufficient. In order to overcome some of the difficulties discussed in [12, 25], one needs to assume thatwhere is a slowly varying function of as , possibly depending on . We use the symbol as follows. We say that for if , . Taking the assumption above, one can derive the extreme value laws for the distance observables discussed in [12]. Let’s now estimate as a function of in the case of generic quadratic tangency between the the hyperplane and the unstable manifold in

In the case of Axiom A systems, since the invariant measure is SRB, we have almost everywhere on the attractor [46, 47], where is the Hausdorff dimension. As discussed in [46], we have that in this case all of the generalized Renyi dimensions have the same value. Moreover, we can conjecture that , where is the Kaplan–Yorke dimension [48]:where the ’s are the Lyapunov exponents of the systems, ordered from the largest to the smallest, is such that is positive and is negative. Following [46], we can also write , where , and are the dimensions of the attractor restricted to the stable, unstable and neutral directions, respectively, at the point . We have that is equal to the number of positive Lyapunov exponents , , for Axiom A systems is unitary, while is given by . Note that we can also express as followswhere indicates the integer part of , while its fractional part , because the last term gives a positive contribution smaller than 1, and and are both integer.

We follow the construction proposed by Holland et al. [30] for low dimensional maps. We derive the result by considering the following heuristic argument. Near , the attractor could be seen as the cartesian product of a multidimensional paraboloid (of dimension and of a fractal set of dimension immersed in . Note that this excludes for example conservative chaotic systems, whose attractor has the same dimension of the phase space, and systems that can be decomposed into a conservative part and a purely contractive part, whose attractor also has integer dimension. The mass of the paraboloid , where is the distance from the minimum. Instead, for each point of the paraboloid, the mass of corresponding fractal set , where is the distance along the the cartesian projection and is a slowly varying function of as , possibly depending on . This is where the non-trivial slowly varying pre-factor appears relevant.

In our case, and : for the former relation we assume a generic relation between the direction of the gradient of and the stable directions, while the latter relation results from the functional form of the paraboloid, see also Fig. 1. Hence, we obtain that whereand is also a slowly varying function of its argument. This construction can be made more formal by considering the disintegration of the SRB measure along the stable and unstable directions of the flow [48].

As a side note, we emphasize that in the case of more general tangencies between the the unstable and neutral manifold and the manifold , we have that , where ( ) gives the number of directions along the unstable (neutral) manifold where a tangency of order with the manifold is found. We obviously have that .

We continue our discussion considering the case of generic tangency. Following the same argument as above, we have that , where . We defineand obtain . Using the definition of slowly varying function and considering Eq. 5, we derive that in the limit :Note that the corresponding cdf is given by . Comparing Eqs. 3 and 9, one obtains that belongs to the GPD family, and that the GPD parameters can be expressed as follows: These results complement and extend what obtained in [30] using the GEV approach.

It is important to remark that Eq. 9 has been obtained in the limit of , and under the assumption that is a regularly varying function of degree D as . When considering a finite range , one should not expect deviations of the empirical distributions of extremes of from what prescribed in Eq. 9, which are intrinsic to the fractal nature of the measure. See also discussion and Fig. 1 in [12]. When finite ranges for are considered, one expects that, in some averaged sense, Eq. 9 fits well the distributions of extremes of and Eqs. 10 and 11 give the value of the two relevant parameters of the GPD, analogously to the idea that the number of points of the attractor at distance smaller than a small but finite from the point scales approximately, on the average as .

Comments

Equation 10 provides a very valuable information,as it shows that the shape parameter of the GPD does not depend on the considered observable, but only on the dimensions of the stable and of the unstable manifold. Moreover, the shape parameter is always negative, which is hardly surprising as we are considering compacts attractor and a well-behaved observable, whose values on the attractor have an upper bound. Note that for Axiom A systems, and are constant almost everywhere in the attractor , so that the information gathered for has a global value. Therefore, the expression for is universal, in the sense that we can gather fundamental properties of the dynamical system by looking at the shape parameters of the extremes of a generic observables with the properties described above. Note also that can be used to put upper and lower bounds to the Kaplan–Yorke dimension of the system, as and , so that . On the other hand, these inequalities can be read as constraints to the shape parameters of the extremes of a general observable for a system for which we know the Kaplan–Yorke dimension: .

On the other hand, the expression we obtain for provides clear support for calling it the scale parameter. We derive, as anticipated, that and we observe that it is proportional to the actual range of values considered as extremes of the observable , by incorporating the difference between the absolute maximum of the observable and the selected threshold . Therefore, if we consider as observable , with and take as threshold for the value , we have that and . In physical terms, changes if we change the unit of measure of the observable, whereas does not. More generally, we can make the following construction. Let’s define . If we select an observable , with , differentiable and positive in a sufficiently wide neighbourhood around so to ensure monotonicity of near , we get and , where .

From the Extremes to the Partial Dimensions Along the Stable and Unstable Directions of the Flow

It is worth considering the following strategy of investigation of the local properties of the invariant measure near , where . By performing statistical inference of the extremes of we can deduce as a result of the data fitting the best estimate of . If, following [12], we select as observable, e.g. , , we have that the extremes of the observable feature as shape parameter and scale parameter [12], where is a constant and is the chosen threshold.

We can then easily derive: where, as discussed above, we can take . Therefore, using rather general classes of observables, we are able to deduce the partial dimensions along the stable and unstable manifolds, just by looking at the properties of extremes related to . It is to be noted that, as clear from the results presented in [12], similar conclusions can be drawn choosing powers of forms for are possible. Note that, more generally, and can be deduced from the knowledge of any pair of values , , , and .

Expressing the Shape Parameter in Terms of the GPD Moments and of the Invariant Measure of the System

We consider the physical observable . We denote bythe density corresponding to the cumulative distribution given in Eq. 3. We can express its first two moments as follows: It is easy to derive thatandwhere we indicate explicitly that we refer to the observable and we have introduced the index of dispersion , the ratio between the variance and the squared first moment of the considered stochastic variable.

We now try to connect the previous formulas to the properties of the invariant measure of the dynamical system. As we know, the GPD is the exact asymptotic model for the extremes of the observable , so that we can express the results in terms of the conditional invariant measure as follows: where is the usual Heaviside distribution and, in general,We then obtain the following expression for the shape and dispersion parameters, respectively:andwhere these result are exact in the limit for . As a check, it is useful to verify that the right hand side of Eq. 22 gives the same general result as given in Eq. 10. By definition we have:We derive the following expression using repeatedly the distributional relation : so thatandBy plugging these expression into Eq. 22, we indeed obtain , which agrees with Eq. 10.

We also note that it is possible to generalize the results given in Eqs. 25–28. Using the fundamental theorem of calculus, it is possible to derive that:Moreover, it is notable that what presented in this subsection can be replicated step by step for the distance observables observable discussed above. We obtain:andwhere the quantities , are constructed analogously to how described in Eq. 21.

We wish to remark that Eqs. 22–23 and Eqs. 29–30 could in fact provide a very viable method for estimating the GPD parameters from data, since moments estimator are in general more stable than maximal likelihood methods, and then deriving the value of and using Eqs. 12–13.

Response Theory for the Extremes of General Observables

We wish to present some ideas on how to use response theory and the specific expressions given in Eqs 10–11 to derive a response theory for extremes of physical and distance observables in Axiom A dynamical systems. Let’s assume that we alter the Axiom A dynamical system under consideration as , where is a small parameter and is a smooth vector field, so that the evolution operator is transformed as and the invariant measure is transformed as . Ruelle’s response theory allows to express the change in the expectation value of a general measurable observable as a perturbative series as , with indicating the order of perturbative expansion, whereis the expectation value of over the perturbed invariant measure anddefines the unperturbed expectation value of . The term corresponding to the perturbative order of expansion is given by , where can be expressed in terms of the time-integral of a suitably defined Green function [33]. At this stage, we limit ourselves to the linear response of the system. We consider the following useful formula:and take into account the case.

Sensitivity of Shape Parameter as Determined by the Changes in the Moments of the Distribution

We wish to propose a linear response formula for the parameter using Eqs. 19–22. We start by considering that in Eq. 22 the shape parameter is expressed for every as a function of actual observables of the system. Unfortunately, in order to apply Ruelle’s response theory, we need the observables to be smooth, which is in contrast with the presence of the in the definition of the terms . Nonetheless, replacing the ’s with a smooth approximation , the Ruelle response theory can be rigorously applied. We now consider a sequence of approximating such that the measure of the support of is smaller than . It is reasonable to expect that as , the effect of the smoothing becomes negligible, because a smaller and smaller portion of the extremes is affected, and the response of the smoothed observable approaches that of . Therefore, we can retain the in the definition of the and define rigorously for every :andBy expanding the derivative in Eq. 31, the previous expression can be decomposed in various contributions entailing the linear response of the system to the perturbation for the observables , , and their values in the unperturbed case for .

We wish to remark the special relevance of the observable , which is normalizing factor in Eqs. 19–20, and, in practice, measures the fraction of above--threshold events. Therefore, once is chosen, the sensitivity of with respect to informs on whether the -perturbation to the vector flow leads to an increase or decrease in the number of extremes. We obtain:where is the derivative of the function with all the caveats discussed above. The result can be interpreted as follows. In the last formula, is the adjoint of the tangent linear of the unperturbed flow, and indicates the partial derivative with respect to the variable . At each instant we consider, in the unperturbed system, all the trajectories starting in the infinite past from points distributed according to the invariant measure such that the observable has value equal to . For each of these trajectories, we can measure whether the presence of the perturbation field would lead to a decrease or increase in at time . Summing over all trajectories, we get whether there is a net positive or negative change in the above threshold events at time . We integrate over and get the final result. Considering the geometrical construction given in Fig. 1, the previous formula can also be approximated as follows:Therefore, Eqs. 31–32 provides recipes for computing the sensitivity of and at for any case of practical interest, where is indeed finite, because in order to collect experimental data or process the output of numerical simulations we need to select a threshold which is high enough for discriminating true extremes and low enough for allowing a sufficient number of samples to be collected for robust data processing. Note that all statistical procedures used in estimating GPD parameters from data are actually based on finding a reasonable value for such that both conditions described above apply by testing that parameters’ estimates do not vary appreciably when changing [4]. So, if when investigating the extremes of for the unperturbed and -perturbed dynamics we find a common value of such that the GPD statistical inference of extremes is satisfactory, Eqs. 31–32 provide correct formulas for the sensitivities.

We wish to underline that apparently formal problems emerge when taking the limit in Eqs.31–32 for higher and higher values of . It is indeed not clear at this stage whetherexists, because we cannot apply the smoothing argument presented above in the limit of vanishing . Moreover, it is not clear whether such limit is equal towhich seems at least as well suited for describing the change of the shape observable given in Eq. 22 due to an -perturbation in the dynamics. Obviously, if the two limits given in Eqs. 35 and 36 exist and are equal, then a rigorous response theory for can be established. Same applies when considering the properties of .

The same derivation and discussion can be repeated for the observables introduced above and we can derive the corresponding formulas for and , where the relevant limit for is .

Let’s try to give a more intuitive interpretation to the results given above. Let’s consider Eq. 17 and assume that, indeed, is differentiable with respect to . We have:which implies that the sensitivity of the shape parameter is half of the opposite of the sensitivity of the inverse of the index of dispersion . Therefore, a positive sensitivity of the index of dispersion (larger relative variability of the extremes of the observable with positive values of implies a larger value (closer to 0) of , and so the possibility that larger and larger extremes are realized. Same interpretation applies for the observables.

Sensitivity of the Shape Parameter as Determined by the Modification of the Geometry

In the previous subsection we have shown that the Ruelle response theory supports the idea that the shape parameters descriptive of the extremes of both the physical observables and the distance observables change with a certain degree of regularity when considering -perturbations to the dynamics.

In this subsection, we wish to look at the sensitivity of extremes with respect to perturbation from another angle, i.e. through the relationship between the shape parameters and and the partial dimension of the attractor along the stable, neutral and unstable manifolds of the underlying dynamical system, see Eqs. 12–13. As long as the -perturbation is small, the modified dynamical system belongs to the Axiom A family, so that the results presented above apply. Therefore, we can write in more general terms: In the following, we introduce somewhat carelessly derivatives with respect to of quantities that are not, a priori, differentiable. The main point we want to make is that if and are differentiable with respect to , then various quantities describing the structure of the attractor are also differentiable. Therefore, the existence of the limits given in Eqs. 35 and 36 (and their equivalent for the observables) would have far-reaching consequences. We will discuss the obtained results at the end of the calculations. Another caveat we need to mention is that Eqs. 38–39 are in general true almost anywhere, so that we may have to interpret the derivatives in this section in some suitable weak form.

It seems relevant to add the additional hypothesis of strong transversality for the unperturbed flow, which is equivalent to invoking structural stability [48]. We take such pragmatic point of view and proceed assuming that derivatives with respect to are well defined. Linearizing the dependence of on around in Eq. 38, we obtain:We have that , as, thanks to structural stability, small perturbations do not alter the qualitative properties of the dynamics, and cannot change in a step-wise way the number of expanding or neutral directions. We now separate the quantity into its integer component and the rest, which is generically non vanishing:Only the last term is different from zero, because, thanks to structural stability, . Using Eq. 7, we obtain:where is defined as in Eq. 6. Since is not altered for infinitesimal perturbations to the dynamical system, we have that, using Eq. 6:Expanding the previous expressions, the final formula reads as follows:This implies that the shape parameter increases, thus attaining a value closer to zero ( is always negative) when the perturbation increases the Kaplan–Yorke dimension of the attractor, so, in qualitative sense, if it favors “forcing” over “dissipation”. This matches quite well, at least qualitatively, with the discussion following Eq. 37.

We have that, when considering a distance observable of the form , following the same steps described above one gets the following result:such result can be easily generalized by considering the class of observables described in [12].

Combining Eq. 31 with Eq. 44, and the derivative with respect to of Eq. 29 with Eq. 45, we can derive two expressions for the derivative of the sensitivity of the Kaplan–Yorke dimension at : where we take the limit for in Eq. 46 and in Eq. 47.

The previous results imply that if one of , or the Kaplan–Yorke dimension of the underlying Axiom A system change smoothly with -perturbations to the dynamics, so do the other two quantities. This may suggest ways to study the regularity of the Kaplan–Yorke dimension by resorting to the analysis of the regularity of a much simpler expressions involving moments of given observables.

This result provides useful insight also not considering the problematic limits discussed above. Taking a more qualitative point of view, this suggests that when considering small perturbation in the dynamics of chaotic systems behaving like Axiom A systems there is a link between the presence (or lack) of regularity of the parameters describing the extremes of a wide class of observables and of the regularity of the Kaplan–Yorke dimension. This matches with the fact that numerical practice with moderate to high dimensional strongly chaotic systems shows that, actually, the parameters describing extremes of energy like observables and the Lyapunov exponents both have—within numerical precision—a smooth dependence on the parameters of the system. See [3, 49] for an extensive discussion in a simplified yet relevant fluid dynamical model. A detailed investigation of the apparent regularity for all practical purposes of the Lyapunov exponents with respect to small perturbations in the dynamics of intermediate complexity to high-dimensional models has been presented in [50]. We also wish to remark that if these regularity hypotheses were not satisfied, the very widespread (and practically successful) procedure of parametric tuning of high-dimensional models of natural, engineered or social phenomena would be absolutely hopeless, and delicate numerical procedures such as those involved in data assimilation of high-dimensional dynamical systems would lack any sort of robustness, contrary to the accumulated experience.

Numerical Experiments

As thoroughly discussed in [23-25], it is far from trivial to devise suitable numerical experiments for studying to what extent the theoretically derived asymptotic extreme value laws for dynamical systems can be detected in finite datasets obtained as outputs of simulations. In this Section we would like to present simple numerical experiments providing some heuristic and preliminary support to the fact that the universal properties for the extremes can be observed when considering specific observables for dynamical systems. More detailed numerical studies, where inference of the geometrical properties of the attractor is performed using the statistics of extremes of suitable observables will be reported elsewhere. Hence, we consider as toy model the smooth map introduced by Hénon [45]: As well known, depending on the value of the two parameters and , the Hénon map can exhibit either regular or chaotic behaviour, where, in the latter case, the invariant measure is supported on a strange attractor [51]. We consider two sets of parameter values for which chaotic behavior is observed, and . In the first case, the largest Lyapunov exponent and the Kaplan–Yorke dimension is estimated as , where and . In the second case, the largest Lyapunov exponent and the Kaplan–Yorke dimension is estimated as . Note that it is reasonable to expect that the Hénon maps considered here do not possess an SRB measure, because the considered pair of values of and do not seem to belong to the Benedicks–Carleson set of parameters. As a consequence, these systems are not exact dimensional and the local dimension does not have the same value almost everywhere on the attractor, so, in rigorous terms, we have no a-priori reasons to expect that our results should necessarily apply. Nonetheless, in order to assess the robustness of our findings, it is interesting to check to what extent our theoretical predictions are met, at least qualitatively, in such a basic model of chaotic dynamics.

We proceed as follows for both pairs of parameters and . We first choose, for sake of simplicity, the observable , where . The initial conditions are selected in the basin of attraction of the strange attractors. We perform long integrations (order of iterations) and select the maximum value of , which we denote as , and define as the unique point belonging to the attractor such that . We then construct the observable , which measures the distance between the orbit and the point . As discussed in [12], the asymptotic properties of the extremes (maxima) of the observable allow to derive easily the local dimension . We then repeat the investigation using, instead, the observable . In all the analyses presented below, we have chosen extremely high thresholds for studying the statistical properties of the extremes of the and observables, in such a way to include only about a fraction of about or less of the total number of points of the orbit. All the results are insensitive to choice of , which suggests that we are well into the asymptotic regime.

The results obtained for the Hénon system featuring are shown in Fig. 2, where we present the complementary cumulative distribution of excesses (see Eqs. 5 and 9) for () and for the corresponding in panel a) (panel b)). The empirical values of for the and observables are shown by the blue and red curves, respectively, and the power law behavior given by the theory (assuming Axiom A properties!) are shown by the black and magenta lines, respectively.

Fig. 2

Results of numerical simulations performed on the Hénon map with parameters’ value and . a Blue curve empirical for the observable , with . Black line power law behavior deduced from the theory. Red curve empirical for the observable , with . Magenta line power law deduced from the theory. b Same as a, for the observable , with and . c Approximation to the attractor with blow-ups of the portions of the invariant measure corresponding to the extremes of the observables ( regions); the vertical lines indicate the thresholds. In both inserts, we consider . See also Fig. 1 (Color figure online)

The error bars on the empirical (estimated by varying the initial conditions of the simulation) are for almost all values of so small that they cannot be graphically reproduced. Instead, the flat region obtained for very low values of results from the finiteness of the sampling and gives the baseline uncertainty. Note that the straight lines are obtained out of the theoretical predictions, without any procedure of optimization or of fit, so that no uncertainties are involved. The empirical and theoretical distributions obey the same normalization.

We first observe that the local dimension in the vicinity of both ’s is extremely close to the , as scales to a very good approximation with an exponent ; compare the red curves and the magenta lines. Note that, considering that the local dimension has rather large variations across the attractor of the Hénon system, such a correspondence was not intentionally pursued. However, these are favorable circumstances to check the theory. We find that the distributions for the observables and also obey accurately the power law scaling with exponent introduced in Eq. 8, compare the blue curves and the black lines. In panel c) we present a simple description of the geometry of the problem, by showing an approximation to the map’s attractor with blow-ups of the portions of the invariant measure corresponding to the extremes of the observables (the regions introduced in Fig. 1). Even if the geometrical properties of the regions of the attractor around the two ’s seem indeed different, when zooming in, the two regions look similar. The presence of many parabolas-like smooth curves stacked according to what looks qualitatively like a Cantor set fits with the comments and calculations given in Sect. 2.

In Fig. 3 we report the corresponding results obtained for the Hénon system featuring . By looking at the empirical of the observables, we note that also in this case the local dimension is close to the value of for both extremal points ’s (compare the red curves and the magenta lines in panels a and b). Nonetheless, the slope of the empirical data is slightly steeper than the theoretical value. The agreement between the predicted value of the power law scaling for the of the observables is not as good as in the case reported in Fig. 2. The predicted scaling exponent seems to overestimate the very large extremes. Nonetheless, a power law scaling is apparent for the empirical . Note that the bias between the theoretical and empirical scalings is of the same sign for both the and observables, suggesting that also for the observables part of the disagreement is due to the discrepancy between the local dimension and the Kaplan–Yorke dimension (there is a shift in the values of the slopes). Also here, panel c provides an approximate representation of the attractor of the system, and, in particular of the regions: by comparing it with panel c of Fig. 3, and considering that they contain the same number of points, one can intuitively grasp that the local dimension is lower in this case.

Fig. 3

Same as in Fig. 2, but for parameters’ value and . In this case in a and , and in b and

We would like to emphasize that in panels a and b for Figs. 2 and 3, we observe deviations of the empirically obtained from the power law behaviour, in the form of fluctuations above and below a straight line in a – plot (this is quite clear in Fig. 3). As discussed in Sect.  2 B, the presence of such modulations across scales result from the fact that gaps are present along the stable manifold containing , with a Cantor set-like structure. See the inserts in Figs. 2c and 3c, where the stable manifold (not shown) is, as opposed to the unstable manifold, not orthogonal to the gradient of (the direction, in this case). So, when we integrate the density of states along the direction of the gradient of the observable starting in in order to obtain and , we get a factor () coming from the (local) paraboloidal form of the unstable manifold discussed in Sect. 2, times a devil’s staircase which can, on the average, be approximated by the power law . The same geometric arguments apply when considering integrations along the spherical shells centered in for constructing the extreme value laws for observables. Smooth approximation to devil’s staircases, appear in a closely related context when using the GEV approach for studying extreme values laws in random dynamical systems whose attractor is the actual Cantor set [25].

Such preliminary results suggest that it is indeed promising to use the combined statistical properties of the extremes of physical and distance observables for determining the geometry of the attractor in terms of its partial dimensions along the stable and the unstable manifold.

Multiple Time Scales

We briefly wish to mention here some additional features which may appear and be extremely relevant at finite time in practical cases, where the dynamics can deviate from Axiom A when certain time scales are considered. Let’s assume that one can to a first approximation partition the (unique) attractor of a chaotic dynamical system into, say, two pieces, so that the system has two time scales, a short one related to the transitive dynamics within each of the two pieces, and a long one corresponding to intermittent jumps from one to the other piece. In this case, if we observe the system for a time scale intermediate between the two, the properties of the extremes will depend only on the properties of the visited portion of the attractor and we will observe a Weibull distribution, as discussed here, as the dynamics may be Axiom A-like. When our observation time nears the long time scale, we might observe extraordinary large events, corresponding to excursions directed towards the other piece of the attractor, until a jump, corresponding to an irreversible (on the short time scale) transition will take place. Such extraordinary events will not fit the Weibull law found on smaller time scales, because they result from the global properties of the attractor, which have not been sampled yet. Therefore, in these intermediate scale, the results proposed here will not be valid. Instead, one may interpret such extraordinary events as Dragon Kings [55], which will manifest as outliers spoiling the Weibull statistics and pushing the statistics of extremes towards an (apparently) unphysical Frechét distribution. Observing extremes over even longer time scales, so that the orbit visits many times both parts of the attractor, we shall recover a Weibull law, which reflects the global properties of the attractor. In a system with these properties, small perturbations to the dynamics might impact substantially the long time scale discussed above, with the result of having a high sensitivity of the statistics of extremes when a fixed time window of observation is considered. The presence of such strong sensitivity has been proposed as a method for detecting precursors of global stability thresholds related to critical transitions [56].

Conclusions

This paper has addressed the problem of studying the EVT for general observables of mixing Axiom A dynamical system. We have set in a common framework the investigation of general properties of distance observables , for which we had derived some basic results in [12], and of physical observables, first discussed in [30]. By physical observables we mean rather general classes of smooth observables (e.g. quadratic, energy-like quantities) which achieve an absolute maximum at a specific, non-critical point (the gradient does not vanish) of the attractor of the system.

We have built up from the recent results of Holland et al. [30], who have studied accurately this problem on discrete maps with specific mixing properties using the GEV framework and have derived, as fundamental result, that the extremes of A indeed obey a GEV extreme value law and that , the shape parameter of the distribution, is equal to , with , where is the partial dimension along the stable manifold and is the partial dimension along the unstable manifold.

In this paper, using the GDP approach, thus considering exceedances above a given threshold , and considering the physically relevant case of mixing Axiom A systems, we derive through direct integration that the shape parameter can also be expressed as . We have framed our results for continuous flows, so , where is the dimension along the neutral direction and is unitary. In the case of discrete maps, we obtain the same results as in [30]. We have also been able to derive the explicit expression for the scale parameter .

It is clear that the parameter is always negative (so that the distribution of extremes is upper limited), reflecting the fact that the observable is smooth and the attractor is a compact set. Moreover, measuring allows us to provide an upper and lower bound for and vice versa, because , or, conversely . In particular, we have that is small and negative if and only if the Kaplan–Yorke dimension of the attractor is large. If we consider a chaotic system with a high dimensional attractor (e.g. in the case of an extensive chaotic system with many degrees of freedom), we derive that . This may well explain why in a multitude of applications in natural sciences such as hydrology, meteorology, oceanography the special member of the GPD family—the exponential model—given in Eq. 3 usually gives a good first guess of the statistics of observed extremes [9].

Alternatively, this result suggests that if we perform a statistical analysis using the POT method (using an empirical threshold ) of the extremes for a high-dimensional chaotic system and obtain as a result of the statistical inference of the collected data for the GPD model a shape parameter or , we should conclude that our sample is not yet suited for an EVT statistical fit. This may depend on the fact that we have selected an insufficiently stringent value for . Obviously, choosing higher values for implies that we need to have longer time series of the observable under investigations.

Interestingly, by combining the expression for obtained in this paper for a physical observable and the expression for of the GPD describing the extremes of observable of the form , where dist is the distance function [12], we can express the partial dimension of the attractor on the stable and unstable dimension as simple functions of and . The straightforward result is that and . The same can be obtained using, instead, the parameters and . This provides further support to the idea that extremes can be used as excellent diagnostic indicators for the detailed dynamical properties of a system. The message seems to be that one can construct observables whose large fluctuations give precise information on the dynamics. While considering various sorts of anisotropic scalings of the neighborhood of a point of the attractor allows to derive its partial dimensions [52], the specific result we obtain here is that choosing an arbitrary physical observable and studying its extremes, we automatically select a special, non ellipsoidal neighborhood, where the degree of anisotropy between the stable and unstable (and neutral) directions is generically universal and given by the factor 1/2.

We wish to make an additional remark. Let’s assume, instead, that the gradient of is vanishing in and that at leading order near , where is a negative definite symmetric matrix and the square brackets indicate the scalar product. It is clear that, apart from a linear change in the coordinates and rescaling, the statistical properties of the extremes of will match those of .

In the second part of the paper we have tackled the problem of studying how the properties of extremes of the observable change when an -perturbation is added to the system. As theoretical framework, we have taken the point of view of Ruelle [33, 34], who has shown that the SRB measure of Axiom A systems is differentiable with respect to -perturbations to the dynamics and has provided explicit formulas for studying how the expectation values of generic observables of Axiom A systems change when the system is subjected to perturbations.

We have used the fact that the GPD is an exact asymptotic model for extreme events in order to find a simple functional relation between and (as well as and ) and the first two moments of the probability of above-threshold exceedance expressed in terms of the invariant measure of the system. These expressions are amenable to direct treatment with Ruelle’s response theory, at least when we do not consider the limit but stick to the practical situation where we need to consider a finite range for the extremes. The differentiability properties of and are hard to ascertain in the limit. We have also found an explicit expression for the sensitivity of the number of extremes—seen an over threshold events—to the -perturbation.

We have then taken into consideration our results on the relationship between and and the partial dimensions of the attractor. Interestingly, it seems that there is an intimate connection between the differentiability with respect to of , and , so that either all of them or none of them is differentiable with respect to . Under the hypothesis of differentiability, we have been able to derive that the sensitivity of and with respect to is proportional to the sensitivity of the Kaplan–York dimension with respect to . Specifically, we obtain that if the perturbation tends to increase the dimensionality of the attractor (thus, in physical terms, favoring forcing over dissipation), the value of becomes closer to zero, so that the occurrence of very large extreme events becomes more likely. The system, in this case, has more freedom to perform large fluctuations.

Taking a more pragmatic point of view, these results at least provide a rationale for the well-known fact that in moderate to high-dimensional strongly chaotic systems the Kaplan–Yorke dimension (and, actually, all the Lyapunov exponents) change smoothly with the intensity of the perturbating vector field, as discussed in [49, 50], and similar behavior is found for the parameters describing the extremes of energy-like quantities [3].

Since it has been shown in [36] that the linear response vanishes for any observable in the case of stochastic forcing of rather general nature, whereas the second order response gives the leading order of perturbation, we expect that adding a moderate noise to a chaotic dynamical system will not alter significantly the shape parameter describing the EVT of both physical and distance observables.

Finally, we have performed a set of simple numerical experiments using the celebrated Hénon map for two different pairs of parameters—(a,b) = (1.4,0.3) and (1.2,0.3). While these maps are definitely not Axiom A, it seemed to us worthwhile to test the robustness of the theory in more general classes of systems and to have an indication of whether the asymptotic properties discussed here are practically observable. One has to keep in mind that, when considering extremes, the approach to asymptotic behavior is far from being trivial to detect in finite datasets [24, 25]. We find encouraging agreement between our theory and the outputs of numerical experiments for both sets of parameters, which suggests that it is worthwhile to study accurately more comprehensive models in order to see whether one can practically derive the geometrical properties of the attractor from the statistics of extremes of distance and physical observables. Further numerical investigations are needed for studying whether it is possible to find satisfactory numerical evidence for the response theory for the extremes developed here and in particular for the relationship between the sensitivity of the EVT’s parameters of physical and distance observables and the sensitivity of the Kaplan Yorke dimension to perturbations to the underlying dynamical system.

In this work we have considered the case where the observable has a unique maximum restricted to in . If and share some symmetries, is not unique, and instead there is a set of points ’s belonging to , finite or infinite, depending of the kind of symmetries involved, where reaches its maximum value restricted to . Let’s consider the relevant case where and share a discrete symmetry, so that the set of the maximal point ’s, has finite cardinality. The results discussed here for the extremes of will nonetheless apply, because we can perform an equivalent geometrical construction as in Fig. 1 for each element of . When we consider an -perturbation to the dynamics which respects the discrete symmetry, it is clear that all the results of the response presented here apply. Finally, one can deduce that if the considered perturbation, instead, breaks the discrete symmetry, the results presented here will still be valid as the break of the degeneracy will make sure that only one of the ’s (or a subset of , if the corresponding perturbed vector flow obeys to a a subgroup of the original symmetry group) still accounts for the extreme events of .

This work may constitute a new theoretical viewpoint for studying extremes in a rather general setting and understanding how they change when the dynamical system is slightly perturbed. The results on how to express the sensitivity of the Kaplan–Yorke dimension with respect to an -perturbation to the dynamics seems also useful. Our findings might find applications in many sectors of physical, engineering and social sciences, and, just to provide a basic example of crucial relevance, in the investigation of the impact of climate change on climate extremes.

Obviously, as in the case of all results pertaining to EVT, it is important to test numerically the practical verification of the findings presented here. In this paper we have extensively discussed the relevance of finite-time effects in the selection of the extremes of physical observables and in the definition of the relevant sensitivities. We need to mention that, recently, some renormalization group methods have been applied for deriving systematic finite size corrections to extreme value laws [53, 54]. These results seem extremely promising and might lead to improved methods for fitting extreme value statistics to given datasets.

We need to remark that such results have been derived using some intuitive geometrical construction and assuming generic relations between the direction of the gradient of at and the stable directions. It is possible to devise special pairs of Axiom A systems and observables such that the strange attractors do not fulfill such generic conditions. One can easily construct a situation where the gradient of is orthogonal also to stable manifold by immersing the attractor in a higher dimensional space and taking observables defined on such a space. In this case, the factor appearing in Eq. 8 will affect also the stable dimensions. Nonetheless, the results that for high-dimensional systems the distribution of extremes is indistinguishable from the Gumbel as the shape parameters tends to zero from below is not be affected by this correction.

We believe that typical combinations of Axiom A systems and observable functions allow for the generic conditions to be obeyed. We still need to understand how to frame consistently such a concept of genericity, which obviously differs from the traditional one, which focuses either on the observables, or on the systems. This should be the subject of theoretical investigation and accurate numerical testing.