Abdelaati Daouia

Laurent Gardes

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

On kernel smoothing for extremal quantile regression

Keywords: Regression, extreme quantile, extreme-value index, kernel smoothing, von-Mises condition, asymptotic normality

Nonparametric regression quantiles obtained by inverting a kernel estimator of the conditional distribution of the response are long established in statistics. Attention has been, however, restricted to ordinary quantiles staying away from the tails of the conditional distribution. The purpose of this paper is to extend their asymptotic theory far enough into the tails. We focus on extremal quantile regression estimators of a response variable given a vector of covariates in the general setting, whether the conditional extreme-value index is positive, negative, or zero. Specifically, we elucidate their limit distributions when they are located in the range of the data or near and even beyond the sample boundary, under technical conditions that link the speed of convergence of their (intermediate or extreme) order with the oscillations of the quantile function and a von-Mises property of the conditional distribution. A simulation experiment and an illustration on real data were presented. The real data are the American electric data where the estimation of conditional extremes is found to be of genuine interest.

Introduction

Quantile regression plays a fundamental role in various statistical applications. It complements the classical regression on the conditional mean by offering a more useful tool for examining how a vector of regressors X ∈ R p influences the entire distribution of a response variable Y ∈ R. The nonparametric regression quantiles obtained by inverting a kernel estimator of the conditional distribution function are used widely in applied work and investigated extensively in theoretical statistics. See, for example [START_REF] Berlinet | Asymptotic normality of convergent estimates of conditional quantiles[END_REF][START_REF] Samanta | Non-parametric estimation of conditional quantiles[END_REF][START_REF] Stone | Consistent nonparametric regression (with discussion)[END_REF][START_REF] Stute | Conditional empirical processes[END_REF], among others. Attention has been, however, restricted to conditional quantiles having a fixed order α ∈ (0, 1). In the following, the order α has to be understood as the conditional probability to be larger than the conditional quantile. In result, the available large sample theory does not apply sufficiently far in the tails.

There are many important applications in ecology, climatology, demography, biostatistics, econometrics, finance, insurance, to name a few, where extending that conventional asymptotic theory further into the tails of the conditional distribution is an especially welcome development. This translates into considering the order α = α n → 0 or α n → 1 as the sample size n goes to infinity. Motivating examples include the study of extreme rainfall as a function of the geographical location [START_REF] Gardes | Conditional extremes from heavy-tailed distributions: An application to the estimation of extreme rainfall return levels[END_REF], the estimation of factors of high risk in finance [START_REF] Tsay | Analysis of financial time series[END_REF], the assessment of the optimal cost of the delivery activity of postal services [START_REF] Daouia | Frontier estimation and extreme value theory[END_REF], the analysis of survival at extreme durations [START_REF] Koenker | Reappraising medfly longevity: A quantile regression survival analysis[END_REF], the edge estimation in image reconstruction [START_REF] Park | On nonparametric estimation of data edges[END_REF], the accurate description of the upper tail of the claim size distribution for reinsurers [START_REF] Beirlant | Statistics of extremes: Theory and applications[END_REF], the analysis of environmental time series with application to trend detection in ground-level ozone [START_REF] Smith | Extreme value analysis of environmental time series: An application to trend detection in ground-level ozone (with discussion)[END_REF], the estimation of autoregressive models with asymmetric innovations [START_REF] Feigin | Linear programming estimators and bootstrapping for heavy tailed phenomena[END_REF], etc.

There have been several efforts to treat the asymptotics of extreme conditional quantile estimators in semi/parametric and other nonparametric regression models. For example, Chernozhukov [START_REF] Chernozhukov | Extremal quantile regression[END_REF] and Jurecková [START_REF] Jurecková | Remark on extreme regression quantile[END_REF] considered the extreme quantiles in the linear regression model and derived their asymptotic distributions under various distributions of errors. Other parametric models are proposed in [START_REF] Davison | Models for exceedances over high thresholds[END_REF][START_REF] Smith | Extreme value analysis of environmental time series: An application to trend detection in ground-level ozone (with discussion)[END_REF], where some extreme-value based techniques are extended to the point-process view of high-level exceedances. A semi-parametric approach to modeling trends in sample extremes, based on local polynomial fitting of the Generalized extreme-value distribution, has been introduced in [START_REF] Davison | Local likelihood smoothing of sample extremes[END_REF]. Hall and Tajvidi [START_REF] Hall | Nonparametric analysis of temporal trend when fitting parametric models to extreme-value data[END_REF] suggested a nonparametric estimation of the temporal trend when fitting parametric models to extreme values. Another semi-parametric method has been developed in [START_REF] Beirlant | Local polynomial maximum likelihood estimation for Pareto-type distributions[END_REF], where the regression is based on a Pareto-type conditional distribution of the response. Fully nonparametric estimators of extreme conditional quantiles have been discussed in [START_REF] Beirlant | Local polynomial maximum likelihood estimation for Pareto-type distributions[END_REF][START_REF] Chavez-Demoulin | Generalized additive modelling of sample extremes[END_REF], where the former approach is based on the technique of local polynomial maximum likelihood estimation, while spline estimators are fitted in the latter by a maximum penalized likelihood method. Recently, [START_REF] Gardes | A moving window approach for nonparametric estimation of the conditional tail index[END_REF][START_REF] Gardes | Functional nonparametric estimation of conditional extreme quantiles[END_REF] proposed, respectively, a moving-window based estimator for the tail index and extreme quantiles of heavy-tailed conditional distributions, and they established their asymptotic properties.

In the context of kernel-smoothing, the asymptotic theory for quantile regression in the tails is relatively unexplored and still in full development. Daouia et al. [START_REF] Daouia | Kernel estimators of extreme level curves[END_REF] have extended the asymptotics further into the tails in the particular setting of a heavy-tailed conditional distribution, while [START_REF] Girard | Frontier estimation via kernel regression on high powertransformed data[END_REF][START_REF] Girard | Central limit theorems for smoothed extreme value estimates of Poisson point processes boundaries[END_REF] have analyzed the case α n = 1/n in the particular situation where the response Y given X = x is uniformly distributed. The purpose of this paper is to develop a unified asymptotic theory for the kernel-smoothed conditional extremes in the general setting where the conditional distribution can be short, light or heavy-tailed. We will focus on the α n → 0 case, which corresponds to the class of large quantiles of the upper conditional tail. Similar considerations evidently apply to the case α n → 1. Specifically, we first obtain the asymptotic normality of the extremal quantile regression under the 'intermediate' order condition nh p α n → ∞ where h = h n → 0 stands for the bandwidth involved in the kernel smoothing estimation. Next we extend the asymptotic normality far enough into the 'most extreme' order-β n regression quantiles with β n /α n → 0, thus providing a conditional analog of modern extreme-value results [START_REF] De Haan | Extreme value theory: An introduction[END_REF]. We also analyze kernel-smoothed Pickands type estimators of the conditional extreme-value index as in the familiar nonregression case [START_REF] Drees | Refined Pickands estimators of the extreme value index[END_REF].

The paper is organized as follows. Section 2 contains the basic notations and assumptions. Section 3 states the main results. Section 4 presents some simulation evidence and practical guidelines. Section 5 provides a motivating example in production theory, and Section 6 collects the proofs.

The setting and assumptions

Let (X i , Y i ), i = 1, . . . , n, be independent copies of a random pair (X, Y ) ∈ R p × R. The conditional survival function (csf) of Y given X = x is denoted by F (y|x) = P(Y > y|X = x)
and the probability density function (pdf) of X is denoted by g. We address the problem of estimating extreme conditional quantiles

q(α n |x) = F ← (α n |x) = inf{t, F (t|x) ≤ α n },
where α n → 0 as n goes to infinity. In the following we denote by y F

(x) = q(0|x) ∈ (-∞, ∞] the endpoint of the conditional distribution of Y given X = x. The kernel estimator of F (y|x) is defined for all (x, y) ∈ R p × R by Fn (y|x) = n i=1 K h (x -X i )I{Y i > y} n i=1 K h (x -X i ), (1)
where I{.} is the indicator function and h = h n is a nonrandom sequence such that h → 0 as n → ∞. We have also introduced K h (t) = K(t/h)/h p where K is a pdf on R p . In this context, h is called the window-width. Similarly, the kernel estimators of conditional quantiles q(α|x) are defined via the generalized inverse of Fn (.|x):

qn (α|x) = F ← n (α|x) = inf{t, Fn (t|x) ≤ α}, (2)
for all α ∈ (0, 1). Many papers are dedicated to the asymptotic properties of this type of estimator for fixed α ∈ (0, 1): weak and strong consistency are proved respectively in [START_REF] Stone | Consistent nonparametric regression (with discussion)[END_REF] and [START_REF] Gannoun | Estimation non paramétrique de la médiane conditionnelle, médianogramme et méthode du noyau[END_REF], asymptotic normality being established in [START_REF] Stute | Conditional empirical processes[END_REF][START_REF] Samanta | Non-parametric estimation of conditional quantiles[END_REF][START_REF] Berlinet | Asymptotic normality of convergent estimates of conditional quantiles[END_REF]. In Theorem 1 below, the asymptotic distribution of ( 2) is investigated when estimating extreme quantiles, i.e when α = α n goes to 0 as the sample size n goes to infinity. The asymptotic behavior of such estimators then depends on the nature of the conditional distribution tail. In this paper, we assume that the csf satisfies the following von-Mises condition, see for instance [19, equation (1.11.30)]:

(A.
F (y|x) F ′′ (y|x) ( F ′ ) 2 (y|x) = γ(x) + 1,
where F ′ (.|x) and F ′′ (.|x) are respectively the first and the second derivatives of F (.|x).

Here, γ(.) is an unknown function of the covariate x referred to as the conditional extreme-value index. Let us consider, for all z ∈ R, the classical K z function defined for all u ∈ R by

K z (u) = u 1 v z-1 dv.
The associated inverse function is denoted by K -1 z . Then, (A.1) implies that there exists a positive auxiliary function a(.|x) such that, lim

y↑y F (x) F (y + t(x)a(y|x)|x) F (y|x) = 1 K -1 γ(x) (t(x)) , (3)
where t(x) ∈ R is such that 1 + t(x)γ(x) > 0. Besides, (3) implies in turn that the conditional distribution of Y given X = x is in the maximum domain of attraction (MDA) of the extreme-value distribution with shape parameter γ(x), see [START_REF] De Haan | Extreme value theory: An introduction[END_REF]Theorem 1.1.8] for a proof. The case γ(x) > 0 corresponds to the Fréchet MDA and F (.|x) is heavy-tailed while the case γ(x) = 0 corresponds to the Gumbel MDA and F (.|x) is light-tailed. The case γ(x) < 0 represents most of the situations where F (.|x) is short-tailed, i.e F (.|x) has a finite endpoint y F (x), this is referred to as the Weibull MDA.

The convergence (3) is also equivalent to

b(t, α|x) := q(tα|x) -q(α|x) a(q(α|x)|x) -K γ(x) (1/t) → 0 (4)
for all t > 0 as α → 0, see [START_REF] De Haan | Extreme value theory: An introduction[END_REF]Theorem 1.1.6]. For all (x, x ′ ) ∈ R p × R p , the Euclidean distance between x and x ′ is denoted by d(x, x ′ ). The following Lipschitz condition is introduced:

(A.2) There exists c g > 0 such that g(x) -g(x ′ ) ≤ c g d(x, x ′ ).
The last assumption is standard in the kernel estimation framework.

(A.3) K is a bounded pdf on R p , with support S included in the unit ball of R p .

Main results

Let B(x, h) be the ball centered at x with radius h. The oscillations of the csf are controlled by

∆ κ (x, α) := sup (x ′ ,β)∈B(x,h)×[κα,α] F (q(β|x)|x ′ ) β -1 ,
where (κ, α) ∈ (0, 1) 2 . Under assumption (A.1), F (.|x) is differentiable and the associated conditional density will be denoted in the sequel by f (.|x). We first establish the asymptotic normality of qn (α n |x).

Theorem 1. Suppose (A.1), (A.2) and

(A.3) hold. Let 0 < τ J < • • • < τ 2 < τ 1 ≤ 1
where J is a positive integer and x ∈ R p such that g(x) > 0. If α n → 0 and there exists κ ∈ (0, τ J ) such that

nh p α n → ∞, nh p α n (h ∨ ∆ κ (x, α n )) 2 → 0,
then, the random vector

f (q(α n |x)|x) nh p α -1 n (q n (τ j α n |x) -q(τ j α n |x)) j=1,...,J
is asymptotically Gaussian, centered, with covariance matrix K 2 2 /g(x)Σ(x) where Σ j,j ′ (x) = (τ j τ j ′ ) -γ(x) τ -1 j∧j ′ for (j, j ′ ) ∈ {1, . . . , J} 2 . Let us remark that, in the particular case where J = 1, τ 1 = 1 and α n = α is fixed in (0, 1), we find back the result of [START_REF] Berlinet | Asymptotic normality of convergent estimates of conditional quantiles[END_REF]Theorem 6.4]. Theorem 1 can be equivalently rewritten as Corollary 1. Under the assumptions of Theorem 1, the random vector

nh p α n q(α n |x) a(q(α n |x)|x) qn (τ j α n |x) q(τ j α n |x) -1 j=1,...,J
is asymptotically Gaussian, centered, with covariance matrix K 2 2 /g(x) Σ(x) where Σj,j ′ (x) = (τ j τ j ′ ) -(γ(x)∧0) τ -1 j∧j ′ for (j, j ′ ) ∈ {1, . . . , J} 2 . Moreover [19, Theorem 1.2.5] and [19, page 33] show that lim

y↑y F (x) a(y|x) y = γ(x) ∨ 0. (5)
Under the assumptions of Theorem 1, and from [START_REF] Chernozhukov | Extremal quantile regression[END_REF], it follows that qn (τ j α n |x)/q(τ j α n |x) P -→ 1 when n → ∞ which can be read as a weak consistency result for the considered estimator. Besides, if γ(x) > 0, then collecting [START_REF] Chernozhukov | Extremal quantile regression[END_REF] and Corollary 1 shows that the random vector

nh p α n qn (τ j α n |x) q(τ j α n |x) -1 j=1,...,J
is asymptotically Gaussian, centered, with covariance matrix K 2 2 γ 2 (x)/g(x) Σ(x) where the coefficients of the covariance matrix can be simplified Σj,j ′ (x) = τ -1 j∧j ′ for (j, j ′ ) ∈ {1, . . . , J} 2 . Our results thus build on and complement the analysis given by [START_REF] Daouia | Kernel estimators of extreme level curves[END_REF]Theorem 2] in the case γ(x) > 0. As pointed out in [START_REF] Daouia | Kernel estimators of extreme level curves[END_REF], the condition nh p α n → ∞ implies α n > log p (n)/n eventually. This condition provides a lower bound on the order of the extreme conditional quantiles for the asymptotic normality of kernel estimators to hold. We now propose a scheme to estimate extreme conditional quantiles without this restriction. Let α n → 0 and β n /α n → 0 as n → ∞. Suppose one has γn (x) and ân (x) two estimators of γ(x) and a(q(α n |x)|x) respectively. Then, starting from the estimator qn (α n |x) of q(α n |x) defined in [START_REF] Beirlant | Statistics of extremes: Theory and applications[END_REF] and making use of (4), it is possible to build an estimator qn (β n |x) of q(β n |x) which is an extreme conditional quantile of higher order than q(α n |x):

qn (β n |x) = qn (α n |x) + K γn(x) (α n /β n )â n (x). (6)
Let us consider, for all z ∈ R, the function defined for all u > 1 by

K ′ z (u) = ∂K z (u) ∂z = u 1 v z-1 log(v)dv.
The following result provides a quantile regression analog of [START_REF] De Haan | Extreme value theory: An introduction[END_REF]Theorem 4.3.1].

Theorem 2. Suppose (A.1) holds and let α n → 0, β n /α n → 0. Let qn (α n |x) be the kernel estimator of q(α n |x) defined in [START_REF] Beirlant | Statistics of extremes: Theory and applications[END_REF]. Let γn (x) and ân (x) be two estimators of γ(x) and a(q(α n |x)|x) respectively such that

Λ -1 n γn (x) -γ(x), ân (x) a(q(α n |x)|x) -1, qn (α n |x) -q(α n |x) a(q(α n |x)|x) t d -→ ζ(x), (7)
where ζ(x) is a non-degenerate R 3 random vector,

Λ n log(α n /β n ) → 0 and Λ -1 n b(β n /α n , α n |x) K ′ γ(x) (α n /β n ) → 0 as n → ∞. Then, Λ -1 n qn (β n |x) -q(β n |x) a(q(α n |x)|x)K ′ γ(x) (α n /β n ) d -→ c(x) t ζ(x),
where c(x

) t = 1, -(γ(x) ∧ 0), (γ(x) ∧ 0) 2 .
As an illustration, for all r ∈ (0, 1), let us consider τ j = r j-1 , j = 1, . . . , J. The following estimators of γ(x) and a(q(α n |x)|x) are introduced

γRP n (x) = 1 log r J-2 j=1 π j log qn (τ j α n |x) -qn (τ j+1 α n |x) qn (τ j+1 α n |x) -qn (τ j+2 α n |x) âRP n (x) = 1 K γRP n (x) (r) J-2 j=1 π j r γRP n (x)j (q n (τ j α n |x) -qn (τ j+1 α n |x)),
where (π j ) is a sequence of weights summing to one. Let us highlight that γRP n (x) is an adaptation to the conditional case of the Refined Pickands estimator introduced in [START_REF] Drees | Refined Pickands estimators of the extreme value index[END_REF]. The joint asymptotic normality of (γ RP n (x), âRP n (x), qn (α n |x)) is established in the next theorem. Theorem 3. Suppose (A.1), (A.2) and (A.3) hold. Let x ∈ R p such that g(x) > 0. If α n → 0 and there exists κ ∈ (0, τ J ) such that

nh p α n → ∞, nh p α n   h ∨ ∆ κ (x, α n ) ∨ J j=1 b(τ j , α n |x)   2 → 0,
as n → ∞, then the random vector

nh p α n γRP n (x) -γ(x), âRP n (x) a(q(α n |x)|x) -1, qn (α n |x) -q(α n |x) a(q(α n |x)|x)
t is asymptotically centered and Gaussian.

The asymptotic covariance matrix is denoted by S(x). It can be explicitly calculated from [START_REF] Ruppert | Semiparametric Regression[END_REF] in the proof of Theorem 3, but the result would be too complicated to be reported here. As a consequence of the two above theorems, one obtains the asymptotic normality of the extreme conditional quantile estimator built on γRP n (x) and âRP n (x):

qRP n (β n |x) := qn (α n |x) + K γRP n (x) (α n /β n )â RP n (x).
Corollary 2. Suppose (A.1), (A.2) and (A.3) hold. Let x ∈ R p such that g(x) > 0. If α n → 0, β n /α n → 0 and there exists κ ∈ (0, τ J ) such that

nh p α n (log(α n /β n )) 2 → ∞, nh p α n   h ∨ ∆ κ (x, α n ) ∨ J j=1 b(τ j , α n |x) ∨ b(β n /α n , α n |x) K ′ γ(x) (α n /β n )   2 → 0, as n → ∞, then nh p α n qRP n (β n |x) -q(β n |x) a(q(α n |x)|x)K ′ γ(x) (α n /β n ) is asymptotically Gaussian, centered with variance c(x) t S(x)c(x).
Finally, two particular cases of γRP n (x) may be considered. First, constant weights

π 1 = . . . π J-2 = 1/(J -2) yield γRP,1 n (x) = 1 (J -2) log r log qn (τ 1 α n |x) -qn (τ 2 α n |x) qn (τ J-1 α n |x) -qn (τ J α n |x)
.

Clearly, when J = 3, this estimator reduces the kernel Pickands estimator introduced and studied in [START_REF] Daouia | Kernel estimators of extreme level curves[END_REF] in the situation where γ(x) > 0. Second, linear weights

π j = 2j/((J -1)(J -2)) for j = 1, . . . , J -2 give rise to a new estimator γRP,2 n (x) = 2 (J -1)(J -2) log r J-2 j=1 log qn (τ j α n |x) -qn (τ j+1 α n |x) qn (τ J-1 α n |x) -qn (τ J α n |x) ,
which can be read as the average of J -1 estimators γRP,1 n (x). These estimators are now compared on finite sample situations.

Some simulation evidence

Monte Carlo experiments

To evaluate finite-sample performance of the conditional extreme-value index and extreme quantile estimators described above we have undertaken some simulation experiments following the model

Y i = G(X i ) + σ(X i ) U i , i = 1, . . . , n.
The local scale factor, σ(x) = (1 + x)/10, is linearly increasing in x, while the local location parameter

G(x) = x(1 -x) sin 2π(1 + 2 -7/5 ) x + 2 -7/5
has been introduced in [27, Section 17.5.1]. The design points X i are generated following a standard uniform distribution. The U i 's are independent and their conditional distribution given X i = x is chosen to be standard Gaussian, Student t k(x) , or Beta(ν(x), ν(x)), with

k(x) = [ν(x)] + 1, ν(x) = 1 10 + sin(πx) 11 10 - 1 2 exp{-64(x -1/2) 2 } -1
, and [ν(x)] being the integer part of ν(x). Let us recall that the Gaussian distribution belongs to the Gumbel MDA, i.e γ(x) = 0, the Student distribution t k(x) belongs to the Fréchet MDA with γ(x) = 1/k(x) > 0 and the Beta distribution belongs to the Weibull MDA with γ(x) = -1/ν(x) < 0. In all cases we have q(β|x) = G(x) + σ(x) F ← U |X (β|x), for β ∈ (0, 1). All the experiments were performed over 400 simulations for n = 200, and the kernel function K was chosen to be the Triweight kernel

K(t) = 35 32 (1 -t 2 ) 3 I{-1 ≤ t ≤ 1}.
Monte Carlo experiments were first devoted to accuracy of the two conditional extreme-value index estimators γRP,1 n (x) and γRP,2 n (x). The measures of efficiency for each simulation used were the mean squared error and the bias

MSE{γ n (.)} = 1 L L ℓ=1 {γ n (x ℓ ) -γ(x ℓ )} 2 , Bias{γ n (.)} = 1 L L ℓ=1 {γ n (x ℓ ) -γ(x ℓ )}
for γn (x) = γRP,1 n (x), γRP,2 n (x), with the x ℓ 's being L = 100 points regularly distributed in [0, 1]. To guarantee a fair comparison among the two estimation methods, we used for each estimator the parameters (α n , h) minimizing its mean squared error, with α n ranging over A = {0.1, 0.15, 0.2, . . . , 0.95} and the bandwidth h ranging over a grid H of 50 points regularly distributed between

h min = max 1≤i<n |X (i+1) -X (i) | and h max = |X (n) -X (1) |/2, where X (1) ≤ • • • ≤ X (n) are the ordered ob- servations.
The resulting values of MSE and bias are averaged on the 400 Monte Carlo replications and reported in Table 1 for J ∈ {3, 4, 5} and r ∈ {1/J, (J -1)/J}.

) and γRP,2 n (x) -Results averaged on 400 simulations with n = 200. The results may not be available for r = 1/J and J = 5 since the numerator {q n (τ j α n |x) -qn (τ j+1 α n |x)} and the denominator {q n (τ J-1 α n |x) -qn (τ J α n |x)} in the definitions of both estimators might be null when n is not large enough.

It does appear that the results for r = 1/J are superior to those for r = (J -1)/J, uniformly in J. For these desirable results, it may be seen that the estimator γRP,1 n (x) performs better than γRP,2 n (x) in the Gaussian error model, whereas the latter is superior to the former in the Student error model. It may be also seen that there is no winner in the Beta error model in terms of both MSE and bias.

Turning to the performance of the extreme conditional quantile estimators, we consider as above the two measures of performance

MSE{q n (β n |•)} = 1 L L ℓ=1 {q n (β n |x ℓ ) -q(β n |x ℓ )} 2 , Bias{q n (β n |•)} = 1 L L ℓ=1 {q n (β n |x ℓ ) -q(β n |x ℓ )} , for q n (β n |x) = qn (β n |x), qRP,1 n (β n |x), qRP,2 n (β n |x).
The averaged MSE and bias of these three estimators of q(β n |x), computed for β n ∈ {0.05, 0.01, 0.005}, J ∈ {3, 4} and r = 1/J, over 400 Monte Carlo simulations are displayed in Table 2. Here also, we used for each estimator the smoothing parameters (α n , h) minimizing its MSE over the grid of values A × H described above. When comparing the estimators qRP,1 4.2 Benchmark nonparametric estimators of γ(x) and q(β n |x)

Alternative modern smoothing techniques were discussed in e.g. [START_REF] Beirlant | Statistics of extremes: Theory and applications[END_REF]Section 7.5]. For comparison, we focus on the prominent local polynomial maximum likelihood estimation. This contribution fits a generalized Pareto (GP) model to the exceedances Z x i = Y j -u x given Y j > u x , for a high threshold u x , where j denotes the original index of the ith exceedance. Let N x be the number of all exceedances over u x and rearrange the indices of the explanatory variable such that X i denotes the covariate observation associated with exceedance Z x i . If g(z; σ, γ) stands for the GP density, then the local polynomial maximum likelihood approach maximizes the kernel weighted log-likelihood function

L Nx (β 1 , β 2 ) = 1 N x Nx i=1 log g   Z x i ; p 1 j=0 β 1j (X i -x) j , p 2 j=0 β 2j (X i -x) j   K h (X i -x) with respect to (β ′ 1 , β ′ 2 ) = (β 10 , • • • , β 1p 1 , β 20 , • • • , β 2p 2 )
to get the estimates σGP n (x) = β10 and γGP n (x) = β20 of the parameter functions σ(x) and γ(x) of the GP distribution fitted to the exceedances over u x . Note that local polynomial fitting also provides estimates of the derivatives of σ(x) and γ(x) up to order p 1 and p 2 , respectively. In order to not overload the estimation procedure, we confine ourselves to p 1 = p 2 = 0. The Monte Carlo results for γGP n (.) are reported in Table 3 (l-h.s). For each simulation, we used the parameters (h, u) that minimize the MSE{γ GP n (.)}, with the bandwidth ranging over the grid H described above and the threshold ranging over the αth sample quantiles of Y , where α ∈ A. The estimator γGP n has clearly smaller MSEs than the γRP n estimators in the Gaussian and Student error models, but it seems to be less efficient in the Beta error model than both γRP,1 n and γRP,2 n for J = 4 and r = 1/J. From a theoretical point of view, it should be clear that the pointwise asymptotic normality of γGP n (x) is proved in [START_REF] Beirlant | Local polynomial maximum likelihood estimation for Pareto-type distributions[END_REF] only in case γ(x) > 0. Moreover, the proof is restricted to the setting where the design points X i are deterministic. On the other hand, as suggested in [2, Section 7.5.2] and [START_REF] Beirlant | Local polynomial maximum likelihood estimation for Pareto-type distributions[END_REF], the extreme conditional quantile q(β n |x) can be estimated by

MSE{γ GP n } Bias{γ GP n } Gaussian 0.
qGP n (β n |x) := u x + σGP n (x) γGP n (x) n ⋆ xh β n k x -γ GP n (x) -1
where n ⋆ xh is the number of observations in [x -h, x + h] and k x is the number of exceedances receiving positive weight. Table 3 (r-h.s) reports the Monte Carlo estimates obtained by using in each simulation the parameters (h, u) that minimize the MSE{q GP n (β n |•)}, where h ∈ H and u ranges over the αth sample quantiles of Y with α ∈ A. In all cases, the regression quantile (RQ) estimator qn (β n |•) do appear to be more efficient than qGP n (β n |•). Compared with the qRP n (β n |•) estimators (for J = 3 and r = 1/J), qGP n (β n |•) seems to be more efficient only in the Gaussian error model for β n = 0.005 = 1/n, but not by much. A typical realization of the experiment in each simulated scenario is shown in Figure 1, where the smoothing parameters of each estimator were chosen in such a way to minimize its MSE. From a theoretical viewpoint, unlike our estimators qn (β n |x) and qRP n (β n |x), the asymptotic distribution of qGP n (β n |x) is not elucidated yet.

Data-driven rules for selecting the parameters h and α n

The use of the 'RQ' estimator qn (β n |x) ≡ qn (β n |x; h), which relies on the inversion of Fn (•|•), requires only the choice of the bandwidth h in an interval H of lower and upper bounds given respectively by, say, h min := max 1≤i<n X (i+1) -X (i) and h max := X (n) -X (1) /4. One way to select this parameter is by employing the cross-validation criterion as in [START_REF] Daouia | Kernel estimators of extreme level curves[END_REF] to obtain

h cv = arg min h∈H n i=1 n j=1 I(Y i ≥ Y j ) -Fn,-i (Y j |X i ) 2 ,
where Fn,-i (

•|•) is the estimator Fn (•|•) computed from the sample {(X j , Y j ), 1 ≤ j ≤ n, j = i}.
The empirical procedure of [START_REF] Yu | Local linear quantile regression[END_REF] could be used to get the alternative data-driven global bandwidth

h yj = h cv β n (1 -β n ) φ(Φ -1 (β n )) 2 1/5
, where φ and Φ stand respectively for the standard normal density and distribution functions. However, the use of the 'RP' estimators qRP,i n (β n |x) and γRP,i n (x), for i = 1, 2, requires in addition the selection of an appropriate order α n . To simplify the discussion we set α n at k/n ⋆ xh , where the integer k varies between 1 and n ⋆ xh -1, for each h ∈ H. We also consider the value J = 3 for which γRP,1 h,k) should one use to compute the estimates in practice could be the automatic ad hoc data driven-rule employed in [START_REF] Daouia | Frontier estimation and extreme value theory[END_REF]. The main idea is to evaluate first the estimates, for each x in a chosen grid of values, and then to select the parameters where the variation of the results is the smallest. This can be achieved in two ways:

Selecting h and k separately.

Step 1. Select a data-driven global bandwidth h, for example, h cv or h yj .

Step , as functions of k may be so unstable that reasonable values of k (which would correspond to the true value of q(β n |x), respectively γ(x)) may be hidden in the graphs. In result, the estimators may exhibit considerable volatility as functions of x itself.

A typical realization is shown in Figure 2 when the bandwidths h cv (left panels) and h yj (right panels) are used in Step 1. It may be seen that the method affords reasonable estimates in both Gaussian and Beta error models regarding the difficult curvature of the extreme quantile regression and the very small sample size n = 200. However, it seems that the method fails in the case of Student noise, where the superiority of qRP Selecting h and k simultaneously.

Step 1. For each h ∈ H, proceed to Step 2 described in the separate parameters' selection. Set the value of k where the standard deviation is minimal to be k xh and calculate the corresponding estimate qRP n (β n |x; h, k xh ).

Step 2. Compute the standard deviation of the estimates qRP n (β n |x; h, k xh ) over a window of (say 10) successive values of h. Select the bandwidth where the standard deviation is minimal and then evaluate the corresponding estimate.

In our simulations, we used a refined grid H of 50 points between min(h cv , h yj -h cv ) and h yj + 2h cv . Any other limit bounds of H could of course be chosen near h cv below and near h yj above. See Figure 3 for a typical realization in each simulated scenario. Here also the method is not without disadvantage as can be seen from the case of Student noise, where good results require a large sample size. selection methods except for the heavy-tailed case, corresponding to generally severe events. A sensible practice would be to verify whether the resulting 'RP', 'GP' and 'RQ' estimators point toward similar conclusions: the hard question of how to pick out the smoothing parameters (h, k) simultaneously in an optimal way might thus become less urgent. In contrast, if the estimators look clearly different, this might diagnose a heavy-tailed conditional distribution with a great variability in severity: thereby our technique might be viewed as an exploratory tool, rather than as a method for final analysis. Doubtlessly, further work to define a concept of selecting appropriate values for the crucial parameters (J, r) in the qRP n (β n |x) estimators will yield new refinements.

Concluding remarks

The case of multiple covariates. We have discussed the asymptotic distributional properties of both

(β n |•) in red with β n = 1/n. Its estimators qn (β n |•) = qn (β n |•, h cv ) in magenta, qRP,1 n (β n |•) ≡ qRP,2 n (β n |•)
in black with (r, J) = (1/3, 3), and qGP n (β n |•) in green. The (X i , Y i )'s are depicted as blue points.

estimators qn (β n |x) and qRP n (β n |x) in detail in Sections 2-3 for multiple regressors X ∈ R p , but our contributions are probably only of a theoretical value in the case p > 1. Indeed, as in the ordinary setting where the quantile order does not depend on the sample size, the kernel-smoothing method suffers from the 'curse of dimensionality'. In our setting of extreme quantile regression, the curse is exacerbated by several degrees of magnitude and drastically increases in higher dimensions. To overcome this vexing defect, one can use dimension reduction techniques such as ADE (Average Derivative Estimator), see for instance [START_REF] Hardle | Investigating smooth multiple regression by the method of average derivatives[END_REF]. Nevertheless, the theoretical properties of such methods are not yet established in the extreme-value framework.

Data example

Data on 123 American electric utility companies were collected and the aim is to investigate the economic efficiency of these companies (see, e.g., [START_REF] Greene | A gamma-distributed stochastic frontier model[END_REF]). A possible way to measure this efficiency is by looking at the maximum level of produced goods which is attainable for a given level of inputs-usage. From a statistical point of view, this problem translates into studying the upper boundary of the set of possible inputs X and outputs Y , the so-called cost/econometric frontier in production theory. Hendricks and Koenker [START_REF] Hendricks | Hierarchical spline for conditional quantiles and the demand for electricity[END_REF] stated: "In the econometric literature on the estimation of production technologies, there has been considerable interest in estimating so called frontier production models that correspond closely to models for extreme quantiles of a stochastic production surface". The present paper may be viewed as the first 'purely' nonparametric work to actually investigate theoretically the idea of Hendricks and Koenker.

For our illustration purposes, we used the measurements on the variable Y = log(Q), with Q being the production output of a firm, and the variable X = log(C), with C being the total cost involved in the production. Figure 4 shows the n = 123 observations, together with estimated extreme conditional quantiles qn (β n |x), qGP n (β n |x) and qRP,1

n (β n |x) = qRP,2
n (β n |x) at r = 1/J with J = 3. Given the small sample size, it was enough to use β n = 1/n in describing the conditional distribution tails. For selecting the window width h and the number k of extremes, we maintained the automatic empirical data-driven rules described above. It appears that the extreme conditional quantile estimates are similar for both simultaneous (top) and separate (bottom) selection methods. Following their evolution, it may be seen that the American electric utility data do not correspond to the situation hoped for by the practitioners of a heavily short-tailed production process. Indeed, one may distinguish between two different behaviors of the extreme regression quantiles: They indicate a short-tailed conditional distribution for companies working at (transformed) input-factors larger than, say, 2. In contrast, the tail distribution for the smallest companies with inputs X i < 2 seems to be moderately heavy. Therefore, the theoretical economic assumption that producers should operate on the upper boundary of the joint support of (X, Y ) rather than on its interior is clearly not fulfilled here, revealing a certain lack of production performance in this sector of activity. The estimated graph of qRP,1 n (β n |x), qn (β n |x) or qGP n (β n |x) might be interpreted as the set of the most efficient firms. It is then clear that the firms achieve significantly lesser output than that predicted by the extremal quantile frontiers. This indicates a relative economic inefficiency especially in the population of the (sparse) smallest companies in terms of inputs.

6 Appendix: Proofs

Preliminary results

We begin with a homogeneous property of the quantile function.

Lemma 1. Suppose (A.1) holds. If α n → 0 as n → ∞, then, lim n→∞ q(ξα n |x) q(α n |x) = ξ -(γ(x)∨0) ,
for all ξ > 0.

Proof. From (4), we have

q(ξα n |x) q(α n |x) = 1 + K γ(x) (1/ξ) a(q(α n |x)|x) q(α n |x) (1 + o(1))
and the conclusion follows using [START_REF] Chernozhukov | Extremal quantile regression[END_REF].

The following lemma states that the convergence in (3) is uniform.

Lemma 2. Under (A.1), if z n (x) ↑ y F (x)
as n → ∞, then for all sequence of functions t n (x) such that t n (x) → t 0 (x) as n → ∞ where t 0 (x) is such that there exists η > 0 for which 1 + γ(x)t 0 (x) ≥ η then, decreasing function, we have:

lim n→∞ F (z n (x) + t n (x)a(z n (x)|x)|x) F (z n (x)|x) = 1 K -1 γ(x) (t 0 (x)) . Proof. Since t n (x) → t 0 (x) as n → ∞, for all ε 1 > 0 such that |γ(x)|ε 1 < η, there exists N 1 ≥ 0 such that for all n ≥ N 1 , t 0 (x) -ε 1 ≤ t n (x) ≤ t 0 (x) + ε 1 . Since a(z n (x)|x) > 0 and F (.|x) is a -3 -2 -1 0 1 2 3 4 5
(z n (x) + (t 0 (x) + ε 1 )a(z n (x)|x)|x) F (z n (x)|x) ≤ F (z n (x) + t n (x)a(z n (x)|x)|x) F (z n (x)|x) ≤ F (z n (x) + (t 0 (x) -ε 1 )a(z n (x)|x)|x) F (z n (x)|x) . Now, since |γ(x)|ε 1 < η, it is easy to check that 1 + γ(x)(t 0 (x) + ε 1 ) ∧ 1 + γ(x)(t 0 (x) -ε 1 ) > 0.
Hence, under (A.1), for all ε 2 > 0, there exists

N 2 ≥ 0 such that for all n ≥ N 2 1 -ε 2 K -1 γ(x) (t 0 (x) + ε 1 ) ≤ F (z n (x) + t n (x)a(z n (x)|x)|x) F (z n (x)|x) ≤ 1 + ε 2 K -1 γ(x) (t 0 (x) -ε 1 )
.

Since ε 1 and ε 2 can be chosen arbitrarily small, this concludes the proof.

Let us remark that the kernel estimator (1) can be rewritten as Fn (y|x) = ψn (y, x)/ĝ n (x) where

ψn (y, x) = 1 n n i=1 K h (x -X i )I{Y i > y}
is an estimator of ψ(y, x) = F (y|x)g(x) and ĝn (x) is the classical kernel estimator of the pdf g(x) defined by: ĝn

(x) = 1 n n i=1 K h (x -X i ).
Lemma 3 gives standard results on the kernel estimator (see [START_REF] Parzen | On estimation of a probability density function and mode[END_REF] for a proof) whereas Lemma 4 is dedicated to the asymptotic properties of ψn (y, x).

Lemma 3. Suppose (A.2), (A.3) hold. If nh p → ∞, then, for all x ∈ R p such that g(x) > 0,

(i) E(ĝ n (x) -g(x)) = O(h), (ii) var(ĝ n (x)) = g(x) K 2 2 nh p (1 + o(1)).
Therefore, under the assumptions of the above lemma, ĝn (x) converges to g(x) in probability.

Let us introduce some further notations. In the following, y n (x) is a sequence such that y n (x) ↑ y F (x) and

y n,j (x) = y n (x) + K γ(x) (1/τ j )a(y n |x)(1 + o(1)) for all j = 1, . . . , K. Recall that 0 < τ J < • • • < τ 2 < τ 1 ≤ 1.
Moreover, the oscillations of the csf are controlled by

ω n (x) := max j=1,...,J sup x ′ ∈B(x,h) F (y n,j (x)|x ′ ) F (y (x)|x) -1 .
Lemma 4. Suppose (A.1), (A.2) and (A.3) hold and let

x ∈ R p such that g(x) > 0. If ω n (x) → 0 and nh p F (y n (x)|x) → ∞ then, (i) E( ψn (y n,j (x), x)) = ψ(y n,j (x), x) (1 + O(ω n (x)) + O(h)), for j = 1, . . . , J.
(ii) The random vector nh p ψ(y n (x), x) ψn (y n,j (x), x) -E( ψn (y n,j (x), x)) ψ(y n,j (x), x) j=1,...,J is asymptotically Gaussian, centered, with covariance matrix K 2 2 V where V j,j ′ = τ -1 j∧j ′ for (j, j ′ ) ∈ {1, . . . , J} 2 .

Proof. (i) Since the (X i , Y i ), i = 1, . . . , n are identically distributed, we have

E( ψn (y n,j (x), x)) = R p K h (x -t) F (y n,j (x)|t)g(t)dt = S K(u) F (y n,j (x)|x -hu)g(x -hu)du, under (A.3). Let us now consider |E( ψn (y n,j (x), x))-ψ(y n,j (x), x)| ≤ F (y n,j (x)|x) S K(u)|g(x -hu) -g(x)|du (8)
+ F (y n,j (x)|x)

S K(u) F (y n,j (x)|x -hu) F (y n,j (x)|x) -1 g(x-hu)du. (9)
Under (A.2), and since g(x) > 0, we have

(8) ≤ F (y n,j (x)|x)c g h S d(u, 0)K(u)du = ψ(y n,j (x), x)O(h), (10)
and, in view of [START_REF] Drees | Refined Pickands estimators of the extreme value index[END_REF],

(9) ≤ F (y n,j (x)|x)ω n (x) S K(u)g(x -hu)du = F (y n,j (x)|x)ω n (x)g(x)(1 + o(1)) = ψ(y n,j (x), x)ω n (x)(1 + o(1)). (11)
Combining [START_REF] Drees | Refined Pickands estimators of the extreme value index[END_REF] and [START_REF] Feigin | Linear programming estimators and bootstrapping for heavy tailed phenomena[END_REF] concludes the first part of the proof.

(ii) Let β = 0 in R J , Λ n (x) = (nh p ψ(y n (x), x)) -1/2 , and consider the random variable

Ψ n = J j=1
β j ψn (y n,j (x), x) -E( ψn (y n,j (x), x))

Λ n (x)ψ(y n,j (x), x)

= n i=1 1 nΛ n (x)    J j=1 β j K h (x -X i )I{Y i ≥ y n,j (x)} ψ(y n,j (x), x) -E   J j=1 β j K h (x -X i )I{Y i ≥ y n,j (x)} ψ(y n,j (x), x)      := n i=1 Z i,n .
Clearly, {Z i,n , i = 1, . . . , n} is a set of centered, independent and identically distributed random variables with variance

var(Z i,n ) = 1 n 2 h 2p Λ 2 n (x) var   J j=1 β j K x -X i h I{Y i ≥ y n,j (x)} ψ(y n,j (x), x)   = 1 n 2 h p Λ 2 n (x) β t Bβ,
where B is the J × J covariance matrix with coefficients defined for (j, j ′ ) ∈ {1, . . . , J} 2 by

B j,j ′ = A j,j ′ ψ(y n,j (x), x)ψ(y n,j ′ (x), x) , A j,j ′ = 1 h p cov K x -X h I{Y ≥ y n,j (x)}, K x -X h I{Y ≥ y n,j ′ (x)} = K 2 2 E 1 h p Q x -X h I{Y ≥ y n,j (x) ∨ y n,j ′ (x)} -h p E(K h (x -X)I{Y ≥ y n,j (x)})E(K h (x -X)I{Y ≥ y n,j ′ (x)}),
with Q(.) := K 2 (.)/ K 2 2 also satisfying assumption (A.3). As a consequence, the three above expectations are of the same nature. Thus, remarking that, for n large enough, y n,j (x) ∨ y n,j ′ (x) = y n,j∨j ′ (x), part (i) of the proof implies

A j,j ′ = K 2 2 ψ(y n,j∨j ′ (x), x) (1 + O(ω n (x)) + O(h)) -h p ψ(y n,j (x), x)ψ(y n,j ′ (x), x) (1 + O(ω n (x)) + O(h)) leading to B j,j ′ = K 2 2 ψ(y n,j∧j ′ (x), x) (1 + O(ω n (x)) + O(h)) -h p (1 + O(ω n (x)) + O(h)) = K 2 2 ψ(y n,j∧j ′ (x), x) (1 + o(1)), since ψ(y n,j∧j ′ (x), x) → 0 as n → ∞. Now, from Lemma 2, lim n→∞ ψ(y n,j∧j ′ (x), x) ψ(y n (x), x) = 1 K -1 γ(x) (K γ(x) (1/τ j∧j ′ )) = τ j∧j ′ entailing B j,j ′ = K 2 2 V j,j ′ ψ(y n (x), x) (1 + o(1)),
and therefore, var(Z i,n ) ∼ K 2 2 β t V β/n, for all i = 1, . . . , n. As a preliminary conclusion, the variance of Ψ n converges to K 2 2 β t V β. Consequently, Lyapounov criteria for the asymptotic normality of sums of triangular arrays reduces to

n i=1 E |Z i,n | 3 = nE |Z 1,n | 3 → 0. Remark that Z 1,n is a bounded random variable: |Z 1,n | ≤ 2 K ∞ J j=1 |β j | nΛ n (x)h p ψ(y n,J , x) = 2 τ J K ∞ J j=1 |β j |Λ n (x)(1 + o(1))
and thus,

nE |Z 1,n | 3 ≤ 2 τ J K ∞ J j=1 n|β j |Λ n (x)var(Z 1,n )(1 + o(1)) = 2 τ J K ∞ K 2 2 J j=1 |β j |β t V βΛ n (x)(1 + o(1)) → 0
as n → ∞. As a conclusion, Ψ n converges in distribution to a centered Gaussian random variable variance K 2 2 β t V β for all β = 0 in R p . The result is proved. Let us now focus on the estimation of small tail probabilities F (y n (x)|x) when y n (x) ↑ y F (x) as n → ∞. The following result provides sufficient conditions for the asymptotic normality of Fn (y n (x)|x). is asymptotically Gaussian, centered, with covariance matrix K 2 2 /g(x)V where V j,j ′ = τ -1 j∧j ′ for (j, j ′ ) ∈ {1, . . . , J} 2 .

Proof of Proposition 1. Keeping in mind the notations of Lemma 4, the following expansion holds

Λ -1 n (x) J j=1 β j Fn (y n,j (x)|x) F (y n,j (x)|x) -1 = T 1,n + T 2,n -T 3,n ĝn (x) , (12)
where

T 1,n = g(x)Λ -1 n (x) J j=1
β j ψn (y n,j (x), x) -E( ψn (y n,j (x), x)) ψ(y n,j (x), x)

T 2,n = g(x)Λ -1 n (x) J j=1
β j E( ψn (y n,j (x), x)) -ψ(y n,j (x), x) ψ(y n,j (x), x)

T 3,n =   J j=1 β j   Λ -1 n (x) (ĝ n (x) -g(x)) .
Let us highlight that assumptions nh p F (y n (x)|x)ω 2 n (x) → 0 and nh p F (y n (x)|x) → ∞ imply that ω n (x) → 0 as n → ∞. Thus, from Lemma 4(ii), the random term T 1,n can be rewritten as

T 1,n = g(x) K 2 β t V βξ n , (13)
where ξ n converges to a standard Gaussian random variable. The nonrandom term T 2,n is controlled with Lemma 4(i):

T 2,n = O(Λ -1 n (x)(h + ∆(y n (x), x)) = O (nh p F (y n (x)|x)) 1/2 (h ∨ ω n (x)) = o(1). (14)
Finally, T 3,n is a classical term in kernel density estimation, which can be bounded by Lemma 3:

T 3,n = O(hΛ -1 n (x)) + O P (Λ -1 n (x)(nh p ) -1/2 ) = O nh p+2 F (y n (x)|x) 1/2 + O P ( F (y n (x)|x)) 1/2 = o P (1). (15)
Collecting ( 12)-( 15), follows that

ĝn (x)Λ -1 n (x) J j=1 β j Fn (y n,j (x)|x) F (y n,j (x)|x) -1 = g(x) K 2 β t V βξ n + o P (1).
Finally, ĝn (x)

P -→ g(x) yields nh p F (y n (x)|x) J j=1 β j Fn (y n,j (x)|x) F (y n,j (x)|x) -1 = K 2 β t V β g(x) ξ n + o P (1)
and the result is proved.

The last lemma establishes that K γn(x) (r n ) inherits from the convergence properties of γn (x).

Lemma 5. Suppose ξ

(γ) n (x) := Λ -1 n (γ n (x) -γ(x)) = O P (1), where Λ n → 0. Let r n ≥ 1 or r n ≤ 1 such that Λ n log(r n ) → 0. Then, Λ -1 n K γn(x) (r n ) -K γ(x) (r n ) K ′ γ(x) (r n ) = ξ (γ) n (x)(1 + o P (1)).
Proof. Since γn (x) P -→ γ(x), a first order Taylor expansion yields

K γn(x) (r n ) = K γ(x) (r n ) + Λ n ξ (γ) n (x)K ′ γn(x) (r n ), where γn (x) = γ(x) + Θ n Λ n ξ (γ) n (x) with Θ n ∈ (0, 1). As a consequence Λ -1 n K γn(x) (r n ) -K γ(x) (r n ) K ′ γ(x) (r n ) = ξ (γ) n (x) K ′ γn(x) (r n ) K ′ γ(x) (r n ) = ξ (γ) n (x) 1 + rn 1 (s γn(x)-γ(x) -1)s γ(x)-1 log(s)ds
|(s γn(x)-γ(x) -1)| ≤ 2(log r n )|γ n (x) -γ(x)| = o P (1). As a conclusion, Λ -1 n K γn(x) (r n ) -K γ(x) (r n ) K ′ γ(x) (r n ) = ξ (γ) (x)(1 + o P (1))
and the result is proved. The case r n ≤ 1 is easily deduced since

K γ(x) (1/r n ) = -K -γ(x) (r n ) and K ′ γ(x) (1/r n ) = K ′ -γ(x) (r n ).

Proofs of main results

Proof of Theorem 1. Let us introduce v n = (nh p α -1 n ) 1/2 , σ n (x) = (v n f (q(α n |x)|x)) -1 and, for all j = 1, . . . , J, W n,j (x) = v n Fn (q(τ j α n |x) + σ n (x)z j |x) -F (q(τ j α n |x) + σ n (x)z j |x) a n,j (x) = v n τ j α n -F (q(α n,j |x) + σ n (x)z j |x) and z j ∈ R. We examine the asymptotic behavior of J-variate function defined by

Φ n (z 1 , . . . , z J ) = P   J j=1 σ -1 n (x)(q n (τ j α n |x) -q(τ j α n |x)) ≤ z j   = P   J j=1 {W n,j (x) ≤ a n,j (x)}   .
Let us first focus on the nonrandom term a n,j (x). For each j ∈ {1, . . . , J} there exists θ n,j ∈ (0, 1) such that

a n,j (x) = v n σ n (x)z j f (q n,j (x)|x) = z j f (q n,j (x)|x) f (q(α n |x)|x) ,
where

q n,j (x) = q(τ j α n |x) + θ n,j σ n (x)z j = q(τ j α n |x) + θ n,j z j τ j (nh p α n ) -1/2 τ j α n f (q(τ j α n |x)|x) f (q(τ j α n |x)|x) f (q(α n |x)|x) = q(τ j α n |x) + θ n,j z j τ γ(x) j (nh p α n ) -1/2 τ j α n f (q(τ j α n |x)|x) (1 + o(1)), since y → f (q(y|x)
(t|x) = F (t|x) f (t|x) (1 + o(1)), (16)
leading to

q n,j (x) = q(τ j α n |x) + θ n,j z j τ γ(x) j (nh p α n ) -1/2 a(q(τ j α n |x)|x)(1 + o(1)). Lemma 2 with z n (x) = q(τ j α n |x), t n (x) = θ n,j z j τ γ(x) j (nh p α n ) -1/2 (1 + o( 1
)) and t 0 (x) = 0 yields F (q n,j (x)|x)

τ j α n → K -1 γ(x) (0) = 1
as n → ∞. Recalling that y → f (q(y|x)|x) is regularly varying, we have

f (q n,j (x)|x) f (q(α n |x)|x) → τ γ(x)+1 j
as n → ∞ and therefore a n,j (x) = z j τ γ(x)+1 j

(1 + o(1)), j = 1, . . . , J.

Let us now turn to the random term W n,j (x). Let us define z n,j (x) = q(τ j α n |x) + σ n (x)z j for j = 1, . . . , J, y n (x) = q(α n |x), and consider the expansion

z n,j (x) -y n (x) a(y n (x)|x) = q(τ j α n |x) -q(α n |x) a(q(α n |x)|x) + σ n (x)z j a(q(α n |x)|x)
.

From (4), we have

lim n→∞ q(τ j α n |x) -q(α n |x) a(q(α n |x)|x) = K γ(x) (1/τ j ), and
lim n→∞ σ n (x)z j a(q(α n |x)|x) = 0, leading to z n,j (x) = y n (x) + K γ(x) (1/τ j )a(y n (x)|x)(1 + o(1)
). Introducing β n,j (x) = F (y n,j (x)|x), the oscillation ω n (x) can be rewritten as

ω n (x) = max j=1,...,J sup x ′ ∈B(x,h) F (q(β n,j (x)|x)|x ′ ) β n,j (x) -1 .
For all κ ∈ (0, τ J ) and j = 1, . . . , J, we eventually have z n,j (x) ∈ [y n (x), z n (x)] where z n (x) :=

y n (x) + K γ(x) (2/
W n,j (x) = F (z n,j (x)|x) α n ξ n,j = τ j ξ n,j (1 + o(1))
where ξ n = (ξ n,1 , . . . , ξ n,j ) t converges to a centered Gaussian random vector with covariance matrix K 2 2 /g(x)V . Taking into account of ( 17), the results follows.

Proof of Corollary 1. Let us remark that, from (16), f (q(α n |x)|x) nh p α -1 n (q n (τ j α n |x) -q(τ j α n |x)) j=1,...,J = q(α n |x) f (q(α n |x)|x) α n q(τ j α n |x) q(α n |x) nh p α n qn (τ j α n |x) q(τ j α n |x) -1 j=1,...,J = q(α n |x) a(q(α n |x)|x) q(τ j α n |x) q(α n |x) nh p α n qn (τ j α n |x) q(τ j α n |x) -1 j=1,...,J Proof of Theorem 2. By definition, q n (β n |x) = q n (α n |x) + (K γ(x) (α n /β n ) + b(β n /α n , α n ))a(q(α n |x)|x)

and thus, the following expansion can be easily established:

Λ -1 n qn (β n |x) -q(β n |x) a(q(α n |x)|x)K ′ γ(x) (α n /β n ) = Λ -1 n qn (α n |x) -q(α n |x) a(q(α n |x)|x)K ′ γ(x) (α n /β n ) + Λ -1 n K γn(x) (α n /β n ) -K γ(x) (α n /β n ) K ′ γ(x) (α n /β n ) ân (x) a(q(α n |x)|x) + Λ -1 n K γ(x) (α n /β n ) K ′ γ(x) (α n /β n ) ân (x) a(q(α n |x)|x) -1 + Λ -1 n b(β n /α n , α n ) K ′ γ(x) (α n /β n ) =: T n,1 + T n,2 + T n,3 + T n,4 .
see [START_REF] De Haan | Extreme value theory: An introduction[END_REF], it follows that qn (τ j α n |x) -qn (τ j+1 α n |x) a(q(α n |x)|x) = ε n (ξ j,n -ξ j+1,n ) + K γ(x) (1/τ j ) -K γ(x) (1/τ j+1 ) + b(τ j , α n |x) -b(τ j+1 , α n |x)

= ε n (ξ j,n -ξ j+1,n ) + K γ(x) (r)r -γ(x)j + O(η n (x)). [START_REF] Jurecková | Remark on extreme regression quantile[END_REF] Replacing in [START_REF] Hendricks | Hierarchical spline for conditional quantiles and the demand for electricity[END_REF], we obtain (log r)γ n,j (x) = log ε n (ξ j,n -ξ j+1,n ) + K γ(x) (r)r -γ(x)j + O(η n (x))

ε n (ξ j+1,n -ξ j+2,n ) + K γ(x) (r)r -γ(x)(j+1) + O(η n (x)) , or equivalently, (log r)(γ n,j (x) -γ(x)) = log 1 + ε n (ξ j,n -ξ j+1,n )r γ(x)j K γ(x) (r) + O(η n (x))

-log 1 + ε n (ξ j+1,n -ξ j+2,n )r γ(x)(j+1) K γ(x) (r) + O(η n (x))

A first order Taylor expansion yields (log r)ε -1 n (γ n,j (x) -γ(x)) = r γ(x)j K γ(x) (r) ξ j,n -(1 + r γ(x) )ξ j+1,n + r γ(x) ξ j+2,n + O ( ε -1 n η n (x)) + o P (1) thus, under the assumption (nh p α n ) 1/2 η n (x) → 0 as n → ∞, nh p α n (γ RP n (x)-γ(x)) = 1 (log r)K γ(x) (r) J-2 j=1 π j r γ(x)j ξ j,n -(1 + r γ(x) )ξ j+1,n + r γ(x) ξ j+2,n +o P (1).

Defining for the sake of simplicity π -1 = π 0 = π J-1 = π J = 0, β

ii) Second, let us now consider a n,j (x) = r γRP n (x)j (q n (τ j α n |x) -qn (τ j+1 α n |x)) K γRP n (x) (r) such that ân (x) = J-2 j=1 π j a n,j (x). From [START_REF] Jurecková | Remark on extreme regression quantile[END_REF], it follows that, for all j = 1, . . . , J, a n,j (x) a(q(α n |x)|x) = r γRP n (x)j K γRP n (x) (r) ε n (ξ j,n -ξ j+1,n ) + K γ(x) (r)r -γ(x)j + O(η n (x)) .

Remarking that γRP n (x) = γ(x) + (nh p α n ) -1/2 ξ (γ) n (x), Lemma 5 yields a n,j (x) a(q(α n |x)|x) = 1 + r γ(x)j K γ(x) (r) ε n (ξ j,n -ξ j+1,n ) + O(η n (x))

1 + A first order Taylor expansion thus gives nh p α n a n,j (x) a(q(α n |x)|x) -1 = ξ (γ) n (x) j log(r) -K ′ γ(x) (r) K γ(x) (r) + r γ(x)j K γ(x) (r) (ξ j,n -ξ j+1,n ) + O nh p α n η n (x) + o P (1).

K ′ γ(x) (r) K γ(x) (r) (nh p α n ) -1/2 ξ
Recalling that π -1 = π 0 = π J-1 = π J = 0 and introducing

E(π) = J j=1 jπ j , β (a) 1
= -r -γ(x) -(r -γ(x) + 1) E(π) -

in view of [START_REF] Koenker | Reappraising medfly longevity: A quantile regression survival analysis[END_REF]. iii) Third, Corollary 1 states that

ξ 1,n = √ nh p α n a(q(α n |x)|x) (q n (α n |x) -q(α n |x)) (26)
is asymptotically Gaussian. Finally, collecting [START_REF] Koenker | Reappraising medfly longevity: A quantile regression survival analysis[END_REF], ( 25) and ( 26

Subsection 4 .

 4 1 provides Monte Carlo evidence that the extreme quantile function estimator qRP,1 n (β n |x) is efficient relative to the version qRP,2 n (β n |x), whether γ(x) is positive, negative or zero, and outperforms the estimator qn (β n |x) for heavy-tailed conditional distributions. Subsection 4.2 provides a comparison with the promising local smoothing approach introduced in [1] and [2, Section 7.5.2]. Practical guidelines for selecting the bandwidth h and the order α n are suggested in Subsection 4.3.

Figure 1 :

 1 Figure 1: Typical realizations for simulated samples of size n = 200. From left to right and from top to bottom, Y |X is Gaussian, Student, Beta. The true quantile function q(β n |•) in red with β n = 1/n. Its estimators qn (β n |•) in magenta, qRP,1 n (β n |•) ≡ qRP,2 n (β n |•) in black with r = 1/J and J = 3, and qGP n (β n |•) in green. The observations (X i , Y i ) are depicted as blue points.

n

  (x) ≡ γRP,2 n (x) := γRP n (x; h, k) and qRP,1 n (β n |x) ≡ qRP,2 n (β n |x) := qRP n (β n |x; h, k), with r = 1/J. An empirical way to decide what values of (

  n (β n |•) over both qGP n (β n |•) and qn (β n |•), demonstrated via the Monte Carlo study, is clearly sacrificed. This failure is probably due to the arbitrary choice (3, 1/3) of the parameters (J, r) in qRP n (β n |•). It might also be seen that, apart from the student error model, the three estimators qRP n (β n |•), qn (β n |•) and qGP n (β n |•) point toward similar results.

Figure 2 : 1 n

 21 Figure 2: Separate parameters' selection using the same simulated samples as in Figure 1. From left to right, the used bandwidth is h cv , h yj . From top to bottom, Y |X is Gaussian, Student, Beta. The true quantile function q(β n |•) in red with β n = 1/n. Its estimators qn (β n |•) in magenta, qRP,1 n (β n |•) ≡ qRP,2 n (β n |•) in black with r = 1/J and J = 3, and qGP n (β n |•) in green. The observations (X i , Y i ) are depicted as blue points.

Figure 3 :

 3 Figure 3: Simultaneous parameters' selection using the same simulated samples as above. From left to right and from top to bottom, Y |X is Gaussian, Student, Beta. The true quantile function q(β n |•) in red with β n = 1/n. Its estimators qn (β n |•) = qn (β n |•, h cv ) in magenta, qRP,1

Figure 4 :

 4 Figure 4: Scatterplot of the American electric utility data (blue points). Results obtained via the simultaneous (top) and separate (bottom) selection methods for β n = 1/n: the estimator qn (β n |•) is drawn in magenta with h yj (top) and h cv (bottom), qGP n (β n |•) in green and qRP,1 n (β n |•) in black with (r, J) = (1/3, 3).

Proposition 1 .

 1 Suppose (A.1), (A.2) and (A.3) hold and letx ∈ R p such that g(x) > 0. If nh p F (y n (x)|x) → ∞ and nh p F (y n (x)|x)(h ∨ ω n (x)) 2 → 0, then, the random vector nh p F (y n (x)|x)Fn (y n,j (x)|x) F (y n,j (x)|x) -1 j=1,...,J

rn 1 s

 1 γ(x)-1 log(s)ds . Suppose for instance r n ≥ 1. The assumptions yield (log r n )(γ n (x) -γ(x)) P -→ 0 and thus, for n large enough, with high probability, sup s∈[1,rn]

( 1 +

 1 o(1)) = q(α n |x) a(q(α n |x)|x) τ -(γ(x)∨0) j nh p α n qn (τ j α n |x) q(τ j α n |x) -1 j=1,...,J (1 + o(1)),in view of Lemma 1. The result follows from Theorem 1.

2 =

 2 r -γ(x) log(r) , we end up withξ (γ) n (x) := nh p α n (γ RP n (x) -γ(x)) j-2 ξ j,n + o P (1).

n

  (x)(1 + o P (1)) exp ξ (γ) n (x)j log(r)(nh p α n ) -1/2 .

2 = 1 =

 21 r -γ(x) E(π) -K ′ γ(x) (r) log(r)K γ(x) (r) it follows that ξ (a) n (x) := nh p α n a n (x) a(q(α n |x)|x) -E(π) log(r) -K ′ γ(x) (r) K γ(x) (r) ξ (γ)(x) n + 1 K γ(x) (r) J j=1 r γ(x)j (π j -r -γ(x) π j-1 )ξ j,n + o P (1) = 1 K γ(x) (r) J j=1 r γ(x)j β (a) 0 π j + β (a) 1 π j-1 + β (a)2 π j-2 ξ j,n + o P[START_REF] Beirlant | Local polynomial maximum likelihood estimation for Pareto-type distributions[END_REF]

2 A 2

 22 ), (ξ(γ) n (x), ξ (a)n (x), ξ 1,n ) t = 1 K γ(x) (r) A(x)ξ n + o P (1),where A(x) is the 3 × J matrix defined byA 1,j (x) = r γ(x)j β ,j (x) = r γ(x)j β j (x) = K γ(x) (r)I{j = 1},for all j = 1, . . . , J. It is thus clear that the random vector (ξ(γ) n (x), ξ

n 2 2g

 2 (x), ξ 1,n ) t converges in distribution to a centered Gaussian random vector with covariance matrixK (x)K 2 γ(x) (r)A(x)Σ(x)A t (x) =: S(x).

Table 1 :

 1 Performance of γRP,1 n

				r = 1/J							r = (J -1)/J
			MSE			Bias			MSE			Bias
		γRP,1 n	(x)	γRP,2 n	(x)	γRP,1 n	(x)	γRP,2 n	(x)	γRP,1 n	(x)	γRP,2 n	(x)	γRP,1 n	(x)	γRP,2 n	(x)
	Gaussian
	J=3	0.2026	0.2026	-0.2415	-0.2415	0.7656	0.7656	-0.3213	-0.3213
	J=4	0.1915	0.2018	-0.3270	-0.3501	0.6730	0.7960	-0.3455	-0.3747
	J=5	NaN	NaN	NaN	NaN	0.7305	0.9128	-0.4104	-0.4107
	Student
	J=3	0.2882	0.2882	-0.2964	-0.2964	1.1109	1.1109	-0.4497	-0.4497
	J=4	0.3350	0.2837	-0.4167	-0.3480	0.9991	1.1997	-0.4384	-0.4597
	J=5	NaN	NaN	NaN	NaN	1.1245	1.3331	-0.5715	-0.5872
	Beta
	J=3	0.1157	0.1157	-0.0730	-0.0730	0.6737	0.6737	-0.2591	-0.2591
	J=4	0.0510	0.0597	-0.0811	-0.0750	0.5861	0.6891	-0.2338	-0.2432
	J=5	NaN	NaN	NaN	NaN	0.6431	0.8167	-0.2185	-0.2757

Table 2 :

 2 n |x) themselves with qn (β n |x), the results (both in terms of MSE and bias) indicate that qRP,2 n |x) in all cases, and that the latter is appreciably more efficient than qn (β n |x) only in the Student error model. It may be also noticed that qn (β n |x) is more efficient but not by much (especially when J = 3) in the Gaussian and Beta error models. Performance of qRP,1

						n	(β n |x) is slightly less efficient than qRP,1
					r = 1/J, βn = 0.05
				MSE					Bias
		qRP,1 n	(βn|x)	qRP,2 n	(βn|x)	qn(βn|x)	qRP,1 n	(βn|x)	qRP,2 n	(βn|x)	qn(βn|x)
	Gaussian
	J=3	0.0110	0.0110	0.0108	0.0001	0.0001	0.0063
	J=4	0.0591	0.0796	0.0108	0.1136	0.1131	0.0063
	Student
	J=3	0.0307	0.0307	0.0771	-0.0134	-0.0134	0.0871
	J=4	0.0532	0.0743	0.0771	0.0792	0.0792	0.0871
	Beta
	J=3	0.0091	0.0091	0.0022	0.0505	0.0505	0.0135
	J=4	0.0745	0.1002	0.0022	0.1746	0.1752	0.0135
					r = 1/J, βn = 0.01
				MSE					Bias
		qRP,1 n	(βn|x)	qRP,2 n	(βn|x)	qn(βn|x)	qRP,1 n	(βn|x)	qRP,2 n	(βn|x)	qn(βn|x)
	Gaussian
	J=3	0.0265	0.0265	0.0161	-0.0776	-0.0776	-0.0360
	J=4	0.0693	0.0926	0.0161	0.1092	0.1225	-0.0360
	Student
	J=3	0.1115	0.1115	0.6825	-0.0895	-0.0895	-0.0959
	J=4	0.1304	0.3992	0.6825	0.0018	0.1089	-0.0959
	Beta
	J=3	0.0143	0.0143	0.0034	0.0523	0.0523	0.0212
	J=4	0.1038	0.1265	0.0034	0.1964	0.2064	0.0212
					r = 1/J, βn = 0.005
				MSE					Bias
		qRP,1 n	(βn|x)	qRP,2 n	(βn|x)	qn(βn|x)	qRP,1 n	(βn|x)	qRP,2 n	(βn|x)	qn(βn|x)
	Gaussian
	J=3	0.0354	0.0354	0.0203	-0.0981	-0.0981	-0.0524
	J=4	0.0719	0.0932	0.0203	0.0982	0.1073	-0.0524
	Student
	J=3	0.2919	0.2919	0.9782	-0.1623	-0.1623	-0.2605
	J=4	0.4569	0.9748	0.9782	-0.1920	0.0280	-0.2605
	Beta
	J=3	0.0155	0.0155	0.0038	0.0536	0.0536	0.0239
	J=4	0.1130	0.1337	0.0038	0.1871	0.2111	0.0239

n (β n |x) and qRP,2 n (β n (β n (β n |x), qRP,2 n (β n |x) and qn (β n |x) with β n = 0.05 (top), β n = 0.01 (middle) and β n = 0.005 (bottom) -Results averaged on 400 simulations with n = 200.

Table 3 :

 3 Performance of γGP

n and qGP n (β n |•) -Results averaged on 400 simulations with n = 200.

  2. Evaluate qRP n (β n |x; h, k) at k = 1, . . . , n ⋆ xh -1. Then compute the standard deviation of the estimates over a 'window' of (say [ n ⋆ xh ]) successive values of k. The value of k where this standard deviation is minimal defines the desired parameter. the Y i 's such that |X i -x| ≤ h. The main difficulty when employing such a separate choice of h and k is that both qRP n (β n |x; h cv , k) and qGP n (β n |x; h cv , k), respectively γRP

	The same considerations evidently apply to γRP n (x) and to the 'benchmark' estimators γGP n (x) :=
	γGP n (x; h, k) and qGP n (β n |x) := qGP n (β n |x; h, k), defined in Subsection 4.2, with the covariate depen-dent threshold being u x := Y x (n ⋆ xh -k) , and Y x (1) ≤ • • • ≤ Y x (n ⋆ xh ) being the sequence of ascending order statistics corresponding to

n (x; h cv , k) and γGP n (x; h cv , k)

  Monte Carlo evidence. The experiments indicate that qn (β n |x) is efficient relative to the modern smoothing estimator qGP n (β n |x) first introduced in [1, 2]. The simulations also indicate that the performance of the alternative estimator qRP,1

					Conditional 0.995 quantiles (Gaussian noise)							Conditional 0.995 quantiles (Gaussian noise)
		2			True quantile									1
					RP estimate
					RQ estimate									0.8
		1.5			GP estimate
														0.6
		1												0.4
														0.2
	Yi	0.5											Yi
														0
		0												-0.2
														-0.4
		-0.5																					True quantile
														-0.6									RP estimate
																							RQ estimate
																							GP estimate
		-1	0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1	-0.8	0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1
								Xi												Xi
					Conditional 0.995 quantiles (Student noise)							Conditional 0.995 quantiles (Student noise)
		3.5			True quantile									3.5			True quantile
		3			RP estimate RQ estimate									3			RP estimate RQ estimate
					GP estimate												GP estimate
		2.5												2.5
		2												2
		1.5												1.5
	Yi	1											Yi	1
		0.5												0.5
		0												0
		-0.5												-0.5
		-1												-1
		-1.5	0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1	-1.5	0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1
								Xi												Xi
						Conditional 0.995 quantiles (Beta noise)							Conditional 0.995 quantiles (Gaussian noise)
		1												1.2
														1
														0.8
		0.5												0.6
														0.4
	Yi												Yi
														0.2
		0												0
														-0.2
											True quantile											True quantile
											RP estimate		-0.4									RP estimate
											RQ estimate											RQ estimate
											GP estimate											GP estimate
		-0.5	0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1	-0.6	0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1
								Xi												Xi
	n (β n |x), at least in terms of MSE. In comparison with qn (β n |x) and qGP (β n |x) is quite remarkable in comparison with its n (β n |x), the variability and the bias of both qRP,1 analog qRP,2 n n (β n |x) and qRP,2 n (β n |x) are quite respectable and it seems that the heavier is the conditional tail, the better the estimators qRP n (β n |x) are. It should be also clear that qRP,1

n (β n |x) and qRP,2 n (β n |x) can be improved appreciably by tuning the choice of the parameters J and r.

Tuning parameters selection in practice. The simulations worked reasonably well for our 'ad hoc'

  κ)a(y n (x)|x) and thusβ n,j (x) ∈ [ F (z n (x)|x), α n ] eventually. Now, Lemma 2 implies that F (z n (x)|x)/α n → κ/2 as n → ∞ and thus, for n large enough, β n,j (x) ∈ [κα n , α n ].Consequently, ω n (x) ≤ ∆ κ (α n , x). Applying Proposition 1 and Lemma 2 yields

Acknowledgements. The authors thank the Editor and both anonymous Associate Editor and reviewer for their valuable feedback on this paper.

Introducing

(ξ (γ) n (x), ξ (a) n (x), ξ (q) n (x)) := Λ -1 n γn (x) -γ(x), ân (x) a(q(α n |x)|x) -1, qn (α n |x) -q(α n |x) a(q(α n |x)|x) , from and remarking that, when u → ∞,

the first term can be rewritten as

Second, Λ n → 0 and (7) entail ân (x)/a(q(α n |x)|x) P -→ 1 and thus

from Lemma 5. From ( 7), [START_REF] Greene | A gamma-distributed stochastic frontier model[END_REF], and in view of

the third term can be rewritten as

Finally, T n,4 = o P (1) by assumption and the conclusion follows from ( 7), [START_REF] De Haan | Extreme value theory: An introduction[END_REF], [START_REF] Hall | Nonparametric analysis of temporal trend when fitting parametric models to extreme-value data[END_REF] and [START_REF] Hardle | Investigating smooth multiple regression by the method of average derivatives[END_REF].

Proof of Theorem 3. The proof consists in deriving asymptotic expansions for the three considered random variables. i) Let us first introduce

such that γRP n (x) = J-2 j=1 π j γ n,j (x). From Theorem 1 and in view of (4), we have, for all j = 1, . . . , J,

n and where the random vector ξ n = (ξ j,n ) j=1,...,J is asymptotically Gaussian, centered, with covariance matrix K


Sofiane Moad
email: sofiane.moad@irisa.fr

Mohammad Abdul Azim

Nizar Bouabdallah
email: nizar.bouabdallah@inria.fr

Rami Langar
email: rami.langar@lip6.fr

CDEEC: A Connectivity Degree-Based Energy Efficient Clustering Protocol for Wireless Sensor Networks

Keywords: Wireless sensor networks, energy efficiency, clustering

In this paper we propose CDEEC, a Connectivity Degree-based and Energy Efficient Clustering protocol, for wireless sensor networks in order to achieve further energy conservation while extending the network lifetime. The main idea behind our proposal is to consider, in addition to the information about the nodes' residual energy, the specific topology characteristics of the network, such as the connectivity degree, in the selection process of the cluster heads. To evaluate the performance of our proposal, we compare it to classical HEED-based (Hybrid, Energy-Efficient Distributed Clustering) networks. Simulation results show that significant energy conservation can be achieved.

I. INTRODUCTION

Wireless sensor networks (WSNs) consist of tiny nodes capable of sensing, computing and communicating [START_REF] Al-Karaki | Routing techniques in wireless sensor networks: a survey[END_REF]. The main role of WSNs is to collect and route data back to a collection point known as the Base Station (BS). Due to large network sizes, limited power supplies and inaccessible remote deployment environments, WSNs introduce new challenges compared to traditional wireless networks [START_REF] Estrin | Connecting the physical world with pervasive networks[END_REF]. One of the major concerns is the short lifespan of such networks. To extend the network lifetime, designing energy-efficient protocols is critical. Exploiting clustering is one of the most interesting ways in achieving such an objective [START_REF] Younis | Node clustering in wireless sensor networks: Recent developments and deployment challenges[END_REF]. A number of clusterbased protocols have been developed to make these networks practical and efficient [START_REF] Abbasi | A survey on clustering algorithms for wireless sensor networks[END_REF].

Cluster-based protocols were originally proposed concerning wireline networks in order to alleviate scalability issues. However, clustering is used in WSNs for efficient energy consumption in order to extend network lifetime. Typically, nodes with relatively higher residual energy can operate as Cluster Heads (CHs). They send aggregated information from the nodes with lower residual energy, which are only used to perform sensing in the proximity of the target. This helps to balance energy consumption inside the network, and thus improves the overall network lifetime. To achieve energy efficiency, a multi-hop routing protocol can be integrated into the cluster-based protocols in two ways. The first way concerns integrating an intra-communication protocol within a cluster [START_REF] Mhatre | Design guidelines for wireless sensor networks: communication, clustering and aggregation[END_REF]. Instead of sending data directly to the CHs, a routing protocol is incorporated within a cluster to achieve such an objective. The second way concerns integrating an inter-communication protocol between clusters [START_REF] Neander | Asymmetric multihop communication in large sensor networks[END_REF]. In this case, CHs organize themselves into a routing architecture to deliver data to the BS, instead of direct data transmission to the BS.

To this regard, various clustering protocols have been proposed in the literature [START_REF] Younis | Node clustering in wireless sensor networks: Recent developments and deployment challenges[END_REF]- [START_REF] Abbasi | A survey on clustering algorithms for wireless sensor networks[END_REF]. In order to reduce the signalling overhead, the proposed protocols usually use heuristic probabilistic techniques, in which a node becomes a CH with a certain probability, as opposed to deterministic techniques, which require sending overhead messages for the CH's election.

In this paper, we chose the HEED protocol [START_REF] Younis | HEED: A hybrid, energy-efficient, distributed clustering approach for ad-hoc sensor networks[END_REF], one of the well-known probabilistic clustering protocol in WSNs, to extend the network lifetime. In HEED, the probability of becoming a CH is based on the nodes' current residual energy. However, in HEED no emphasis is made on the specific network topology characteristics, such as the degree of connectivity in the selection of CHs. The main idea behind our proposed CDEEC protocol is then to consider, in addition to the information about the nodes' residual energy, the degree of connectivity in CH selection. The main objective of our proposal is to allow nodes with a high degree of connectivity to become CHs and hence elect a lesser number of, and more efficient, CHs. CDEEC exploits the aggregation operation at CH level and incorporates a routing protocol inside clusters to conserve further energy, thereby extending the network lifetime.

The remainder of this paper is organized as follows. Section II presents the related work. In Section III, we expose the network model. Section IV describes the proposed clustering protocol, followed by the performance evaluation in section V. Finally, section VI contains our concluding remarks.

II. RELATED WORK

Recently there has been increased interest in studying energy-efficient clustering algorithms in the context of both ad hoc and sensor networks. The main aim of clustering protocols in ad hoc networks is to generate the minimum number of clusters while maintaining network connectivity. In these algorithms the election of CHs is based mainly on the identity of nodes [START_REF] Baker | The architectural organization of a mobile radio network via a distributed algorithm[END_REF], the degree of connectivity [START_REF] Gerla | Multicluster, mobile, multimedia radio network[END_REF] or the connected dominating set [START_REF] Wan | Distributed construction of connected dominating set in wireless ad hoc networks[END_REF]. These techniques are discussed in depth in [START_REF] Yu | A survey of Clustering schemes for Mobile Ad hoc Networks[END_REF].

In the case of WSNs, the main objective of clustering protocols is to minimize energy consumption by the network in order to extend the network lifetime (surveys dealing with WSN clustering protocols can be found in [START_REF] Younis | Node clustering in wireless sensor networks: Recent developments and deployment challenges[END_REF] and [START_REF] Abbasi | A survey on clustering algorithms for wireless sensor networks[END_REF]).

We classify WSN clustering protocols into two categories: probabilistic, in which a node becomes a CH with a certain probability, and deterministic, which requires an exchange of overhead messages for the CH's election. Here, we review some of the protocols proposed in the literature for each class: PEGASIS [START_REF] Lindsey | PEGASIS: Power-efficient gathering in sensor information systems[END_REF], DWEHC [START_REF] Ding | Distributed energy-efficient hierarchical clustering for wireless sensor networks[END_REF], and TASC [START_REF] Virrankoski | TASC: topology adaptive spatial clustering for sensor networks[END_REF] in the deterministic class, and EEHC [START_REF] Bandyopadhyay | An energy efficient hierarchical clustering algorithm for wireless sensor networks[END_REF], EECS [START_REF] Ye | EECS: an energy efficient clustering scheme in wireless sensor networks[END_REF], and HEED [START_REF] Younis | HEED: A hybrid, energy-efficient, distributed clustering approach for ad-hoc sensor networks[END_REF] in the probabilistic class.

Power-Efficient GAthering in Sensor Information Systems (PEGASIS) [START_REF] Lindsey | PEGASIS: Power-efficient gathering in sensor information systems[END_REF] is composed of two phases: 1) steady and 2) gathering. The steady phase consists of a formation of chains instead of clusters. In the chain formation, the BS and sensor nodes are connected via a chain using a greedy algorithm. One of the nodes, in the chain, is selected in turns to represent the head. In the data-gathering phase, each node delivers the sensing data to the nearest neighboring node until the data reaches the head node which aggregates and delivers the sensing data to the BS.

Distributed Weight-Based Energy-Efficient Hierarchical Clustering (DWEHC) [START_REF] Ding | Distributed energy-efficient hierarchical clustering for wireless sensor networks[END_REF] is another deterministic clustering protocol, which aims at high energy efficiency by generating balanced cluster sizes and optimizing the intra-cluster topology. To become a CH, in DWEHC, each sensor node calculates its weight after exchanging its geographic location with its neighbors. The weight is a function of the sensor's residual energy and the distance to its neighbors. In a neighborhood, the node with the highest weight would be elected as a CH and the remaining nodes become CMs. The CMs at this stage are considered as first level because they are attached directly to the CH. Then, the CMs progressively adjust their level to reach a CH using the minimum amount of energy. Specifically, if a CM finds an intermediary node with the minimum energy cost to reach its CH, then it can choose this intermediary node as its parent.

Similar to the idea in [START_REF] Ding | Distributed energy-efficient hierarchical clustering for wireless sensor networks[END_REF], Topology Adaptive Spatial Clustering (TASC) [START_REF] Virrankoski | TASC: topology adaptive spatial clustering for sensor networks[END_REF] decomposes large non-uniform networks into smaller locally-uniform clusters. TASC assumes nodes are aware of 2-hop neighborhood information and also that they know the distance to their neighborhood. This is achieved by deriving a set of weights that includes distance, connectivity and density information within the locality of each node. The derived weights form the terrain for holding a CH election procedure in which each node selects the node closer to the center of mass of its 2-hop neighborhood to becoming its CH.

Energy-Efficient Hierarchical Clustering (EEHC) [START_REF] Bandyopadhyay | An energy efficient hierarchical clustering algorithm for wireless sensor networks[END_REF] is a probabilistic clustering algorithm. The basic operation of the EEHC algorithm consists of electing CHs with probability p and each CH announces its election to the k-hop away neighboring nodes. Any node that receives such a CH election announcement, if it is not itself a CH, becomes a member of the closest cluster. In addition, if the election announcement does not reach a node within a specific time interval, the node becomes a forced CH. However EEHC is extended to a corresponding multi-level architecture. In fact, the basic clustering process is recursively triggered at different steps. At each step i, EEHC elects levels of CHs with a corresponding election probability p i , and a forwarding parameter k i . The communication flow in EEHC is as follows. Ordinary sensor nodes, nodes that are not elected as CHs, transmit their collected data to the corresponding first-level (level 1) CHs, the CHs of the first-level clusters transmit the aggregated data to the second-level CHs and so on, till the highest level, h, of the clustering hierarchy is reached. Then the CHs of those h level clusters transmit their final aggregated data reports to the BS.

In [START_REF] Ye | EECS: an energy efficient clustering scheme in wireless sensor networks[END_REF], the authors proposed the Energy-Efficient Clustering Scheme (EECS) protocol. In this protocol, CH candidates compete for the ability to elevate to a CH with a certain probability. This competition involves candidates broadcasting their residual energy to neighboring candidates. If a given node does not find a node with more residual energy, it becomes a CH. EECS extends this algorithm by the dynamic sizing of clusters based on cluster distance from the BS. The result is an algorithm that addresses the problem that clusters at a greater range from the BS require more energy for transmission than those that are closer. Ultimately, this improves energy distribution throughout the network, resulting in better resource usage and extended network lifetime.

Hybrid, Energy-Efficient Distributed Clustering (HEED) [START_REF] Younis | HEED: A hybrid, energy-efficient, distributed clustering approach for ad-hoc sensor networks[END_REF] is one of the well-known clustering protocols in WSNs, and on which CDEEC is built. The HEED protocol operates in two main phases: (1) the set-up phase where clusters are formed and (2) the steady phase where the sensor nodes transmit their data using Time Division Multiple Access (TDMA) frames. The choice of the CHs with HEED is made in an iterative way. The aim is to achieve a better CH distribution within the WSN, but this is made at the cost of more complexity and increased CH overhead within the WSN.

The first sub-phase is the initialization. Nodes exchange hello messages to discover their neighborhoods. The second sub-phase consists of a competition process and allows the election of CHs in the network. The third sub-phase is the finalization and allows nodes to join their corresponding CH based on the degree of connectivity.

Before a node starts executing HEED, it calculates its probability of becoming a CH, CH prob , as follows:

CH prob = C prob × E residual E max (1)
where C prob is an initial percentage of CHs among all nodes, E residual is the current node residual energy and E max is a reference maximum energy.

The clustering process requires a number of iterations at each node. During iteration i, every uncovered node that does not hear any CH notification from its neighbors, elects itself as a CH with a probability CH prob . After iteration i, the set of potential CHs, S CH , is set to: {CHs after iteration i-1 new CH in iteration i}. A node n j joins the CHs with the lowest cost in S CH . S CH may include n j itself if it is selected as a tentative CH. Then, every node doubles its CH prob and goes to the next iteration.

If a node elects itself as a CH, it sends a tentative message including node ID, status and cost. If a node completes HEED execution without being assigned to a final CH, it considers itself uncovered and announces itself as a CH with final state.

However, the election of CHs in HEED does not take into a consideration the degree of connectivity, which can increase the efficiency of the CH election. We propose a mechanism that allows nodes with a high degree of connectivity to become CHs and, hence benefit from the aggregation operation at CH level to conserve more energy. In addition, as HEED does not implement a routing scheme inside the clusters, we extend the HEED's intra-cluster communication by incorporating a routing protocol within a cluster.

III. NETWORK MODEL

We consider a WSN consisting of N sensors deployed over a vast field to continually monitor the environment. We denote the i-th sensor node by n i . We use the following assumptions concerning sensor nodes and the underlying network model:

• Sensor nodes and the BS are stationary after deployment.

The BS is far away from the square sensing field and can be reached by sensor nodes under a single high transmission range T R2.

• Nodes are dispersed in the A × A area. Nodes are homogenous and have the same capabilities. A unique identifier ID is assigned to each node.

• Each node n i can reach its neighbors n j (n j ∈ Set neighbor i , where Set neighbor i denotes the set of n i ' neighbors) and the BS with transmission ranges T R1 and T R2, respectively (T R1< T R2). We refer here to T R1 and T R2 as the minimum and the maximum energy level available at the sensor node, respectively, in which nodes can choose to reach the required destination.

• Links are symmetric, i.e., if n i ∈ Set neighbor j , then n j ∈ Set neighbor i .

• Aggregation is performed only at CH level. CHs use perfect aggregation to eliminate data redundancy and reduce communication load i.e., any number of packets can be aggregated to one packet [START_REF] Heinzelman | An Application-Specific Protocol Architecture for Wireless Microsensor Networks[END_REF]. Other aggregation techniques, such as those proposed in our previous work [START_REF] Abdul Azim | SAG: Smart Aggregation Technique for Continuous-Monitoring in Wireless Sensor Networks[END_REF], are still possible. In order to estimate the mote's residual energy, we use the energy model from Heinzelman et al [START_REF] Heinzelman | An Application-Specific Protocol Architecture for Wireless Microsensor Networks[END_REF]. We assume that only transmission, reception and aggregation consume energy. No energy is consumed during sleep mode. According to this model, a sensor spends E elec = 50nJ to run the transmitter receiver circuitry. To transmit a packet of size k over a distance d, the amount of energy consumed E T x is given by the following equation.

E T x = (E elec * k) + E f s * k * d 2 , if d < d 0 (E elec * k) + E mp * * d 4 , otherwise. (2)
where E mp and E f s depends on the transmitter amplifier model to use, d 0 is the distance threshold between the transmitter and the receiver over which the multi-path fading channel model is used. To receive a message of size k, the energy required by the receiver is given by: E Rx = (E elec * k). The energy consumed in aggregating the data is E DA .

IV. PROPOSED PROTOCOL

This section describes our proposed protocol. Firstly, we study the case of a simple topology using HEED, and then we provide some hints on our proposed CDEEC protocol used to alleviate issues that may arise within HEED. An additional mechanism is also provided in order to enhance the performance of our proposal in general topologies.

A. A case study

To illustrate the operation of HEED we use the simple topology presented in Fig. 1. Here, we iterate sub-phase two (i.e., CH selection sub-phase) of the HEED set-up phase. Tables I and II draw the different iterations for nodes n 1 and n 4 , respectively. We note that nodes n 1 , n 2 , n 5 and n 6 behave similarly since they are symmetric. Likewise for nodes n 3 and n 4 .

In our analysis we use the following notations. CH previous denotes the probability of becoming a CH at the previous iteration i.e., iteration i -1. We set CH previous to 0.2 for all the nodes at iteration 0.

In our example (see tables I and II), we assume that only nodes n 1 and n 5 access the medium and send tentative messages for their neighbors. Following each iteration the set S CH (n j ) of CHs, that are one hop away from node n j (j = 1, ..., N ), is updated. After each iteration i, node n j doubles its CH prob . Once CH previous of a node n j reaches 1, it joins the least costly CH of the set S CH (n j ). Otherwise (i.e., S CH (n j ) = ∅), node n j forces itself to become a CH.

In our example, at the end of sub-phase two, nodes n 1 and n 5 become CHs. At the finalizing sub-phase (sub-phase three), nodes n 3 and n 4 join their CHs nodes n 1 and n 5 , respectively. On the other hand, nodes n 2 and n 6 force themselves to become CHs since S CH (2) = ∅ and S CH (6) = ∅.

It is worth noting that with HEED the cost of a CH is considered to be its degree of connectivity (i.e., the count of neighbors it has). Hence, nodes with a lower degree of connectivity are preferred to become CHs.

It is also important to highlight that with HEED, once a node n j sends a tentative message at iteration i to possibly become a CH, it prevents its neighbors henceforth from sending tentative messages (i.e., from competing to become CHs). As a result the first node among its neighborhood that succeeds in sending a tentative message becomes the CH for its neighbors. If many nodes send tentative messages at the same iteration, the least costly nodes become CHs.

Intuitively, with HEED, nodes with a high degrees of connectivity have a low probability of becoming CHs. This is due to two main reasons:

• Firstly, there is a high probability that one of its neighbors sends a message before the highly-connected node. This probability increases, the higher the nodes' degree of connectivity.

• Secondly, if a highly-connected node succeeds in sending a tentative message at the same iteration as some of its neighbors, the least-costly node (i.e., the node with the lowest degree of connectivity) will be elected as a CH. To illustrate this again, let us consider the following star network example presented in Fig. 2. Suppose that all the N nodes have the same CH prob at the beginning denoted by p. It is easy to see that once one of the leaf nodes sends a tentative message, it prevents the node at the center n 0 from becoming a CH. Specifically, the probability that node n 0 is elected as a CH among its neighbors can be written as follows:

P rob{A|B} = p(1 -p) N -1 1 -(1 -p) N (3)
where the events A and B are:

A = {n 0 = a CH}. B =
{there is at least one node that has sent a tentative message}, and thus the probability that leaf nodes become CHs is:

P rob{ Ā|B} = 1 - p(1 -p) N -1 1 -(1 -p) N (4)
We can see that P rob{A|B} ≪ P rob{ Ā|B}. In this case, all the leaf nodes become CHs, and node n 0 will join one of them. Then, each of the leaf nodes has to send its report directly to the distant sink node, which results in excessive energy consumption. However if the node in the center, n 0 , was elected as a CH, it would aggregate all the data from the leaf node and send the aggregate packet at once to the sink node. As a result, high energy communication with the sink node is replaced by low energy and short distance communication with the CH.

This example illustrates the limitations of the HEED, which may lead to excessive energy consumption. In what follows we propose our CDEEC protocol to cope with the aforementioned issue.

To alleviate the aforementioned problem and to allow the election of efficient CHs that maximally exploit the aggregation operation, we propose to block nodes with a lower degrees of connectivity (i.e., with a degree of connectivity < M , where M is a pre-specified threshold) from competing to become CHs. This proposed technique is inspired by [START_REF] Ye | EECS: An energy Efficient Clustering Scheme in wireless Seneor networks[END_REF], in which nodes are blocked from CH competition with a predefined probability. By doing so, we aim to block nodes that are not suitable to become CHs in advance and to, therefore allow only targeted nodes to CH competition.

Obviously, a blocked node will force itself to become a CH at sub-phase three since its set of CHs S CH = ∅. To illustrate this let us revisit the example of Fig. 1. In this example we block nodes with a degree of connectivity < 2 from participating in the CH selection process i.e., we set M = 2. In doing so nodes n 1 , n 2 , n 5 , n 6 do not participate in the CH selection process. Only nodes n 3 and n 4 compete to become CHs. Suppose that node n 1 is elected as a CH. In this case, nodes n 3 , n 2 and n 4 will join CH n 3 as Cluster Members (CMs). On the other hand, the sets of CHs associated with nodes n 5 and n 6 are null. As such, nodes n 5 and n 6 force themselves to become CHs. We can see that, with CDEEC, three nodes (either {n 3 , n 5 , n 6 } or {n 1 , n 2 , n 4 }) are elected as CHs. On the other hand, using HEED, nodes {n 1 , n 2 , n 5 , n 6 } (in addition to nodes n 3 and n 4 ) will compete from the beginning in the CH selection process. As such, the following four nodes {n 1 , n 2 , n 5 , n 6 } are likely to become CHs. This results in increasing energy consumption in the network since it increases the number of high energy-consuming communications between the sensor nodes and the distant sink node (four CHs instead of three). Using our method, however, replaces one high energy consuming communication with a local and lower-consuming energy communication.

B. Additional mechanism

To further reduce the number of CHs, we also propose to permanently block nodes with degrees of connectivity < M from becoming CHs. To illustrate this, let us revisit Fig. 1. Recall that in this example only node n 3 is selected as a CH at the end of sub-phase two of the clustering process. Moreover nodes n 1 , n 2 and n 4 join CH n 3 as CMs. On the other hand, the sets of CHs associated with nodes n 5 and n 6 are null. Instead of forcing these nodes to become CHs at the end of the clustering process, we propose to let them join CH n 3 by using a multi-hop routing. In other words, nodes n 5 and n 6 join CH n 3 as CMs. Specifically, both nodes send their reports to the CH n 3 through n 4 .

In doing so, we replace the high energy-consuming communications between sensor nodes n 5 and n 6 with the distant sink node by local and short multi-hop communications with the CH. As a result, the CH n 3 is the only one responsible for communicating the aggregated data to the distant sink node. Compared to HEED, we can see that energy conservation can be realized. Based on this simple example, only one CH is elected instead of four by enabling our scheme. In other words, we increase the aggregation capacity of the elected CH as, instead of four packets, only one is sent to the sink node.

Generally, in order to avoid isolating blocked nodes from the networks a blocked node,n i , can decide to unblock itself, following two conditions:

1) If all of n i 's neighbors are not participating in the clustering (due to dead nodes or is the blocking process), and 2) If n i has the highest degree of connectivity among its neighbors (in case of the same degree of connectivity, they all decide to participate in the clustering process). By doing so, we ensure at least one unblocked neighbor for any n i . As a result, the blocked node chooses one of its neighbors as a next-hop node as follows:

• If n i hears final messages from CHs, it registers itself as a CM and joins the node with the highest cost among the CH announcements.

• If n i does not hear any final message, it chooses one of the unblocked neighbors as the next-hop.

V. PERFORMANCE EVALUATION

In this section, we describe the results of simulations carried out. To achieve this, we use TOSSIM under TinyOS [START_REF] Levis | TOSSIM: Accurate and Scalable Simulation of Entire TinyOS Applications[END_REF]. The application we use in our study is a continuous monitoring application, where data is generated continuously with a predefined reporting frequency. Nodes can use two different levels of power depending on the type of transmissions. On one hand, CMs use transmission range T R1 to reach CHs. On the other hand, CHs use transmission range T R2 to reach the BS. The clustering process is triggered every round time in order to balance the energy consumption of CHs. Therefore, for every round time, new CHs and CMs are elected. The network is configured with the parameters listed in Table III.

A. Simulation results

Firstly, we show the results concerning the case study that has already been described in subsection IV-A. Then, we show the results concerning the additional mechanism generalized into random topologies.

1) Case study: Fig. 3(a) and Fig. 3(b) show the probability of becoming a CH for each node for both HEED and CDEEC protocols. As explained before, the probability of nodes n 3 and n 4 in CDEEC is higher than that of nodes in HEED. This is because of blocking nodes n 1 , n 2 , n 5 and n 6 from participating in the cluster process. The probability of node n 3 is higher than that of node n 4 since the priority of becoming CHs is given to nodes with a smaller node ID, in case of those with the same degree of connectivity. Fig. 4 shows the average network residual energy over time. The result clearly demonstrates that CDEEC protocol conserves more energy, when compared to HEED.

The small number of reports received at the BS (see Fig. 5) in CDEEC compared to HEED can be explained by the fact that only a small number of CHs are elected, thereby reducing the amount of data to be transmitted and thus saving energy. This does not mean that less information is reported to the BS. Rather it means that CDEEC benefits from aggregation operation. The small difference in the number of reports issued using the two protocols is explained by the fact that a small number of nodes are used in the network, and thus a small number of CHs are elected. Finally, Fig. 6 depicts the average energy consumed for both protocols. As we can see, CDEEC saves further energy at the expense of high energy consumption of nodes with a high degree of connectivity. Recall that energy balancing is achieved by reclustering the network in each round. Therefore, nodes that are elected at any given round will have little chance of becoming CHs in the next round as CH election depends highly on the residual energy of nodes. 2) Additional mechanism: The results in this section concern our CDEEC protocol while considering the generalized additional mechanism using random topologies.

To validate the performance of CDEEC we used random topologies rather than specific topologies. Three random topologies of 20 nodes were considered. These topologies differ based on the ACD parameter, defined as the Average Connectivity Degree of the network. As the proposed protocol is tightly based on the degree of connectivity, considering random topologies with different ACD parameter is sufficient to validate the results. Therefore we use different values of T R1: 20m, 35m and 55m to create three random topologies of different connectivities: small (i.e., ACD= 2), moderate (i.e., ACD= 6.2) and highly-connected networks (i.e., ACD= 10), respectively.

We set the pre-specified threshold M to ACD. Fig. 7 and Fig. 8 show the effect of M on both protocols. We can see that HEED is suitable for high-connected topologies. The reason is that more CHs are elected with small ACD value in HEED and therefore consume excessive energy. In contrast, CDEEC outperforms HEED since it restricts the percentage of CHs. We can also see from these figures that CDEEC performs better in the topology with ACD = 6.2 than in the topology with ACD = 10. Though in these situations almost the same percentage of CHs are formed, the better performance of the ACD = 6.2 case is achieved due to the fact that nodes use a smaller transmission range compared to the other topology. The choice of the threshold M limits the number of nodes that will participate in the clustering process, which has a direct effect on the multi-hop routing process. Indeed, with a small M value, HEED is preferred as the majority of nodes will participate in the clustering process and tend not to use a multihop routing inside the clusters. However, choosing a high value for M decreases the number of nodes competing in clustering but at the same time makes the multi-hop process difficult to accomplish. This difficulty is related to whether the blocked nodes will be surrounded by at least one unblocked neighbor. Choosing an appropriate M threshold requires knowledge of the network, such as the network structure (to specify which nodes are to be blocked or not), or the ACD parameter, as we have already seen.

VI. CONCLUSION

This paper aims at reducing energy consumption in wireless sensor networks. To achieve this, we considered the node's degree of connectivity, along with their residual energy, in the cluster head selection process. Specifically, we increased the probability of becoming CHs for nodes with a high degree of connectivity. By doing so, we benefit from the aggregation capacity of such nodes, thus enabling additional energy conservation. To evaluate the performance of our proposed CDEEC protocol, we implemented CDEEC in TinyOS and compared it with HEED. Results showed that significant energy conservation can be achieved.

2 Fig. 1 .

 21 Fig. 1. Network topology.

Fig. 2 .

 2 Fig. 2. Star topology.

Fig. 3 .

 3 Fig. 3. Probability of becoming a CH vs node ID.

Fig. 4 .

 4 Fig. 4. Average residual energy of the network.

Fig. 5 .

 5 Fig. 5. Total reports received at the BS Vs time.

Fig. 6 .

 6 Fig. 6. Average energy consumed Vs Node degree.

Fig. 7 .

 7 Fig. 7. Average energy consumed in the network.

Fig. 8 .

 8 Fig. 8. Percentage of CHs in the network.