Laurent Gardes

Samuel Maistre

Nonparametric asymptotic confidence intervals for extreme quantiles

Keywords: Extreme quantiles, confidence interval, heavy-tailed distribution. is

In this paper, we propose new asymptotic confidence intervals for extreme quantiles i.e., for quantiles located outside the range of the available data.

We restrict ourselves to the situation where the underlying distribution is heavy-tailed. While asymptotic confidence intervals are mostly constructed around a pivotal quantity, we consider here an alternative approach based on the distribution of order statistics sampled from a uniform distribution.

The convergence of the coverage probability to the nominal one is established under a classical second-order condition. The finite sample behavior is also examined and our methodology is applied to a real dataset.

Introduction

Let X be a random variable defined on a probability space (Ω, A, P). The distribution of X is entirely characterized by its cumulative distribution function F (•) := P(X ≤ •) or equivalently by its quantile function given by 1 F ← (•) = inf{x ∈ S | F (x) ≥ •}, where S is the support of the distribution. In what follows, we assume that F is a continuous and strictly increasing function so that F ← = F -1 is the inverse function of F .

The question of the estimation of a large quantile F ← (1 -β) =: Q(β) for a level β close to 0 is crucial in actuarial science or in finance in order to be able to compute risk measures such as the Value-at-Risk or the Tail-Value-at-Risk, see Embrechts et al. [START_REF] Embrechts | Modelling Extremal Events for Insurance and Finance[END_REF] and Konstantinides [START_REF] Konstantinides | Risk theory. A Heavy Tail Approach[END_REF]. Large quantiles are also useful in climate science to calculate return levels of extreme precipitation or temperature, see for instance Katz et al. [START_REF] Katz | Statistics of extreme in hydrology[END_REF] and Naveau et al. [START_REF] Naveau | Statistical methods for extreme event attribution in climate science[END_REF].

There is a large literature dedicated to quantile estimation when n independent copies X 1 , . . . , X n of X are observed. When the level β = β n is not too small in the sense that β n → 0 with nβ n → c ∈ [1, ∞], the quantile Q(β n ) is usually estimated by inverting the empirical cumulative distribution function that is to say by Q (E) n (β n ) := X n-nβn ,n where X 1,n ≤ . . . ≤ X n,n is the sample arranged in ascending order. When the level β n is eventually smaller than 1/n, i.e., when nβ n → c ∈ [0, 1), the quantile Q(β n ) is located outside the range of the sample {X 1 , . . . , X n }. This kind of quantile is referred to as extreme quantile and its estimation requires a priori knowledge of the shape of the tail-distribution. For instance, one can assume that the distribution is heavy-tailed or equivalently that the quantile function Q is a regularly varying function at 0 with index -ξ < 0 that is to say that for all t > 0,

lim α→0 Q(tα) Q(α) = t -ξ .
Under this assumption, Weissman [START_REF] Weissman | Estimation of parameters and large quantiles based on the k largest observations[END_REF] proposed to estimate the quantile Q(β n ) by the statistic

Q (W) n (β n ) = Q (E) n (α n ) α n β n ξn ,
where α n is a sequence converging to 0 and such that nα n → ∞ and ξ n is a consistent estimator of the tail-index ξ.

Despite its practical interest, little attention has been paid to the construction of asymptotic confidence intervals for Q(β n ) when β n → 0. Recall that for a preselected probability γ ∈ (0, 1), an asymptotic confidence interval for Q(β n )

where F beta (•; a, b) is the distribution function of a Beta distribution with parameters a > 0 and b > 0. Since β n ∈ (0, 1) is a sequence converging to 0 with nβ n → ∞, the sets S inf (β n , n, δ) and S sup (β n , n, δ) are non empty for all δ ∈ (0, 1) (see Lemma 1, i)) allowing us to introduce the integers

L γ (β n ) := max S inf β n , n, 1 -γ 2 and R γ (β n ) := min S sup β n , n, 1 -γ 2 .
The confidence interval proposed in Gardes [START_REF] Gardes | Nonparametric confidence interval for conditional quantiles with large-dimensional covariates[END_REF] is given by [X Lγ (βn),n , X Rγ (βn),n ] and we have the following result.

Proposition 1 Let β n ∈ (0, 1) be a sequence converging to 0 with nβ n → ∞.

If, as n → ∞,

τ n := ln(n) (nβ n ) 1/2 → 0, one has P [X Lγ (βn),n , X Rγ (βn),n ] Q(β n ) = γ + o(τ n ).
Unfortunately, this confidence interval cannot be used for an extreme quantile.

Indeed, if nβ n → c ∈ [0, 1), lim n→∞ ln(n) (nβ n ) 1/2 = ∞,
and thus Proposition 1 does not hold for the extreme quantile Q(β n ). Moreover, for n large enough, S sup (β n , n, (1 -γ)/2) = ∞ for all γ > 1 -2 exp(-c), see Lemma 1,ii). It thus appears that an extrapolation of the confidence interval bounds beyond the range of the data is required to deal with extreme quantiles.

The aim of this paper is to adapt the construction procedure of quantile confidence intervals introduced in Gardes [START_REF] Gardes | Nonparametric confidence interval for conditional quantiles with large-dimensional covariates[END_REF] to the extreme case ). The distribution is assumed to be heavy-tailed and the extrapolation method is inspired by the construction of Weissman's estimator.

nβ n → c ∈ [0, 1
The rest of the paper is organized as follows. The new confidence interval for an extreme quantile is given in Section 2 and the convergence of its coverage probability to γ is established. A procedure to reduce the bias of the coverage probability is described in Section 3. The finite sample performance is examined in Section 4 and an application to a real dataset is presented in Section 5. All the proofs are gathered in Appendix.

Construction of the confidence interval

The following model is considered for the distribution of the random variable X of interest.

(A.1) There exist positive constants d and ξ such that α ξ Q(α) → d as α → 0. Under (A.1), the quantile function Q is a regularly varying function at 0 with index -ξ < 0. More specificaly, Q(α) = α -ξ (α -1 ) where the function converges to d at infinity. Note that is a slowly varying function at infinity i.e., such that for all t > 0,

lim x→∞ (tx) (x) = 1. (1)
Condition (A.1) holds for a large number of heavy-tailed distributions such as Fréchet, Burr, student, among others.

Given an estimator ξ n of the tail-index ξ, the proposed confidence interval for the extreme quantile Q(β n ) is defined below.

Definition 1 Let X 1 , . . . , X n be n independent random variables with common quantile function Q satisfying (A.1). For a preselected probability γ ∈ (0, 1), a sequence β n such that nβ n → c ∈ [0, 1) and a sequence α n converging to 0 with nα n → ∞, the two-sided confidence interval for the extreme quantile

Q(β n ) is I(γ, β n , ξ n ) := X n-nαn ,n (t n,L (γ)) ξn , X n-nαn ,n (t n,R (γ)) ξn ,
where, denoting by F ← beta (•; a, b) the quantile function of a Beta distribution of parameters a > 0 and b > 0,

t n,L (γ) = β -1 n F ← beta 1 -γ 2 , nα n + 1, n -nα n and t n,R (γ) = β -1 n F ← beta 1 + γ 2 , nα n + 1, n -nα n .
The idea leading to this confidence interval is given hereafter. We start by remarking that 1-F (X n-nαn ,n ) follows a Beta distribution with parameters nα n + 1 and n -nα n . Consequently,

P X n-nαn ,n (t n,L (γ)) ξ > Q(β n ) = F beta (z n (γ), nα n + 1, n -nα n ) , ( 2
)
where [START_REF] De Haan | Extreme Value Theory: An Introduction[END_REF] in the proof of Theorem 1), one has that

z n (γ) = 1 -F ((t n,L (γ)) -ξ Q(β n )). Under (A.1), x 1/ξ (1 -F (x)) → d 1/ξ as x → ∞. Replacing x by Q(u), we obtain that u ξ Q(u) → d as u → 0. Since t n,L (γ) ∼ α n /β n (see equation
(t n,L (γ)) -ξ Q(β n ) ∼ dα -ξ n → ∞. Hence, as n → ∞, z n (γ) ∼ d 1/ξ t n,L (γ)[Q(β n )] -1/ξ ∼ t n,L (γ)β n = F ← beta 1 -γ 2 , nα n + 1, n -nα n . Replacing z n (γ) in (2) by F ← beta 1-γ 2 , nα n + 1, n -nα n , we finally obtain P X n-nαn ,n (t n,L (γ)) ξ ≤ Q(β n ) ≈ 1 - 1 -γ 2 = 1 + γ 2 . ( 3
)
A similar approximation can be obtained for t n,R (γ). Replacing ξ by a consistent estimator conducts us to the definition of the confidence interval I(γ, β n , ξ n ). It is worth noting that even if z n (γ) ∼ α n , the sequence F beta (z n (γ), nα n + 1, n -nα n ) still depends on γ when n goes to infinity.

We will prove hereafter that the approximation (3) is theoretically justified and more specifically that the coverage probability of the random interval I(γ, β n , ξ n ) converges to γ. To establish this result, one needs to control the rate of convergence in (A.1).

(A.2) There exist a function A converging to 0 at infinity and ρ < 0 such that lim

α→0 1 A(α -1 ) Q(tα) Q(α) -t -ξ = t -ξ t -1 1 u ρ-1 du. Condition (A.
→ 0 as n → ∞ such that σ -1 n ξ n -ξ ξ d -→ D. If (nα n ) 1/2 max A(α -1 n ); σ n ln(α n /β n ) → 0, then lim n→∞ P(I(γ, β n , ξ n ) Q(β n )) = γ.
Note that under the conditions of Theorem 1,

P X n-nαn ,n (t n,L (γ)) ξn ≤ Q(β n ) → 1 + γ 2 .
Hence, as a by-product of Theorem 1, we obtain a one-sided confidence interval for Q(β n ) given by X n-nαn ,n (t n,L (2γ -1)) ξn , ∞ .

To estimate ξ, the natural way is to use the classic Hill estimator [START_REF] Hill | A simple general approach to inference about the tail of a distribution[END_REF] defined for a sequence k n ∈ {1, . . . , , n} with k n → ∞ and k n /n → 0, by

ξ (H) n (k n ) := 1 k n kn i=1 log X n-i+1,n X n-kn,n . ( 4
)
According to [16, Theorem 3.2.5], under (A.2) and if k

1/2 n A(n/k n ) → 0, k 1/2 n ξ (H) n (k n ) -ξ ξ d -→ N (0, 1) . ( 5
)
We thus have the following result which is a direct consequence of Theorem 1.

Corollary 1 Assume that conditions (A.1) and (A.2) hold. Let α n ∈ (0, 1), k n ∈ {1, . . . , n} and β n ∈ (0, 1) be sequences such that

α n , k n /n → 0, k n , nα n → ∞ and nβ n → c ∈ [0, 1). If (nα n ) 1/2 max A(α -1 n ); (k n ) -1/2 ln(α n /β n ) → 0 and k 1/2 n A(n/k n ) → 0 then lim n→∞ P(I(γ, β n , ξ (H) n (k n )) Q(β n )) = γ.
The choice of the sample fraction k n in the Hill estimator is a difficult task

and remains an open problem. Several procedures of selection can be found in the literature. Dekkers et de Haan [START_REF] Dekkers | Optimal choice of sample fraction in extreme-value estimation[END_REF] propose to minimize the mean squared-

error of ξ (H) n (k n ) -ξ.
Adaptation of boostrap methods are investigated by Danielsson et al. [START_REF] Danielsson | Using a bootstrap method to choose the sample fraction in tail index estimation[END_REF] and Gomes and Oliveira [START_REF] Gomes | The bootstrap methodology in statistics of extremes -Choice of the optimal sample fraction[END_REF]. Graphical based methods are considered by Sousa and Michailids [START_REF] De Sousa | A diagnostic plot for estimating the tail index of a distribution[END_REF] and Beirlant et al. [START_REF] Beirlant | Tail index estimation, Pareto quantile plots and regression diagnostics[END_REF][START_REF] Beirlant | On exponential representations of log-spacings of extreme order statistics[END_REF]. In the simulation study (see Section 4), we propose an ad-hoc procedure to select the optimal number k n of order statistics.

Finally, note that we can obtain similar results to Corollary 1 by using other tail index estimators such as for instance the bias-reduced Hill estimator proposed by Cairo et al. [START_REF] Cairo | Direct reduction of bias of the classical Hill estimator[END_REF]. Because of its simplicity, we choose to focus on the classical Hill estimator which provides good finite sample performances (see Section 4).

A bias reduction procedure

We propose a heuristic method to reduce the bias in the coverage probability of the confidence interval I(γ, β n , ξ n ). For i ∈ {1, . . . , n}, let

U i = 1 -F (X i ).
Since F is a continuous and strictly increasing function, U 1 , . . . , U n are independent uniform random variables and for all k ∈ {1, . . . , n} the order statistic U k,n follows a beta distribution with parameters k and n -k + 1. In the proof of Theorem 1, it is shown that

P X n-nαn ,n (t n,L (γ)) ξn ≤ Q(β n ) = P U nαn +1,n ≥ β n t n,L (γ)W n,L (γ) ,
and

P X n-nαn ,n (t n,R (γ)) ξn ≥ Q(β n ) = P U nαn +1,n ≤ β n t n,R (γ)W n,R (γ) , with W n,• (γ) = exp ξ n -ξ ξ ln t n,• (γ) ×    1 -F Q(β n )(t n,• (γ)) -ξn β n (t n,• (γ)) ξn/ξ    .
The construction of the confidence interval is based on the fact that, for n large enough, W n,• (γ) ≈ 1. Indeed, since

β n t n,L (γ) = F ← beta 1 -γ 2 , nα n + 1, n -nα n ,
we obtain the approximation

P X n-nαn ,n (t n,L (γ)) ξn ≤ Q(β n ) ≈ P U nαn +1,n ≥ β n t n,L (γ) = 1 + γ 2 ,
and similarly for the right-bound of the interval.

To reduce the bias in the coverage probability of the confidence interval, we propose to find a more accurate approximation w n,• (γ) of W n,• (γ). In Lemma 4, we show that under the conditions of Theorem 1,

W n,• (γ) = 1 + ξ n -ξ ξ ln t n,• (γ)(1 + o P (1)), (6)
where o P (1) is a notation for a random variable converging to 0 in probability.

It appears that the error in the approximation of W n,• (γ) by 1 is mainly due to the estimation of the tail index. In view of ( 6), we decide to approximate could consists in replacing ξ n by a bias-corrected estimator of the tail index, see e.g., Cairo et al. [START_REF] Cairo | Direct reduction of bias of the classical Hill estimator[END_REF]. This requires the estimation of the second order parameter which is a quite difficult task in practice. In this paper, we propose an alternative solution.

From ( 6), we have Proposition 2 Let ξ n and ξn be two estimators of ξ such that

W n,• (γ) ≈ 1 + sgn( ξ n -ξ)| ξ n -ξ|/ξ ln t n,• (γ)
σ -1 n ξ n -ξ ξ d -→ D and ς -1 n ξn -ξ ξ d -→ D,
for some non-degenerate distributions D and D and where σ n and ς n are two sequences converging to 0 and such that σ n /ς n → ∞. Then, the estimator

s n := ( ξ n -ξn )/| ξ n -ξn | is such that lim n→∞ P s n = sgn( ξ n -ξ) = 1.
We finally propose the following bias-corrected confidence interval.

Definition 2 Let X 1 , . . . , X n be n independent random variables with quantile function Q satisfying (A.1). For a preselected probability γ ∈ (0, 1), a sequence β n such that nβ n → c ∈ [0, 1) and a sequence α n converging to 0 with nα n → ∞, the bias-reduced confidence interval for the extreme quantile

Q(β n ) is I BR (γ, β n , ξ n ) := X n-nαn ,n (t n,L ( γ n,L )) ξn , X n-nαn ,n (t n,R ( γ n,R )) ξn ,
where γ n,• is the solution (if it exists) of the equation in u ∈ (0, 1)

t n,• (u) (1 + r s n σ n ln t n,• (u)) = t n,• (γ). (7)
Conditions for the existence and the uniqueness of the solution of ( 7) are provided in the particular case where the Hill estimator is used to estimate the tail index, see hereafter.

Application with the Hill estimator

Let us give the expression of the bias-reduced confidence interval when the Hill estimator ξ

(H) n (k n ) is used
for the estimation of the tail-index ξ. Recall that under the second-order condition (A.2) and for sequences k n and kn converging to infinity such that

k 1/2 n A(n/k n ) → 0 and k1/2 n A(n/ kn ) → 0, one has k 1/2 n ξ (H) n (k n ) -ξ ξ d -→ N (0, 1) and k1/2 n ξ (H) n ( kn ) -ξ ξ d -→ N (0, 1).
If the sequence kn is such that k n / kn → 0, one can use the result of Proposition 2 to estimate the sign of ξ

(H) n (k n ) -ξ.
This estimator is given by

s (H) n := ξ (H) n (k n ) -ξ (H) n ( kn ) | ξ (H) n (k n ) -ξ (H) n ( kn )| . ( 8
)
An example of sequences k n and kn satisfying the above conditions is given by k n = d(ln(n)) 2 for some positive constant d and kn = k n ln(ln(n)). If, for instance, we take for r the expectation of the folded normal distribution |N (0, 1)|, we obtain the bias-reduced confidence interval

X n-nαn ,n t n,L ( γ (H) n,L ) ξ (H) n (kn) , X n-nαn ,n t n,R ( γ (H) n,R ) ξ (H) n (kn) , (9)
where γ

(H) n,• is the solution (if it exists) of the equation in u ∈ (0, 1) t n,• (u) 1 + 2 k n π s (H) n ln t n,• (u) = t n,• (γ). ( 10
)
We close this section by a result dedicated to the existence and the uniqueness of the solution of [START_REF] Embrechts | Modelling Extremal Events for Insurance and Finance[END_REF].

Proposition 3 If the sequences k n , α n and β n involved in [START_REF] Embrechts | Modelling Extremal Events for Insurance and Finance[END_REF] are chosen in such a way that t n,R (0

) = t n,L (0) =: t n (0) > 1, k n > 2 π max (ln(β n )) 2 (1 -β n t n,R (γ)) 2 ; [1 + ln(β -1 n )] 2 , ( 11
)
and

k n > 2(ln(t n (0)) 2 π max 1 - t n,L (γ) t n (0) -2 ; 1 - t n,R (γ) t n (0) -2 , ( 12
)
then the solution of [START_REF] Embrechts | Modelling Extremal Events for Insurance and Finance[END_REF] exists and is unique.

As an example of sequences satisfying the conditions of Proposition 3, one can take

β n = c/n with c > 0, α n = m n /n with m n → ∞ and m n (ln(m n )) 2 = o(n)
and

k n = d max((ln(n)) 2 , m n (ln(m n ))
2 ) for some d > 2/π, see Lemma 5 for a proof. In particular, a concrete choice is given by

k n = d(ln(n)) 2 for some d > 2/π, β n = c/n with c > 0 and α n = (ln(k n )) a
/n for some a > 0. This choice for the sequence α n is considered in the simulation study, see Section 4.

Simulation study

Let X 1 , . . . , X n be n independent copies of a random variable X. We propose to look at the finite sample behavior of the confidence interval for Q(β n ) with

β n = 1/n.
Four models are considered for the quantile function of X.

• Model 1: Burr distribution. For κ > 0 and c > 0, the quantile function is given for all u ∈ (0, 1) by

Q(u) = u -1/κ -1 1/c
. The tail-index is ξ = 1/(cκ) and the second-order parameter is ρ = -1/κ.

• Model 2: Fréchet distribution. The quantile function is given for all u ∈ (0, 1) by

Q(u) = [-ln(1 -u)] -ξ .
The second-order parameter is ρ = -1.

• Model 3: Absolute Student distribution. For µ > 0, the quantile function is given for all u ∈ (0, 1) by Q(u) = qt(1 -u/2; µ), where qt(•; µ) is the quantile function of a Student distribution with µ degree of freedom.

The tail-index is ξ = 1/µ and the second-order parameter is ρ = -1.

• Model 4: Log-Gamma distribution with shape parameter a = 2 and rate parameter r = 2. The tail-index is given here by ξ = 1/r = 1/2 and the second-order parameter is ρ = 0.

Note that condition (A.2) is satisfied by models 1 to 3 but not for model 4

for which ρ = 0. For each model, we consider three sample sizes : n ∈ {200, 500, 1000}. The tail index is estimated by the Hill estimator ξ

(H) n (k n )
for some k n ∈ {1, . . . , n -1}. The sequence α n is taken such as

nα n = max 3; (ln(k n )) 0.85 .
This choice is motivated by the remark made after Proposition 3 since for k n large enough, nα n ∼ (ln(k n )) a with a = 0.85. Note that our choice ensures that nα n ≥ 3 whatever the value of k n . For the bias-corrected version of the confidence interval given in ( 9), we take s

(H)
n as in [START_REF] Danielsson | Using a bootstrap method to choose the sample fraction in tail index estimation[END_REF] with the sequence kn = min ( k n ln(ln(n)) ; n -1) .

This choice corresponds to the example given after the definition of s When the uniqueness is not guaranteed, the bias-reduced interval is replaced by the standard one, as presented in Definition 1.

The coverage probability is approximated by generating N = 10000 independent samples of size n for the four above mentioned models. This approximation is denoted γN (k n ). In the same way, we approximate the coverage probabilities We compare our intervals' performances with the "new normal" ridge confidence interval presented in Buitendag et al. [START_REF] Buitendag | Confidence intervals for extreme Pareto-type quantiles[END_REF]. This interval is based on a ridge regression estimator for the extreme value index, as developed in Buitendag et al. [START_REF] Buitendag | Ridge regression estimators for the extreme value index[END_REF], as well as on a new asymptotically normal quantity.

P X n-nαn ,n (t n,L (γ)) ξn ≤ Q(β n ) and P X n-nαn ,n (t n,R (γ)) ξn ≥ Q(β n ) ,
We choose this interval as a competitor among the 12 presented in Buitendag et al. [START_REF] Buitendag | Confidence intervals for extreme Pareto-type quantiles[END_REF] since it ensured accurate coverage in most of their simulations. This new normal ridge confidence interval is based on an estimator of the second order parameter ρ from Fraga Alves, Gomes, and de Haan [START_REF] Alves | A new class of semi-parametric estimators of the second order parameter[END_REF].

However, this estimator is not always well-defined. When such problem occurs, we replace the new normal ridge confidence interval by the new normal Hill confidence interval, which is also presented in Buitendag et al. [START_REF] Buitendag | Confidence intervals for extreme Pareto-type quantiles[END_REF] and does not rely on the estimation of ρ. The coverage comparison of our two confidence intervals (standard and bias reduced) with the new normal ridge is given in Figure 1.

Our confidence intervals provide good coverage probabilities for the distributions we considered, as we can see on the left column of Figure 1. m * . We finally take for k n the integer

k (opt) n := 1 |J (opt) m * | j∈J (opt) m * j = 1 2 max J (opt) m * + min J (opt) m * .
This procedure is applied with M = 5 and the obtained results for the twosided confidence intervals are shown in Tables 1 to 6. We also give in Tables 7 to 18 (see Appendix), the empirical coverage probabilities for the one-sided confidence intervals (left and right bound).

Since our selection procedure for k (opt) n

is not tailored for the new normal ridge estimator, we also use a "rule of thumb" with k = n/5 . We notice that our bias reduced confidence interval nearly always performs well when applied for

k (opt) n
. The worst case is attained for a sample size of n = 200 with a Student distribution where ξ = 1 and ρ = -1 for which the estimated coverage is 0.9381, see Table 2. The bias reduction is mostly beneficial.

On the contrary, the performances of the new normal ridge confidence interval are often better using n/5 , but it only outperforms our confidence intervals when n = 200 for the Student and the Fréchet distributions.

Application to a real dataset

We apply our procedure to the norwegian fire insurance claims dataset as studied in Beirlant et al. [START_REF] Beirlant | Practical Analysis of Extreme Values[END_REF]. Data have been corrected for inflation and the claims are expressed in millions of 2012 Norwegian Krone (NOK). The heavy tailness of this data set has already been validated in Gardes and Girard [START_REF] Gardes | On the estimation of the variability in the distribution tail[END_REF]. as described in Section 4. An illustration of the procedure is given in Figure 2 for year 1992, for which k

(opt) n = 195.
On Figure 3, we check that the procedure is consistent, showing the rough alignment of points

-log i/k (opt) n , log X n-i,n /X n-k (opt) n ,n , 1 ≤ i < k (opt) n , on a line of slope ξ n (k (opt) n
) passing through the origin.

The resulting bias corrected intervals are illustrated in Figure 4.

Note that the larger intervals are obtained for the higher centers of the intervals, which gives an insight on the quantile estimate. This is not surprising since a large quantile is associated with a large tail-index and a high variability in the tail of the distribution. The width of the interval provides us an indication on the tail-variability and in that sense, our results are consistent with those exposed in Gardes and Girard [START_REF] Gardes | On the estimation of the variability in the distribution tail[END_REF]. Note that the confidence in-

-log i/k (opt) n log X n-i,n /X n-k (opt) n ,n 20
Hence, if 1 -δ > 1 -exp(-c) then, for n large enough, S sup (α n , n, δ) = ∅.

In all what follows, U 1 , . . . , U n are n independent standard uniform random variables. The associated order statistics are denoted U 1,n ≤ . . . ≤ U n,n .

Lemma 2 Let α n ∈ (0, 1) be a sequence converging to 0 with nα n → ∞ as n → ∞. If there exist δ ∈ (0, 1) and k n ∈ {1, . . . , n} such that P(U kn,n ≤

α n ) → δ then k n /(nα n ) → 1.
Proof -Let us first show that necessarily, k n → ∞. Suppose that this is not the case. There exists A > 0 such that for all N ∈ N \ {0}, there exists n ≥ N for which k n < A. Now, let E 1 , . . . be independent standard exponential random variables. From Rényi's representation of ordered standard uniform random variables,

P(U kn,n ≤ α n ) = P nT kn T n+1 ≤ nα n ,
where for i ∈ N \ {0},

T i = E 1 + . . . + E i .
By the law of large numbers, n/T n+1 converges in probability to 1 and for all N ∈ N \ {0}, there exists n ≥ N such that

T kn ≤ A +1 i=1 E i .
Hence, for all N ∈ N \ {0}, there exists n ≥ N such that

P(U kn,n ≤ α n ) ≥ P   n T n+1 A +1 i=1 E i ≤ nα n   → 1, since nα n → ∞. Hence, lim sup n→∞ P(U kn,n ≤ α n ) = 1,
which is in contradiction with the fact that

P(U kn,n ≤ α n ) → δ ∈ (0, 1). Now, for all ε > 0, let A n,ε = {nU kn,n /k n ∈ [1 -ε, 1 + ε]}. One has P(U kn,n ≤ α n ) = P({U kn,n ≤ α n } ∩ A n,ε ) + P({U kn,n ≤ α n } ∩ A C n,ε ),
where A C n,ε is the complement of the set A n,ε . Since k n → ∞, Rényi's representation of the order statistic U kn,n , and the law of large number entail that nU kn,n /k n converges in probability to 1. Hence,

P({U kn,n ≤ α n } ∩ A C n,ε ) ≤ P(A C n,ε ) → 0.
Moreover,

P({U kn,n ≤ α n } ∩ A n,ε ) =            0 if nα n /k n < 1 -ε, P(A n,ε ) if nα n /k n > 1 + ε.
Since P(A n,ε ) → 1, we thus have proven that for all ε > 0,

P(U kn,n ≤ α n ) → 0 if nα n /k n < 1 -ε and P(U kn,n ≤ α n ) → 1 if nα n /k n > 1 + ε. Since by assumption, P(U kn,n ≤ α n ) → δ ∈ (0, 1
), it appears that for all ε > 0, for n large enough, 1 -ε ≤ nα n /k n ≤ 1 + ε and the proof is complete.

Let us now give a technical result which is a consequence of the second-order condition.

Lemma 3 Under conditions (A.1) and (A.2), if α n ∈ (0, 1) and β n ∈ (0, 1)

are sequences converging to 0 with α n /β n → ∞ then, for every sequence t n such that t n ∼ (α n /β n ) -ξ , one has

lim n→∞ 1 |A(α -1 n )| 1 -F (Q(β n )t n ) t -1/ξ n β n -1 = - 1 ξρ .
Proof -Condition (A.2) and [START_REF] De Haan | Extreme Value Theory: An Introduction[END_REF]Theorem 2.3.9] entail that for all t > 0,

lim x→∞ 1 A((1 -F (x)) -1 ) 1 -F (tx) 1 -F (x) -t -1/ξ = ξ -2 t -1/ξ t 1 u ρ/ξ-1 du.
We can then use [START_REF] De Haan | Extreme Value Theory: An Introduction[END_REF]Remark B.3.15] to obtain

lim x→∞ 1 A((1 -F (x)) -1 ) t 1/ξ x 1 -F (xt x ) 1 -F (x) -1 = - 1 ρξ , where t x → ∞ as x → ∞. Applying this result with x = Q(β n )t n ∼ dα -ξ n → ∞ and t x = t -1 n ∼ (α n /β n ) ξ → ∞ as n → ∞ leads to β n t -1/ξ n 1 -F (Q(β n )t n ) 1 A((1 -F (Q(β n )t n ) -1 ) × 1 -F (Q(β n )t n ) t -1/ξ n β n -1 → 1 ξρ ,
as n → ∞. From (A.1), x 1/ξ (1 -F (x)) → d 1/ξ as x → ∞ and thus, 1 -

F (Q(β n )t n ) ∼ α n .
Hence, since |A| is a regularly varying function, |A((1 -

F (Q(β n )t n ) -1 )| ∼ |A(α -1 n )| and thus lim n→∞ 1 |A(α -1 n )| 1 -F (Q(β n )t n ) t -1/ξ n β n -1 = - 1 ξρ . W n,• (γ) = exp ξ n -ξ ξ ln t n,• (γ) ×    1 -F Q(β n )(t n,• (γ)) -ξn β n (t n,• (γ)) ξn/ξ    , is established.
Lemma 4 Assume that conditions (A.1) and (A.2) hold. Let α n ∈ (0, 1) be a sequence converging to 0 with nα n → ∞ and β n ∈ (0, 1) a sequence such that nβ n → c ∈ [0, 1). Let also ξ n be an estimator of the tail-index ξ for which there exist a non-degenerate distribution D and a sequence σ n → 0 as n → ∞ Taylor expansion of the exponential function leads to

such that σ -1 n ξ n -ξ /ξ d -→ D. If (nα n ) 1/2 max A(α -1 n ); σ n ln(α n /β n ) → 0, then W n,• (γ) = 1 + ( ξ n -ξ)/ξ ln t n,• (γ)(1 + o P (1)). Proof -Recall that W n,• (γ) = F 1,n × F 2,n where F 1,n = exp(Z n ) with Z n := ln t n,• (γ)( ξ n -ξ)/ξ and F 2,n :=    1 -F Q(β n )(t n,• (γ)) -ξn β n (t n,• (γ)) ξn/ξ    .

Let us first focus on the factor

F 1,n = 1 + Z n (1 + o P (1)). (13)
Let us now deal with the factor F 2,n . Applying Lemma 3 with t n = (α n /β n ) -ξn leads to

F 2,n = 1 - 1 ξρ A(α -1 n )(1 + o(1)), ( 14
)
where o( 1) is the generic notation for a sequence converging to 0. Multiplying ( 13) and ( 14), we obtain

W n,• (γ) = 1 + Z n 1 - 1 ξρ Z -1 n A(α -1 n )(1 + o(1)) + o P (1) .
Using again the facts that ln t n,• ∼ ln(α n /β n ) and σ -1

n ( ξ n -ξ)/ξ d -→ D, one has that Z -1 n = O P (σ -1 n (ln(α n /β n )) -1 ).
Remarking that

A(α -1 n )/σ n ln(α n /β n ) → 0, we have Z -1 n A(α -1 n ) = o P (1)
and the proof is complete.

In the last lemma, we propose an example of sequences α n , β n and k n satisfying the conditions of Proposition 3 and that can be used to construct the bias-reduced confidence interval [START_REF] Dekkers | Optimal choice of sample fraction in extreme-value estimation[END_REF].

n = d max((ln(n)) 2 , m n (ln(m n )) 2 ) for some d > 2/π. Proof -First note that α n → 0, nα n = m n → ∞ and that β n /α n → ∞. Since t n (0) ∼ α n /β n → ∞, one has for n large enough that t n (0) > 1. Next, since ln(β n ) ∼ ln(n), [1 -ln(β n )] 2 ∼ (ln(n)) 2 and 1 -β n t n,R (γ) → 1, it is readily seen that condition (11) is satisfied since k n max (ln(β n )) 2 (1 -β n t n,R (γ)) 2 ; [1 + ln(β -1 n )] 2 ∼ k n (ln(n)) 2 ≥ d > 2 π .
Let us now focus on condition [START_REF] Alves | A new class of semi-parametric estimators of the second order parameter[END_REF]. One has that ln

t n (0) ∼ ln(α n /β n ) = ln(m n ) -ln(c) ∼ ln(m n ).
Moreover, a Taylor expansion together with the use of [13, Lemma 1] entail that for n large enough, there are constants 0 < c 1 < c 2 such that

m 1/2 n t n,• (γ) t n (0) -1 ∈ [c 1 , c 2 ].
Hence, for n large enough, there exists constants 0 < c 3 < c 4 such that 12) is checked and the proof is complete.

m n (ln(m n )) 2 ln t n (0) 1 -t n,• (γ)/t n (0) ∈ [c 3 , c 4 ]. Since k n /[m n (ln(m n )) 2 ] ≥ d > 2/π, condition (

Proof of main results

Proof of Theorem 1 -Let us prove that

P X n-nαn ,n (t n,L (γ)) ξn ≤ Q(β n ) → 1 + γ 2 , the proof of P X n-nαn ,n (t n,R (γ)) ξn ≥ Q(β n ) → 1 -γ 2 ,
being similar. Let U i := 1 -F (X i ) for i ∈ {1, . . . , n}. Since F is a continuous and strictly increasing function,

P X n-nαn ,n (t n,L (γ)) ξn ≤ Q(β n ) = P U nαn +1,n ≥ Z n ( ξ n ) ,
where

Z n ( ξ n ) = 1 -F Q(β n ) (t n,L (γ)) -ξn .
By assumption, σ -1 n ( ξ n -ξ) = O P (1) that is to say that for all ε > 0, there exists c ε > 0 such that for n large enough, P(σ -1 n | ξ n -ξ| > c ε ) < ε. We have the decomposition

P U nαn +1,n ≥ Z n ( ξ n ) = P {U nαn +1,n ≥ Z n ( ξ n )} ∩ {| ξ n -ξ| > σ n c ε } + P {U nαn +1,n ≥ Z n ( ξ n )} ∩ {| ξ n -ξ| ≤ σ n c ε } =: T 1,n + T 2,n .
For n large enough,

0 < T 1,n < ε. ( 15
)
To deal with the second term, the definition of t n,L (γ) and the fact that U nαn +1,n is distributed as a Beta distribution with parameters nα n + 1

and n -nα n , entail that

P(U nαn +1,n ≥ β n t n,L (γ)) = 1 -γ 2 .
Hence, from Lemma 2,

t n,L (γ) ∼ α n /β n , (16)
and thus, for n large enough,

(1 -ε) ln α n β n ≤ ln (t n,L (γ)) ≤ (1 + ε) ln α n β n .
As a consequence, when the event

{| ξ n -ξ| ≤ σ n c ε } occurs, exp -c ε (1 -ε)σ n ln α n β n ≤ (t n,L (γ)) ξ-ξn ≤ exp c ε (1 + ε)σ n ln α n β n .
Since by assumption, σ n ln(α n /β n ) → 0, one has for n large enough

1 -2c ε σ n ln α n β n ≤ (t n,L (γ)) ξ-ξn ≤ 1 + 2c ε σ n ln α n β n , and thus 1 -F (Q(β n )t n,+ (ξ)) ≤ Z n ( ξ n ) ≤ 1 -F (Q(β n )t n,-(ξ)) , where t n,± (ξ) := (t n,L (γ)) -ξ 1 ± 2c ε σ n ln α n β n .
Applying Lemma 3 with t n = t n,± (ξ) leads to

lim n→∞ 1 |A(α -1 n )| 1 -F (Q(β n )t n,± (ξ)) (t n,± (ξ)) -1/ξ β n -1 = - 1 ξρ .
Using this convergence, when the event {| ξ n -ξ| ≤ σ n c ε } occurs, one has for n large enough that

β n [t n,+ (ξ)] -1/ξ 1 + |A(α -1 n )| 1 + ε ρξ ≤ Z n ( ξ n ) ≤ β n [t n,-(ξ)] -1/ξ 1 -|A(α -1 n )| 1 + ε ρξ .
As a consequence, for n large enough,

p n,--ε ≤ T 2,n ≤ p n,+ , (17)
where

p n,± := P U nαn +1,n ≥ β n [t n,± (ξ)] -1/ξ 1 ± |A(α -1 n )| 1 + ε ρξ .
Since U nαn +1,n is distributed as a Beta distribution with parameters nα n +

) 18
Gathering ( 17) and ( 18) entail that for n large enough,

1 + γ 2 -2ε ≤ T 2,n ≤ 1 + γ 2 + ε. ( 19
)
Collecting ( 15) and ( 19) conclude the proof. 2/π and remark that when u ∈ (0, 1), t n,R (u) ∈ (t n (0), β -1 n ). If we prove the existence and the uniqueness of the solution of the equation

Proof of

f s (H) n ,dn (x) = t n,R (γ), (20)
when x ∈ (t n (0), β -1 n ), then necessarily equation [START_REF] Embrechts | Modelling Extremal Events for Insurance and Finance[END_REF] will also admit an unique solution given by γ The inequality

k n > 2 π
(ln(t n (0)) 2 (1 -t n,R (γ)/t n (0)) 2 , entails that f 1,dn (t n (0)) < t n,R (γ). Moreover, since t n (0) > 1, on has n )] 2 yields to β -1 n < exp(d -1 n -1). Hence, the function f -1,dn is increasing on the interval (t n (0), β -1 n ). Moreover, since t n (0) > 1, f -1,dn (t n (0)) < t n (0) < t n,R (γ) and finally, the inequality

f 1,dn (β -1 n ) > β -1 n > t n,
k n > 2 π
(ln(β n )) 2 (1 -β n t n,R (γ)) 2 , entails that f -1,dn (β -1 n ) > t n,R (γ) and the proof is complete.

  where sgn is the sign function. Using the fact that σ -1 n | ξ n -ξ|/ξ d -→ |D| for some sequence σ n → 0, we propose to replace | ξ n -ξ|/ξ by rσ n where r > 0 is any constant in the support of the distribution of |D|. For instance, r can be the expectation of |D| or a random value generated from |D|. To approximate W n,• (γ), we also need to estimate sgn( ξ n -ξ). An estimator s n is proposed in the next result.

  ) with the constraint kn ≤ n -1. This bias-corrected version is only used when the uniqueness of a solution to (20) has been shown. This corresponds to k ≥ 26 for n = 200, k ≥ 34 for n = 500 and k ≥ 40 for n = 1000.

  by γN,L (k n ) and γN,R (k n ). Of course, it is expected to obtain for γN (k n ) a value close to γ and for γN,L (k n ) and γN,R (k n ), values close to (1 + γ)/2.

  The bias reduction enlarges the range of values of k for which the interval has a coverage rate close to the objective. The new normal ridge confidence interval performs especially well for the Student distribution (ξ = 1, ρ = -1) and the Fréchet distribution (ξ = 0.5, ρ = -1). However, it is not very accurate for the log-gamma distribution (ξ = 0.5, ρ = 0), while our bias reduced interval is always close to the nominal level of 0.95 even if the second-order condition (A.2) is not valid in this case.For the left bound, in the central column of Figure1, our intervals coverage rates are very close to the nominal level of 0.975, as long as k does not exceed 160. As for the right bound, in the right column of Figure1, the nominal level is only hit once, for a value of k below 200, depending on the distribution and the interval considered (bias reduced or not).However, coverage curves only indicate approximately what will happen if wechoose a fixed k. Usually, we select a value of k which seems adequate with regard to the Hill plot, as it will be presented in the next paragraph.Selection of the sequencek n -To select the number of observations used to estimate the tail-index, we propose for each of the N samples to pick k n in a stability region of the plot {(j, ξ (H) n (j)); j ∈ J} where J := { 0.05n , . . . , 0.5n }. To detect this region, we use the following procedure. Let J (H) := min j∈J ξ (H) n (j); max j∈J ξ (H) n (j) ,

Figure 1 :

 1 Figure 1: Coverage rate of the interval and of the left and right bounds (from left to right) for absolute Student, Fréchet, Burr and log-gamma distributions and sample size n = 500.

  terval for year 1985[177, 850] does not match confidence intervals for years1990 and 1991 (respectively [46, 158] and[32, 116]).

Figure 2 :

 2 Figure 2: Selection procedure for k (opt) n

Figure 3 :

 3 Figure 3: Q-Q plot for year for the k (opt) n

Figure 4 :

 4 Figure 4: Application to Norwegian Fire Insurance data set for years 1985 to 1992.

  F 1,n . Since ln t n,• ∼ ln(α n /β n ) and, by assumption, σ -1 n ( ξ n -ξ)/ξ d -→ D, one has that Z n /[σ n ln(α n /β n )] = O P (1), where O P (1) represents a stochastically bounded random variable. Since σ n ln(α n /β n ) → 0, we thus have that Z n = O P (σ n ln(α n /β n )) = o P (1). A

Lemma 5

 5 Let m n be a positive sequence such that m n → ∞ and m n (ln(m n )) 2 = o(n). The conditions of Proposition 3 are satisfied by taking β n = c/n with c > 0, α n = δ n /n and k

Proposition 2 -Proof of Proposition 3 -

 23 Let us introduce the notations ζ 1,n := σ -1 n ξ n -ξ and ζ 2,n := ς -1 n ξn -ξ .One has sn = sgn( ξ n -ξ)sgn(1 -∆ n ), where ∆ n := ς n σ n ζ 2,n ζ 1,n .As a consequence,P s n = sgn( ξ n -ξ) = P (sgn(1 -∆ n ) = 1) = P(∆ n < 1).The conclusion is straightforward since, under the assumptions of Proposi-We focus on the equation (10) with • = R. The case • = L is very similar and is omitted for the sake of conciseness. For s ∈ {-1, 1} and d > 0, let us introduce the function f s,d : [0, ∞) → R defined for all x ≥ 0 by f s,d (x) = x(1 + sd ln x). Let us also introduce the sequence d n := k -1/2 n

  = 2F beta (β n x * n , nα n + 1, n -nα n ) -1, where x *n is the solution of[START_REF] Kocherginsky | Practical confidence intervals for regression quantiles[END_REF].Let us start with the study of the function f s,d . When s = 1, the functionf s,d is decreasing on (0, exp(-(1 + d -1))) and increasing elsewhere. We alsoremark that f 1,d (x) ≤ 0 for all x ≤ exp(-(1 + d -1 )) and that f s,d (x) → +∞ as x → ∞. When s = -1, the function is increasing on (0, exp(d -1 -1))and decreasing elsewhere. Moreover, 0 < f s,d (x) ≤ d exp(d -1 -1) forx ≤ exp(d -1 -1) and f s,d (x) → -∞ as x → ∞.First, let us consider the situation when s (H) n = 1 in equation[START_REF] Kocherginsky | Practical confidence intervals for regression quantiles[END_REF].

= - 1 .

 1 R (γ). Hence, in view of the above study of the function f 1,d , the existence and uniqueness of the solution is proved. Now, let us focus of the case s (H) n The inequality k n > (2/π)[1 + ln(β -1

Table 1 :

 1

	Coverage	Gardes	Gardes and Maistre, Buitendag et al.,
	n
	probability and Maistre bias reduced	new normal ridge
	γN (k (opt) n	)	0.9355	0.9463	0.9732
	200
	γN ( n/5 ) 0.9236	0.9580	0.9531
	γN (k (opt) n	)	0.9446	0.9520	0.9764
	500
	γN ( n/5 ) 0.9520	0.9653	0.9604
	γN (k (opt) n	)	0.9510	0.9535	0.9717
	1000
	γN ( n/5 ) 0.9531	0.9530	0.9713
	Coverage	Gardes	Gardes and Maistre, Buitendag et al.,
	n
	probability and Maistre bias reduced	new normal ridge
	γN (k (opt) n	)	0.9066	0.9381	0.9506
	200
	γN ( n/5 ) 0.8663	0.9147	0.9484
	γN (k (opt) n	)	0.9240	0.9427	0.9571
	500
	γN ( n/5 ) 0.9037	0.9353	0.9415
	γN (k (opt) n	)	0.9400	0.9512	0.9563
	1000
	γN ( n/5 ) 0.9191	0.9399	0.9294

Empirical coverage probabilities for a theoretical coverage probability of γ = 0.95 and a Student distribution (ξ = 0.5 and ρ = -1).

Table 2 :

 2

	Coverage	Gardes	Gardes and Maistre, Buitendag et al.,
	n
	probability and Maistre bias reduced	new normal ridge
	γN (k (opt) n	)	0.9203	0.9458	0.9543
	200
	γN ( n/5 ) 0.8823	0.9241	0.9495
	γN (k (opt) n	)	0.9386	0.9500	0.9374
	500
	γN ( n/5 ) 0.9171	0.9453	0.9426
	γN (k (opt) n	)	0.9501	0.9556	0.9300
	1000
	γN ( n/5 ) 0.9306	0.9483	0.9260

Empirical coverage probabilities for a theoretical coverage probability of γ = 0.95 and a Student distribution (ξ = 1 and ρ = -1).

Table 3 :

 3

	Coverage	Gardes	Gardes and Maistre, Buitendag et al.,
	n
	probability and Maistre bias reduced	new normal ridge
	γN (k (opt) n	)	0.9429	0.9519	0.9656
	200
	γN ( n/5 ) 0.9308	0.9613	0.9467
	γN (k (opt) n	)	0.9465	0.9475	0.9684
	500
	γN ( n/5 ) 0.9588	0.9681	0.9433
	γN (k (opt) n	)	0.9493	0.9463	0.9707
	1000
	γN ( n/5 ) 0.9535	0.9498	0.9565

Empirical coverage probabilities for a theoretical coverage probability of γ = 0.95 and a Fréchet distribution (ξ = 0.5 and ρ = -1).

Table 4 :

 4

	Coverage	Gardes	Gardes and Maistre, Buitendag et al.,
	n
	probability and Maistre bias reduced	new normal ridge
	γN (k (opt) n	)	0.9319	0.9496	0.9611
	200
	γN ( n/5 ) 0.9061	0.9460	0.9493
	γN (k (opt) n	)	0.9410	0.9507	0.9670
	500
	γN ( n/5 ) 0.9384	0.9578	0.9441
	γN (k (opt) n	)	0.9518	0.9548	0.9696
	1000
	γN ( n/5 ) 0.9458	0.9530	0.9444

Empirical coverage probabilities for a theoretical coverage probability of γ = 0.95 and a Burr distribution (ξ = 0.5 and ρ = -√ 0.5).

Table 5 :

 5

	Coverage	Gardes	Gardes and Maistre, Buitendag et al.,
	n
	probability and Maistre bias reduced	new normal ridge
	γN (k (opt) n	)	0.9285	0.9503	0.9550
	200
	γN ( n/5 ) 0.8951	0.9316	0.9514
	γN (k (opt) n	)	0.9464	0.9541	0.9200
	500
	γN ( n/5 ) 0.9336	0.9530	0.9301
	γN (k (opt) n	)	0.9550	0.9557	0.8846
	1000
	γN ( n/5 ) 0.9441	0.9527	0.8938

Empirical coverage probabilities for a theoretical coverage probability of γ = 0.95 and a Burr distribution (ξ = 0.5 and ρ = -1).

Table 6 :

 6 We only use data from years 1985 to 1992, for which the sample sizes are similar, ranging from 607 in 1985 to 827 in 1988. For the quantile level, we take β = 1/638 where 638 is the median value of the sample sizes from years 1985 to 1992. For years with a sample size larger than 638, the level β can thus be considered as extreme. For each year, we select k

	(opt)
	n

Empirical coverage probabilities for a theoretical coverage probability of γ = 0.95 and a log-Gamma distribution (ξ = 0.5 and ρ = 0).

  + 1, n -nα n (1 + δ n,± ); nα n + 1, n -nα n , , nα n + 1, n -nα n (1 + θδ n,± ) ∼ α n , with θ ∈ (0, 1). From [13, Lemma 1] f beta (ζ n , nα n + 1, n -nα n ) = O n 1/2 -γ)/2, nα n + 1, n -nα n ) ∼ α n , we have p n,± = (1+γ)/2+ O δ n,± (nα n ) 1/2 . Since by assumption δ n,± (nα n ) 1/2 = O (nα n ) 1/2 max A(α -1 n ); σ n ln

									α	n 1/2	.
	Since F ← beta ((1 α n β n	= o(1),
	we finally obtained
						p n,± =	1 + γ 2	+ o(1).	(
						n )|	1 + ε ρξ
	= F ← beta	1 -γ 2	, nα n + 1, n -nα n		1 ± 2c ε σ n ln	α n β n	1 ± |A(α -1 n )|	1 + ε ρξ	,
	one has
	p n,± = 1-F beta F ← beta , nα n where 1 -γ 2
	δ n,± = 1 ± 2c ε σ n ln	α n β n	1 ± |A(α -1 n )|	1 + ε ρξ	-1 = O σ n log	α n β n	+O A(α -1 n ) .
	The notation O(1) stands for an asymptotically bounded sequence. From the
	mean value theorem,
	p n,± =	1 + γ 2	+ δ n,± F ←

1 and n -nα n and since

β n [t n,± (ξ)] -1/ξ 1 ± |A(α -1 beta 1 -γ 2 , nα n + 1, n -nα n × f beta (ζ n , nα n + 1, n -nα n ) ,

where

ζ n = F ← beta 1 -γ 2

Appendix -Proofs and additional results

Preliminaries results

Lemma 1

i) If α n is a sequence converging to 0 with nα n → ∞ then for all δ ∈ (0, 1) and for n large enough, S inf (α n , n, δ) = ∅ and S sup (α n , n, δ) = ∅.

ii) If β n ∈ (0, 1) is a sequence such that nβ n → c ∈ [0, 1) then for n large enough, S sup (α n , n, δ) = ∅ for all δ < exp(-c).

Proofi) The function j → F beta (α n ; n -j + 1, j) is increasing. Hence, for all j ∈ {1, . . . , n},

Since α n n → 0 and 1 -(1 -α n ) n → 1, the conclusion is straightforward. ii) As shown before, for all j ∈ {1, . . . , n},

Additional results

Additional results on the one-sided confidence intervals are given below. 18: Empirical coverage probabilities for the right bound (with an objective of (1 + γ)/2 = 0.975) and a log-Gamma distribution (ξ = 0.5 and ρ = 0).


Gabriel Peyré
email: gabriel.peyre@polytechnique.fr

RANDOM SENSING OF GEOMETRIC IMAGES

Keywords: Compressed sensing, best basis, geometry of images

This paper proposes an extension of compressed sensing that allows to express the sparsity prior in a dictionary of bases. This enables the use of the universal sampling strategy of compressed sensing together with an adaptive recovery process that adapts the basis to the structure of the sensed signal. A fast greedy scheme is used during reconstruction to estimate the best basis using an iterative refinement. Numerical experiments on geometrical images show that adaptivity is indeed crucial to capture the structures of complex natural signals.

Introduction

Classical Sampling vs. Compressed Sensing

The classical sampling theory of Shannon is based on uniform smoothness assumptions (low frequency spectral content). Under band limited condition, finely enough sampled functions can be recovered from a set of n pointwise measurements.

However, being bandlimited is not a good model for natural images and one usually assume that such a function f has a decomposition with few elements on some fixed orthogonal basis B. This kind of assumption is at the heart of the compression of an image using a wavelet expansion. Under such sparseness assumption, one can hope to use a much smaller number n < N of measurements, which are linear projections Φf = { f, φ i } n i=1 on a set of fixed vectors φ i ∈ R N . The price to pay for this compressed sampling strategy is a non-linear reconstruction procedure to recover f from the compressed representation Φf . This theory of compressed acquisition of data has been pushed forward during last few years conjointly by Candès and Tao [START_REF] Candès | Robust uncertainty principles: Exact signal reconstruction from highly incomplete frequency information[END_REF] and Donoho [START_REF] Donoho | Compressed sensing[END_REF].

In order for this recovery to be effective, one needs sensing vectors φ i that are incoherent with the vectors of B. A convenient way to achieve this property is to use random vectors φ i , which cannot be sparsely represented in basis B.

Application in imaging.

Compressed sensing acquisition of data could have an important impact for the design of imaging devices where data acquisition is expensive. For instance in seismic or magnetic imaging one could hope to use few random projections of the object to acquire together with a high precision reconstruction.

Analogies in physiology.

This compressed sampling strategy could potentially lead to interesting models for various sensing operations performed biologically. Skarda and Freeman [START_REF] Skarda | Does the brain make chaos in order to make sense of the world[END_REF] have proposed a non-linear chaotic dynamic to explain the analysis of sensory inputs. This chaotic state of the brain ensures robustness toward unknown events and unreliable measurements, without using too many computing resources. While the theory of compressed sensing is presented here as a random acquisition process, its extension to deterministic or dynamic settings is a fascinating area for future research in signal processing.

The Best Basis Approach

Frames vs. dictionary of bases. Fixed orthogonal bases are not flexible enough to capture the complex redundancy of natural images. For instance the orthogonal wavelet transform [START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF] lacks of translation and rotation invariance and is not efficient to compress geometric images [START_REF] Pennec | Bandelet Image Approximation and Compression[END_REF]. It is thus useful to consider families of vectors that are redundant but offer a stable decomposition. For instance, frames of translation invariant wavelets have been used for image denoising and frames of rotation-invariant Gabor functions are useful to characterize textures [START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF]. However, to capture the complex structure of natural images, one needs a very large number of such elementary atoms. Frame theory suffers from both theoretical difficulties (lack of stability) and technical problems (computational complexity) when the number of basis vectors increases too much. To cope with these problems, one can consider a dictionary D = {B λ } λ∈Λ of orthogonal bases B λ . Choosing an optimal basis in such a dictionary allows to adapt the approximation to the complex content of a specific signal.

Images and geometry.

Cartoon images is a simple model that captures the sketch content of natural images. Figure 1,(a), shows such a geometrically regular image, which contains smooth areas surrounded by regular curves. The curvelet frame of Candès and Donoho [START_REF] Candès | New tight frames of curvelets and optimal representations of objects with piecewise C 2 singularities[END_REF] can deal with such a regularity and enjoys a better approximation rate than traditional isotropic wavelets. This results can be enhanced using a dictionary of locally elongated functions that follow the image geometry. Bandelets bases of Le Pennec and Mallat [START_REF] Pennec | Bandelet Image Approximation and Compression[END_REF] provide such a geometric dictionary together with a fast optimization procedure to compute a basis adapted to a given image.

Adaptive biological computation. Hubel and Wiesel

have shown that low level computation done in area V1 of the visual cortex are well approximated by multiscale oriented linear projections [START_REF] Hubel | Receptive fields and functional architecture of monkey striate cortex[END_REF]. Olshausen and Field proposed in [START_REF] Olshausen | Emergence of simple-cell receptive-field properties by learning a sparse code for natural images[END_REF] that redundancy is important to account for sparse representation of natural inputs. However further non-linear processings are done by the cortex to remove high order geometrical correlations present in natural images. Such computations are thought to perform long range groupings over the first layer of linear responses [START_REF] Lee | Computations in the early visual cortex[END_REF] and thus correspond to an adaptive modification of the overall neuronal response. The best-basis coding strategy could thus offer a signal-processing counterpart to this neuronal adaptivity. In this paper we restrict ourselves to orthogonal best basis search. Future extensions of this approach to dictionaries of redundant transforms are likely to improve numerical results and better cope with biological computations.

Best Basis Computation

In the following we use the norms

||f || p ℓ p = i |f [i]| p and ||f || ℓ 0 = # {i \ f [i] = 0 }. A dictionary D Λ = {B λ } λ∈Λ is a set of orthogonal bases B λ = {ψ λ m } m of R N .
We associate a cost pen(λ) which is a prior complexity measure on each basis B λ that satisfies λ∈Λ 2 -pen(λ) = 1. A fixed weight pen(λ) = log 2 (M ) can be used if the size M of D is finite. The parameter pen(λ) can be interpreted as the number of bits needed to specify a basis B λ . Following the construction of Coifman and Wickerhauser [START_REF] Coifman | Entropy-based algorithms for best basis selection[END_REF], a best basis B λ * adapted to a signal minimizes a Lagrangian L

λ * = argmin λ∈Λ L(f, λ, t) (1)
where

L(f, λ, t) = ||Ψ λ f || ℓ 1 + C 0 t pen(λ),
where

Ψ λ = [ψ λ 0 , . . . , ψ λ N -1 ]
T is the transform matrix defined by B λ . This best basis B λ * is thus the one that gives the sparsest description of f as measured by the ℓ 1 norm. The Lagrange multiplier C 0 t weights the penalization pen(λ) associated to the complexity of a basis B λ . The parameter t is the level of deterministic noise caused by not having a perfectly k-sparse function but rather a function whose coefficients in basis B λ * decay fast. This parameter is not available since it depends on the decay speed and on k (which in turn depends on the number of sensed vectors). However the recovery algorithm presented in section 2.3 does not need the exact value of t since an upper bound t s is estimated during iterations.

For a dictionary D that enjoys a multiscale structure, the optimization of L is carried out with a fast procedure, see [START_REF] Coifman | Entropy-based algorithms for best basis selection[END_REF][START_REF] Pennec | Bandelet Image Approximation and Compression[END_REF] for practical examples of this process.

Compressed Sensing Reconstruction

In this paper, the sampling matrix Φ = [φ 0 , . . . , φ n-1 ] T is defined by random points {φ i } i of unit length although other random sensing schemes can be used, see [START_REF] Candès | Robust uncertainty principles: Exact signal reconstruction from highly incomplete frequency information[END_REF].

Basis Pursuit Formulation

Searching for the sparsest signal f * in some basis B = {ψ m } m that matches the sensed values y = Φf leads to consider a penalized variational problem

f * = argmin g∈R N 1 2 ||Φg -y|| 2 ℓ 2 + t||Ψg|| ℓ 1 . ( 2
)
where

Ψ = [ψ 0 , . . . , ψ N -1 ]
T is the transform matrix defined by B. The Lagrange multiplier t accounts both for stabilisation against noise and for approximate sparsity, which is common in practical applications.

Candès and Tao in [START_REF] Candès | Robust uncertainty principles: Exact signal reconstruction from highly incomplete frequency information[END_REF] and Donoho in [START_REF] Donoho | Compressed sensing[END_REF] have shown that, if f is sparse enough in some basis B, f can be recovered from the sensed data y = Φf . More precisely, they show that it exists a constant

C such that if ||Ψf || ℓ 0 ≤ k then if n ≥ C log(N ) k one recovers f * = f .

Best Basis Pursuit

The compressed sensing machinery is extended to a dictionary of bases D Λ by imposing that the recovered signal is sparse in at least one basis of D Λ . To avoid using too complex basis the recovery process from noisy measurements takes into account a complexity pen(λ) of the optimal basis B λ . The original recovery procedure (2) is replaced by the following minimization

f * = argmin g∈R N min λ∈Λ 1 2 ||Φg -y|| 2 ℓ 2 + t ||Ψ λ g|| ℓ 1 + C 0 t 2 pen(λ) , (3)
where the penalization C 0 t 2 pen(λ) is the same as in equation (1).

Iterative Thresholding For Sparsity Minimization

The recovery procedure suggested by equation (3) corresponds to the inversion of the operator Φ under sparsity constraints on the observed signal f . In this paper we extends the algorithm of Daubechies et al. [START_REF] Daubechies | An iterative thresholding algorithm for linear inverse problems with a sparsity constraint[END_REF] to the setting of a dictionnary of bases. Searching in the whole dictionary D for the best basis that minimizes formulation [START_REF] Coifman | Entropy-based algorithms for best basis selection[END_REF] is not feasible for large dictionaries, which typically contain of the order of 2 N bases. Instead we propose a greedy search for the best basis during the recovery process. We make repeated use of the soft thresholding operator defined in an orthogonal basis B = {ψ m } m by S t (g) = m s t ( g, ψ m ) ψ m where s t (x) = sign(x)(|x| -t) + in the following algorithm • Initialization. Set s = 0, f 0 = 0 and choose λ 0 ∈ Λ at random or using some default choice (such as a DCT basis in 1D or a wavelet basis in 2D).

Best Bandelet Basis Compressed Sensing

Adapted Bandelet Transform

The bandelet bases dictionary was introduced by Le Pennec and Mallat [START_REF] Pennec | Bandelet Image Approximation and Compression[END_REF] to perform adaptive approximation of images with geometric singularities, such as the cartoon image in figure 1,(a). We present a simple implementation of the bandelet transform inspired from [START_REF] Peyré | Surface compression with geometric bandelets[END_REF].

A bandelet basis B λ is parameterized by λ = (Q, {θ S } S∈Q ), where Q is a quadtree segmentation of the pixels locations and θ S ∈ [0, π[ ∪ Ξ is an orientation (or the special token Ξ) defined over each square S of the segmentation, see figure 1,(a). The bandelet transform corresponding to this basis applies independently over each square S of the image either • if θ S = Ξ: a 2D isotropic wavelet transform, • if θ S = Ξ: a 1D wavelet transform along the direction defined by the angle θ S . We now detail the latter transform. The position of a pixel x = (x 1 , x 2 ) ∈ S with respect to the direction θ S is p x = sin(θ S )x 1 -cos(θ S )x 2 . The m pixels {x (i) } in S are ranked according to the 1D ordering p x (0) ≤ p x (1) ≤ . . . ≤ p x (m-1) . This ordering allows to turn the image {f 1,(c). The bandelet transform of the image f inside S is defined as the 1D Haar transform of the signal f 1D , see figure 1,(d). This process is both orthogonal and easily invertible, since one only needs to compute the inverse Haar transform and pack the retrieved coefficients at the original pixels locations. Keeping only a few bandelet coefficients and setting the others to zero performs an approximation of the original image that follows the local direction θ S , see figure 1 (f).

[x]} x∈S defined over S into a 1D signal f 1D [i] = f [x (i) ], see fig- ure
In order to restrict the number of tested geometries θ for a square S ∈ Q containing #S pixels, we follow [START_REF] Peyré | Surface compression with geometric bandelets[END_REF] and use the set of directions that pass through two pixels of S. The number of such directions is of the order of (#S) 2 . For this bandelet dictionary, the penalization of a basis B λ where λ = (Q, {θ S } S ) is defined as pen(λ) = #Q + S∈Q 2 log 2 (#S), where #Q is the number of leaves in Q. A fast best basis search, described in [START_REF] Peyré | Surface compression with geometric bandelets[END_REF], allows to define a segmentation Q and a set of directions {θ S } S adapted to a given image f by minimizing (1). This process segments the image into squares S on which f is smooth, thus setting θ S = Ξ and squares containing an edge, where θ S closely matches the direction of this singularity.

Numerical Results

The geometric image depicted in figure 2, (a) is used to compare the performance of the original compressed sensing algorithm in a wavelet basis to the adaptative algorithm in a best bandelet basis. Since the wavelet basis is not adapted to the geometric singularities of such an image, reconstruction (b) has strong ringing artifacts. The adapted reconstruction (c) exhibits fewer such artifacts since the bandelet basis functions are elongated and follow the geometry. The segmentation is depicted after the last iteration, together with the chosen direction θ S which closely matches the real geometry. On figure 2,(d/e/f), one can see a comparison for a natural image containing complex geometric structures such as edges, junctions and sharp line features. The best bandelet process is able to resolve these features efficiently.

Conclusion

Using a dictionary of bases decouples the approximation process from the redundancy needed for adaptivity and requires the design of a penalization cost on the set of bases. This framework is not restricted to orthogonal bases, although it is a convenient mathematical way to ensure the compressed sensing recovery condition. This best basis approach to sensing and recovery is also a promising avenue for interactions between biological processing, where a deterministic or chaotic process is highly probable and signal processing, where randomization has proven useful to provide universal coding strategies.

Figure 1 .

 1 Figure 1. (a) a geometric image together with some adapted dyadic segmentation Q ; (b) a square S together with some adapted direction θ S ; (c) the 1D signal f 1D obtained by mapping the pixels values f (x (i) ) on a 1D axis ; (d) the 1D Haar coefficients of f 1D ; (e) the 1D approximation obtained by reconstruction from the 20 largest Haar coefficients ; (f) the corresponding square approximated in bandelet.

Figure 2 .

 2 Figure 2. (a/d) original image ; (b/e) compressed sensing reconstruction using the wavelet basis, n = N/6 (b: PSNR=22.1dB, b: PSNR=23.2dB) ; (c/f) reconstruction using iteration in a best bandelet basis (c: PSNR=24.3dB, f: PSNR=25.1dB).

• Step 1: Updating the estimate.

  Enforce condition y = Φf s+1 using f s+1 = f s + Φ T (y -Φ f s ). • Step 2: Denoising the estimate. Compute the value of the current threshold t s = 3σ s using the estimator of the noise level σ s = median(|Ψf s+1 |)/0.6745. Perform the denoising f s+1 = S ts (f s+1 ) where S ts is the threshold operator at t s in the basis B λs . • Step 3: Update best basis. Update the best basis using λ s+1 = argmin λ L(f s+1 , λ, t s ). For typical dictionaries such as the ones considered in this paper, this minimization is carried out with a fast procedure, as seen in subsection 1.3. • Stopping criterion. If s < s max , go to step 1, otherwise stop iterations. In all our experiments, the number of iterations is set to s max = 20. The stochastic noise level t s is due to random acquisition corrupting the current estimate f s . It is computed though the median estimator and its value decays during iterations toward the deterministic noise value t that is due to approximate sparsity.