Milan Stehlik
email: milan.stehlik@jku.at

Pedro Aguirre
email: pablo.aguirre@usm.cl

Stéphane Girard
email: stephane.girard@inria.fr

Pavlina Jordanova
email: pavlina_kj@abv.bg

Jozef Kiseľák
email: jozef.kiselak@upjs.sk

Sebastian Torres-Leiva
email: sebastian.torresle@alumnos.usm.cl

Z Sadovsky

Andres Rivera
email: arivera@cecs.cl

M Stehlík

P. c , Girard, S. d , Jordanova, P. e , Kise ľák, J. f , Torres Leiva, S. c , Sadovský, Z. g Aguirre

On ecosystems dynamics

Keywords: stochasticity, determinism, chaos, autonomous system, Guanaco glacier, snow extremes

We show how a dynamical system given by a t-score function for some class of monotonic data transformations generates consistent extreme value estimators. The variation of their values increases the uncertainty of proper assessment of climate change. Two important examples illustrate the methodology: mass balance measurements on Guanaco glacier, Chile, and extreme snow loads in Slovakia. We experience singular learning of the transitions in ecosystems.

Introduction

In the past several decades scientific effort has been focused on studying and understanding global climate changes. The effect of climatic changes has become more and more visible and in many regions of the world these changes are represented by increasing of weather extremes [START_REF] Chan | Typhoon Haiyan and beyond[END_REF], [START_REF] Coumou | A decade of weather extremes[END_REF], [START_REF] Tank | Trends in indices of daily temperature and precipitation extremes in Europe, 1946-99[END_REF].

All ecosystems (Methan [START_REF] Sabolova | Entropy based statistical inference for methane emissions released from wetland[END_REF], Guanaco Glacier [START_REF] Jordanova | Weak properties and robustness of t-Hill estimators[END_REF], Snow extremes [START_REF] Stehlík | Statistical Analysis related to Exceptional Snow Loads[END_REF]) are oscillating. Decomposition to deterministic, stochastic and chaotic part have been studied by [START_REF] Stehlík | Missing chaos in global climate change data interpreting?[END_REF]. We can understand contributions to oscillations in at least three following ways:

1) Extreme Value Index (EVI) ξ oscillates around 0 (it can have positive, negative or zero values). As [START_REF] Penalva | Extreme value analysis -A brief overview with an application to flow discharge rate data in a hydrometric station in the north of Portugal[END_REF] pointed out, difficulties may rise with the "regularity conditions" for the maximum likelihood (ML) estimation ( [START_REF] Smith | Maximum likelihood estimation in a class of nonregular cases[END_REF]), it is shown that the usual property of asymptotic normality holds provided the extreme value parameter ξ is larger than -0.5. For all environments we can consider ξ > -1. ( [START_REF] Penalva | Extreme value analysis -A brief overview with an application to flow discharge rate data in a hydrometric station in the north of Portugal[END_REF]) Recently, Zhou ([46,[START_REF] Zhou | The extent of the maximum likelihood estimator for the extreme value index[END_REF]) showed that the ML estimators verify the property of asymptotic normality for ξ > -1. The Second Order Regularity condition (SOC) can be difficult to be checked (or even satisfied) in practical application. E.g. if the observed random variable (r.v.) is power of Uniform or has power law behavior at the finite right end point (see Example 3.3.15 and 3.3.16, page 137, [START_REF] Embrechts | Modelling Extremal Events for Insurance and Finance[END_REF]), there is not unique SOC parameter ρ.

2) Aside of 1), the consistent estimators of tail parameters can be build up upon t-scores ( [START_REF] Jordanova | Weak properties and robustness of t-Hill estimators[END_REF]). The parameters of Harmonic mean estimators (HME) are consisting dynamical system which can surprisingly always find a monotonic representing data function (t-score function) η. This process contributes to deterministic dynamics of [START_REF] Stehlík | Missing chaos in global climate change data interpreting?[END_REF].

3) The use of Negative t-Hill(n-t-Hill) for estimation of the EVI index ξ < 0) can give several limiting behaviors, however, limits can be given by symmetric (normal) or classical (Weibull) distributions, which both are special cases of generalized gamma distribution (ggd), see [START_REF]Stehlík Homogeneity and Scale Testing of Generalized Gamma Distribution[END_REF].

The paper is organized as follows. In the next session we study autonomous system of t-score functions. In section 3 we study mass balance measurements from Guanaco glacier and we show that both negative and positive EVIs are obtained. In section 4 we study the extremal snow loads in Slovakia, again receiving both negative and positive EVIs. To maintain the readability of manuscript we put technicalities to Appendix.

Dynamical systems of t-score functions

The transformation-based score ( [START_REF] Fabián | Induced cores and their use in robust parametric estimation[END_REF], [START_REF] Stehlík | On the favourable estimation of fitting heavy tailed data[END_REF]) or shortly the t-score for the density f is defined as

T η (x; θ) = - 1 f (x; θ) d dx 1 η (x) f (x; θ) .
It expresses a relative change of a basic component of the density, i.e., density divided by the Jacobian of mapping η. The t-score is a suitable function for using the generalized moment method for the estimation of parameters of heavy-tailed distributions. Let X 1 , . . . , X n be independent identically distributed (i.i.d) sample from F with probability density function (p.d.f.) f . The parametric version of the so-called t-mean, which can be considered as a measure of central tendency of distributions, yields the moment estimation equations for θ in the form

1 n n i=1 T (x i ; θ) = 0.
The solution θ is strongly consistent and asymptotically normal (see [START_REF] Fabián | Induced cores and their use in robust parametric estimation[END_REF]). For t-Hill estimator [START_REF] Fabián | On robust and distribution sensitive hill like method[END_REF], we have bounded score

S(x; θ) = T η (x; θ) = θ 1 - θ + 1 θx ,
and for generalized t-Hill estimator [START_REF] Beran | The harmonic moment tail index estimator: asymptotic distribution and robustness[END_REF] (Pareto distribution and η(x) = ln(x -1), x > 1), we have the score

S(x; θ, β) = θ 1 -θ+β-1 θx β-1
, for β = 1,

θ 1 θ -ln x , for β = 1. ( 1
)
where β > 0 is tuning parameter. For β = 2 we receive t-Hill, with "typical" transformation of the support of the distribution (1, ∞) to the whole real line (-∞, ∞) is η(x) = ln(x -1), x > 1. Here an important inverse problem arises. For a given score S, does there exist one or several sufficiently smooth functions η such that equation

T η = S (2)
holds? Which qualitative properties do they posses? Consider now the Pareto distribution with the probability density function (p.d.f.)

f (x, θ) = θx -θ-1 , x > 1
where θ > 0 is a shape parameter (the tail index). Let us modify the equation ( 2) by multiplying by f > 0 in order to receive exact 2nd-order differential equation in the form

h(x) + d dx f (x) η (x) = 0, (3)
where h(x) = S(x; θ, β) f (x). Now, integrate equation ( 3) to obtain an equation, which is solvable by quadrature, of the form

H(x) + f (x) η (x) = C,
where H(x) is an antiderivative of h. Its form (under the condition β = 1θ) is:

H(x) = θ 2 1 - θ + β -1 θx β-1 x -θ-1 dx = θx -θ x 1-β -1 + C 1 .
This yields several classes1 of solutions expressible in general in the form of special functions (a nonelementary antiderivatives). But this is an obstacle, since they can be hardly jointly analyzed because of their transcendental nature. This difficulties motivate us to study equation [START_REF] Balakrishnan | The Exponential Distribution: Theory, Methods, and Applications[END_REF], by a different approach, applicable for general density f and score function S. In order to analyse it is more convenient to define some extra variables w = (x, y, z) := (t + a, η, η ), a ∈ supp(f ) = {x ∈ R, : f (x) = 0}. Under the assumption η = 0 equation ( 2) is equivalent to the system ẇ = W(x, y, z), where W(x, y, z) = (1, z, Ψ(x, z)),

Ψ(x, z) = z 2 S + z d dx ln(f (x)) and (x, y, z) ∈ D 0 := [a, ∞) × [a, ∞) × R \ {0}.
We use this approach in details for [START_REF] Beran | The harmonic moment tail index estimator: asymptotic distribution and robustness[END_REF], where a = 1, x ≥ 1 is the independent variable, η(x) ≥ 1 is the unknown function with η (x) = 0 and (β, θ) ∈ R + are parameters. Here (x, y, z) := (t + 1, η, η ). In this way, (3) is equivalent to the following set of autonomous ordinary differential equations:

     ẋ = 1, ẏ = z, ż = ϕ(x, z), (4)
where

ϕ(x, z) = - θ + 1 x z + θ 1 - θ + β -1 θx β-1 z 2 ,
and (x, y, z) ∈ D 0 .

In our setting any initial condition (x 0 , y 0 , z 0 ) ∈ D 0 defines a unique smooth solution of (4)and, hence, a unique differentiable solution y = η(x) of (3). Each solution of (4) can be represented as a smooth orbit {(x(t), y(t), z(t))} in R 3 parameterised by t ∈ R; see [START_REF] Guckenheimer | Nonlinear Oscillations, Dynamical Systems, and Bifurcations of Vector Fields[END_REF] for more details.

The (unique) orbit through a given point (x, y, z) ∈ D 0 is tangent to the vector (1, z, ϕ(x, z)) at the point (x, y, z). Hence, an orbit always flows forward in the direction of x and never "comes back" near any point already visited in the same orbit. More precisely, there is no dense orbit of ( 4) in any open region of the phase space R 3 . Hence, there cannot be topological mixing, which is one of the necessary ingredients of chaotic dynamics [START_REF] Guckenheimer | Nonlinear Oscillations, Dynamical Systems, and Bifurcations of Vector Fields[END_REF][START_REF] Hasselblatt | A First Course in Dynamics: With a Panorama of Recent Developments[END_REF].

For the fixed initial condition we are able to obtain monotonic solution for t-score for almost all possible cases of parameters. The t-score defines consistent estimator of tail parameter θ. The choice of parameter β is an issue of experience for the statistician/ecologists. Experienced choice of parameter β brings a proper trade-off between robustness and efficiency (see [START_REF] Beran | The harmonic moment tail index estimator: asymptotic distribution and robustness[END_REF]). The graph of y = η(x) in the (x, y)-plane is determined by the initial condition (x 0 , y 0 , z 0 ) at t = t 0 . In particular, from (4) it follows that z 0 is the initial slope η (x 0 ) = ẏ(t 0 ) = z 0 of this solution. By the continuity of the solutions of (4), since z = η (x) = 0, then sign(z 0 ) = sign (z(t)) for all t > 0; namely, the sign of z 0 determines the (constant) sign of ẏ = z. Hence, y(t) is a monotone function of t, and, hence, any solution η(x) of ( 3) is either an increasing or decreasing function for every x. This a is very important property, since monotone transformations of data are statistically optimal.

2.1.2. The set ϕ -1 (0)

In spite of η(x) being a monotone function, its graph could still have a number of different shapes, it could be bounded or unbounded, etc. Statistical reasons for having η a monotonic function are as follows: if we have a differentiable inverse of η we can compute the induced distribution after the transformation. From the ecological point of view, it is as parsimonious view on ecological dynamics that one can consider.

In particular, the level set ϕ -1 (0) in the (x, z)-plane determines the values of x where d 2 η(x) dx 2 = 0 and, hence, where η(x) has an inflection point. Moreover, the domain D 0 of system (4) can be continuously extended in order to include the plane z = 0. Hence, in what follows we consider a continuous extension of the domain D 0 given by

D 1 := [1, ∞) × [1, ∞) × R so that the set M 0 := {(x, z) ∈ [1, ∞) × R : z = 0} ⊂ ϕ -1 (0).
In order to study ϕ -1 (0), let us define the function

g(x) := θ + 1 x θ -(θ + β -1)x 1-β , (5)
and let gr(g) be the graph of z = g(x) in the domain [1, ∞) × R. Furthermore, let us define the following sets in the parameter space (β, θ):

Ω 1 = {(β, θ) ∈ R + : 0 < β < 1, θ + β -1 < 0}; (6)
Ω 2 = {(β, θ) ∈ R + : 0 < β < 1, θ + β -1 > 0}; (7)
Ω 3 = {(β, θ) ∈ R + : 1 < β < 2}; (8)
Ω 4 = {(β, θ) ∈ R + : 2 < β}; (9)
and their boundaries:

T 1 = {(β, θ) ∈ R + : θ + β -1 = 0}; (10)
T 2 = {(β, θ) ∈ R + : β = 1}; (11)
T 2 = {(β, θ) ∈ R + : β = 2}. ( 12
)
The sets Ω k , k = 1, 2, 3, 4 and T j , j = 1, 2, 3 are shown in panel (e) of Figure 1. From a statistical point of view, by choosing of system (4) we decided for a specific form of dynamical system, driven by autonomous system for t-scores of Pareto distribution. Here we consider statistical learning based on t-score function with monotonous transformation η. From an ecological point of view, this was a convenient and parsimonious approach to model the underlying dynamics for extreme value estimators under the statistical constraints of Pareto tail, and a monotonous smooth η. The usefulness of Lemma 1 and Figure 1 is that we can use the information on the nature of the set ϕ -1 (0) -as a function of parameters β and θ-to give some geometric insight into the possible shapes of η(x) in terms of slope and inflection points. This will be greatly illustrated in the next subsection.

Fact:

The following statements hold, see Appendix C.1:

Lemma 1.

1. If (β, θ) ∈ Ω 1 , the set ϕ -1 (0) = {M 0 , M 1 } consists of two branches, where M 1 = gr(g). The set ϕ -1 (0) is qualitatively as in the sketch of Figure 1 3. If (β, θ) ∈ Ω 2 , the set ϕ -1 (0) = {M 0 , M 1 , M 2 } consists of three branches, where M 1 ∪ M 2 = gr(g). The set ϕ -1 (0) is qualitatively as in the sketch of Figure 1(c). 4. If (β, θ) ∈ T 2 , the set ϕ -1 (0) = {M 0 } consists of the single branch M 0 . The set ϕ -1 (0) is qualitatively as in the sketch of Figure 1

(d). 5. If (β, θ) ∈ Ω 3 , the set ϕ -1 (0) = {M 0 , M 1 , M 2 } consists of three branches where M 1 ∪ M 2 =
gr(g). The set ϕ -1 (0) is qualitatively as in the sketch of Figure 1

(f ). 6. If (β, θ) ∈ T 3 , the set ϕ -1 (0) = {M 0 , M 1 , M 2 } consists of three branches where M 1 ∪ M 2 is the graph of z = x(θ + 1) θx 2 -(θ + 1)x
. The set ϕ -1 (0) is qualitatively as in the sketch of Figure

1(g). 7. If (β, θ) ∈ Ω 4 , the set ϕ -1 (0) = {M 0 , M 1 , M 2 } consists of three branches where M 1 ∪ M 2 =
gr(g). The set ϕ -1 (0) is qualitatively as in the sketch of Figure 1(h).

0 0 1 1 0 0 1 1 0 0 0 1 1 1 0 0 0 1 1 1 M 0 M 1 0 x z (a) Ω 1 0 0 1 1 0 0 1 1 0 0 0 1 1 1 0 0 0 1 1 1 M 0 M 1 0 x z (b) T 1 0 0 1 1 0 0 1 1 0 0 0 1 1 1 0 0 0 1 1 1 M 0 M 1 M 2 x a 0 z (c) Ω 2 0 0 0 1 1 1 0 0 0 1 1 1 0 0 1 1 0 0 1 1 M 0 0 x z (d) T 2 T 1 Ω 1 T 2 T 3 Ω 2 Ω 3 Ω 4 θ 0 1 1 2
(e)

β 0 0 0 1 1 1 0 0 0 1 1 1 0 0 1 1 0 0 1 1 M 0 M 1 M 2 0 x a x z (f) Ω 3 0 0 0 1 1 1 0 0 0 1 1 1 0 0 1 1 0 0 1 1 M 0 M 1 M 2 0 1
x a x z (g)

T 3 0 0 0 1 1 1 0 0 0 1 1 1 0 0 1 1 0 0 1 1 M 0 M 1 M 2 1 0 x a x z (h) Ω 4

The shapes of η(x)

In this section we describe the geometric properties of the function y = η(x). We compute different solutions of (4) for fixed -and representative-values of (β, θ) in each parameter regime.

The initial conditions are chosen of the form (x 0 , y 0 , z 0 ) = (1, 1, z (k) 0 ), with z (k) 0 = 0. For each value of z (k) 0 , the graphs of y = η k (x) and z = η k (x) correspond to the projection of the solution (x k (t), y k (t), z k (t)) of (4) onto the (x, y) and (x, z)-planes, respectively.

Rather than computing each solution as mere integration from a given initial condition, we obtain each desired orbit (x k (t), y k (t), z k (t)) with high accuracy as an element of a family of solutions of a well-posed boundary value problem -which is solved by continuation in Auto [START_REF] Doedel | AUTO-07p: Continuation and bifurcation software for ordinary differential equations[END_REF]; see also [START_REF] Doedel | Lecture notes on numerical analysis of nonlinear equations[END_REF]. This numerical procedure is explained in Appendix C.5.

Figures 2 and3 show the graphs of a set of selected solutions y = η k (x) and z = η k (x), where the index k > 0 if z (k) 0 > 0, and k < 0 if z (k) 0 < 0. Panels (a1) and (a2) of Figure 2 show the projection of each solution to the (x, y) and (x, z)-plane, respectively, for (β, θ) = (0.5, 0.4) ∈ Ω 1 . Also shown in Figure 2(a2) is the set ϕ -1 (0) = {M 0 , M 1 } (grey lines); compare with Figure 1 Let us now describe in more details the properties of y = η(x) in each scenario. In Figure 2(a1), for (β, θ) ∈ Ω 1 , every solution y = η k (x) with positive initial slope z (k) 0 > 0 is a monotone increasing function. From Figure 2(a2), the corresponding derivative functions z = η k (x) grow unbounded at certain finite values

x ∞ k < ∞. Hence, lim x→x ∞ k η k (x) = ∞.
On the other hand, in Figure 2(a1) every solution y = η k (x) with negative initial slope z (k) 0 < 0 is a monotone decreasing function with bounded negative derivative -and, hence, η k (x) eventually changes sign and takes negative values-; in fact, from Figure 2(a2), the corresponding derivatives z = η k (x) tend to M 0 as x → ∞.

In Figure 2(b1), for (β, θ) ∈ Ω 2 , there exists a threshold initial slope z * 0 > 0 such that if z (k) 0 > z * 0 , the solution y = η k (x) grows unbounded much as in the previous case for (β, θ) ∈ Ω 1 ; see the curve η 1 for instance. However, if z (k) 0 < z * 0 , the derivatives z = η k (x) in panel (b2) start as increasing functions until they cross the curve M 1 near x = x as -hence, attaining a maximum at the intersection point with M 1 -and decrease towards M 0 . As a consequence, the corresponding functions y = η k (x) have an inflection point and remain bounded; furthermore, each solution y = η k (x) tends to a horizontal asymptote y = y ∞ k . On the other hand, the behaviour of solutions for z (k) 0 < 0 for (β, θ) ∈ Ω 2 is qualitatively as in Ω 1 . In Figure 3

(a2), for (β, θ) ∈ Ω 3 , if z (k) 0
> 0, the derivatives z = η k (x) decay monotonously to M 0 and fail to intersect the curve M 2 . On the other hand, if z (k) 0 < 0, the derivatives z = η k (x) may not be monotonic functions -initially some of them are decreasing functions that cross the curve M 1 to start to increase; see for instance η -1 , η -2 and η -3 -, but eventually all of them tend monotonously towards M 0 . Therefore, in Figure 3(a1), every curve y = η k (x) tends to a horizontal asymptote y = y ∞ k ; in particular, η -1 , η -2 and η -3 have an inflection point at certain value x < x as where z = η k (x) crosses the curve M 1 . For (β, θ) ∈ Ω 4 , in Figure 3(b1), virtually the same qualitative behaviour of y = η k (x) is observed. In this last scenario, the only minor difference is that all the solutions y = η k (x) with negative initial slope z 0 k < 0 have an inflection point as η k crosses the curve M 1 ; see also Figure 3(b2).

Finally, it remains to explore the shape of y = η k (x) when parameters (β, θ) cross from one open region to another in Figure 1(e). However, notice that, if (β, θ) ∈ T 1 , the set ϕ -1 (0) is qualitatively as in the case for (β, θ) ∈ Ω 1 ; one can expect y = η(x) and z = η (x) to be qualitatively as in Figure 2(a1) and (a2), respectively. Similarly, if (β, θ) ∈ T 3 , the set ϕ -1 (0) is qualitatively as in the case for (β, θ) ∈ Ω 4 for x ≥ 1; hence, the functions y = η(x) and z = η (x) are qualitatively as

η 1 η 2 η 3 η 4 η 5 η 6 η -1 η -2 ✄ ✄ ✄ ✎ η -3 ✄ ✄ ✄ ✎ η -4 η -5 x y (a1) Ω 1 5
η ′ 1 η ′ 2 η ′ 3 η ′ 4 η ′ 5 η ′ 6 M 1 M 0 ◗ ◗ ◗ ❦ x z (a2)
η ′ 1 η ′ 2 η ′ 3 ✠ η ′ 4 ✑ ✑ ✰ η ′ 5 ✘ ✘ ✾ η ′ 6 M 1 M 0 ✁ ✁ ✁ ☛ M 2 x z (b2) Ω 2 Figure 2:
The graph of the functions y = η k (x) and z = η k (x) -for different initial conditions-for (β, θ) = (0.5, 0.4) ∈ Ω 1 in panels (a1) and (a2), respectively, and for (β, θ) ∈ Ω 2 in panels (b1) and (b2), respectively.

in Figure 3(b1) and (b2), respectively. Therefore, it only remains to see the case when (β, θ) ∈ T 2 .

Figure 4 shows the graphs of a set of solutions y = η k (x) and

z = η k (x) for (β, θ) = (2, 0.7) ∈ T 2 . In panel (b), all the functions z = η k (x) are monotonously decreasing (if k > 0) or increasing (if k < 0) towards M 0 as x → ∞. Thus, in panel (a), the functions y = η k (x) increase (if k > 0)
or decrease (if k < 0) monotonously, and converge to a finite horizontal asymptote as x → ∞. In Figures 2,3 and4 the evolutions of transformations η(x) and their derivatives z = η(x) are plotted. They, as a set form the statistical model "per se" in the sense of [START_REF] Mccullagh | What is a statistical model? (with discussion)[END_REF]. Here we consider that to each transformation η a specific estimator and inference function (T-score) are given (see also [START_REF] Jordanova | Weak properties and robustness of t-Hill estimators[END_REF]).

η ′ -1 η ′ -2 ❍ ❍ ( η ′ -3 ✐ η ′ -4 ❳ ❳ ❳ ② η ′ -5 ✠ M 0 ✁ ✁ ☛ M 1 M 2 x z (a2)
η ′ -1 η ′ -2 ❍ ❍ ( η ′ -3 ❍ ❍ ❍ ❍ ( η ′ -4 ✐ η ′ -5 ❳ ❳ ❳ ❳ ❳ ② M 0 M 1 M 2 x z (b2) Ω 4

Application 1: Guanaco Glacier

A glacier is a volume of snow and ice higher than 0.01 km 2 with evidences of flow, that is fed by solid water (snow, hail or hoarfrost), that is transformed into ice by densification and compaction, loosing mass mainly by melting, sublimation or calving. The balance between the gains and losses during a period of time (normally the hydrological year between April 1 and March 31 of the following year) is defined as Glacier Mass Balance [START_REF] Kaser | A manual for monitoring the mass balance of mountain glaciers[END_REF]. For measuring this glacier mass balance, a number of stakes or poles are installed over a glacier (location determined mainly by the altimetry of the zone), where their protruding height is repeatedly measured during each season together with the density of the intervening mass. The changes are converted into water equivalent by multiplying the poles heights by the density of the snow or ice gained or losses during the surveyed period. The mass balance per stake is then aggregated relative to the area of influence of each one, as in the following equation [START_REF] Kaser | A manual for monitoring the mass balance of mountain glaciers[END_REF]

η ′ -1 η ′ -2 ✐ η ′ -3 ✐ η ′ -4 ❳ ❳ ❳ ❳ ❳ ② η ′ -5 ✛ M 0 η ′ 1 η ′ 2 ✟ ✟ ✟ ✙ η ′ 3 ✟ ✟ ✟ ✙ η ′ 4 ✟ ✟ ✟ ✙ η ′ 5 ✟ ✟ ✟ ✙ η ′ 6 ✘ ✘ ✘ ✾ x z (b) T 2
B = i b i s i S .
Guanaco glacier is located in the semi-arid region of Chile (latitude 29S) at an altitude range between 5,000 and 5,300 m above sea level, with a surface area of 1:61 [km 2 ] in 2015 and a maximum thickness of 120[m] [START_REF] Rabatel | Glacier changes in the pascua-lama region, chilean andes (29s): recent mass balance and 50yr surface area variations[END_REF]. This glacier, together with other minor glaciers and glacierets (Figure 5), is located close by Pascua Lama [START_REF] Gascoin | Glacier contribution to streamflow in two headwaters of the huasco river, dry andes of chile[END_REF], a gold mining project that was interrupted in 2015. This mining project was monitoring these glaciers for environmental purposes since 2002 with the main aim of assessing possible impacts of their mining activities on the glacier mass balance among several other variables.

Estimators of the EVI

It is of interest to consider the statistical behavior of maxima X n:n = max{X 1 , ..., X n }, where the sequence of independent random variables X 1 , ..., X n has the common cumulative distribution function F , that is X i ∼ F . The Fisher-Tippett-Gnedenko theorem [START_REF] De Haan | Extreme Value Theory: An Introduction[END_REF] shows us that if there exists a sequence of constants {a n > 0} and b n such that

P X n:n -b n a n ≤ z → G(z) as n → ∞,
where G(z) is a non-degenerate distribution function, then G corresponds to a Generalized Extreme Value distribution (GEV)

G(z) = exp -1 + ξ z -µ σ -1/ξ . ( 13
)
This result quite resembles the Central Limit Theorem, but now with the maximum of the sequence of random variables as the statistic of interest, not the average. The proof of this theorem can be found in [START_REF] Leadbetter | Extremes and related properties of random sequences and processes[END_REF], while a summary of some results with weaker hypothesis (such as loss of independency or stationarity) can be found in [START_REF] Coles | An introduction to statistical modeling of extreme values[END_REF]. This limiting cumulative distribution function (c.d.f.) G(z) can be classified into three types according to its shape parameter ξ (called also EVI), these are often called Weibull (ξ < 0), Fréchet (ξ > 0) and Gumbell (ξ = 0) types. If, for example, F is such that G is a Weibull type, then it is said that F belongs to the Weibull Domain of Attraction of G, usually denoted as F ∈ DA(G α ). The shapes of these distributions for different parameters are given on Figure 6 . The practical relevant range of the EVI is [-1 2 , 1 2 ]. In particular, several practical applications can be found in the case of the Weibull domain of attraction, for example, the ultimate world record in a specific athletic event given today's state of the art [START_REF] Einmahl | Records in Athletics Through Extreme-Value Theory[END_REF], the estimation of the efficiency frontier in economics [START_REF] Farrell | The measurement of productive efficiency[END_REF] or the limit behavior of the distance of two random points over a convex set [START_REF] Mayer | Limit theorems for the diameter of a random sample in the unit ball[END_REF]. In general, the estimation of the finite right endpoint is linked to estimate ξ < 0.

If the previous theorem holds, then the conditional probability for large values of u, has a Generalized Pareto distribution (GP)

P (X > u + y | X > u) ,
F (y) = 1 -1 + ξ y σ -1/ξ , (14)
where σ = σ + ξ(uµ), this limit c.d.f. can be also classified into three types, often called, Beta (ξ < 0), Pareto (ξ > 0) and Exponential (ξ = 0). This approach is known as Peaks over threshold (POT) [START_REF] Coles | An introduction to statistical modeling of extreme values[END_REF] and suits better in a case when the whole data set is available and not only the maximum value of a predefined block. However it is necessary to choose an appropriate threshold. Several estimators were took into account for this work. Here X i:n denotes the i-th order statistic of the sample X 1 , ..., X n , where X n:n = max {X 1 , ..., X n } These estimators depend on the number of order statistics (k) which are used:

• An estimator based on the method of moments [START_REF] Hosking | Parameter and quantile estimation for the generalized pareto distribution[END_REF], labeled as MOM:

ξMOM (k) = 1 2 1 - µ 2 0 (k) σ 2 (k)
,

where:

µ 0 (k) = 1 k k i=1 [X n-i+1:n -X n-k:n ] , σ 2 (k) = 1 k -1 k i=1 [X n-i+1:n -X n-k:n -µ 0 (k)] 2 .
• The method of probability weighted moments estimator [START_REF] Hosking | Parameter and quantile estimation for the generalized pareto distribution[END_REF], labeled as PWM, enjoys good properties in the small sample context [START_REF] Furrer | Probability weighted moments properties for small samples[END_REF]:

ξPWM (k) = 2 - µ 0 (k) µ 0 (k) -2µ 1 (k) ,
where µ 0 (k) is defined as in the previous estimator, and:

µ 1 (k) = 1 k(k -1) k i=1 (i -1) [X n-i+1:n -X n-k:n ] .
• Pickands' estimator [START_REF] Pickands | Statistical inference using extreme order statistics[END_REF], labeled as PICK, is computed because it only needs four distinct order statistics (which makes it a robust estimator), and is location and scale invariant:

ξPICK (k) = 1 log 2 log X n-r k +1:n -X n-2r k +1:n X n-2r k +1:n -X n-4r k +1:n ,
where r k = k/4 with • the floor function.

• A moment estimator [START_REF] Arnold | On the estimation of the extreme-value index and large quantile estimation[END_REF], labeled as m-Hill:

ξm-Hill = M (1) n + 1 - 1 2   1 - M (1) n 2 M (2) n    -1 , M (j) n = 1 k k-1 i=0 (log X n-i:n -log X n-k:n ) j .
• Falk's estimator [START_REF] Falk | Some best parameter estimates for distributions with finite endpoint[END_REF], labeled as n-Hill, defined as:

ξn-Hill = 1 k k-1 i=1 [log (X n:n -X n-i:n ) -log (X n:n -X n-k:n )] , k ∈ {1, ..., n -1}.

Negative t-Hill estimator

The t-Hill estimator [START_REF] Fabián | On robust and distribution sensitive hill like method[END_REF] is a robust EVI estimator which works in the Fréchet domain of attraction, defined as:

ξ = -1 + 1 k k i=1 X k+1:n X i:n -1
, where: k ∈ {1, ..., n -1}.

This one is a particular case of the HME [START_REF] Beran | The harmonic moment tail index estimator: asymptotic distribution and robustness[END_REF], which relates to the studied score function (1):

ξ = H (β) n:k :=              1 β -1   1 k k i=1 U 1-β ik -1 -1   , if β > 1 1 k k i=1 log (U ik ) , if β = 1
where

U ik := X n-i+1:n X n-k:n .
In order to obtain a negative version of t-Hill estimator we applied a procedure similar to the one introduced by Falk [START_REF] Falk | Some best parameter estimates for distributions with finite endpoint[END_REF][START_REF] De Haan | Extreme Value Theory: An Introduction[END_REF], resulting in the following estimator for the EVI, labeled as n-t-Hill:

ξn-t-Hill = 1 - 1 k k i=1 X n:n -X n-i+1:n X n:n -X n-k:n -1
, where: k ∈ {2, ..., n -1}.

In more detail, we have been considering the fact that if X belongs to the Weibull domain of attraction with ξ < 0, therefore, Y defined as:

Y = 1 x * -X (15)
belongs to the Fréchet domain of attraction, i.e., Y ∈ DA(G -ξ ) (because -ξ > 0). In [START_REF] Arnold | On the estimation of the extreme-value index and large quantile estimation[END_REF], x * denotes the unknown value of the finite right endpoint, as an estimator of this value, the maximum observation of the sample was used, which is recommended for ξ < -1 2 [START_REF] De Haan | Extreme Value Theory: An Introduction[END_REF]. The following asymptotic results are detailed in Section Appendix A and proved in the Appendix:

ξn-t-Hill P -----→ k→∞ ξ, √ k ξn-t-Hill -ξ d -----→ k→∞ N 0, ξ 2 (1 -ξ) 2 (1 -2ξ) , ξ < - 1 2 , k -ξ ξn-t-Hill -ξ d -----→ k→∞ ξ(1 -ξ)Weibull 1, - 1 ξ , 0 > ξ > - 1 2 .

Computational Results

For this section, the chosen software was R [START_REF] Core | R: A Language and Environment for Statistical Computing[END_REF], the code with the implementation of the estimation of a partially smooth c.d.f. was cordially shared by S. Müller and implemented in this work. All previous estimators were implemented plus their smoothed versions [START_REF] Müller | Partially smooth tail-index estimation for small samples[END_REF], labeled with an "s" at the end of their names.

Small sample: Generalized Pareto distribution

In order to better contrast the results of [START_REF] Müller | Partially smooth tail-index estimation for small samples[END_REF], the same setup was studied, that is, n random samples were generated from a Generalized Pareto (GP) distribution of parameters ξ = {-1, -0.75, -0.5, -0.25, -0.1}, µ = 0, σ = 1 with n = {16, 32, 64} and a fixed seed for the random number generator (RNG) algorithm of value 200905 the same that was used in [START_REF] Müller | Partially smooth tail-index estimation for small samples[END_REF]. The results are summarized in It can be seen a similar behavior among different values of n, but there are more differences among them than in [START_REF] Müller | Partially smooth tail-index estimation for small samples[END_REF], therefore, results for n = {16, 32, 64} are shown and not only for n = 32. Smoothed estimators show great advantage over their non-smoothed counterparts, however, the lack of theoretical results about their order of convergence or their behavior over more general distributions, makes difficult to put smoothed estimators over non-smoothed ones in a more general view.

In the large sample case, the estimation is way more reliable than in the small one. One can observe here that all estimators are converging to the original value of the simulation, however, it can also be seen, the slower rate of convergence of both Negative t-Hill (n-t-Hill) and Falk's estimator (n-Hill) when ξ < -0.5, this result is expected given the asymptotic results obtained in this work.

Robustness testing for Pareto tails

The robust properties of t-Hill estimator have already been studied [START_REF] Beran | The harmonic moment tail index estimator: asymptotic distribution and robustness[END_REF] they are explained by the fact that this estimator uses the Harmonic mean instead of the Arithmetic one on the data. For this experiment, the c.d.f. of the contaminated samples is defined as:

F (x) = (1 -α)F 1 (x) + αF 2 (x),
where F 1 is the c.d.f. of a transformed Pareto distribution (see equation ( 15)) of parameter -ξ with ξ < 0, while F 2 is similar but with parameter -ξ/2.

Table 2 shows the relative efficiency of the estimators among themselves, for each value of α = {0, 0.05, 0.1, 0.15, 0.2} and ξ = {-0.1, -0.25, -0.5, -0.75}. The MSE and the percentage of contribution of each estimator to the global MSE is computed. The m-Hill estimator was left out due to its large MSE results. A starting value of k = 10 was chosen for the computation of the percentage, given the large instability of the estimation when almost all order statistics are took into account. It can be stated from the previous table that both n-Hill and n-t-Hill estimators are the only ones to perform better while the contamination level grows when ξ > -0.75. Also, it can be seen how all smoothed estimators are performing better than their non-smoothed counterpart when ξ > -0.5.

In table 3 it can be found a similar analysis, but now taking the average of 20 order statistics around the optimal value of k for each estimator, the optimal was chosen given its MSE. A notion of relative stability can be obtained from these results if we make a comparison with table 2, because a decrease in the percentage across the tables indicates a larger zone of reliable estimation of the EVI, this topic is quite important, because the optimal value of k is a topic of active research.

After a contamination both Negative t-Hill (n-t-Hill) and Falk's estimator (n-Hill) are almost invariant in their estimation, which clearly shows the robust properties of these estimators in this case of contamination. Estimation of the EVI on mass balance data of several glaciers from Pascualama region, from 2002 to 2014 can be found in Figure 7. We can see from the values of several EVI estimators that for a complete modelling of the EVI on mass balance data from Pascualama region we need to consider several values of EVI, not only one.

Application 2: Extremal snow loads

Understanding of snow extremes plays an important role for both climatology and civil engineering. An approach combining engineering and climatology to assess accidental snow loads on structures is suggested in [START_REF] Sadovský | Assessment of accidental snow loads for the design of structures[END_REF]. We use the data of collection and analysis of snow loads in Slovakia carried out recently [START_REF] Sadovský | Collection and analysis of climatic measurements for the assessment of snow loads on structures[END_REF] and [START_REF] Sadovský | Special features of the collection and analysis of snow loads[END_REF]. The long-term weekly measurements of snow water equivalent (SWE) of snow cover at rain-gauge stations are employed. Out of the rain-gauge stations, meteorological stations were selected at which daily SWE values have been recalculated using other climatological measurements, like depth of snow cover etc. (see [START_REF] Sadovský | Special features of the collection and analysis of snow loads[END_REF]). Preliminary statistical analysis has been made in [START_REF] Sadovský | Assessment of accidental snow loads for the design of structures[END_REF] and in [START_REF] Stehlík | Statistical Analysis related to Exceptional Snow Loads[END_REF].

The SWE records of winter seasons are well suited for the assessment of the characteristic snow load on the ground, which is defined as 98 % quantile of a suitable extreme value distribution fitted to the yearly snow load maxima. It is assumed that the maximum is a member of the same population, however, with a mean return period of say about 1000 years and more. Following [START_REF] Sanpaolesi | Phase 1 Final Report to the European Commission, Scientific Support Activity in the Field of Structural Stability of Civil Engineering Works: Snow Loads[END_REF] the largest snow load value is exceptional if the ratio k of the load to the characteristic snow load determined without that value is greater than 1.5. The snow loads identified as exceptional should be treated in accidental design situations as accidental actions (loads), cf. Eurocode EN 1990 (see [START_REF] Sanpaolesi | Phase 1 Final Report to the European Commission, Scientific Support Activity in the Field of Structural Stability of Civil Engineering Works: Snow Loads[END_REF]).

The novelty of the approach for the assessment of the accidental snow loads by [START_REF] Sadovský | Assessment of accidental snow loads for the design of structures[END_REF] can be briefly described. First the k values in excess of 1.5 are identified. Then by the expertise of climatologists based on the geomorphology of Slovakia, regions of similar climate conditions for the occurrence of accidental snow loads are determined (see Figure 8). Within a given region, the values of the empirical distribution function F restricted to the N ordered k values in excess of 1.5 is calculated as where N R is the sum of winter seasons over all stations in the region and i ∈ (N R -N + 1, ..., N R ).

F (k i ) = i N R + 1 , (16)
The obtained empirical upper tail for k ratios is approximated, e.g. by nonlinear regression analysis using Pareto, exponential and Generalized Extreme Value (GEV) distributions. The extremes of 0.999 and 0.9999 quantiles of the distributions are of particular interest. Under Generalised Pareto Distribution (GPD) with parameters µ, ξ and δ we understand the distribution with c.d.f.

1 -1 + ξ x -µ δ -1/ξ , x > µ - δ ξ . (17)
Under Generalized Extreme Value (GEV) distribution with parameters ξ, µ and δ we understand the one with c.d.f.

exp -1 + ξ x -µ δ -1/ξ , x > µ - δ ξ . ( 18
)
We fit GPD to data from the companion paper (Sadovský et al. [START_REF] Sadovský | Assessment of accidental snow loads for the design of structures[END_REF]). Particularly, we study regions of Slovakia separately. We refine study of Regions 2 and the composite Region 4, within which the mountain basins are considered as one region. For the exceptional snow loads and their corresponding k values in Region 1, treated in Sadovský et al. [START_REF] Sadovský | Exceptional snowfalls and the assessment of accidental snow loads for structural design[END_REF], a statistical dependence on the altitude is studied. The idea is to check the anticipated low dependence of k values on the altitude inferred from their definition, which comprises the altitude dependence already in the characteristic values.

Region 1

Here we consider 63 data, maximal values of the exceptional snow loads in Region 1, described in [START_REF] Stehlík | Statistical Analysis related to Exceptional Snow Loads[END_REF], where a subset of these data is considered. We assume that the observations are independent. Their mean excess plot shows that the observed random variable has c.d.f. with Pareto tail. In order to fit the GPD to the exceedances of the threshold we use maximum likelihood approach, implemented in function gpd.f it in R. The threshold that we choose is µ = 97. The number of exceedances is 28. The estimated values of the parameter ξ and δ of the c.d.f. ( 17) are given in Table 4. The GPD pp-plot of the exceedances is given on Figure 9. It proves that we have made a good fit of the distribution of the data over the threshold.

In the above estimation we used only the upper order statistics of the observations in order to estimate the parameters. In this way we lose information about the observed values below the threshold 97. If we can estimate the c.d.f. in all range of the data it always would be better. Therefore now we make GEV fit of the data. The maximum likelihood approach is implemented in function gev.f it in R. The estimated values of the parameters µ, ξ and δ of the c.d.f. ( 18) are given in Table 5. Looking on the pp-plot on Figure 10 we can conclude that this GEV [START_REF] Fabián | On robust and distribution sensitive hill like method[END_REF] fit with the parameters, given in Table 5

Fitting the distribution of k in Region 1

Here we consider component k -the ratio of the load to the characteristic snow load determined without the largest snow load value for Region 1. The number of observations is again 63. The mean excess plot shows that the observed random variable has GPD with negative parameter ξ. In order to fit this distribution we use the maximum likelihood approach, implemented in function gpd.f it in R. For the threshold we use the minimal observed value µ = 1.5. The number of exceedances is 39. The estimated values of other parameters ξ and δ of the c.d.f. ( 18) are given in Table 6.

The pp-plot of the exceedances on Figure 13 shows that our GPD fit of the distribution of the data is very good. The same conclusion could be made having in mind the similarity between the theoretical and empirical c.d.fs.

We can compare the Negative Hill estimator, proposed by [START_REF] Falk | Some best parameter estimates for distributions with finite endpoint[END_REF] and the Negative t-Hill introduced in Section 3.3. Both, the Negative Hill like plot on Figure 14 and the Negative t-Hill like plot on Figure 15 confirmed the above values of the EVI.

Region 4

In contrast to the Region 1, where we have relatively much amount of data, in Region 4 we have only 17 observations. Therefore instead of using the functions gev.git and gpd.f it in R, we use the regression in order to estimate the coefficients in these distributions. This approach is considered in [START_REF] Stehlík | Statistical Analysis related to Exceptional Snow Loads[END_REF] with respect to the distribution of k, where these data are considered. In order to simplify the regression model, instead of ( 17) we use the following definition for Pareto distribution Under exponential c.d.f. we understand

F (x) = 1 - λ x 1/ξ , x > λ. (19)
F (x) = 1 -e -x-µ a , x > µ. (20)
The mean excess plot of the maximal values of the exceptional snow loads in Region 4 shows that the observed random variable has c.d.f. with Pareto tail. Here we compare three types of distributions for modelling these data: exponential type (20), Pareto type [START_REF] Falk | Some best parameter estimates for distributions with finite endpoint[END_REF] and Generalized extreme value type [START_REF] Fabián | On robust and distribution sensitive hill like method[END_REF]. For the estimation of the coefficients we used lm function in R and the results are compared with the corresponding Hill and t-Hill estimators.

We start with fitting of c.d.f. [START_REF] Falk | Some best parameter estimates for distributions with finite endpoint[END_REF]. Due to the fact that we have relatively small amount of data we will use for the threshold the minimal observed value. Having in mind that ( 20) is equivalent to

log(1 -F (x)) = log λ ξ - log x ξ , x > λ
we make the following transformation of the data F 1 = log(1cdf ) and X1 = log sort(M ), where cdf are the values in the empirical c.d.f. of the data and M denoted the Maximal snow load. The scatter plot of F 1 and X1 is given on Figure 16. Therefore we use the regression model

F 1 = log λ ξ -ξ -1 X1 + ε, (21)
where ε is the random error. The estimators of the coefficients are given in Table 7. Both are statistically significant. Residual standard error of the model in 0.1841 on 14 degrees of freedom.

We have the multiple R 2 = 0.9492 and adjusted R 2 = 0.9456. The empirical value of the Fisher characteristic is 261.7. Its p-value is 1.866e-10. The degrees of freedom are correspondingly 1 and 14. The last means that this model is adequate. The normal qq-plot of the error terms on Figure 17 shows that the distribution of the error term in the last regression is not a normal one. The GEV qq-plot on Figure 18 shows that the observed random variable is GEV distributed. The estimators of the coefficients in [START_REF] Fabián | On robust and distribution sensitive hill like method[END_REF] are given in Table 8. Now we come back to the initial coefficients and obtain that ξ = 0.3070782 and λ = 113.4926. The qq-plot on Figure 19 show the quality of this fit. The confidence intervals are wide therefore we also apply GPD and GEV models. We can use Hill and t-Hill estimators in order to obtain the estimated value of ξ. Depending on the number of upper order statistics that are included in the estimators, the Hill and t-Hill plots are given correspondingly in Figure 20 and Figure 21 Table 9: Estimators of the coefficients in GPD of the error term in [START_REF] Embrechts | Modelling Extremal Events for Insurance and Finance[END_REF].

Now we will use the fact that the Pareto distribution appears mainly in exceedances over high threshold. We use the function gpd.f it in R, over the threshold 140 and obtain the parameter µ in ( 17) is equal to 140, and the estimators of ξ and δ are given in Table 9. The qq-plot on Figure 22 shows that this fit could be useful. Now we use the same approach including all the values. Although the sample size is only 17, we use the function gpd.f it in R, over the threshold 105.99 and we obtain μ = 105.99 in [START_REF] Embrechts | Modelling Extremal Events for Insurance and Finance[END_REF], ξ and δ are given in Table 10. In this case we have very wide confidence intervals. The qq-plot is given on Figure 23.

The differences in Pareto predictions show again the well known fact that if we use a small amount of data and the distribution of the observed variable is not exact Pareto, the Pareto fit is not stable. Therefore we make GEV [START_REF] Fabián | On robust and distribution sensitive hill like method[END_REF] fit of the data. The estimated values of the parameters µ, ξ and δ of the c.d.f. are given in Table 11. ξ still have wide confidence intervals, however µ and δ are statistically significant. The pp-plot on Figure 24 shows that this fit is relatively good.

We make also exponential fit with c.d.f. [START_REF] Farrell | The measurement of productive efficiency[END_REF] [2] this is the estimator for µ. Again by [START_REF] Balakrishnan | The Exponential Distribution: Theory, Methods, and Applications[END_REF], the best estimator for scale a is a = 17(X n -106) 16 = 17(165.9 -106) 16 = 63.64375.

The corresponding exponential pp-plot, given on Figure 25 shows that again we observe a good fit.

Region 2

In Region 2 we have only 9 observations. Here we proceed with fitting the distributions of the maximal snow loads and k in this region.

The maximal values of the exceptional snow loads in Region 2 are given on Figure 26. In sense of reproducible research we present the data m 144 86 88 110 84 101 80 89 97

Although we have very small amount of data their mean excess plot on Figure 27, shows that the observed random variable is possibly to have a c.d.f. with Pareto tail. In order to fit the Generalised Pareto Distribution (GPD) to the exceedances of the threshold we use maximum likelihood approach, implemented in function gpd.f it in R. The threshold that we choose is the minimal observation that is 85. The number of exceedances is 7. The estimated value of the parameter ξ and δ of the c.d.f.

1 -1 + ξ x -85 δ -1/ξ , x > 85 - δ ξ (22)
are given in Table 12. Due to the small number of observations we have wide confidence intervals. The pp-plot of the exceedances is given on Figure 28 the corresponding c.d.f. are in Figure 29. It proves that we have made a relatively good fit of the distribution of the data over the threshold.

As we have already mentioned if we can estimate the c.d.f. in all range of the data it always would be better. Therefore now we make GEV fit of the data. More precisely we estimate the parameters of the c.d.f. [START_REF] Fabián | On robust and distribution sensitive hill like method[END_REF]. We use the function gev.f it in R and obtain the estimated values of the parameters µ, ξ and δ of the c.d.f. [START_REF] Fabián | On robust and distribution sensitive hill like method[END_REF]. They are given in Table 13. Looking on Figure 31 we can compare the theoretical GEV c.d.f [START_REF] Fabián | On robust and distribution sensitive hill like method[END_REF] with the parameters, given in Table 13, with the empirical c.d.f. of the observed data. This together with the pp-plot on Figure 30 show that this fit is better than the GPD fit.

Fitting the distribution of k in Region 2

Here we consider component k -the ratio of the load to the characteristic snow load determined without the largest snow load value for Region 2. In sense of reproducible research we present the data k 1.72 1.60 1. [START_REF] Zhou | Existence and consistency of the maximum likelihood estimator for the extreme value index[END_REF] The mean excess plot of the data shows that the observed random variable has GEV distribution therefore we estimate the parameters of the c.d.f. [START_REF] Fabián | On robust and distribution sensitive hill like method[END_REF]. The estimated values of the parameters µ, ξ and δ of the c.d.f. ( 18) are given in Table 14. Looking on Figure 33 we can compare the theoretical GEV c.d.f [START_REF] Fabián | On robust and distribution sensitive hill like method[END_REF] with the parameters, given in Table 14, with the empirical c.d.f. of the observed data. This together with the pp-plot on Figure 32 shows that although the confidence intervals of the coefficients are wide this fit is good.

Conclusions

It is clear that oscillations of natural ecological systems are measured imprecisely. To this imprecision contribute both non-chaotical and chaotical dynamical systems. In the present paper we explained non-chaotical dynamics of autonomous system of t-scores, which underline statistical estimates of entropy. Other source of contribution is switching between negative and positive extreme value indexes. This is well visible from both applications, Guanaco glacier in Chile and extreme snow loads in Slovakia. These observations provide new illustrations of the decomposition in deterministic, stochastic and chaotic parts introduced in [START_REF] Stehlík | Missing chaos in global climate change data interpreting?[END_REF]. Therein we studied the methane emission example and we outlined the problem of deceptiveness regarding a complete certainty of the climate change. Several other contributions to imprecision can be made by fusion of several p-values ( [START_REF] Baker | Statisticians issue warning over misuse of P values[END_REF]).

More care should be taken in interpretations and the use of various observations of the climate change as well as in their methodologies and approaches.

Thus based on these new observations, we doubt the possibilities of a 100% judgment (a 100% paradigm of expectations) that the climate change has been caused solely by human effect. The currently obtained consensus of 95-98% to justify a human impact on the global climate change is probably the maximal possible threshold. effect. The currently reached consensus of 95-98% to justify a human impact on the global climate change is probably the maximal attainable threshold. where ζ is a standard Gaussian random variable, see [START_REF] De Haan | Extreme Value Theory: An Introduction[END_REF], Theorem 3.6.4. It is thus possible to compare the asymptotic behavior of ξ and ξ for ξ ∈ (-1, -1/2). First, it is clear that ξ has a larger asymptotic variance (given by ξ 2 (1ξ) 2 /(1 -2ξ)) than ξ (which is ξ 2 ). It is possible to show that ξ has a smaller asymptotic bias (given by (1ξ)/(1ρξ)) than ξ (which is ξ/(ρ(1ρ)(1 + ξ))) for all ξ ≤ ξ 0 where ξ 0 =

F θ (x) = 1 -exp(-x θ ), x > 0.
1-ρ- √ (1-ρ)(1-5ρ+4ρ 2 -4ρ 3 ) 2(ρ 2 -ρ-1)
. In particular, the asymptotic bias of ξ explodes when ξ approaches -1 while the asymptotic bias of ξ remains bounded. Finally, note that, in case (ii) even when λ = 0, ξ has a negative asymptotic bias given by (1ξ)ξΓ(1ξ)k ξ . Nevertheless, this bias can easily be estimated and corrected.

Lemma 2. Under (A.1), the n-t-Hill estimator can be rewritten as

1 1 - ξ - 1 1 -ξ = T 1,n -1 1-ξ 1 -T 2,n + ξ 1 -ξ T 2,n 1 -T 2,n , with T 2,n = U (Y n-k,n ) U (Y n,n ) , T 1,n = 1 k k-1 j=0 U (Y n-k,n ) U (Y n-j,n ) ,
and where Y 1,n ≤ • • • ≤ Y n,n are ordered statistics associated with n independent standard Pareto random variables.

The next lemma provides some consequences of the second condition (A.2).

Lemma 3. Under (A.2), there exists a function A 0 asymptotically equivalent to A such that, for all ε > 0 there exists t 0 > 0 such that, for all t ≥ t 0 and x ≥ 1,

x ξ U (tx) U (t) = 1 + R 1 (t, x)A 0 (t), (C.3) We first establish the asymptotic distribution of T 1,n . in view of [START_REF] De Haan | Extreme Value Theory: An Introduction[END_REF], page 75. Second, using the same arguments as previously yields

U (t) U (tx) = x ξ - x ξ+ρ -x ξ ρ A 0 (t) -x ξ R 0 (t,
1 k k-1 j=0 1 ρ U (Y n-j:n ) U (Y n-k:n ) ξ+ρ - U (Y n-j:n ) U (Y n-k:n ) ξ d = 1 k k-1 j=0 1 ρ (Y * j ) ξ+ρ -(Y * j ) ξ .
The laws of large numbers shows that this quantity converges in probability to τ (ξ, ρ) and therefore

2. 1 .

 1 The qualitative behaviour of the solutions 2.1.1. The function y = η(x) is monotone

(a). 2 .

 2 If (β, θ) ∈ T 1 , the set ϕ -1 (0) = {M 0 , M 1 } consists of two branches, where M 1 is the graph of z = θ + 1 θx . The set ϕ -1 (0) is qualitatively as in the sketch of Figure 1(b).

Figure 1 :

 1 Figure 1: Sketches of the set ϕ -1 (0) in the (x, z)-plane (in panels (a)-(d) and (f)-(h)) for different parameter values (β, θ) in the regions shown in panel (e).

  (a). Similarly, Figures 2(b1)-(b2), Figures 3(a1)-(a2) and Figures 3(b1)-(b2) show the corresponding graphs of y = η k (x) and z = η k (x) (and the set ϕ -1 (0)) for (β, θ) = (0.5, 0.7) ∈ Ω 2 , (β, θ) = (1.5, 0.7) ∈ Ω 3 and (β, θ) = (3, 0.7) ∈ Ω 4 , respectively.

Figure 3 :

 3 Figure 3: The graph of the functions y = η k (x) and z = η k (x) -for different initial conditions-for (β, θ) = (1.5, 0.7) ∈ Ω 3 in panels (a1) and (a2), respectively, and for (β, θ) = (0.5, 0.7) = (3, 0.7) ∈ Ω 4 in panels (b1) and (b2), respectively.

Figure 4 :

 4 Figure 4: The graph of the functions y = η k (x) and z = η k (x) -for different initial conditions-for (β, θ) = (2, 0.7) ∈ T 2 in panels (a) and (b), respectively.

Figure 5 :

 5 Figure 5: Geographical location of glaciers nearby Guanaco glacier, the three largest are Guanaco, Estrecho and Ortigas 1; accounting for almost the 95% of the area of the 7 studied ice bodies.

Figure 6 :

 6 Figure 6: Different behavior of the Generalized Extreme Value distribution according to the value of its shape parameter, the EVI.

Figure 7 :

 7 Figure 7: EVI computed over the mass balance data of several glaciers from Pascualama region, from 2002 to 2014.

Figure 8 :

 8 Figure 8: Regions of exceptional snow loads (Sadovský et al. [36]). Parameter Estimation Standard error ξ 0.3325008 0.2979518 δ 36.8023369 12.7791203

Figure 9 :Figure 10 :

 910 Figure 9: GPD pp-plot of the maximal values of the exceptional snow loads in Region 1.

Figure 11 :

 11 Figure 11: Hill plot of the Maximal values of the exceptional snow loads in Region 1.

Figure 12 :Table 6 :

 126 Figure 12: t-Hill plot of the Maximal values of the exceptional snow loads in Region 1.

Figure 13 :

 13 Figure 13: GPD pp-plot of k in Region 1.

Figure 14 :

 14 Figure 14: Hill plot of k in Region 1.

Figure 15 :

 15 Figure 15: t-Hill plot of k in Region 1.

  .

F1Figure 16 :

 16 Figure 16: Maximal values of the exceptional snow loads in Region 4.

Figure 17 :

 17 Figure 17: Normal qq-plot of the error terms in model (21).

Figure 18 :

 18 Figure 18: GEV qq-plot of the error terms in model (21).

Figure 19 :

 19 Figure 19: Exponential (see (19)) qq-plot of the maximal snow loads in Region 4.

Figure 20 :

 20 Figure 20: Hill plot of the Maximal values of the exceptional snow loads in Region 4.

Figure 21 :

 21 Figure 21: t-Hill plot of the Maximal values of the exceptional snow loads in Region 4.

Figure 22 :

 22 Figure 22: GPD pp-plot of the maximal values of the exceptional snow loads in Region 4, threshold 140.

Figure 23 :

 23 Figure 23: GPD pp-plot of the maximal values of the exceptional snow loads in Region 4, threshold 105.99.

Figure 24 :

 24 Figure 24: GEV pp-plot of the maximal values of the exceptional snow loads in Region 4.

Figure 25 :

 25 Figure 25: Exonential pp-plot of the maximal values of the exceptional snow loads in Region 4.

Figure 26 :Figure 27 :

 2627 Figure 26: Maximal values of the exceptional snow loads in Region 2.

Figure 28 :

 28 Figure 28: GPD pp-plot of the maximal values of the exceptional snow loads in Region 2.

Figure 29 :

 29 Figure 29: Empirical c.d.f. and the corresponding estimated GPD c.d.f.

Figure 30 :

 30 Figure 30: GEV pp-plot of the maximal values of the exceptional snow loads in Region 2.

Figure 31 :

 31 Figure 31: Empirical c.d.f. and the corresponding estimated GEV c.d.f.

Figure 32 :

 32 Figure 32: GPD pp-plot of the values of k in Region 2.

Figure 33 :

 33 Figure 33: Empirical c.d.f. and the corresponding estimated GEV c.d.f.

Theorem 1 .Corollary 1 .

 11 Suppose (A.2) holds and let θ = min(1/2, -ξ). Let k → ∞ such as k/n → 0 and k θ A(n/k) → λ ∈ R as n → ∞. Then, the following asymptotic expansion holdsk θ ( ξξ) = k θ-1/2 ξ(1ξ) √ 1 -2ξ ζ (1 + o P (1)) + k θ+ξ (1ξ)ξ ζ ξ (1 + o P (1)) + λ(1ξ) 1ρξ (1 + o P (1)),where ζ is a standard Gaussian random variable and ζ ξ is a random variable following the Weibull distribution W (-ξ).As a consequence, two different cases appear: Assume that the assumptions of Theorem 1 hold with λ = 0.(i) If ξ < -1/2, then k 1/2 ( ξξ) is asymptotically Gaussian centred with variance ξ 2 (1ξ) 2 /(1 -2ξ). (ii) If ξ > -1/2, then k -ξ ( ξξ) converges in distribution to (1ξ) 2 W (-1/ξ).As a comparison, let us stress that the asymptotic distribution of the n-Hill estimator defined by ξ n:n -X n-i+1:n )log(X n:n -X n-k:n ) (see for instance[START_REF] De Haan | Extreme Value Theory: An Introduction[END_REF], paragraph 3.6.2) has been established only for ξ ∈ (-1, -1/2) whereas Theorem 1 holds for all ξ < 0. Under the assumptions of Theorem 1 and assuming ξ ∈ (-1, -1/2), it has been shown that k 1/2 ( ξξ) = ξ ζ (1 + o P (1)) + λξ ρ(1ρ)(1 + ξ) (1 + o P (1)),

  x)A 0 (t) + R 2 x)A 2 0 (t), (C.4) with |R 0 (t, x)| ≤ εx ρ+ε , |R 1 (t, x)| ≤ ε -1/ρ, |R 2 (t, x)| ≤ 4(ε 2 + 1/ρ 2 )x ξ .

Lemma 4 .Lemma 5 . 1 ( 1 5 .R 0 R 2 3

 45115023 Let k → ∞ such that k/n → 0 as n → ∞. Then, under (A.2), k -ξ T 1,n d → ζ ξ where ζ ξ is a random variable following the Weibull distribution W (-1/ξ).Second, we focus on the asymptotic distribution of T2,n . Let k → ∞ such that k/n → 0 as n → ∞. Then, under (A.2), + σ(ξ)k -1/2 ζ(1 + o P (1)) + τ (ξ, ρ)A(n/k)(1 + o P (1)), where ζ is a standard Gaussian random variable, τ (ξ, ρ) = -ξ)(1-ρ-ξ) and σ(ξ) Applying expansion (C.4) of Lemma 3 with t = Y n-k:n P → ∞ (see[START_REF] De Haan | Extreme Value Theory: An Introduction[END_REF], Corollary 2.2.2) and x = Y n-j:n /Y n-k:n ≥ 1, j = 0, . . . , k -1 yields the expansionT 2,n = T 3,n -T 4,n -T 5,n + T 6,n (Y n-k:n , Y n-j:n Y n-k:n , T 6,n = A 2 0 (Y n-k:n ) Y n-k:n , Y n-j:n Y n-k:n .The four terms are studied separately. In view of Rényi representation, (Y n-j:n /Y n-k:n ) j=0,...,k-1 d = (Y * k-j:k ) j=0,...,k-1 where Y * 1 , . . . , Y * k is a sample of independent random variables from the standard Pareto distribution and Y * 1:k ≤ • • • ≤ Y * k:k are the associated ordered statistics. Consequently, T 3,n can be rewritten as T Remarking that E((Y * 1 ) ξ ) = 1/(1ξ) and var((Y * 1 ) ξ ) = σ 2 (ξ), the central limit theorem entails that T 3,n = 1 1ξ + k -1/2 σ(ξ)ζ(1 + o P (1)), (C.6) where ζ is a standard Gaussian random variable. Let us now focus on T 4,n . First, remark that A 0 (Y n-k:n ) = A(Y n-k:n )(1 + o P (1)) = A(n/k)(1 + o P (1)),

T 4 ,

 4 n = τ (ξ, ρ)A(n/k)(1 + o P (1)). (C.7)

Table 1 ,

 1 which tries to resemble[30, table 3].

		n=16		n=32		n=64
	ξ	Best estimator(s) Range(s) of k Best estimator(s) Range(s) of k Best estimator(s) Range(s) of k
	-0.1	n-Hills MOMs	k < 4 k ≥ 4	n-Hills MOMs	k < 4 k ≥ 4	n-Hills MOMs	k < 4 k ≥ 4
	-0.25	n-Hills	All k	n-Hills	All k	n-Hills	All k
	-0.5 -0.75	n-Hills n-t-Hills	All k All k	n-Hills n-t-Hills n-t-Hills	k < 8 k ≥ 8 All k	n-Hills n-t-Hills n-t-Hills n-Hill	k < 5 k ≥ 5 k ≤ 15 k > 15
	-1	n-t-Hills n-Hills	k ≤ 10 k > 10	n-t-Hills n-Hill	k < 7 k ≥ 7	n-t-Hills n-Hill	k < 6 k ≥ 6

Table 1 :

 1 Best estimator according to MSE for different sample sizes and EVI. Estimators ending with "s" represent smoothed versions according to[START_REF] Müller | Partially smooth tail-index estimation for small samples[END_REF].

Table 2 :

 2

	ξ	α	n-Hill n-Hills n-t-Hill n-t-Hills MOM MOMs PWM PWMs PICK PICKs
		0	13.62 8.96	13.38	9.56	2.63	1.56	3.62	2.08	31.89 12.69
		0.05 13.55 8.96	13.33	9.53	2.68	1.64	3.68	2.19	31.58 12.87
	-0.1	0.1	13.02 8.58	12.74	9.13	2.74	1.72	3.84	2.37	32.47 13.39
		0.15 12.53 8.38	12.32	8.90	2.82	1.88	4.05	2.64	32.58 13.91
		0.2	12.01 7.94	11.73	8.43	3.04	2.11	4.36	3.00	33.01 14.38
		0	9.50	2.48	9.71	3.30	4.34	2.21	5.91	3.14	43.01 16.39
		0.05 9.11	2.33	9.25	3.13	4.36	2.38	5.97	3.35	43.34 16.79
	-0.25	0.1	8.56	2.15	8.68	2.91	4.57	2.65	6.23	3.70	43.32 17.23
		0.15 7.68	1.92	7.84	2.65	4.71	3.10	6.53	4.30	43.15 18.13
		0.2	6.85	1.68	6.98	2.37	5.04	3.66	7.01	5.02	42.56 18.85
		0	3.74	5.12	4.80	2.94	7.64	4.62	8.91	5.13	44.80 12.30
		0.05 3.27	5.37	4.26	3.26	7.38	4.88	8.78	5.45	44.56 12.78
	-0.5	0.1	2.84	5.42	3.80	3.50	7.47	5.32	8.88	6.01	43.16 13.60
		0.15 2.44	5.30	3.38	3.64	7.58	5.91	9.11	6.72	41.34 14.59
		0.2	2.12	5.18	3.02	3.78	7.92	6.57	9.48	7.51	38.96 15.45
		0	1.31	22.60	2.18	16.12	6.47	8.30	6.50	7.59	23.11 5.81
		0.05 1.19	22.41	2.04	16.21	6.39	8.37	6.50	7.75	22.84 6.30
	-0.75	0.1	1.21	21.75	2.04	15.99	6.47	8.51	6.65	8.01	22.40 6.98
		0.15 1.25	20.54	2.09	15.38	6.74	8.70	7.01	8.34	22.08 7.87
		0.2	1.36	19.17	2.23	14.64	7.17	8.94	7.46	8.73	21.47 8.84

Percentages of relative efficiency (the lesser the better) of each estimator for different levels of contamination and EVI. The minimum value can be seen in bold.

Table 3 :

 3

	ξ	α	n-Hill n-Hills n-t-Hill n-t-Hills MOM MOMs PWM PWMs PICK PICKs
		0	19.35 14.58	17.34	13.58	1.22	1.24	1.85	1.90	16.57 12.37
		0.05 18.41 13.84	16.31	12.74	1.40	1.55	2.23	2.21	17.40 13.91
	-0.1	0.1	16.30 12.23	14.27	11.15	2.06	2.02	3.24	2.94	19.32 16.47
		0.15 14.15 10.67	12.29	9.63	2.53	2.44	4.10	3.62	21.91 18.68
		0.2	12.33 9.25	10.56	8.26	3.05	2.91	4.98	4.34	23.76 20.56
		0	13.95 4.71	11.51	4.35	2.63	2.70	3.95	4.04	30.41 21.76
		0.05 11.18 3.50	8.78	3.05	3.22	3.48	4.94	4.83	31.78 25.25
	-0.25	0.1	8.11	2.33	6.02	1.91	4.43	4.45	6.70	6.05	32.83 27.17
		0.15 5.31	1.40	3.72	1.08	5.06	5.09	7.70	6.85	34.90 28.89
		0.2	3.43	0.83	2.22	0.61	5.50	5.49	8.33	7.47	36.31 29.80
		0	3.55	7.35	2.73	5.13	5.12	5.27	6.65	6.72	34.31 23.17
		0.05 1.63	8.67	1.10	6.89	5.74	5.95	7.44	7.21	31.45 23.92
	-0.5	0.1	0.81	8.44	0.74	6.95	6.76	6.49	8.52	7.68	29.70 23.91
		0.15 0.43	4.67	0.46	44.40	4.10	4.07	5.38	4.79	17.79 13.91
		0.2	0.49	4.81	0.52	39.06	4.63	4.64	6.04	5.44	19.62 14.74
		0	1.00	31.47	1.14	19.10	4.52	4.43	5.05	4.77	17.83 10.68
		0.05 0.84	29.52	1.09	20.05	4.87	5.00	5.50	5.30	16.86 10.97
	-0.75	0.1	1.06	26.15	1.32	19.38	6.00	5.42	6.84	5.68	17.74 10.41
		0.15 1.13	23.42	1.41	18.60	6.26	6.18	7.13	6.47	18.37 11.03
		0.2	1.18	21.69	1.44	18.14	6.45	6.69	7.39	6.94	18.99 11.09

Percentages of relative efficiency (the lesser the better) of each estimator for different levels of contamination and EVI, taking into account 20 order statistics around the optimal MSE value. The minimum value can be seen in bold.

Table 4 :

 4

Maximum likelihood estimators of the parameters of GPD

[START_REF] Embrechts | Modelling Extremal Events for Insurance and Finance[END_REF]

.

Table 5 :

 5 is relatively good. Both Hill on Figure11and t-Hill plot on Figure12confirmed the above values of the EVI, which is positive in this region. Maximum likelihood estimators of the parameters of GEV[START_REF] Fabián | On robust and distribution sensitive hill like method[END_REF].

	Parameter Estimation Standard error
	µ	85.5114413	3.34024537
	ξ	0.3345636	0.09937401
	δ	23.8166815	2.87468837

Table 7 :

 7 The estimators of the coefficients in[START_REF] Furrer | Probability weighted moments properties for small samples[END_REF].

	Parameter Estimation Standard error t value log λ/ξ 15.4089 1.0102 15.25 4.08e -10 Pr(> |t|) -1/ξ -3.2565 0.2013 -16.18 1.87e -10
	Parameter	Estimation Standard error
	µ	0.01233436	1.999988e-06
	ξ	0.19224841	2.150387e-03
	δ	-1.02813582	1.999988e-06

Table 8 :

 8 Estimators of the coefficients in GEV of the error term in[START_REF] Furrer | Probability weighted moments properties for small samples[END_REF].

Table 10 :

 10 . The minimal observed value is 106. According to Estimators of the coefficients in GPD of the error term in[START_REF] Embrechts | Modelling Extremal Events for Insurance and Finance[END_REF].

	Parameter	Estimation Standard error
	ξ	-0.02649701	0.2509478
	δ	61.52754876	21.4716669

Table 11 :

 11 The estimators of the parameters in the GEV (18) model.

Table 12 :

 12 The estimated value of the parameter ξ and δ of the c.d.f.[START_REF] Gascoin | Glacier contribution to streamflow in two headwaters of the huasco river, dry andes of chile[END_REF].

  1.42 1.32 1.29 1.24 1.22 1.17

	Parameter Estimation Standard error
	µ	87.6472012	3.2619087
	ξ	0.4756571	0.3795467
	δ	8.2649512	3.0653720

Table 13 :

 13 The estimated values of the parameters µ, ξ and δ of the c.d.f.[START_REF] Fabián | On robust and distribution sensitive hill like method[END_REF].

Table 14 :

 14 The estimated values of the parameters µ, ξ and δ of the c.d.f.[START_REF] Fabián | On robust and distribution sensitive hill like method[END_REF].

E.g. for β = 1 (Hill or MLE estimator) η(x) = -θ ln ln x+const., x > 1 is the example of η which can be expressed in terms of elementary functions.

Acknowledgement Milan Stehlík acknowledges the support of Fondecyt Proyecto Regular No. 1151441. Pablo Aguirre was partially funded by Fondecyt grant 11150306 and Proyecto Basal CMM Universidad de Chile. This work was also supported by the Slovak Research and Development Agency under the contract No. SK-AT-2015-0019. Jozef Kise ľák was partially supported by grant VEGA M Š SR 1/0344/14.

Appendix A. General information

Let X, X 1 , . . . , X n be independent random variables with common c.d.f. F . In the following, we assume that:

(A.1) F belongs to the maximum domain of attraction of the Weibull distribution.

This assumption entails that F has a finite right endpoint x * and that F has a negative EVI ξ < 0. Our goal is to estimate ξ. To this end, the so-called n-t-Hill estimator is considered:

, where X 1,n ≤ • • • ≤ X n,n are the ordered statistics associated with X 1 , . . . , X n and k ∈ {1, . . . , n}.

Introducing the auxiliary random variable Z := (x * -X) -1 and denoting by G its c.d.f., extremevalue theory shows that G belongs to the maximum domain of attraction of the Fréchet distribution. Under (A.1), G has an infinite right endpoint and a positive EVI given by -ξ. Additionally, the tail quantile function of Z defined by

is regularly varying with index -ξ, i.e. U (tx)/U (t) → x -ξ as t → ∞ for all x > 0.

(A.1)

See [START_REF] Bingham | Regular Variation[END_REF] for a general account on regular variation. In extreme-value theory, the second order condition aims at quantifying the rate of convergence in (A.1):

(A.2) There exist ρ < 0 and some positive or negative function

See for instance [START_REF] De Haan | Extreme Value Theory: An Introduction[END_REF] p. 74.

Appendix B. Main properties

Our main result establishes an asymptotic representation for the n-t-Hill estimator in terms of Weibull and Gaussian random variables. To this end, recall that the Weibull distribution with parameter θ > 0 denoted by W (θ) is defined by the c.d.f.

Appendix C. Proofs

Appendix C.1. Reasoning of lemma 1

From the Implicit Function Theorem, for every z 0 = 0, x 0 ≥ 1, the equation ϕ(x, z) = 0 implicitly defines the function z = g(x) near (x 0 , z 0 ), provided θ -(θ + β -1)x 1-β 0 = 0. Hence, the graph of g(x) determines branches of the set ϕ -1 (0).

If (β, θ) ∈ Ω 1 , let u = x 1-β . Then, the graph of q(u) := θ -(θ + β -1)u is a straight line with positive slope -(θ + β -1). Since q(0) = θ > 0, then q(u) > 0 for every u ≥ 0. It follows that θ -(θ + β -1)x 1-β > 0 and, hence, g(x) > 0 and its graph on the positive (x, z)-plane consists of a single branch M 1 . Moreover, one can readily check that lim x→0 + g(x) = ∞ and lim x→∞ g(x) = 0. In addition,

where

Since F (x) > 0 for every (β, θ) ∈ Ω 1 , then g(x) is a strictly monotone decreasing function. In this way, the set ϕ -1 (0) consists of the two branches M 0 and M 1 , where M 1 is defined as the graph of g(x), that is,

Hence, (adopting the same notation of Ω 1 ) the set T -1 1 (0) = {M 1 } where M 1 is the hyperbola given by z = (θ + 1)/(θx); see the qualitative sketch in Figure 1(b).

If (β, θ) ∈ Ω 2 , from (5) it follows that the line

is a vertical asymptote of ϕ -1 (0). Hence, the graph of z = g(x) consists of two branches. Let M 1 and M 2 be the branches defined for x < x as and x > x as , respectively. Then, it is easy to see that g(x) > 0 for x < x as , and g(x) < 0 for x > x as . Moreover, lim x→0 + g(x) = ∞ and lim x→∞ g(x) = 0. Furthermore, from (C.1) and (C.2) it follows that g(x) has a local minimum at

x as < x as . On the other hand, g (x) > 0 for every x > x as ; hence, M 2 is defined by a strictly increasing function. As a consequence, the graph of ϕ -1 (0) if (β, θ) ∈ Ω 2 is qualitatively as in Figure 1(c).

The analysis for the other cases is similar.

Appendix C.2. Auxiliary lemmas

The first lemma provides an expansion of the t-Hill estimator in terms of two random variables T 1,n and T 2,n which derive its asymptotic distribution.

Appendix C.3. Proofs of auxiliary lemmas Proof of Lemma 2. Introducing Z j = (x * -X j ) -1 for j = 1, . . . , n, the n-t-Hill estimator can be rewritten as 1

Besides, {Z j } j=1,...,n d = {U (Y j )} j=1,...,n where Y 1 , . . . , Y n is a sample of independent random variables from a standard Pareto distribution. Therefore,

and the conclusion follows.

Proof of Lemma 3. From [START_REF] De Haan | Extreme Value Theory: An Introduction[END_REF], Theorem B.2.18, it is possible to control the rest in the convergence (A.2): There exists a function A 0 asymptotically equivalent to A such that, for all ε > 0 there exists t 0 > 0 such that, for all t ≥ t 0 and x ≥ 1,

Letting R 1 (t, x) := (x ρ -1)/ρ + R 0 (t, x), it follows that

It straightforwardly follows that

.

Since, for t large enough, |A 0 (t)| ≤ 1/(2ε -2/ρ), the second part (C.4) of the lemma follows.

Proof of Lemma 4.. Applying expansion (C.3) of Lemma 3 with t = Y n-k:n P → ∞ (see [START_REF] De Haan | Extreme Value Theory: An Introduction[END_REF], Corollary 2.2.2) and x = Y n:n /Y n-k:n ≥ 1 yields k:k /k converges in distribution to the extreme-value distribution Φ 1 (see for instance [START_REF] Embrechts | Modelling Extremal Events for Insurance and Finance[END_REF], Table 3

From the definition of T 5,n , it follows that, for all ε ∈ (0, 1ξρ),

.

Similarly, from the definition of T 5,n , it follows that, for all ε ∈ (0, 1ξρ),

in view of Lemma 4 and 5. Let θ = min(1/2, -ξ). It follows that

under the assumption k θ A(n/k) → λ as n → ∞. As a consequence, ξ → P ξ and therefore

and the result is proved.

Proof of Corollary 1. Two cases arise. If ξ < -1/2 then θ = 1/2 and thus k 1/2 ( ξξ) is asymptotically Gaussian with mean λ(1ξ) 2 τ (ξ, ρ) and variance σ 2 (ξ)(1ξ) 4 . Conversely, if ξ > -1/2, then θ = -ξ and k -ξ ( ξξ) converges in distribution to λ(1ξ) 2 τ (ξ, ρ) + (1ξ)ξW (-1/ξ) where W (-1/ξ) is the Weibull distribution with shape parameter -1/ξ.

Appendix C.5. Explanation for continuation in Auto [START_REF] Doedel | AUTO: A program for the automatic bifurcation analysis of autonomous systems[END_REF] In order to compute a family of solutions of a three-dimensional system

we consider a function u : [0, 1] → R 3 satisfying a rescaled version of (C.11) given by the differential equation u = T f (u(t)), (C.12)

where f stands for the vector field defined by the system (4) and T > 0 is the integration time (also known as the "period") of an orbit segment of f . Note that in (C.12), the period T appears as an explicit parameter and the actual integration time over an orbit segment is always 1. Geometrically, the function u represents a unique orbit segment {u(t) = (x(t), y(t), z(t)) ∈ R 3 | 0 ≤ t ≤ 1} provided that suitable boundary conditions are posed at one or both end points u(0) and u(1). In our case, we consider

where µ, ν are dummy parameters that determine the 'initial' coordinates y 0 and z 0 , respectively, of a given solution.

The boundary value problem (C.12)-(C.13) defines a (µ, ν, T )-dependent family of orbit segments. For any fixed (µ, T ) = (µ 0 , T 0 ) we have a uniquely defined one-parameter family of orbit segments (parameterised by ν) with fixed integration time T 0 and fixed initial coordinates x(0) = 1, y(0) = µ 0 . In order to compute this ν-family by continuation in Auto [START_REF] Doedel | AUTO: A program for the automatic bifurcation analysis of autonomous systems[END_REF][START_REF] Doedel | AUTO-07p: Continuation and bifurcation software for ordinary differential equations[END_REF], we need to specify an initial orbit segment u ν0 that satisfies (C.12)-(C.13) for some fixed ν = ν 0 . To this end, a possible choice is to take advantage of the continuous extension of (4) to the line z = 0 and consider the constant solution -i.e., a trivial orbit segment-u(t) ≡ (1, 0, 0) of (C.12) with T = 0; continuation in T for fixed u(0) = (1, 0, 0) up to T = T 0 yields the desired initial orbit segment u 0 satisfying (C.12)-(C.13) for µ 0 = 0 and ν 0 = 0. A collection of orbits segments is then obtained by fixing T = T 0 and continuing u 0 in µ up to a user-defined value; in our case, the stopping condition is µ = 1. Finally, the desired family of orbit segments is obtained by allowing ν to vary while keeping T = T 0 and µ = 1 fixed. In this way, each orbit segment in this family corresponds to a solution {x(t) = (x(t), y(t), z(t)) ∈ R 3 | 0 ≤ t ≤ T 0 } of (4) with initial condition (x 0 , y 0 , z 0 ) = (1, 1, ν).