Pierre Ribereau

Esterina Masiello

Philippe Naveau

Skew Generalized Extreme Value Distribution: Probability Weighted Moments Estimation and Application to Block Maxima Procedure

Keywords: AMS Subject classifications: 62F10, 62E17, 62P12 Skew distributions, Generalized Extreme Value distribution, Extreme Value Theory, Probability Weighted Moments, Return Level Estimation

) on the skew normal distribution, we propose an extension of the Generalized Extreme Value (GEV) distribution, the SGEV. This new distribution allows for a better fit of maxima and can be interpreted as both the distribution of maxima when maxima are taken on dependent data and when maxima are taken over a random block size. We propose to estimate the parameters of the SGEV distribution via the Probability Weighted Moments method. A simulation study is presented to provide an application of the SGEV on block maxima procedure and return level estimation. The proposed method is also implemented on a real-life data.

Skew distributions and Extreme Value Theory

The idea of modifying a symmetric probability density function of a random variable (r.v.) in a multiplicative fashion, introducing skewness, has been implemented for the univariate normal distribution by Azzalini ([2], [START_REF] Azzalini | Further results on a class of distributions which includes the normal ones[END_REF]) yielding the so-called skew-normal distribution having the probability density function of the form 2φ(x)Φ(λx), x, λ ∈ R (where φ(•) and Φ(•) denote respectively the probability distribution function (pdf) and cumulative distribution function (cdf) of the standard normal distribution). This distribution presents the advantage to look like a normal one in the core data but allowing some asymmetry. The model has since been extended to the gaussian related distributions (multivariate normal or Student distribution) and its extensions are provided in the book edited by Genton [START_REF] Genton | Skew-Elliptical Distributions and Their Applications: A Journey Beyond Normality[END_REF] and in the review paper [START_REF] Azzalini | The skew-normal distribution and related multivariate families[END_REF] and the monograph by Azzalini [START_REF] Azzalini | The Skew-Normal and Related Families[END_REF] who regularly updates list of references on his web page. The skew-normal distribution has been also applied in many fields like climatology (Flecher et al. [START_REF] Flecher | A Stochastic Daily Weather Generator for Skewed Data[END_REF]) or insurance (Vernic [START_REF] Vernic | Multivariate skew-normal distributions with applications in insurance[END_REF]) and adapted to classical statistical tools like Space-state models (Naveau et al. [START_REF] Naveau | A skewed Kalman filter[END_REF]) or random fields (Allard and Naveau [START_REF] Allard | A New Spatial Skew-Normal Random Field Model[END_REF]). Since then, several authors have tried to generalize these results to skewing arbitrary symmetric (or not) pdf's. See Gupta et al. [START_REF] Gupta | Some skewsymmetric models[END_REF] for a most detailed discussion of skew-symmetric models based on the normal, Student's t, Cauchy, Laplace, logistic and uniform distributions. In that case, the main feature of these models is that a new parameter λ is introduced to control skewness and kurtosis. Nadarajah and Kotz ( [START_REF] Nadarajah | Skew distributions generated from different families[END_REF], [START_REF] Nadarajah | Skew Models I[END_REF] and [START_REF] Nadarajah | Skew Models II[END_REF]) proposed some skew models based on Cauchy, Laplace, Logistic, Bessel and Pearson (among others). Some relations between Extreme Value Theory and Skew Distribution have already been studied by Chang and Genton [START_REF] Chang | Extreme Value Distributions for the Skew-Symmetric Family of Distributions[END_REF] in the case of skewsymmetric family of distributions and by Lysenko et al. [START_REF] Lysenko | Multivariate extremes of generalized skew-normal distributions[END_REF] in the case of multivariate generalized skew-normal distribution.

In the same spirit, we propose here an extension of the Generalized Extreme Value (GEV) distribution. This distribution, introduced by Jenkinson [START_REF] Jenkinson | The frequency distribution of the annual maximum (or minimum) values of meteorological elements[END_REF], is the limiting distribution for extremes, both normalized maxima and minima of a random sample. In a more formal way, it is known that the only possible limiting form of a normalized maximum of a random sample when a non-degenerate limit exists is captured by the GEV

G(x; σ, γ, µ) =        exp -(1 + γ x-µ σ ) -1/γ , if 1 + γ x-µ σ > 0, γ = 0, exp -exp(-x-µ σ ) , if x ∈ R, γ = 0,
where µ ∈ R, σ > 0 and γ are called the location, scale and shape parameters, respectively.

The GEV model has been applied to environmental sciences to fit, for example, yearly maximal temperatures or yearly maximal river discharges. However, the limiting distribution does not always provide a good fit in the finite samples that arise in practice, since the extreme value theory is based on asymptotic arguments.

Our purpose, here, is to generalize the GEV distribution by addition of an extra parameter, λ, introduced to model flexibility in the data and so better adapt to real situations. The main features of this distribution will be discussed together with two methods of estimation of its parameters.

The paper is organized as follows. In Section 2, we will define the distribution, study some of its characteristics and present some possible interpretations. In Section 3, the parameters of the distribution will be estimated by the classical methods of maximum likelihood and probability-weighted moments. Finally, a simulation study and an application to real data will be performed (Section 4 and 5) to investigate the behavior of the estimator in small and moderate samples.

The skew generalized extreme value distribution

Definition 1. A random variable X is said to have a Skew Generalized Extreme Value (SGEV) distribution if its probability density function is given by f

(x) = (λ + 1)g(x)G λ (x) (1)
where λ > -1, g(.) and G(.) are respectively the density function and the cumulative distribution function of a GEV distribution with location parameter µ, scale parameter σ and shape parameter γ.

The cumulative distribution fucntion is given by

F (x) = G λ+1 (x).
The distribution will be denoted for short by SGEV (µ, σ, γ, λ). It follows from (1) that the pdf of X can be written as

f (x) = λ+1 σ 1 + γ x-µ σ -1/γ-1 exp -(1 + λ) 1 + γ x-µ σ -1/γ if 1 + γ x-µ σ > 0, γ = 0, ( λ+1 σ ) exp -(λ + 1) exp -x-µ σ exp -x-µ σ if x ∈ R, γ = 0.
When λ = 0, equation ( 1) reduces to the GEV pdf. The parameter λ allows for a greater degree of flexibility and we can expect this to be useful in many more practical situations.

The quantile function of the cdf F is given by

Q(u) =      µ + σ γ -log u λ+1 -γ -1 if γ = 0, µ -σ log -log u λ+1 if γ = 0
and allows to compute quantiles in an easy way.

We also find that the expectation, the second-order moment and the variance of X are respectively given by

E(X) = µ - σ γ [1 -(λ + 1) γ Γ(1 -γ)] , E(X 2 ) = µ - σ γ + σ γ(λ + 1) -γ Γ(1 -γ) 2 , V (X) = σ 2 γ 2 (λ + 1) -2γ Γ(1 -2γ) -Γ 2 (1 -γ) , if λ > -1
, where Γ denotes the Gamma function Γ(t) = ∞ 0 e -u u t-1 du, t > 0. Simple modifications of Chang and Genton [START_REF] Chang | Extreme Value Distributions for the Skew-Symmetric Family of Distributions[END_REF] leads to the conclusion that the extreme value index of the SGEV distribution is γ.

The SGEV distribution offers some interesting practical interpretations which are going to be presented in the following remarks.

Remark 1. Let X 1 , X 2 , . . . be i.i.d. r.v.'s with cdf F and let N be a Poisson r.v. with intensity λ * > 0. Moreover, let N be independent of {X i } i≥1 . Define the following maximum with a random block size

M = max(X 1 , . . . , X N ) if N ≥ 1, 0 if N = 0.
Then, we have

P (M ≤ x) = +∞ n=0 F n (x) λ * n e -λ * n! = exp {-λ * [1 -F (x)]} = exp -F (x) λ * . Now, denote G(x) = exp[-F (x)]. It follows P (M ≤ x) = G λ * (x) (2)
which corresponds to SGEVs with strictly positive parameters λ * = λ + 1.

Define the cdf H(x; σ, γ, µ) by

H(x; σ, γ, µ) =      1 -1 + γ( x-µ σ ) -1/γ , if γ = 0, 1 -exp -( x-µ σ ) , if γ = 0,
where x ≥ µ if γ ≥ 0, µ ≤ x ≤ µ -σ/γ if γ < 0 and µ ∈ R is the location parameter, σ > 0 the scale parameter and γ ∈ R the shape parameter. This distribution, which appears as the limit distribution of scaled excesses over high thresholds, is called a Generalized Pareto Distribution (GPD) (Pickands [25]).

If F is a GPD tail, then the cdf of M is the cdf of a SGEV on the interval [µ, ∞[, as it is shown in Figure 1. If x ≤ µ, then P(M ≤ x) = P(N = 0) = exp(-λ * ). The equality (2) offers then a first interpretation for the SGEV distribution: the cdf G λ * (x), λ * = λ + 1, on the interval [µ, ∞[ can be viewed as the cdf of a maximum taken over exceedances in a random Poisson block size.

Remark 2. The SGEV is also the distribution of the maxima when they are taken on dependent data. In that case, λ is between -1 and 0. The extremal index θ (see, for example, Beirlant et al. [START_REF] Beirlant | Statistics of Extremes: Theory and Applications[END_REF]) is then equal to λ + 1. The SGEV is more general since the value of λ needs only to be greater than -1. So, the family of SGEV includes the classical GEV one, but also the GEV one based on dependent data. This is quite interesting because this property could be used to define a new plug-in estimator of the extremal index, θ = λ + 1, where λ is the estimator of λ obtained by the PWM method.

Parameter estimation

Let (X 1 , . . . , X n ) be a sample of SGEV r.v.s. In this section, the parameters (µ, σ, γ, λ) of the SGEV distribution will be estimated by both the method of maximum likelihood (ML) and the probability weighted moments (PWM) method.

Maximum likelihood

In case γ = 0, the log-likelihood function for a sample (X 1 , . . . , X n ) of SGEV r.v.s is given by

logL(µ, σ, γ, λ) = n log(λ + 1) -n log σ - 1 γ + 1 n i=1 log 1 + γ x i -µ σ -(λ + 1) n i=1 1 + γ x i -µ σ -1/γ . If γ = 0, the log-likelihood function reduces to logL(µ, σ, λ) = n log(λ + 1) -n log σ -(λ + 1) n i=1 exp - x i -µ σ - n i=1 x i -µ σ .
The ML estimator (μ, σ, γ, λ) is obtained by maximizing the log-likelihood function. Clearly, no explicit solution exists to these equations so that numerical procedures are called for. For a discussion about the properties of the MLE see Smith [START_REF] Smith | Maximum likelihood estimation in a class of nonregular cases[END_REF]. In the paper by Smith [START_REF] Smith | Maximum likelihood estimation in a class of nonregular cases[END_REF], the ML estimators of the parameters of the GEV distribution are shown to be asymptotically normal whenever γ > -1/2. Since the SGEV distribution is based on the GEV distribution, it is easy to check that the classical properties of the ML estimator of (µ, σ, γ, λ) hold whenever γ > -1/2. Remark 3. The ML method usually gives poor results when it has to estimate a skew parameter: for example, for moderate sample sizes, the ML estimator is infinite with positive probability when dealing with skew normal distribution (Sartori [26]). The issue mentioned here has a close connection with the Fisher singularity problem discussed in Hallin andLey (2012a, 2012b). In Section 4.1 and in Tables 1 to 4, the results obtained by estimating the parameters by Maximum Likelihood are disappointing. This is the reason why, when estimating the return level in Section 4.2, we will avoid the ML method and prefer the Probability Weighted Moments one, which will be described in the next section.

The method of probability weighted moments

The probability weighted moments, introduced by Greenwood et al. [START_REF] Greenwood | Probability weighted moments: definition and relation to parameters of several distributions expressable in inverse form[END_REF], are a simple generalization of the usual moments of a probability distribution. They are defined in terms of the cdf H of a r.v. Y as

E{Y p [H(Y )] r [1 -H(Y )] s }
with p, r and s real numbers. Particularly useful special cases are the probability-weighted moments

β r = E{Y [H(Y )] r }
which have been used by [START_REF] Hosking | Estimation of the Generalized Extreme-Value distribution by the method of probabilityweighted moments[END_REF] to estimate the parameters of the GEV distribution. For the GEV distribution, they equal

β r = (r + 1) -1 µ - σ γ 1 - Γ(1 -γ) (1 + r) -γ γ < 1, γ = 0. ( 3
)
This method of moments has been also applied to the truncated skewnormal distribution (Flecher et al. [START_REF] Flecher | Estimating the Closed Skew-Normal distributions parameters using weighted moments[END_REF]). It can be equally applied to estimate the parameters of the SGEV distribution. Let X be a SGEV (µ, σ, γ, λ) r.v. with cdf F and let Y be a GEV (µ, σ, γ) r.v. with cdf G. It is easy to check that the probability weighted moments of a SGEV distribution, denoted by ν a , are special probability weighted moments of a GEV distribution:

ν a = E[XF a (X)] = xF a (x)f (x)dx = xG (λ+1)a (x)(λ + 1)g(x)G λ (x)dx = (λ + 1) xG (λ+1)a+λ (x)g(x)dx = (λ + 1)E[Y G λ(a+1)+a (Y )] = (λ + 1)β λ(a+1)+a
By taking into account (3), we obtain

ν a = E[XF a (X)] = 1 a + 1 µ - σ γ 1 - Γ(1 -γ) ((λ + 1)(a + 1)) -γ .
In practice, we just need to choose four different values of a and so set up the necessary number of moment equations. Here, we choose a = 0, 1, 2, 3, because they are trivial to estimate. The estimators of γ, σ, µ and λ are then obtained by solving the following system of equations:

                         ν 0 = µ -σ γ [1 -(λ + 1) γ Γ(1 -γ)] , 2ν 1 = µ -σ γ [1 -2 γ (λ + 1) γ Γ(1 -γ)] , 3ν 2 = µ -σ γ [1 -3 γ (λ + 1) γ Γ(1 -γ)] , 4ν 3 = µ -σ γ [1 -4 γ (λ + 1) γ Γ(1 -γ)] , (4)
after replacing ν a by its empirical version νa :

νa = 1 n n i=1 X (i) i n a where X (1) ≤ X (2) ≤ ... ≤ X (n)
represent the ordered sample. This system is quiet hard to solve numerically.

A promising avenue consists of first estimating γ by solving the following equation 3ν

2 -ν0 2ν 1 -ν0 = 3 γ -1 2 γ -1 . (5)
Once γ given, find the estimations of λ and σ by solving the following system of equations

         2ν 1 -ν0 = σ γ Γ(1 -γ)( λ + 1) γ (2 γ -1) , 4ν 3 -3ν 2 = σ γ Γ(1 -γ)( λ + 1) γ (4 γ -3 γ ) , (6)
and finally, σ, γ and λ given, obtain the estimation of µ, that is

μ = ν0 + σ γ 1 -( λ + 1) γ Γ(1 -γ) .
The case γ = 0 can also be covered by this method: a system of three equations with three unknown parameters will then be obtained by taking the first three equations in the system (4) and letting γ go to 0. Nevertheless, in the case γ = 0, all moments exist and then, we do not need a new method based on moments to get good results since the classical moments method will be satisfactory.

To derive the limiting distribution of (μ, σ, γ, λ), we need to study the limiting behavior of νa . By direct application of the central limit theorem, one can easily prove that νa is asymptotically normally distributed with mean ν a and asymptotic variance

(λ + 1) 2 n σ 2 γ 2 [(a + 1)(λ + 1)] 2γ Γ(1 -2γ)K a(λ + 1) + λ (a + 1)(λ + 1) -Γ 2 (1 -γ)
where K denotes the hypergeometric function

K(x) = 1 F 2 (k, 2k; 1 + k; -x).
A similar argument applies to any linear combination of the νa (a = 0, 1, 2, . . .) and it follows that the νa are asymptotically jointly normal.

Theorem 1. Let (X 1 , . . . , X n ) be a sample of SGEV. If γ < 0.5, then the rescaled PWM estimator vector defined by (4) and denoted by

√ n     ν0 ν1 ν2 ν3     converges in distribution, as n → ∞, towards N         ν 0 ν 1 ν 2 ν 3     , Ω    
where Ω denotes the covariance matrix.

The elements of the covariance matrix Ω are given in Appendix A.

Corollary 2. The estimators of µ, σ, γ and λ obtained by the probability weighted moments method are asymptotically normally distributed and are asymptotically unbiased if γ < 0.5.

The proof of the corollary is a direct application of Theorem 1 and the Delta Method.

Remark 4. Moments ν 0 , ..., ν 3 have been chosen in the estimation procedure of the parameters of the SGEV distribution. The same moments could be used to obtain different moment methods. Equations similar to (5) could be written down with any of the 4 sets of 3 moments, among ν 0 , ..., ν 3 . Similarly, any pair of moments among the 6 possible pairs can lead to a system of two equations similar to that in [START_REF] Beirlant | Statistics of Extremes: Theory and Applications[END_REF]. Finally, µ can also be estimated using any of the 4 moments. The estimation variance of the selected combination of moments is the lowest among all possible combinations, since the moments of lowest degree have been selected. This procedure can then be improved using all combinations simultaneously. Consider for example the estimation of γ. Four equations similar to (5) can be written using each triplet of moments. An estimator γ can be defined as the value minimizing the mean squared errors between the left hand side and the right hand side of these equations. A similar approach can then be implemented for estimating σ and λ given γ, and finally for estimating µ given all other estimates.

A simulation study 4.1 Parameters estimation

We perform a simulation study in order to compare the performances of the ML and PWM estimators of the parameters of the SGEV distribution.

Simulations were performed for sample size n = 100 with γ taking values -0.5, 0, 0.5, 1 and 1.5. Without loss of generality, the location and scale parameters were set to µ = 0 and σ = 1. We set also the skew parameter λ to 0.2, 0.5, 1 and 1.5. For each combination of values of γ and λ, we generated 100 random samples from the SGEV distribution and for each sample, we estimate the parameters µ, σ, γ and λ by ML and PWM method. The results are listed in Tables 1 to 4.

Remark that some chosen values for γ are not in the range of validity of the estimation methods. In practice, on real data, the value of γ is unknown, so values which are not in the validity range were used to discuss the robustness of the method. For example, when γ is greater than 0.5 but lesser than 1, moments still exist but the estimators are no more asymptotically normal. In this case, the PWM estimator gives acceptable results when estimating µ and σ (sometimes λ). On the other side, when γ ≥ 1 the PWM method is not satisfactory anymore.

As previously mentioned, the Maximum Likelihood method does not perform well especially when estimating the skew parameter λ. Indeed, in the case γ < 0.5, we can notice the superiority of the PWM method for estimating the skew parameter λ, while ML and PWM provides good and comparable estimations for the three other parameters. In particular, because of the quite large sample size (n = 100) in our simulation study, the ML estimator is preferable to estimate γ (when γ < 0.5). On the other side, both methods do not perform well when γ gets greater (if γ > 0.5), especially for the skew parameter. For the PWM method, this is due to the usual regularity conditions underlying the asymptotic properties of the estimators (γ < 0.5) while for the ML method this may be due to some algorithmic problems. This shows that both methods are not really robust when used outside their range of validity. For a more detailed discussion about the performance of both ML and PWM estimators, see Hosking and Wallis [START_REF] Hosking | Estimation of the Generalized Extreme-Value distribution by the method of probabilityweighted moments[END_REF]. Remark 5. When looking at the results of the simulation study, it appears that the parameter γ is always fairly well estimated by ML. This could suggest to consider a mixed estimation procedure which will consist in estimating γ with the ML method and the other parameters with the PWM method. This kind of procedure would lead to a vector of estimators whose asymptotic law could be quite hard to obtain. An alternative could be to consider the estimation obtained by the PWM method as a starting value for the ML method.

Return level estimation

Let Z 1 ,...,Z n be a sample of i.i.d. random variables from a distribution function H. We are now interested in estimating z p := H -1 (1 -p). The method of Block Maxima consists in replacing the original sample Z 1 , ..., Z n by Y 1 , ..., Y m where Y i is the maximum on the i th block of size k i.e. Y i = max{Z (i-1)k+1 , ...., Z ik } with n = m × k. Since the r.v.'s Y i are asymptotically GEV distributed, the return level can be estimated using (see Beirlant et al. [START_REF] Beirlant | Statistics of Extremes: Theory and Applications[END_REF])

ẑ * p =      µ Y + σ Y γ Y -log(1 -p) k -γ Y -1 , if γ Y = 0, µ Y -σ Y log -log(1 -p) k , if γ Y = 0, (7)
leading to the GEV estimator ẑGEV If Z 1 ,...,Z n are i.i.d. with gaussian distribution N (0, 1) we also have (see Embrechts et al. [START_REF] Embrechts | Modelling Extremal Events for Insurance and Finance[END_REF])

P Z (n,n) -d N n c N n ≤ z → GEV (µ = 0, σ = 1, γ = 0)
where Z (n,n) is the maximum of the sample,

c N n = (2 log n) -1/2 ,
and

d N n = 2 log n - log(4π) + log log n 2(2 log n) 1/2 .
The distribution of Y i can be approximated by a Obviously, the theoretical estimation z BM p will get closer to z p as the block size increases. In practice, we need to choose a block size large enough to ensure a good approximation of the GEV distribution but allowing also enough maxima observations to obtain a good parameter estimation. In that way, the flexibility of the SGEV distribution can give a good compromise. Once the sample of maxima obtained, the return level can be estimated by

GEV (µ = d N k , σ = c N k , γ = 0)
ẑSGEV p =      µ s + σs γs -log(1 -p) k/(λs+1) -γs -1 , if γ s = 0, µ s -σ s log -log(1 -p) k/(λs+1) , if γ s = 0, (8)
where µ s , σ s , γ s et λ s are the parameter estimators obtained after fitting a SGEV on the sample of maxima Y 1 , . . . , Y m .

Similar results stand for the Log-gamma distribution LΓ(α, β) with

c LΓ n = (Γ(β)) -1 (log n) β-1 n 1/α and d n = 0.
We want to answer the two following questions : in the case of slow convergence of the maximum towards the GEV distribution (e.g. the Gaussian case), how the SGEV behaves and how does the SGEV for heavy tailed distributions? Figures 2 and3 represent the estimation of the quantile as a function of the block size. For each chosen distribution (N (0, 1) in Figure 2 and Loggamma(2,2) in Figure 3), we simulated 500 samples of size n = 500 or n = 5000. For each combination, we compute the theoretical value z BM p (red line) to show how good the approximation by a GEV distribution is, the estimation of the quantile based on the maximum by block using both GEV (ẑ GEV p , grey line) and SGEV (ẑ SGEV p , black line) and the 95% empirical confidence intervals (dashed lines). The true value of the quantile is represented by the horizontal line. We also represent, on the bottom row of each figure, the value of the skew parameter estimation λ s , together with a 95% empirical confidence interval. By a first look at these simulations, we observe that the return level estimation is not so good when dealing with Gaussian distributions (cf Figure 2). This is not surprising since it is well known that the normal distribution converges very slowly to the GEV distribution, but it is also an indicator of the limitations of the Extreme Value Theory. Although the SGEV approximation is not great, it seems to propose an improvement in terms of bias and variance.

When dealing with a heavy tailed distribution (cf Figure 3), the quality of the return level estimation using the SGEV approximation gets better even if the approximation of the sample of maxima seems good (see the theoretical z BM p ).

In view of these simulations, the SGEV proposes a complementary method to the classical Extreme Value Theory by improving the return level estimation both in the case of a slow convergence or in the heavy tailed case.

An application to real data

We apply the proposed procedure to the Fort Collins precipitation data. The data consist of the annual daily precipitation maxima (in mm) recorded by a weather station from 1948 to 2001 in Fort Collins (Colorado, USA) (cf Guillou et al. [START_REF] Guillou | Return level bounds for discrete and continuous random variables[END_REF] or Diebolt et al. [START_REF] Diebolt | Improving Probability-Weighted Moments Methods for the Generalized Extreme Value Distribution[END_REF] for more details and explanation on this data set).

The values of the estimated parameters, obtained by ML and PWM methods, are listed in Table 5 when fitting a SGEV distribution or a GEV distribution. For the GEV distribution, the estimation of the return level associated to a period of 100 years is equal to 573.72 (when adopting the PWM method) and 653.33 (when adopting the ML method). By fitting the SGEV distribution, we obtained an estimation of the return level period equal to 608.83 (with the PWM method) and 653.21 (with the ML method). This application was selected to illustrate a number of practical points. First of all, we remark that if we look at the same estimation method (either ML or PWM method) the return level estimation does not vary so much if we take a GEV or a SGEV tail approximation. This is probably due to the fact that the annual maximum is taken on a sample of important size (equal to 365) and so the GEV approximation is good enough to provide interesting results. Indeed, the values of λ are near 0 indicating that the GEV distribution may be a good choice of fit for this sample. In addition, we remark that there is no loss using SGEV instead of GEV distribution. So, in our opinion, the first one should be preferred since it englobes more cases than the classical GEV distribution. Indeed, the SGEV distribution could be used to define the right block size. One can try different block size and choose the one which gives the value of the parameter blocks. Two interpretations of the SGEV distribution can be given. First, the cdf can be viewed as the cdf of a maximum over a random block size when the maximum is taken over exceedances (see Remark 1). Secondly, the SGEV distribution can be seen as the limiting distribution of the maxima when the original data are not independent. In that case, the parameter λ is just θ -1, where θ is the extremal index. This means that the estimator of λ, obtained by the PWM method, can be looked at as a new and efficient estimator of the extremal index θ.

Moreover, the SGEV could be used to find the right block size, the one such that λ is nearest to zero.

Authors would like to thank two anonymous referees for their helpful comments that led to an improvement of the present paper (we especially thank for Remarks 4 and 5).

A Asymptotic covariance matrix

The asymptotic covariance matrix of the estimator vector obtained by the probability weighted moments method (see Theorem 1) can be written as

Ω = (λ + 1) 2 σ 2 γ 2     Ω 0,0 Ω 0,1 Ω 0,2 Ω 0,3 Ω 0,1 Ω 1,1 Ω 1,2 Ω 1,3 Ω 0,2 Ω 1,2 Ω 2,2 Ω 2,3 Ω 0,3 Ω 1,3 Ω 2,3 Ω 3,3    
and the elements on the main diagonal of Ω can be specified as follows -(λ(a + 2) + a + 1) γ Γ(1 -2γ)K (λ + 1)(a + 1) λ(a + 2) + a + 1 +2 ((λ + 1)(a + 1)) γ × (λ(a + 2) + a + 1) γ ((λ + 1)(a + 2)) γ Γ 2 (1 -γ) , Ω a,a+2 = 1 2 ((a + 3)(λ + 1)) 2γ Γ(1 -2γ)K λ(a + 1) + a (λ + 1)(a + 3)

Ω 0,0 = (λ + 1) 2γ Γ(1 -2γ)K λ λ + 1 -Γ 2 (1 -γ) ,
-(λ(a + 3) + a + 2) γ Γ(1 -2γ)K (λ + 1)(a + 1) λ(a + 3) + a + 2 +2 ((λ + 1)(a + 1)) γ × (λ(a + 3) + a + 2) γ ((λ + 1)(a + 3)) γ Γ 2 (1 -γ) , Ω a,a+3 = 1 2 ((a + 4)(λ + 1)) 2γ Γ(1 -2γ)K λ(a + 1) + a (λ + 1)(a + 4)

-(λ(a + 4) + a + 3) γ Γ(1 -2γ)K (λ + 1)(a + 1) λ(a + 4) + a + 3 +2 ((λ + 1)(a + 1)) γ × (λ(a + 4) + a + 3) γ ((λ + 1)(a + 4)) γ Γ 2 (1 -γ) .

p

  after replacing µ Y , σ Y and γ Y by μY , σY and γY , the ML or PWM estimators of the GEV parameters based on the sample Y 1 , ..., Y m .

  and the return level of probability p can be approximated by z BM p after replacing µ Y and σ Y by d N k and c N k , and γ Y by 0 in Equation (7) (in the case γ Y = 0).

Figure 1 :

 1 Figure 1: Cumulative distribution function of M if X i ∼ GP D(µ = 10, σ = 1, γ = 0) and λ = 1 (full line) of the Remark 1 and cumulative distribution function of the GEV distribution (dashed line).

Figure 3 :

 3 Figure 3: Left hand : The same as left hand in Figure 2 but with a Loggamma(2,2) distribution (γ = 1/4). Right hand : Same as left hand but with p = 1/2500.

Ω 1 , 1 =

 11 (2λ + 2) 2γ Γ(1 -2γ)K 2λ + 1 2λ + 2 -Γ 2 (1 -γ) , Ω 2,2 = (3λ + 3) 2γ Γ(1 -2γ)K 3λ + 2 3λ + 3 -Γ 2 (1 -γ) , Ω 3,3 = (4λ + 4) 2γ Γ(1 -2γ)K 4λ + 3 4λ + 4 -Γ 2 (1 -γ) ,where K and Γ are respectively the hypergeometric and gamma functions. The remaining elements of the matrix areΩ a,a+1 = 1 2 ((a + 2)(λ + 1)) 2γ Γ(1 -2γ)K λ(a + 1) + a (λ + 1)(a + 2)

Table 1 :

 1 For each combination of µ, σ and γ, we generate 100 samples of SGEV(µ,σ,γ,λ = 0.2) and calculate the PWM and ML estimations, together with the respective mean and standard deviation (sd).

	Parameter True value Mean for ML Mean for PWM sd for ML sd for PWM
	µ	0	0.184	-0.014	0.452	0.103
	σ	1	0.908	0.953	0.288	0.080
	γ	-0.5	-0.517	-0.433	0.065	0.075
	λ	0.2	0.038	0.213	0.620	0.028
	µ	0	0.097	0.009	0.075	0.105
	σ	1	1.020	0.955	0.086	0.089
	γ	0	0.0168	0.039	0.051	0.072
	λ	0.2	0.066	0.196	0.086	0.007
	µ	0	-0.131	0.018	0.086	0.093
	σ	1	0.921	1.002	0.103	0.126
	γ	0.5	0.511	0.467	0.107	0.104
	λ	0.2	0.383	0.192	0.128	0.069
	µ	0	-0.317	-0.024	0.068	0.264
	σ	1	0.666	1.248	0.096	0.331
	γ	1	1.047	0.772	0.138	0.109
	λ	0.2	0.773	0.364	0.191	0.316
	µ	0	-0.399	-0.528	0.058	0.805
	σ	1	0.388	1.780	0.087	0.877
	γ	1.5	1.492	0.897	0.159	0.085
	λ	0.2	1.317	0.916	0.253	0.937
	λ nearest to zero.
	6 Conclusions
	A new class of distributions, called SGEV distributions, has been proposed
	as a generalization of the classical GEV distribution, thanks to the intro-
	duction of a new parameter λ. Estimators of the parameters and quantiles
	of the SGEV distribution have been derived by ML and PWM methods.
	As mentioned in the above sections, the SGEV distribution is more flexible
	than the classical GEV distribution and so, can fit better maxima when deal-
	ing with small maxima samples or when the maxima come from small size

Table 2 :

 2 For each combination of µ, σ and γ, we generate 100 samples of SGEV(µ,σ,γ,λ = 0.5) and calculate the PWM and ML estimations, together with the respective mean and standard deviation (sd).

	Parameter True value Mean for ML Mean for PWM sd for ML sd for PWM
	µ	0	0.022	0.006	0.445	0.094
	σ	1	1.060	0.951	0.353	0.128
	γ	-0.5	-0.349	-0.439	0.786	0.094
	λ	1	0.124	1.005	1.330	0.027
	µ	0	0.241	0.009	0.249	0.089
	σ	1	1.048	0.965	0.177	0.100
	γ	0	0.196	0.031	0.424	0.069
	λ	1	0.383	0.997	0.331	0.004
	µ	0	0.157	-0.014	0.107	0.105
	σ	1	1.081	1.043	0.116	0.173
	γ	0.5	0.513	0.454	0.128	0.117
	λ	1	0.702	1.004	0.186	0.044
	µ	0	0.011	-0.273	0.074	0.429
	σ	1	0.988	1.326	0.106	0.444
	γ	1	1.028	0.763	0.132	0.119
	λ	1	0.958	1.408	0.162	0.830
	µ	0	-0.100	-2.445	0.085	1.914
	σ	1	0.827	2.334	0.166	1.825
	γ	1.5	1.514	0.894	0.169	0.077
	λ	1	1.262	3.123	0.219	2.416

Table 3 :

 3 For each combination of µ, σ and γ, we generate 100 samples of SGEV(µ,σ,γ,λ = 1) and calculate the PWM and ML estimations, together with the respective mean and standard deviation (sd).

	Parameter True value Mean for ML Mean for PWM sd for ML sd for PWM
	µ	0	0.137	-0.006	0.272	0.098
	σ	1	0.936	0.935	0.222	0.119
	γ	-0.5	-0.523	-0.421	1.005	0.085
	λ	1.5	0.185	1.507	0.808	0.0198
	µ	0	0.179	0.009	0.345	0.100
	σ	1	1.029	0.956	0.192	0.126
	γ	0	0.591	0.027	0.746	0.069
	λ	1.5	0.596	1.497	0.755	0.005
	µ	0	0.185	0.011	0.262	0.124
	σ	1	1.114	1.021	0.180	0.206
	γ	0.5	0.694	0.471	0.407	0.123
	λ	1.5	1.089	1.508	1.991	0.058
	µ	0	0.138	-0.461	0.103	0.676
	σ	1	1.104	1.389	0.154	0.522
	γ	1	1.041	0.770	0.158	0.121
	λ	1.5	1.222	2.219	0.188	1.473
	µ	0	0.006	-2.905	0.070	2.740
	σ	1	0.969	2.732	0.134	1.403
	γ	1.5	1.509	0.904	0.171	0.079
	λ	1.5	1.501	3.127	0.193	2.492

Table 4 :

 4 For each combination of µ, σ and γ, we compute 100 samples of SGEV(µ,σ,γ,λ = 1.5) and calculate the PWM and ML estimations, together with the respective mean and standard deviation (sd).

	Model Method	µ	σ	γ	λ
	SGEV	ML	119.41 50.22 0.35 -0.15
		PWM 107.23 49.99 0.30 0.05
	GEV	ML	111.31 47.39 0.35	-
		PWM 112.47 50.57 0.27	-

Table 5 :

 5 Estimations of µ, σ, γ and λ by ML and PWM methods for the Fort Collins data set, obtained by fitting a SGEV or a GEV distribution.

Acknowledgements

This work has been partially supported by the MIRACCLE-GICC project, the French Research National Agency (ANR) under the reference ANR-08-BLAN-0314-01, the EU-FP7 "ACQWA" Project (www.acqwa.ch) under contract Nr 212250, the ANR-MOPERA project, the PEPER-GIS project, the ANR-McSim project and the research chair Actuariat responsable sponsored by Generali.

realisations by a N (0, 1) distribution, of the estimations of a 1 -p quantile obtained by fitting a GEV (grey line) or a SGEV (dark line) with 95% empirical confidence bounds, as a function of the block size (x abscisse). The horizontal line represents the true quantile and the red line the quantile obtained by the block maxima method, z BM p . Bottom row: estimation of the skew parameter λ; p = 1/500. Right hand : Same as left hand but based on 500 samples of 5000 realizations of a N (0, 1) distribution and estimating a p quantile with p = 1/5000.


Christophe Profeta

Thomas Simon

Windings of the stable Kolmogorov process

Keywords: 2010 Mathematics Subject Classification. 60F99, 60G52, 60J50 Integrated process -Harmonic measure -Hitting time -Stable Lévy process -Winding

come

Introduction and statement of the results

A celebrated theorem by F. Spitzer [START_REF] Spitzer | Some theorems concerning 2-dimensional Brownian motion[END_REF] states that the angular part {ω(t), t ≥ 0} of a two-dimensional Brownian motion starting away from the origin satisfies the following limit theorem 2 ω(t)

log t d -→ C as t → +∞,
where C denotes the standard Cauchy law. An analogue of this result for isotropic stable Lévy processes was given in [START_REF] Bertoin | Stable windings[END_REF], with a slower speed in √ log t and a centered Gaussian limit law. Notice that both these results can be obtained as functional limit theorems with respect to the Skorohod topology. We refer to [START_REF] Doney | Windings of planar stable processes[END_REF] for a recent paper revisiting these problems, with further results and an updated bibliography.

In a different direction, McKean [START_REF] Mckean | A winding problem for a resonator driven by a white noise[END_REF] had observed that the windings of the Kolmogorov diffusion, which is the two-dimensional process Z having a linear Brownian motion as second coordinate and its running integral as first coordinate, obey an almost sure limit theorem. More precisely, if {ω(t), t ≥ 0} denotes the angular part of the process Z starting away from the origin, it is shown in Section 4.3 of [START_REF] Mckean | A winding problem for a resonator driven by a white noise[END_REF] that

ω(t) log t a.s. -→ - √ 3 2 as t → +∞
(the constant which is given in [START_REF] Mckean | A winding problem for a resonator driven by a white noise[END_REF] is actually -√ 3/8, but it will be observed below that the evaluation of the relevant improper integral in [START_REF] Mckean | A winding problem for a resonator driven by a white noise[END_REF] was slightly erroneous). Of course, the degeneracy of the Kolmogorov diffusion makes it wind in a very particular way, since this process visits a.s. alternatively and clockwise the left and right half-planes. The regularity of this behaviour, which contrasts sharply with the complexity of planar Brownian motion, makes it possible to use the law of large numbers and to get an almost sure limit theorem.

The first aim of this paper is to obtain an analogue of McKean's result in replacing Brownian motion by a strictly α-stable Lévy process L = {L t , t ≥ 0}. Without loss of generality, we choose the following normalization for the characteristic exponent

Ψ(λ) = log(E[e iλL 1 ]) = -(iλ) α e -iπαρ sgn(λ) , λ ∈ R, (1.1)
where α ∈ (0, 2] is the self-similarity parameter and ρ = P[L 1 ≥ 0] is the positivity parameter. We refer to [START_REF] Samorodnitsky | Stable Non-Gaussian random processes[END_REF][START_REF] Zolotarev | One-dimensional stable distributions[END_REF] for accounts on stable laws and processes, and to the introduction of our previous paper [START_REF] Profeta | Persistence of integrated stable processes[END_REF] for a discussion on this specific parametrization. Recall that if α = 2, then necessarily ρ = 1/2 and L = { √ 2B t , t ≥ 0} is a rescaled Brownian motion. Introduce the primitive process

X t = t 0 L s ds, t ≥ 0,
and denote by P (x,y) the law of the strong Markov process Z = (X, L) started from (x, y). By analogy with the classical Kolmogorov diffusion [START_REF] Kolmogorov | Zufällige Bewegungen (zur Theorie der Brownschen Bewegung)[END_REF], this process may and will be called the stable Kolmogorov process. When (x, y) = (0, 0), it can be shown without much difficulty -see Lemma 3 below -that under P (x,y) , the process Z never hits (0, 0). Filling in the gaps made by the jumps of L by vertical lines -see the figure below -and reasoning exactly as in [START_REF] Bertoin | Stable windings[END_REF] p.1270 it is possible to define the algebraic angle

ω(t) = (Z 0 , Z t )
measured in the trigonometric orientation.

Figure 1. One path of ((X t , L t ), t ≤ 100) starting from X 0 = -1 and L 0 = 0.

If ρ = 1 resp. ρ = 0, then |L| is a stable subordinator and it is easy to see that Z stays for large times within the positive resp. the negative quadrant with a.s. X t /L t → +∞, so that ω(t) converges a.s. to a finite limit which is

(Z 0 , Ox) resp. (-Z 0 , Ox).
When ρ ∈ (0, 1) and (x, y) = (0, 0), the Lévy process L oscillates and the Kolmogorov process Z winds clockwise and infinitely often around the origin as soon as (x, y) = (0, 0). Indeed, considering the partition R 2 \{(0, 0)} = P -∪ P + with

P -= {x < 0} ∪ {x = 0, y < 0} and P + = {x > 0} ∪ {x = 0, y > 0},
we see that if (x, y) ∈ P -the continuous process X visits alternatively the negative and positive half-lines, starting negative, and that its speed when it hits zero is alternatively positive and negative, starting positive. When (x, y) ∈ P + the same alternating scheme occurs, with opposite signs. In particular, the function ω(t) is a.s. negative for all t large enough. In order to state our first result, which computes the a.s. limit velocity of ω(t), let us finally introduce the parameters

γ = ρα 1 + α ∈ (0, 1/2) and γ = (1 -ρ)α 1 + α ∈ (0, 1/2).
Theorem A. Assume ρ ∈ (0, 1) and (x, y) = (0, 0). Then, under P (x,y) , one has

ω(t) log t a.s.
-→ -2 sin(πγ) sin(πγ) α sin(π(γ + γ)) as t → +∞.

Note that in the Brownian motion case α = 2, we have ρ = 1/2 and γ = γ = 1/3, so that 2 sin(πγ) sin(πγ)

α sin(π(γ + γ)) = √ 3 2 •
The constant -√ 3/8 which is given in Section 4.3 of [START_REF] Mckean | A winding problem for a resonator driven by a white noise[END_REF] is not the right one because of the erroneous evaluation of the integral in (3.8.a) therein: this integral equals actually π/ √ 3, as can be checked by an appropriate contour integration. The proof of Theorem A goes basically along the same lines as in [START_REF] Mckean | A winding problem for a resonator driven by a white noise[END_REF]. We consider the successive hitting times of 0 for the integrated process X :

T (1) 0 = T 0 = inf{t > 0, X t = 0}
and

T (n) 0 = inf{t > T (n-1) 0 , X t = 0},
which can be viewed as the half-winding times of Z. We first check that n → T (n) 0 increases a.s. to +∞ as soon as (x, y) = (0, 0). The exact exponential rate of escape of T (n) 0 , which yields the exact winding velocity, is computed thanks to an elementary large deviation argument involving the law of L T 0 under P (x,0) , a certain transform of the half-Cauchy distribution as observed in [START_REF] Profeta | Persistence of integrated stable processes[END_REF]. Notice that contrary to [START_REF] Mckean | A winding problem for a resonator driven by a white noise[END_REF] where the proof is only sketched, we provide here an argument with complete details.

In the Brownian case α = 2, an expression of the law of the bivariate random variable (T

(n) 0 , |L T (n) 0
|) under P (0,y) has been given in Theorem 1 of [START_REF]Les temps de passage successifs de l'intégrale du mouvement brownien[END_REF], in terms of the modified Bessel function of the first kind. This expression becomes very complicated under P (x,y) when x = 0, even for n = 1 -see Formula (2) p.4 in [START_REF]Les temps de passage successifs de l'intégrale du mouvement brownien[END_REF]. In all cases, this expression is not informative enough to evaluate the upper tails of T (n) 0 . In [START_REF] Profeta | Some limiting laws associated with the integrated Brownian motion[END_REF] it was shown that

P (x,y) [T (n) 0 ≥ t] t -1/4 (log t) n-1 as t → +∞
where, here and throughout, the notation f (t) g(t) means that there exist two constants

0 < κ 1 ≤ κ 2 < +∞ such that κ 1 f (t) ≤ g(t) ≤ κ 2 f (t) as t → +∞.
On the other hand, Theorem A in our previous paper [START_REF] Profeta | Persistence of integrated stable processes[END_REF] shows the non-trivial asymptotics

P (x,y) [T 0 > t] t -θ as t → +∞ for (x, y) ∈ P -, with θ = ρ/(1 + α(1 -ρ))
. By symmetry, the latter result also shows that

P (x,y) [T 0 > t] t -θ as t → +∞ for (x, y) ∈ P + , with θ = (1 -ρ)/(1 + αρ).
The second main result of this paper connects the two above estimates.

Theorem B. Assume that ρ ∈ (0, 1) and (x, y) ∈ P -. For every n ≥ 2, the following asymptotics hold as t → +∞ :

         P (x,y) [T (n) 0 > t] t -θ (log t) [ n-1 2 ] if ρ < 1/2, P (x,y) [T (n) 0 > t] t -θ (log t) [ n 2 ]-1 if ρ > 1/2, P (x,y) [T (n) 0 > t] t -θ (log t) n-1 if ρ = 1/2.
By symmetry, the same result holds for (x, y) ∈ P + with θ and θ switched. In the above statement, the separation of cases is intuitively clear, since for ρ < 1/2 resp. ρ > 1/2 the negative resp. positive excursions below 0 will tend to prevail. The main difficulty in the proof of Theorem B is to show that the required asymptotic behaviour does not depend on the starting point (x, y) ∈ P -. This is handled thanks to a uniform estimate on the Mellin transform of the harmonic measure P (x,y) [L T 0 ∈ .], and a general estimate on the upper tails of the product of two positive independent random variables. These two estimates have both independent interest.

Proofs

Preliminary results.

As mentioned before, we first establish some estimates on the harmonic measure of the left half-plane with respect to the stable Kolmogorov process. When starting from (x, y) ∈ P -, the process (X, L) ends up in exiting P -on the positive vertical axis, and its exit distribution is given by the law of L T 0 under P (x,y) . This distribution is called the harmonic measure since by the generalized Poisson formula, it allows to construct harmonic functions with respect to the degenerate operator L α,ρ y + y ∂ ∂x on the half-plane, where L α,ρ y is the generator of the stable Lévy process L. However, we shall not pursue these lines of research here. Lemma 1. Assume that (x, y) ∈ P -. The Mellin transform s → E (x,y) [L s-1

T 0 ] is real-analytic on (1/(γ -1), 1/(1 -γ)), with two simple poles at 1/(γ -1) and 1/(1 -γ). In particular, the random variable L T 0 has a smooth density f 0

x,y under P (x,y) , and there exist c 1 , c 2 > 0 such that

f 0 x,y (z) ∼ z→0 c 1 z αθ/γ and f 0 x,y (z) ∼ z→+∞ c 2 z -αθ-1 .
Proof. Observe first that the smoothness and the asymptotic behaviour of the density function of L T 0 are a direct consequence of the statement on the Mellin transform, thanks to the converse mapping theorem stated e.g. as Theorem 4 in [START_REF] Flajolet | Mellin transforms and asymptotics: Harmonic sums[END_REF]. This latter statement is also a direct consequence of Theorem B in [START_REF] Profeta | Persistence of integrated stable processes[END_REF] when either x = 0 or y = 0. From now on we shall therefore assume that xy = 0. By Proposition 2 (i) and Equation (3.2) in [START_REF] Profeta | Persistence of integrated stable processes[END_REF] we have

E (x,y) L s-1 T 0 = π +∞ 0 E (x,y) X -ν t 1 {Xt>0} dt (1 + α) 1-ν (Γ (1 -ν)) 2 Γ(1 -s) sin(πs(1 -γ)) (2.1) with s = (1 -ν)(1 + α) ∈ (0, 1
). However, it does not seem easy to study the poles of the right-hand side directly since the integral is not expressed in closed form, and for this reason we shall perform a further Mellin transformation in space. First, we know from Proposition 1 in [START_REF] Profeta | Persistence of integrated stable processes[END_REF] that

E (x,y) [X -ν t 1 {Xt>0} ] = Γ(1 -ν) π ∞ 0
λ ν-1 e -cα,ρλ α t α+1 sin(λ(x + yt) + s α,ρ λ α t α+1 + πν/2) dλ for every ν ∈ (0, 1), with

s α,ρ = sin(πα(ρ -1/2)) α + 1 ∈ (-1, 1) and c α,ρ = cos(πα(ρ -1/2)) α + 1 ∈ (0, 1).
For every β ∈ (0, ν) this yields

0 -∞ |x| β-1 E (x,y) [X -ν t 1 {Xt>0} ]dx = Γ(1 -ν) π +∞ 0 λ ν-1 e -cα,ρλ α t α+1 0 -∞ |x| β-1 sin(λ(x + yt) + s α,ρ λ α t α+1 + πν/2) dx dλ = Γ(1 -ν)Γ(β) π +∞ 0 λ ν-β-1 e -cα,ρλ α t α+1 sin(λyt + s α,ρ λ α t α+1 + π(ν -β)/2) dλ = Γ(1 -ν)Γ(β) Γ(1 -ν + β) E (0,y) X -ν+β t 1 {Xt>0}
where the switching of the first equality is justified exactly as in Lemma 1 of [START_REF] Profeta | Persistence of integrated stable processes[END_REF], and the second equality follows from trigonometry and generalized Fresnel integrals -see (2.1) and

(2.2) in [START_REF] Profeta | Persistence of integrated stable processes[END_REF]. Assume first that y < 0. From Proposition 2 (ii) in [START_REF] Profeta | Persistence of integrated stable processes[END_REF], we obtain

+∞ 0 0 -∞ |x| β-1 E (x,y) [X -ν t 1 {Xt>0} ] dx dt = Γ(1 -ν)Γ(β)Γ(1 -ν + β) π (α+1) 1-ν+β Γ(1-s-β(1+α)) sin(πραβ +πγs) |y| s+β(1+α)-1
for every β ∈ (0, (1 -s)/(α + 1)). Putting this together with (2.1), we finally deduce

0 -∞ |x| β-1 E (x,y) L s-1 T 0 dx = (α + 1) β Γ(β) sin(πραβ + πγs)Γ(1 -ν + β)Γ(1 -s -β(1 + α)) Γ(1 -ν)Γ(1 -s) sin(πs(1 -γ)) |y| s+β(1+α)-1 . (2.2)
We shall now invert this Mellin transform in the variable β in order to get a suitable integral expression for

E (x,y) [L s-1 T 0 ]. Fix β ∈ (0, (1 -s)/(α + 1)
). On the one hand, since ρα < 1, we have

Γ(β) cos((πραβ + πγs)/2) = ∞ 0 x β-1 e -x cos(πρα/2) cos(x sin(πρα/2) + πγs/2) dx and Γ(1 -ν + β) sin((πραβ + πγs)/2) = Γ(1 -ν + β) sin πρα 2 (β + 1 -ν) = ∞ 0
x β-1 x 1-ν e -x cos(πρα/2) sin(x sin(πρα/2)) dx.

On the other hand, a change of variable in the definition of the Gamma function shows that

Γ(1 -s -β(1 + α)) = 1 1 + α +∞ 0 x β-1 x s-1 1+α e -x -1/(1+α) dx. Setting K s (ξ) = +∞ 0 z s 1+α -1 e -cos(πρα/2)(ξ/z+z) cos ξ z sin(πρα/2) + πγs/2 sin(z sin(πρα/2)) dz,
we can now invert (2.2) and obtain, applying Fubini's theorem and using the notation x = x(1 + α)|y| 1+α , a new expression for the Mellin transform of L T 0 :

E ( x,y) L s-1 T 0 = 2|y| s-1 |x| s-1 1+α (1 -s) Γ(s/(1 + α))Γ(2 -s) sin(πs(1 -γ)) +∞ 0 ξ 1-s 1+α -1 e -( ξ |x| ) 1/(1+α) K s (ξ) dξ.
Since 1 1-γ = 1 + αρ α+(1-αρ) < 1 + ρ < 2, it remains to prove that the function (cos(πρα/2))

H s (x) = (1 -s) ∞ 0 ξ 1-s 1+α -1 e -( ξ |x| ) 1/(1+α) K s (ξ) dξ
s 1+α +1 Γ s 1 + α + 1 .
As a consequence, the function H s has an analytic continuation on (-1-α, 1) ⊃ [1/(γ -1), 1) since 1 γ-1 = -1+α 1+α(1-ρ) > -1 -α. Next for every s ∈ (0, 1) an integration by parts shows that

H s (x) = (1 + α) ∞ 0 ξ 1-s 1+α d dξ e -( ξ |x| ) 1/(1+α) K s (ξ) dξ = (1 + α) ∞ 0 ξ 1-s 1+α e -( ξ |x| ) 1 1+α K s (ξ) dξ - 1 + α |x| 1/(1+α) +∞ 0 ξ 1-s-α 1+α e -( ξ |x| ) 1 1+α K s (ξ) dξ,
where K s is well-defined on [0, +∞) for any s > 0, and bounded by

|K s (ξ)| ≤ 2 sin(πρα/2) ∞ 0 z s 1+α -1 e -cos(πρα/2)z dz = 2 sin(πρα/2) (cos(πρα/2)) s/(1+α) Γ s 1 + α .
Consequently, the function H s also admits an analytic continuation on (0, 2)

⊃ (0, 1/(1 -γ)].
This completes the proof in the case y < 0. The case y > 0 may be dealt with in an entirely similar way, and we leave the details to the reader.

Our second preliminary result is elementary, but we could not find any reference in the literature and we hence provide a proof. Lemma 2. Let µ ≥ ν > 0 and n, p ∈ N. Assume that X and Y are two independent positive random variables such that :

P[X ≥ z] z→+∞ z -ν (log z) n and P[Y ≥ z] z→+∞ z -µ (log z) p . Then    P[XY ≥ z] z→+∞ z -ν (log z) n+p+1 if µ = ν, P[XY ≥ z] z→+∞ z -ν (log z) n if µ > ν.
Proof. We first decompose the product as

P[XY ≥ z] = ∞ 0 P[X ≥ zy -1 ] P[Y ∈ dy].
Therefore, for z > A large enough,

P[XY ≥ z] ≥ √ z A P[X ≥ zy -1 ] P[Y ∈ dy] ≥ κ 1 z ν √ z A y ν (log(zy -1 )) n P[Y ∈ dy] ≥ κ 1 (log(z)) n 2 n z ν √ z A y ν P[Y ∈ dy].
Then, integrating by parts,

√ z A y ν P[Y ∈ dy] = A ν P[Y ≥ A] -z ν/2 P[Y ≥ √ z] + ν √ z A y ν-1 P[Y ≥ y]dy.
Now, if ν < µ, this expression remains bounded as z → +∞. Assume therefore that µ = ν.

In this case, we have :

√ z A y ν P[Y ∈ dy] ≥ A ν P[Y ≥ A] - κ 2 2 p (log z) p + νκ 1 √ z A (log y) p y dy ∼ z→+∞ νκ 1 (log z) p+1 2 p+1 (p + 1) ,
which gives the lower bound. To obtain the upper bound, we separate the integral in three parts and proceed similarly, with ε small enough:

1/ε 0 P[X ≥ zy -1 ]P[Y ∈ dy] + εz 1/ε P[X ≥ zy -1 ]P[Y ∈ dy] + ∞ εz P[X ≥ zy -1 ]P[Y ∈ dy] ≤ P[X ≥ εz] + κ 2 (log(εz)) n z ν εz 1/ε y ν P[Y ∈ dy] + P[Y ≥ εz]
and the proof is concluded as before, using an integration by parts and looking separately at both cases ν < µ and ν = µ.

Proof of Theorem A.

By symmetry, it is enough to show Theorem A for (x, y) ∈ P -.

Consider the sequence

T (n) 0 , |L T (n) 0 | n≥1
and set {F n , n ≥ 1} for its natural completed filtration. It is easy to see from the strong Markov and scaling properties of Z that this sequence is Markovian. To be more precise, starting from P -and taking into account the possible asymmetry of the process L, we have the following identities for all p ≥ 1.

T (2p) 0 , |L T (2p) 0 | d = T (2p-1) 0 + |L T (2p-1) 0 | α τ + , |L T (2p-1) 0 | + with (τ + , + ) ⊥ F 2p-1 distributed as (T 0 , |L T 0 |) under P (0,1) , and
T (2p+1) 0 , |L T (2p+1) 0 | d = T (2p) 0 + |L T (2p) 0 | α τ -, |L T (2p) 0 | - with (τ -, -) ⊥ F 2p distributed as (T 0 , |L T 0 |) under P (0,-1)
. The starting term (T 0 , |L T 0 |) has the same law as (τ -, -) if x = 0 and y = -1. By induction we deduce the identities

|L T (2p) 0 | d = |L T 0 | × p-1 k=1 - k × p k=1 + k and |L T (2p+1) 0 | d = |L T 0 | × p k=1 - k × p k=1 + k
where, here and throughout, (τ ± k , ± k ) k≥1 are two i.i.d. sequences distributed as (τ ± , ± ), and all products are assumed independent. From Theorem B (i) in [START_REF] Profeta | Persistence of integrated stable processes[END_REF] and its symmetric version, the Mellin tranforms of ± are given by

E ( -) s-1 = sin(πγs) sin(π(1 -γ)s) and E ( + ) s-1 = sin(πγs) sin(π(1 -γ)s) (2.3)
for each real s in the respective domain of definition, which is in both cases an open interval containing 1. This entails that

E[| log( ± )|] < +∞, with E log( -) = π cot(πγ) > 0 and E log( + ) = π cot(πγ) > 0. (2.4)
The following lemma is intuitively obvious.

Lemma 3. Assume (x, y) ∈ P -. Then one has T

(n) 0

→ +∞ and the process Z never hits the origin, a.s. under P (x,y) .

Proof. To prove the first statement, it is enough to show that

S n = T (n) 0 -T (n-1) 0 → +∞ a.s. as n → ∞. Set κ α,ρ = πα 2 (cot(πγ) + cot(πγ)) = πα sin(π(γ + γ)) 2 sin(πγ) sin(πγ) > 0.
From the above discussion, we have

S 2p d = |L T 0 | α × τ + × p-1 k=1 - k × + k α (2.5)
for every p ≥ 2, with independent products on the right-hand side. For every ε ∈ (0, κ α,ρ ), this entails P (x,y) S 2p ≤ e 2(p-1)(κα,ρ-ε) ≤ P (x,y) |L T 0 | ≤ e -ε(p-1)/2α + P τ + ≤ e -ε(p-1)/2

+ P 1 p -1 p-1 k=1 log( - k ) ≤ π cot(πγ) -ε/2 + P 1 p -1 p-1 k=1 log( + k ) ≤ π cot(πγ) -ε/2 .
From Lemma 1, there exists θ 1 (ε) > 0 such that P (x,y) |L T 0 | ≤ e -ε(p-1)/2α < e -pθ 1 (ε) for p large enough. On the other hand, we have P τ + ≤ e -ε(p-1)/2 ≤ P (0,1) inf{L t , t ≤ e -ε(p-1)/2 } < 0 = P (0,0) sup{ Lt , t ≤ 1} > e ε(p-1)/2α ≤ e -pθ 2 (ε)

for some θ 2 (ε) > 0 and all p large enough, where we have set L = -L, the equality following from translation invariance and self-similarity, and the second inequality from the general estimate of Theorem 12.6.1 in [START_REF] Samorodnitsky | Stable Non-Gaussian random processes[END_REF]. Last, the existence of some θ 3 (ε) > 0 such that both remaining terms can be bounded from above by e -pθ 3 (ε) for p large enough is a standard consequence of (2.3), (2.4) and Cramér's theorem -see e.g. Theorem 1.4 in [START_REF] Hollander | Large deviations. Fields institute monographs[END_REF], recalling that the assumption (I.5) can be replaced by (I.17) therein. We can finally appeal to the Borel-Cantelli lemma to deduce, having let ε → 0, lim inf p→∞ 1 2p log(S 2p ) ≥ κ α,ρ > 0 a.s. (2.6) This shows that S 2p → +∞ a.s. and an entirely similar argument yields S 2p+1 → +∞ a.s. This concludes the proof of the first part of the lemma. The second part is easier. If α ≤ 1, it is well-known that L never hits zero, so that Z never hits the origin. If α > 1, we see from Lemma 1 that L T 0 has no atom at zero under P (x,y) and because ± are absolutely continuous, all L T (n) 0 's have no atom at zero. We finally get P (x,y) [Z visits the origin] = P (x,y)

n≥1 L T (n) 0 = 0 = 0
where the first identification comes from the fact that T

(n) 0

→ +∞ a.s.

We can now finish the proof of Theorem A. Set θ 0 = Z 0 Z T 0 ∈ (-π, 0) a.s. Observing as in [START_REF] Mckean | A winding problem for a resonator driven by a white noise[END_REF] the a.s. identifications

{ω(t) ≥ -(n -1)π + θ 0 } = {T (n) 0 ≥ t} and {ω(t) ≤ -(n -2)π + θ 0 } = {T (n-1) 0 ≤ t}, we see that Theorem A amounts to show that 1 n log(T (n) 0 )
a.s.

-→ πα sin(π(γ + γ)) 2 sin(πγ) sin(πγ) = κ α,ρ as n → +∞.

Firstly, with the above notation, we have a.s. under P (x,y)

lim inf n→∞ 1 n log(T (n) 0 ) ≥ lim inf n→∞ 1 n log(S n ) ≥ κ α,ρ ,
where the second inequality comes from (2.6) and its analogue for n odd. To obtain the upper bound, we will proceed as in the above Lemma 3. Fixing ε > 0, we have

P (x,y) T (n) 0 ≥ e n(κα,ρ+ε) ≤ n k=1 P (x,y) S k ≥ n -1 e n(κα,ρ+ε) ≤ n k=1 P (x,y) S k ≥ e n(κα,ρ+ε/2)
for n large enough, with the above notation for S k and having set S 1 = T 0 . Recalling (2.5) we have for every k = 2p ≤ n P (x,y) S 2p ≥ e n(κα,ρ+ε/2) ≤ P (x,y) |L T 0 | ≥ e nε/8α + P τ + ≥ e nε/8

+ P 2 n p-1 k=1 log( - k ) ≥ π cot(πγ) + ε/4 + P 2 n p-1 k=1 log( + k ) ≥ π cot(πγ) + ε/4 .
Again from Lemma 1, there exists θ 4 (ε) > 0 such that P (x,y) |L T 0 | ≥ e nε/8α < e -nθ 4 (ε) for n large enough, whereas P τ + ≥ e nε/8 ≤ P (0,0) sup{ Lt , t ≤ 1} < e -nε/8α ≤ e -nθ 5 (ε)

for some θ 5 (ε) > 0 and n large enough, the second inequality following e.g. from Proposition VIII.2 in [START_REF] Bertoin | Lévy processes[END_REF]. To handle the third term, we separate according as p ≤ √ n or p > √ n. In the first case, we have the upper bound

P 2 n p-1 k=1 log( - k ) ≥ π cot(πγ) + ε/4 ≤ √ n k=1 P 2 √ n log( - k ) ≥ π cot(πγ) + ε/4 ≤ e -θ 6 (ε) √ n
for some θ 6 (ε) > 0 and n large enough, using Lemma 1 for the second inequality. In the second case, applying Cramér's theorem exactly as in Lemma 3 gives the upper bound

P 2 n p-1 k=1 log( - k ) ≥ π cot(πγ) + ε/4 ≤ e -θ 7 (ε) √ n
for some θ 7 (ε) > 0 and n large enough. The fourth term is estimated in the same way and we finally get the existence of some θ(ε) > 0 such that P (x,y) S 2p ≥ e n(κα,ρ+ε/2) ≤ e -θ(ε) √ n for n large enough. An analogous estimate is obtained for P (x,y) S 2p+1 ≥ e n(κα,ρ+ε/2) and we can apply as usual the Borel-Cantelli lemma to show the required upper bound lim sup and to control the sum separately according as ρ < 1/2, ρ > 1/2 and ρ = 1/2. We leave the details to the reader.

  admits an analytic continuation on [1/(γ-1), 1/(1-γ)]. Observe first that for any s > -1-α, the function K s is uniformly bounded on [0, +∞) by |K s (ξ)| ≤ sin(πρα/2) cos(πρα/2)z dz = sin(πρα/2)

2 . 3 . 1 d===

 231 κ α,ρ + ε, for all ε > 0, a.s. under P (x,y) . Proof of Theorem B. Recall the decomposition T(n) 0 = S 1 + • • • + S n with S |L T 0 | α × τ + × |L T 0 | α × τ -× notation for (τ

  ± , ± ). Let us first investigate the upper tails of the distribution of each S k under P (x,y) . We know thatP (x,y) [T 0 > t] P (x,y) [|L T 0 | α > t] +∞for every k ≥ 1. Suppose then ρ < 1/2 viz. θ < θ, we obtain in a similar wayP (x,y) [S 2p > t] t -θ (log t) p-1 and P (x,y) [S 2p+1 > t] t -θ (log t) p as t → +∞, for every p ≥ 1. Last, if ρ > 1/2 we find P (x,y) [S 2p > t] t -θ (log t) p-1 and P (x,y) [S 2p+1 > t] t -θ (log t) p-1 as t → +∞, for every p ≥ 1. All in all, for all k ≥ 2, this shows thatP (x,y) [S k > t] t -θ (log t) [ k-1 2 ] if ρ < 1/2 P (x,y) [S k > t] t -θ (log t) [ k 2 ]-1 if ρ > 1/2 P (x,y) [S k > t] t -θ (log t) k-1 if ρ = 1/2, and we also know that P (x,y) [S 1 > t] t -θ . The immediate estimate P (x,y) [T P (x,y) [S n > t] yields the required lower bound. To get the upper bounds, it suffices to write P (x,y) [T

	
	   
	   
				(n) 0	> t] ≥
			n
	(n) 0	> t] ≤	P (x,y) [S k > t/n]
			k=1
			P[τ -> t]	P[( -) α > t]	t -θ
	and
	P[τ + > t]	P[( + ) α > t]	t -θ

as t → +∞. Supposing first ρ = 1/2 viz. θ = θ, a successive application of Lemma 2 shows that

P (x,y) [S k > t]

t -θ (log t) k-1 as t →

Acknowledgement. Ce travail a bénéficié d'une aide de la Chaire Marchés en Mutation, Fédération Bancaire Française.