Universal behavior of extreme value statistics for selected observables of dynamical systems
Abstract
The main results of the extreme value theory developed for the investigation of the observables of dynamical systems rely, up to now, on the Gnedenko approach. In this framework, extremes are basically identified with the block maxima of the time series of the chosen observable, in the limit of infinitely long blocks. It has been proved that, assuming suitable mixing conditions for the underlying dynamical systems, the extremes of a specific class of observables are distributed according to the so called Generalized Extreme Value (GEV) distribution. Direct calculations show that in the case of quasiperiodic dynamics the block maxima are not distributed according to the GEV distribution. In this paper we show that, in order to obtain a universal behaviour of the extremes, the requirement of a mixing dynamics can be relaxed if the Pareto approach is used, based upon considering the exceedances over a given threshold. Requiring that the invariant measure locally scales with a well defined exponent  the local dimension , we show that the limiting distribution for the exceedances of the observables previously studied with the Gnedenko approach is a Generalized Pareto distribution where the parameters depends only on the local dimensions and the value of the threshold. This result allows to extend the extreme value theory for dynamical systems to the case of regular motions. We also provide connections with the results obtained with the Gnedenko approach. In order to provide further support to our findings, we present the results of numerical experiments carried out considering the wellknown Chirikov standard map.
pacs:
Valid PACS appear hereI Introduction
Extreme value Theory was originally introduced by Fisher and Tippett (1928) and formalised by Gnedenko (1943), who showed that the distribution of the maxima of a sample of independent identically distributed (i.i.d) stochastic variables converges under very general conditions to a member of the socalled Generalised Extreme Value (GEV) distribution. The attention of the scientific community to the problem of understanding extreme values theory is growing, also because this theory is crucial in a wide class of applications for defining 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 decreases with their magnitude, the damage that they may bring increases rapidly with the magnitude as does the cost of protection against them. From a theoretical point of view, extreme values of observables are related to large fluctuations of the corresponding underlying system. An extensive account of recent results and relevant applications is given in (Ghil et al., 2011).
The traditional (’Gnedenko’) approach for the statistical inference of extremes is related to the original results by Gnedenko (1943): we partition the experimental time series into bins of fixed length, we extract the maximum of each bin, and fit the selected data to the GEV distribution family using methods such as maximum likelihood estimation (MLE) or Lmoments. See (Felici et al., 2007) for a detailed account of this methodology. The selection of just one maximum in a fixed period may lead to the loss of relevant information on the large fluctuations of the system, especially when there are many large values in a given period (Castillo and Hadi, 1997). This problem can be taken care of by considering several of the largest order statistics instead of just the largest one. For such maxima distributions we expect convergence to the Generalized Pareto Distribution (GPD) introduced by Pickands III (1975) and Balkema and De Haan (1974) to model the exceedances over a given threshold. We call this the ’Pareto’ approach. Also this approach has been widely adopted for studying empirically natural extreme phenomena such as those related to waves, winds, temperatures, earthquakes and floods (Simiu and Heckert, 1996; Bayliss and Jones, 1993; Pisarenko and Sornette, 2003).
Both the Gnedenko and Pareto approaches were originally designed to study extreme values for series of i.i.d variables. In this case it is well known that a strong connections exists between the two methodologies, as we have that if block maxima obey the GEV distribution, then exceedances over some high threshold will have an associated GPD. Moreover, the shape parameter of the GPD and that of the corresponding GEV distribution are identical (Leadbetter et al., 1983). As a result, several practical methods (e.g. Hill’s and Pickands’ estimators) developed for estimating the shape parameter of the GEV distribution of the extremes of a given time series are actually based upon comparing the GPD fits at various thresholds (Embrechts et al., 1997; Coles, 2001). In practical terms, it appears that, while the Gnedenko and Pareto approaches provide equivalent information in the asymptotic limit of infinitely long time series, the GPD statistics is more robust when realistic, finite time series are considered (see, e.g., (Ding et al., 2008)).
In recent years, especially under the influence of the rapid development of numerical modelling in the geophysical sciences and of its applications for the investigation of the socioeconomic impacts of extreme events, it has become of great relevance to understand whether it is possible to apply the extreme value theory on the time series of observables of deterministic dynamical systems. Carefully devised 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 on the chosen climatic variable of interest (Kharin et al., 2005; Vannitsem, 2007; Felici et al., 2007; Vitolo et al., 2009a).
Several papers have addressed this issue at a more general level. A first important result is that when a dynamical system has a regular (periodic of quasiperiodic) behaviour, we do not expect, in general, to find convergence to GEV distributions for the extremes of any observable. These results have been presented by Balakrishnan et al. (1995), and more recently, by Nicolis et al. (2006) and by Haiman (2003).
A different mathematical approach to extreme value theory in dynamical systems was proposed in the landmark paper by Collet (2001), which has paved the way for the recent results obtained in the last few years (Freitas and Freitas, 2008; Freitas et al., 2009, 2010; Gupta et al., 2009). The starting point of all of these investigations has been to associate to the stationary stochastic process given by the dynamical system, a new stationary independent sequence which obeys one of the classical three extreme value laws introduced by Gnedenko (1943). The assumptions which are necessary to observe a GEV distribution in dynamical systems rely on the choice of suitable observables (specific functions of the distance between the orbit and the initial condition, chosen to be on the attractor) and the fulfillment of particular mixing conditions that guarantee the independence of subsequent maxima. Recent studies have shown that the resulting parameters of the GEV distributions can be expressed as simple functions of the local (around the initial condition) dimension of the attractor, and detailed numerical investigations have clarified the conditions under which convergence to the theoretical GEV distributions can be satisfactorily achieved when considering finite time series (Faranda et al., 2011a, b, c).
In this paper, we wish to attempt a unification of these two lines of work by using the Pareto rather than the Gnedenko approach. We choose the same class of observables presented in (Freitas and Freitas, 2008; Freitas et al., 2009, 2010; Gupta et al., 2009; Faranda et al., 2011a, b, c) and show that, assuming only that the local measures scales with the local dimension (Bandt, 2007), it is possible to obtain by direct integration a GPD for the threshold exceedences when considering a generic orbit of a dynamical systems, without requiring any special mixing properties. The parameters will depend only on the choice of the threshold and, more importantly, on the local dimension. Note that Castillo and Hadi (1997) had already pointed out that in the case of periodic or quasiperiodic motion the Gnedenko approach to the evaluation of the extreme value statistics is inefficient, basically because in the limit of very large blocks, we tend to observe always the same maximum in all bins. To support our analytical results we provide numerical experiments that we carry out considering the classic Chirikov standard map (Chirikov and Shepelyansky, 2008). This paper is organised as follows. In section 2 we recapitulate the extreme value theory for dynamical systems obtained using the Gnedenko approach. In section 3 we present our general results obtained using the Pareto approach. In section 4 we provide support to our investigation by examining the results of the numerical simulations performed on the standard map. In section 5 we present our final remarks and future scientific perspectives.
Ii Gnedenko Approach: Generalized Extreme Value distributions in dynamical systems
Gnedenko (1943) studied the convergence of maxima of i.i.d. variables
with cumulative distribution function (cdf) where and are normalizing sequences and . Under general hypothesis on the nature of the parent distribution of data, Gnedenko (1943) 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:
(1) 
where
(2) 
This expression holds for , using (location parameter) and (scale parameter) as scaling constants in place of , and (Pickands III, 1968), in particular, in Faranda et al. (2011b) we have shown that and , where 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).
Let us consider a dynamical system , where is the invariant set in some manifold, usually , is the Borel algebra, is a measurable map and an invariant Borel measure. In order to adapt the extreme value theory to dynamical systems, following (Freitas and Freitas, 2008; Freitas et al., 2009, 2010; Gupta et al., 2009), we consider the stationary stochastic process given by:
(3) 
where ’dist’ is a distance in the ambient space , is a given point and is an observable function. The partial maximum in the Gnedenko approach is defined as:
(4) 
Defining , we consider the three classes of observables :
(5)  
(6)  
(7) 
where is a constant and . Using the observable we obtain convergence of the statistics of the block maxima of their time series obtained by evolving the dynamical system to the Type distribution if one can prove two sufficient conditions called and , which basically imply a sort of independence of the series of extremes resulting from the mixing of the underlying dynamics (Freitas and Freitas, 2008). The conditions cannot be simply related to the usual concepts of strong or weak mixing, but are indeed not obeyed by observables of systems featuring a regular dynamics or powerlaw decay of correlations.
A connection also exists between the existence of extreme value laws and the statistics of first return and hitting times, which provide information on how fast the point starting from the initial condition comes back to a neighborhood of , as shown by Freitas et al. (2009) and Freitas et al. (2011). In particular, they proved that for dynamical systems possessing an invariant measure , the existence of an exponential hitting time statistics on balls around almost any point implies the existence of extreme value laws for one of the observables of type described above. The converse is also true, namely if we have an extreme value law which applies to the observables of type achieving a maximum at , then we have exponential hitting time statistics to balls with center . Recently these results have been generalized to local returns around balls centered at periodic points (Freitas et al., 2010).
In Faranda et al. (2011b, c, a) we analised both from an analytical and numerical point of view the extreme value distribution in a wide class of low dimensional maps. We divided the time series of length of the observables into bins each containing the same number of observations, and selected the maximum (or the minimum) value in each of them (Coles et al., 1999). We showed that at leading order (the formulas are asymptotically correct for ), the GEV parameters in mixing maps can be written in terms of (or equivalently ) and the local dimension of the attractor . We have:

type observable:
(8) 
type observable:
(9) 
type observable:
(10)
Moreover, we clearly showed that other kind of distributions not belonging to the GEV family are observed for quasiperiodic and periodic motions.
Iii The Pareto approach: Generalized Pareto Distributions in Dynamical Systems
We define an exceedance as , which measures by how much exceeds the 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 (Leadbetter et al., 1983):
(11) 
where the range of is if and if . We consider the same set up described in the previous section and take into account an observables , such that achieves a maximum for (finite or infinite) and is monotonically decreasing. We study the exceedance above a threshold defined as . We obtain an exceedence every time the distance between the orbit of the dynamical system and is smaller than . Therefore, we define the exceedances . By the Bayes’ theorem, we have that . In terms of invariant measure of the system, we have that the probability of observing an exceedance of at least given that an exceedence occurs is given by:
(12) 
Obviously, the value of the previous expression is 1 if . In agreement with the conditions given on , the expression contained in Eq. (12) monotonically decreases with and vanishes when the radius is given by . Note that the corresponding cdf is given by . In order to address the problem of extremes, we have to consider small radii. At this regard we will invoke, and assume, the existence of the following limit
(13) 
where is the local dimension of the attractor (Bandt, 2007). Therefore, we rewrite the following expression for the tail probability of exceedance:
(14) 
where we have dropped the dependence of to simplify the notation. By substituting with specific observable we are considering, we obtain explicitly the corresponding extreme value distribution law.
By choosing an observable of the form given by either , , or , we derive as extreme value distribution law one member of the Generalised Pareto Distribution family given in Eq. (11). Results are detailed below:

type observable:
(15) 
type observable:
(16) 
type observable:
(17)
The previous expressions show that there is a simple algebraic link between the parameters of the GPD and the local dimension of the attractor around the point . This implies that the statistics of extremes provides us with a new algorithmic tool for estimating the local fine structure of the attractor. These results show that it is possible to derive general properties for the extreme values of the observables , , or independently on the qualitative properties of the underlying dynamics, be the system periodic, quasiperiodic, or chaotic. Therefore, by taking the Pareto instead of the Gnedenko approach, we are able to overcome the mixing conditions (or the requirements on the properties of the hitting time statistics) needed to derive a general extreme value theory for dynamical systems, as proposed in (Collet, 2001; Freitas and Freitas, 2008; Freitas et al., 2009, 2010; Gupta et al., 2009). In (Faranda et al., 2011b, c, a) we had proposed that the reason why a link between the extreme value theory and the local properties of the invariant measure in the vicinity of the point can be explained by the fact that selecting the extremes of the observables , , or amounts to performing a zoom around . In the case of the Gnedenko approach, such a picture is accurate only if the dynamics is mixing (time and spatial selection criteria are equivalent). Instead, in the case of the Pareto approach, this is literally what we are doing when writing Eq. (14), as we are remapping the radius of the ball in a monotonic fashion though the inverse of the functions.
iii.1 Relationship between the Gnedenko and Pareto approaches
The relation between GEV and GPD parameters have been already discussed in literature in case of i.i.d variables (Katz et al., 2005; Coles, 2001; Ding et al., 2008; Malevergne et al., 2006). Coles (2001) and Katz et al. (2005) have proven that the cdf of the GEV defined as can be asymptotically written as that of GPD under a high enough threshold as follows:
(18)  
where , , and , with . In the present case, we have to compare Eqs. (8)(10) for GEV with Eqs. (15)(17) for GPD, keeping in mind that the GEV results hold only under the mixing conditions discussed before. While it is immediate to check that , the other relationships are valid in the limit of large , as expected.
Iv Numerical Investigation
The standard map (Chirikov, 1969) is an areapreserving chaotic map defined on the bidimensional torus, and it is one of the most widelystudied examples of dynamical chaos in physics. The corresponding mechanical system is usually called a kicked rotator. It is defined as:
(19) 
The dynamics of the map given in Eq. (19) can be regular or chaotic. For the motion follows quasi periodic orbits for all initial conditions, whereas if the motion turns to be chaotic and irregular. An interesting behavior is achieved when : in this case we have coexistence of regular and chaotic motions depending on the chosen initial conditions (Ott, 2002).
We perform for various values of ranging from up to an ensemble of 200 simulations, each characterised by a different initial condition randomly taken on the bidimensional torus, and we compute for each orbit the observables , . In each case, the map is iterated until obtaining a statistics consisting exceedances, where the threshold and . We have carefully checked that all the results are indeed robust with respect to the choice of the threshold and of the value of . For each orbit, we fit the statistics of the exceedances values of the observables to a GPD distribution, using a MLE estimation (Castillo and Hadi, 1997) implemented in the MATLAB function gpdfit (Martinez and Martinez, 2002). The results are shown in Fig. 1 for the inferred values of and and should be compared with Eqs. (15)(17). When , we obtain that the estimates of and are compatible with a dimension for all the initial conditions: we have that the ensemble spread is negligible. Similarly, for , the estimates for and agree remarkably well with having a local dimension for all the initial conditions. In the transition regime, which occurs for , the ensemble spread is much higher, because the scaling properties of the measure is different among the various initial conditions. As expected, the ensemble averages of the parameters change monotonically from the value pertaining to the regular regime to that pertaining to the chaotic regime with increasing values of . Basically, this measures the fact that the socalled regular islands shrink with . Note that in the case of the observable , the estimate of the is robust in all regimes, even if, as expected, in the transition between low and high values of the ensemble spread is larger. These results can also be compared with the analysis presented in Faranda et al. (2011a), where we used the Gnedenko approach. In that case, the values obtained in the regular regions were inconsistent with the GEV findings, the very reason being that the dynamics was indeed not mixing. Here, it is clear that the statistics can be computed in all cases, and we have a powerful method for discriminating regular from chaotic behaviors through the analysis of the inferred local dimension.
V Conclusions
The growing attention of the scientific community in understanding the behavior of extreme values have led, in the recent past, to the development of an extreme value theory for dynamical systems. In this framework, it has been shown that the statistics of extreme value can be linked to the statistics of return in a neighborhood of a certain initial conditions by choosing special observables that depend on the distance between the iterated trajectory and its starting point. Until now, rigorous results have been obtained assuming the existence of an invariant measure for the dynamical systems and the fulfillment of independence requirement on the series of maxima achieved by imposing and mixing conditions, or, alternatively, assuming an exponential hitting time statistics (Freitas and Freitas, 2008; Freitas et al., 2009, 2010; Gupta et al., 2009). The parameters of the GEV distribution obtained choosing as observables the function , defined above depend on the local dimension of the attractor and numerical algorithms to perform statistical inference can be set up for mixing systems having both absolutely continuous and singular invariant measures (Faranda et al., 2011b, c; Vitolo et al., 2009b; Holland et al., 2011). Instead, when considering systems with regular dynamics, the statistics of the block maxima of any observable does not converge to the GEV family (Balakrishnan et al., 1995; Nicolis et al., 2006).
Taking a complementary point of view, in this paper we have studied the statistics of exceedances for the same class of observables and derived the limiting distributions assuming only the existence of an invariant measure and the possibility to define a local dimension around the point of interest. To prove that the limiting distribution is a GPD we did not use any further conditions. In particular no assumptions on the mixing nature of the maxima sequence have been made. This means that a GPD limiting distribution holds for the statistics of exceedance of every kind of dynamical systems and it depends only on the threshold value and on the local dimension once we choose the observables , . Other functions can converge to the limiting behaviour of the GPD family if they asymptotically behave like the ’s (compare the discussion in (Freitas et al., 2009)). Nonetheless, this requires, analytically, to perform separately the limit for the threshold going to and that for the radius of the ball going to zero. In practical terms, this requires, potentially, much stricter selection criteria for the exceedances when finite time series are considered.
We note that, as the parameter is inversely proportional to , one can expect that each time we analyse systems of intermediate or high dimensionality, the distributions for and observables will be virtually indistinguishable from what obtained considering the observable: the is in this sense an attracting member of the GPD family. This may also explain, at least qualitatively, why the Gumbel ( for the GEV family) distribution is so efficient in describing the extremes of a large variety of natural phenomena (Ghil et al., 2011).
The universality of this approach allows to resolve the debate on whether there exists or not a general way to obtain information about extreme values for quasiperiodic motions raised in Balakrishnan et al. (1995) and Nicolis et al. (2006). This is due to the fact that considering several of the largest order statistics instead of just the largest one we can study orbits where numerous exceedance are observed in a given block, as it happens for systems with periodic or quasiperiodic behaviors. As an example, one may consider a system with multiple commensurable frequencies: choosing a block length larger or equal to the smallest common period, we select always the same value for any considered observable. On the other hand for mixing maps we find, as expected, an asymptotic equivalence of the results obtained via the Gnedenko and via the Pareto approach.
The Pareto approach provides a way to reconstruct local properties of invariant measures: once a threshold is chosen and a suitable exceedance statistics is recorded, we can compute the local dimension for different initial conditions taken on the attractor. This is true also in the opposite direction: if the knowledge of the exact value for the local dimension is available, once we chose a small enough radius (threshold), it is possible to compute a priori the properties of extremes without doing any further computations. In fact, the expression for the parameters Eqs. (15)(17) do not contain any dependence on the properties of the dynamics except the local dimension.
Besides the analytical results, we have proved that Pareto approach is easily accessible for numerical investigations. The algorithm used to perform numerical simulations is versatile and computationally accessible: unlike the GEV algorithm that requires a very high number of iterations to obtain unbiased statistics, using the Pareto approach we can fix a priori a value for the threshold and the number of maxima necessary to construct the statistics. With the simulations carried out on the standard maps, we obtain meaningful results with a much smaller statistics with respect to what observed when considering the Gnedenko approach.
We hope that the present contribution may provide a tool that is not only useful for the analysis of extreme events itself, but also for characterising the dynamical structure of attractors by giving a robust way to compute the local dimensions, with the new possibility of embracing also the case of quasiperiodic motions.
Future investigations will include the systematic study of the impact on the extreme value statistics of adding stochastic noise to regular and chaotic deterministic dynamical systems, and the use of the Ruelle response theory (Ruelle, 1998, 2009) to study the modulation of the statistics of extremes due to changes in the internal or external parameters of the system, especially in view of potentially relevant applications in geophysics such as in the case of climate studies (Abramov and Majda, 2007; Lucarini and Sarno, 2011).
Acknowledgements.
The authors acknowledge various useful exchanges with R. Blender and K. Fraedrich, and the financial support of the EUERC project NAMASTEThermodynamics of the Climate System.References
 Fisher and Tippett (1928) R. Fisher and L. Tippett, in Proceedings of the Cambridge philosophical society, Vol. 24 (1928) p. 180.
 Gnedenko (1943) B. Gnedenko, The Annals of Mathematics 44, 423 (1943).
 Ghil et al. (2011) M. Ghil et al., Nonlinear Process in Geophysics 18, 295 (2011).
 Felici et al. (2007) M. Felici, V. Lucarini, A. Speranza, and R. Vitolo, Journal of Atmospheric Science 64, 2137 (2007).
 Castillo and Hadi (1997) E. Castillo and A. Hadi, Journal of the American Statistical Association , 1609 (1997).
 Pickands III (1975) J. Pickands III, the Annals of Statistics , 119 (1975).
 Balkema and De Haan (1974) A. Balkema and L. De Haan, The Annals of Probability , 792 (1974).
 Simiu and Heckert (1996) E. Simiu and N. Heckert, Journal of Structural Engineering 122, 539 (1996).
 Bayliss and Jones (1993) A. Bayliss and R. Jones, Peaksoverthreshold flood database (Institute of Hydrology, 1993).
 Pisarenko and Sornette (2003) V. Pisarenko and D. Sornette, Pure and Applied Geophysics 160, 2343 (2003).
 Leadbetter et al. (1983) M. Leadbetter, G. Lindgren, and H. Rootzen, Extremes and related properties of random sequences and processes. (Springer, New York, 1983).
 Embrechts et al. (1997) P. Embrechts, C. Kluppelberg, and T. Mikosh, “Modelling Extremal Events,” (1997).
 Coles (2001) S. Coles, An introduction to statistical modeling of extreme values (Springer Verlag, 2001).
 Ding et al. (2008) Y. Ding, B. Cheng, and Z. Jiang, Advances in Atmospheric Sciences 25, 507 (2008).
 Kharin et al. (2005) V. Kharin, F. Zwiers, and X. Zhang, Journal of Climate 18, 5201 (2005).
 Vannitsem (2007) S. Vannitsem, Tellus A 59, 80 (2007).
 Vitolo et al. (2009a) R. Vitolo, P. Ruti, A. Dell’Aquila, M. Felici, V. Lucarini, and A. Speranza, Tellus A 61, 35 (2009a).
 Balakrishnan et al. (1995) V. Balakrishnan, C. Nicolis, and G. Nicolis, Journal of Statistical Physics 80, 307 (1995).
 Nicolis et al. (2006) C. Nicolis, V. Balakrishnan, and G. Nicolis, Physical review letters 97, 210602 (2006).
 Haiman (2003) G. Haiman, Statistics & Probability Letters 65, 451 (2003).
 Collet (2001) P. Collet, Ergodic Theory and Dynamical Systems 21, 401 (2001).
 Freitas and Freitas (2008) A. Freitas and J. Freitas, Statistics & Probability Letters 78, 1088 (2008).
 Freitas et al. (2009) A. Freitas, J. Freitas, and M. Todd, Probability Theory and Related Fields , 1 (2009).
 Freitas et al. (2010) A. Freitas, J. Freitas, and M. Todd, Arxiv preprint arXiv:1008.1350 (2010).
 Gupta et al. (2009) C. Gupta, M. Holland, and M. Nicol, Lozilike maps, and Lorenzlike maps, preprint (2009).
 Faranda et al. (2011a) D. Faranda, V. Lucarini, G. Turchetti, and S. Vaienti, Arxiv preprint arXiv:1107.5972 (2011a).
 Faranda et al. (2011b) D. Faranda, V. Lucarini, G. Turchetti, and S. Vaienti, Accepted for publication in Journal of Statistical Physics (2011b).
 Faranda et al. (2011c) D. Faranda, V. Lucarini, G. Turchetti, and S. Vaienti, Arxiv preprint arXiv:1106.2299 (2011c).
 Bandt (2007) C. Bandt, in Physics and theoretical computer science: from numbers and languages to (quantum) cryptography security, edited by J. Gazeau, J. Nešetřil, and B. Rovan (IOS Press, Amsterdam, 2007) pp. 91–112.
 Chirikov and Shepelyansky (2008) B. Chirikov and D. Shepelyansky, Scholarpedia 3, 3550 (2008).
 Pickands III (1968) J. Pickands III, The Annals of Mathematical Statistics 39, 881 (1968).
 Freitas et al. (2011) A. Freitas, J. Freitas, and M. Todd, (2011).
 Coles et al. (1999) S. Coles, J. Heffernan, and J. Tawn, Extremes 2, 339 (1999).
 Katz et al. (2005) R. Katz, G. Brush, and M. Parlange, Ecology 86, 1124 (2005).
 Malevergne et al. (2006) Y. Malevergne, V. Pisarenko, and D. Sornette, Applied Financial Economics 16, 271 (2006).
 Chirikov (1969) B. Chirikov, Institute of Nuclear Physics, Novosibirsk (1969).
 Ott (2002) E. Ott, Chaos in Dynamical Systems (Cambridge University Press, New York, 2002).
 Martinez and Martinez (2002) W. Martinez and A. Martinez, Computational statistics handbook with MATLAB (CRC Press, 2002).
 Vitolo et al. (2009b) R. Vitolo, M. Holland, and C. Ferro, Chaos: An Interdisciplinary Journal of Nonlinear Science 19, 043127 (2009b).
 Holland et al. (2011) M. Holland, R. Vitolo, P. Rabassa, A. Sterk, and H. Broer, Arxiv preprint arXiv:1107.5673 (2011).
 Ruelle (1998) D. Ruelle, Phys. Letters A 245, 220 (1998).
 Ruelle (2009) D. Ruelle, Nonlinearity 22, 855 (2009).
 Abramov and Majda (2007) R. Abramov and A. Majda, Nonlinearity 20, 2793 (2007).
 Lucarini and Sarno (2011) V. Lucarini and S. Sarno, Nonlin. Processes Geophys 18, 7 (2011).