Diego Arenas
email: diego.arenas@ifsttar.fr

Rémy Chevrier
email: remy.chevrier@ifsttar.fr

Clarisse Dhaenens
email: clarisse.dhaenens@lifl.fr

Said HANAFI c Joaquin Rodriguez
email: joaquin.rodriguez@ifsttar.fr

Application of a Co-evolutionary Genetic Algorithm to solve the Periodic Railway Timetabling Problem ⋆

Keywords: Genetic Algorithms, PESP, Railway, Scheduling, Timetabling

In train operations, a timetable is used to establish the departure and arrival times for the trains at the stations or other relevant locations in the rail network or a subset of this network. The elaboration of a timetable responds to the commercial needs of the customers, for both passenger and freight traffic, but also, it must respect some security and capacity constraints. The combination of these requirements and constraints makes the preparation of a yearly timetable a complex process that usually takes months to be fully completed. This paper addresses the problem of generating periodic timetables, which means that the trains concerned are operated on a recurrent pattern, e.g., trains of the same line will run every 30 minutes, we present a suitable constraint-based model of the problem. Furthermore, we propose a dedicated genetic algorithm, based on a co-evolutionary scheme with two populations, to create feasible and quality periodic timetables in short periods of time. Finally, two case studies are discussed, both of them representing a subset of the Netherlands railway network.

Introduction

The train timetabling problem involves mainly two procedures: First, to define the departure and arrival times for the trains at the stations or other relevant locations in the rail network or a subset of this network. And second, to establish the platforms and route plans, i.e., the assignment of each train to an appropriate platform and corresponding inbound and outbound routes in every station. For a timetable to be viable or feasible, it must respect some constraints that may be hard constraints: related to safety and capacity issues, e.g., two trains cannot be in the same track segment at the same time, or soft constraints: related to specific requirements for the timetable, e.g, to establish particular passenger connections between two lines of trains.

It is desired that the generation of the arrival and departure times and the selection of the routes through the stations are carried out simultaneously. However, real-size problems are too large and have thousands of constraints. Even modern approaches, using the current available technology are not able to produce a result within an acceptable amount of time. Therefore the two parts of the problem are treated separately: first the timetable is generated and only when a feasible timetable is found, the routings through the stations are determined. However, at this point a feasible routing plan for the given timetable may not be possible, then it is necessary to come back to the previous step and generate another viable train timetable. Many iterations may be necessary in order to finally produce both a viable timetable and routing plan.

In this paper, we focus on the generation of a feasible train timetable, therefore the platforming and routing plans are not discussed. Additionally, besides the timetabling problem, there are two other procedures that are closely related to the railway timetable design problem, see Figure 1. First: the rolling stock plan: which specifies among other things, which engines and carriages will be used for each train line. And second, the crew schedule, which defines the working plan for each person involved on train trips, such as drivers, controllers, etc. On each of the steps described above, an optimization process can be carried out, which consists, for example, in making minor modifications to a feasible timetable in order to optimize it, e.g., to reduce the total number of trains, or minimize the overall waiting times for passengers during connections, or reduce the number of platforms needed in some station at some given time.

Traditionally, the timetable planning problems were solved based on the experience and expertise of the involved railway planners. But at the early nineties it was recognized within the railway community that the application of models and techniques for supporting the solution process of these problems may be truly advantageous. It would not only lead to better solutions, but it also conducts to a considerable reduction of the time needed for the planning process.

There are mainly two types of timetables, periodic and non-periodic. The difference between them is that: Non-periodic timetables are designed based only on the set of requests for trips using the network, i.e., try to allocate the departure and arrival times corresponding to the requested departure and arrival times of the trip. While periodic timetables are built based on a time period T , a timetable pattern is created for this period where all trips are scheduled. Then this pattern is repeated to cover the whole day. This means that each line of passenger trains is operated in a regular recurrent way, e.g., the trains of the same line run every 30, 60 or 120 minutes. See Table 1 for an example of a periodic timetable. One advantage of a periodic timetable is that passengers can easily keep in mind the departure time of their train at their station. However, a drawback is that these kinds of timetables are not very flexible, and it is difficult to offer a large number of direct connections. Another disadvantage, is the fact that a completely periodic timetable may be rather inefficient: trains might have to be operated even at some specific times of a day with only a small number of passengers. Therefore, in practice there are usually exceptions to the completely periodic timetable. For example, there may be some additional trains during rush hours, and frequencies may be reduced during late evening hours [START_REF] Kroon | Railway timetabling from an operations research[END_REF]. In this paper we discuss only the generation of periodic timetables.

Train stations may have many platforms and complex track layouts, which could also be used for resolving timetable conflicts and soft constraints. In this work, we assume that each station has enough number of platforms to accommodate the trains that will finally stop on it.

Finally, in this paper, a line plan is assumed to be known a priori, i.e., the list of lines, and the stop schedule for each line, as well as the desired connections between lines.

The rest of this paper follows the next structure: Section 2 reviews some popular approaches that address this problem. Section 3 presents the formulation used in our approach. Section 4 describes how we use genetic algorithms to solve the train periodic timetable generation problem. Section 5 presents the case studies and discusses the obtained results using our technique.Finally, in Section 6, we conclude this paper and present some perspectives for further works.

Brief Literature Review

There are numerous authors that developed different techniques to address scheduling problems and specifically the train timetabling problem, in this section we will briefly discuss some of the most relevant works that deal with this topic.

Serafini et al.

[?] introduce the Periodic Event Scheduling Problem (PESP), which is a framework for scheduling periodic activities. The PESP is proven to be an exceptionally rich model that can be used in many applications, e.g., airline scheduling, traffic light scheduling and of course, train scheduling, among others.

Odijk [START_REF] Odijk | A constraint generation algorithm for the construction of periodic railway timetables[END_REF] proposes the PCG (PESP Cut Generation) algorithm to generate timetables for train operations using a mathematical model consisting of periodic time window constraints as input. His model is based on the definition of constraints applied on a mathematical formulation to solve the problem.

Kroon et al. [START_REF] Kroon | Railway timetabling from an operations research[END_REF] separate the problem in two parts: the generation of the departure/arrival times and the selection of the routes through the stations. The first part describes a mathematical formulation based on the PESP to solve the scheduling problem, and since the second part is out of scope of this paper it will not be described. On the first part, he defines a set of different types of constraints and then use them as input to the PESP solver, obtaining then a timetable that will be used as input to the second part.

Liebchen's [?] work is also based on the PESP, however two natural optimization variants are presented and discussed: The Max-T-Pesp, and a heuristic method called Cut-Heuristic. Then these two techniques are applied to solve an optimization problem on a real scenario: the subway system of Berlin. The cut-heuristic algorithm is based on two optimization-based observations for the PESP.

Caprara et al. [START_REF] Caprara | Modeling and solving the train timetabling problem[END_REF] concentrate on the problem of a single, one-way track linking two major stations, with a number of intermediate stations in between. They use a directed multi-graph in which nodes correspond to departures and arrivals at a certain station at a given time instant. Then they use this formulation to derive an integer linear programming model which is relaxed in a Lagrangian way. The objective is to maximize the sum of the "profits" of the scheduled trains, the profit achieved for each train depends on: the "shift" defined as the absolute difference between the departure times from a given station in the ideal and actual timetables; and on the "stretch" defined as the non-negative difference between the running times in the actual and ideal timetables. Although Caprara et al. model presented in this work is not applied to create a periodic timetable, it is highly flexible and it can be extended to the PESP model and then used to generate periodic timetables.

Up to this point, the authors base their approaches on the resolution of complex linear programming models. Other authors use evolutionary techniques such as genetic algorithms to completely or partially address the problem. As discussed in [?], Genetic Algorithms are proven to be a solid method to deal with the periodic timetabling problem, and their behavior is more stable than other techniques when used both in small and big instances of the problem.

As an instance of partial use of genetic algorithms, Semet's et al. [START_REF] Semet | An efficient memetic, permutation-based evolutionary algorithm for real-world train timetabling[END_REF] approach focus on the reconstruction of a timetable following a small perturbation. They try to minimize the total accumulated delay by adapting the departure and arrival times for each train and allocation of resources. They use a permutation-based evolutionary algorithm to gradually reconstruct the schedule. The railway network can be seen as a graph where nodes are stations or switches and where interconnecting edges eventually hold several tracks for trains to use. Their objective function consists in minimizing the total accumulated delay (for all trains at all nodes). The algorithm is a hybrid, or memetic algorithm, in the sense that it combines an evolutionary engine with a mathematical programming tool, called ILOG CPLEX. The evolutionary part of the algorithm is used to quickly obtain a good but suboptimal solution, which is fed to CPLEX as an initial solution which will search for the global optimum.

An example of complete use of genetic algorithms to address the problem, Kwan and Mistry [START_REF] Kwan | A co-evolutionary algorithm for train timetabling[END_REF] use a coevolutionary algorithm for the automatic generation train timetables, essentially at the early stages of the process, trying to allocate as many capacity requests as possible, and discover the conflicts that have to be solved by the train operating companies. They assign some penalty weights to each type of constraint violation thus the Objective function is to Minimize the weighted sum of violations. They use three types of populations which are evolved one after another: "departing times", "running times" and "resources options". The combination of three individuals i.e one of each population, result on a timetable. On each step of the evolution a timetable is generated and evaluated, once the termination condition is achieved, the algorithm is stopped and the best timetable is given as result.

Problem Formulation

In this section we describe the model for our solution, first we shall define the data, the variables and the objective function. And then, we describe the constraints that our model takes into account.

Data:

T : Period for one cycle. i, j : Trains. s, t : Stations. L, U : Lower and upper time limits for some constraints.

Variables: D i,s : Departure time of train i from station s A i,s : Arrival time of train i to station s Q e,e ′ : Binary value. If the events e and e ′ take place on the same period T the value of Q e,e ′ is 0, or 1 otherwise.

Objective Function: The objective function is defined as the minimization of the weighted sum of constraint violations. The weight for each type of constraint violation depends on the importance of the constraint, e.g., violating a headway constraint is a lot more serious than violating a connection constraint.:

M inimize[(W r × V r ) + (W h × V h ) + (W s × V s ) + (W c × V c )] M inimize (W x × V x ) (1)
Where: W x : Is the weight value corresponding to a violation of a constraint type x. V x : Represents the constraint violation count of type x x: Type of constraint, can acquire one of the following values: r,h,s,c. Corresponding to: Running-time, Headway, Single-track and Connection.

Constraints:

Running Time: Basic constraint, it will ensure that the travel time for some train i departing from station

• The first vector is a set of integers that represents the departure time of each train line at each station, e.g., if the instance of the problem have 3 lines and each line have 2 intermediate stop stations, the size of this vector would be 9: the first 3 integers are the departure times from the first station, the departure time from the first intermediate station and the departure time from the second intermediate station, the same is for the other 2 train lines. • The second vector is a set of binary values. Each value represents whether the lower or higher running/dwelling time is currently selected, for all journeys and stops of each line. A "0" indicates that the lower (fastest) time is selected, and "1" represents the higher (slowest). These two values are defined a priori. For example the vector [011] represents that the train will make the first journey at its nominal speed (fastest), then it will dwell at the first intermediate station for a "longer" period of time, and will make the second journey at its slower allowed speed. (2) The co-evolutionary Algorithm: The main algorithm is described in Algorithm 1. The essential components of the algorithm are briefly described next:

• Population: The algorithm works with two populations: DT : Departure times and RT : Running/Dwelling times. The individuals of each population are represented by the vectors described above. By combining one individual of each population, a timetable can be built. • Variation Operator: The main operator we use at the evolution stage is mutation for both populations, additional operators can be easily added if necessary. • Individual Selection: Used to select an individual from the populations in order to build a timetable, two selection criteria were implemented:

• Discrete: A random individual is picked from the population.

• By Tournament: Three individuals are randomly picked, the one with better fitness value is selected. • Fitness Function: In order to evaluate one individual from one population, another individual from the other population is selected using one of the previously described selection methods. With these two individuals a new timetable is generated, then it is evaluated using the objective function described in Section 3, therefore, the fitness function consists in minimizing the total weight of the timetable.

Additionally we added a small optimization criteria: we count the total number of physical trains that requires the generated timetable, more trains means heavier weight. For experimenting and validating our work, we use two case studies, next a brief description of each one, they where extracted from [START_REF] Kroon | Railway timetabling from an operations research[END_REF] and [START_REF] Meester | Stochastic delay propagation in railway networks and phase-type distributions[END_REF] respectively. Then the results are presented.

The specifications of the computer used for both case studies are: Processor Intel Core i7 @2.9 GHz x4, 3.5 GB of RAM memory and using Ubuntu-64 v12.10 as operational system. Evolve DT population using the variation operator.

5:

for each individual of DT do 6:

Select an individual from RT population using an individual selection technique.

7:

Combine the two individuals to generate a timetable.

8:

Evaluate the timetable using the fitness function. Evolve RT population using the variation operator.

11:

for each individual of RT do 12:

Select an individual from DT population using an individual selection technique.

13:

Combine the two individuals to generate a timetable.

14:

Evaluate the timetable using the fitness function. In Table 2, the details of the line plan and the minimum running times are reported. Note that in practice the minimum running times may depend on the applied rolling stock type. Here "IC" stands for Intercity and "stop" stands for local train or stop-train. The columns in the table should be read as follows: the IC Amf -Dv dwells in Apd with minimum running times of 24 minutes (Amf -Apd) and 10 minutes (Apd -Dv), respectively. It runs once every half an hour, that represents our cycle time. In the remainder of this section, we denote this train with number 1a in the direction of Dv and 1b in the opposite direction. Stations that are not mentioned in Table 2 are not relevant. That is why for the stop-train Amf -Zl, only the stop in Hd is listed in the table. The dwell times for the other stops are included in the minimum running times. The IC Amf -Zl does not stop in between Amf and Zl. However, to construct the timetable one should know the time it passes Hd. The minimum running times from Amf and Zl to Hd are 15 and 20 minutes, respectively.

Case Study 2

In Figure 4 a network with four lines is depicted; it is part of the Netherlands railway network, it is in fact an extension of the previous case study. However, and for sake of simplicity we are not detailing in this paper the input values for this instance, should the reader is interested, refer to [START_REF] Meester | Stochastic delay propagation in railway networks and phase-type distributions[END_REF]. 5.2 Results

Case Study 1

Designing genetic algorithms requires to set some parameters, one of them is the size of the population which represents the total number of individuals or potential solutions. Other important parameter is the stop criteria, essentially, it establishes when the algorithm should stop and return a solution (even if it is sub-optimal, i.e., some constraints are violated).

Experimentation: In order to analyze and compare the impact of the population size in the resulting timetable, we perform 50 executions with a constant number of evaluations using the following permutation of parameters:

Total number of evaluations: (1 000 , 5 000 , 15 000) with a population of (40 , 80 , 160) individuals. The results are displayed in Figure 5, as we can see, our solution can find a feasible train timetable (zero constraints violated) in 0.64 seconds (by this time, at least 11.33% of the generated timetables are feasible). After 0.87 seconds the probability to obtain a feasible timetable is 95.33%. Finally, after 0.89 seconds (or 15 000 evaluations), we can be sure that all timetables generated (100%) are feasible. By comparing the results with different size of populations, we can see that by increasing the size of the population we are actually reducing the number of generations that are evaluated (because of course each generation would need more evaluations) and thus, reducing the performance of the algorithm. We can conclude that for this sized instances of the problem, the number of generations offers greater improvements than the size of the population.

Case Study 2

As well as in the previous example, we analyze and compare the impact of the population size in the resulting timetable, we perform 40 executions with a constant number of evaluations using the following permutation of parameters:

Total number of evaluations: (0.8 M, 3 M, 5 M) with a population of (700 , 800 , 900) individuals. The results are shown in Figure 6, as we can see, our solution is able to generate a feasible timetable, for this instance of the problem, in less than 25 seconds. By this time, 27.5% of the generated timetables do not violate any constraint. Additionally, in average less than one constraint is violated, and only soft constraints (connections) are involved. After 118.5 seconds of execution, 40% of the generated timetables are feasible. In this case study, and unlike the precedent, the population size seems to affect the final solution in a different way. This can be explained because the size of this instance is larger than the precedent, a greater population offers a richer diversity of the potential solutions increasing the effectiveness of the genetic algorithm. Finally, during our experimentations we did not arrive to achieve a 100% of generated feasible timetables within the max evaluation count of 5 M. We believe that this was because the algorithm converged to local optimal solutions, and the value of the mutation operator was not sufficient to solve the problem.

Conclusion and Perspectives

The periodic train timetable design represents a complex but very interesting and exciting research subject.

There are many different ways to address, model and solve the problem and the results of these may improve the current train timetable design system om the countries exploiting periodic timetables.

The results of our implementation were rather convenient for solving the periodic train timetable generation problem. Indeed, even though it represents a hard and complex problem, the performance of the genetic algorithm solving it, is remarkable. Our solution not only can rapidly find a feasible timetable for the problem, but it also offers a flexible platform in which we can easily add, remove or modify the constraints in order to satisfy some other criteria. For example, it can be easily turned in a solution that focuses on the reduction of the overall waiting time for passengers. Also, by modifying the weight values of the constraints violations we can focus in the improvement of a particular type of constraints, for example by optimizing the connections between lines in the rail network.

However, our current solution does not deal with a very important issue in train timetable design, which is the planning of resources such as the track layout inside the stations, the number of platforms or the available rolling stock. Thus, further improvements to our approach should extend our model allowing it to handle this kind of resources in order to offer a more realistic solution.

Furthermore, in the experimentation process of our solution we have only tested the genetic algorithm with one kind of operator to evolve both populations: mutation; it would be necessary to test it with other kinds of operators in order to determine which one has better results with our data representation. This also concerns the individual selection process: tournament, as there are other techniques that may lead to a faster converging process of the algorithm. Fortunately, our flexible implementation allows a straightforward incorporation of those.

Additionally, we could improve our implementation according to a multi-objective paradigm, that is, deeply studying and analyzing how some constraints affect the final timetable. We could produce different kinds of options that would specialize in some particular interest, like optimization or finding a better connection plan to ensure the desired connectivity on the network. We believe that it can be done by simply adding new kinds of constraints or modifying the existing ones.

Finally, we can also combine our genetic algorithm with other optimization technique, such as mathematical programming, to create a hybrid solution. This can be made, not only to overcome the limitations of the genetic algorithms, e.g., a tendency to converge towards local optimal solutions, but also to take advantage of the strengths of both methods.

Fig. 1 .

 1 Fig. 1. Timetable Planning Problem

Fig. 2 .

 2 Fig. 2. Our approach design schema

1

 1 Description of the Instances 5.1.1 Case Study 1 Considers the triangle Amersfoort (Amf) -Deventer (Dv) -Zwolle (Zl) in the Netherlands, as depicted in Figure 3. The other relevant stations are Apeldoorn (Apd), Harderwijk (Hd) and Olst (Ost).

Fig. 3 .

 3 Fig.3. Case Study 1: The triangle Amf -Dv -Zl[START_REF] Kroon | Railway timetabling from an operations research[END_REF]

Fig. 4 .

 4 Fig. 4. Case Study 2: Subset of the Netherlands railway network[START_REF] Meester | Stochastic delay propagation in railway networks and phase-type distributions[END_REF]

Fig. 5 .

 5 Fig. 5. Resuts of experimentation on Case Study 1

Fig. 6 .

 6 Fig. 6. Resuts of experimentation on Case Study 2

Table 1

 1 Periodic timetable for two TGV (High-Speed Train) lines connecting the cities of Paris (Paris-Nord train station) and Lille (Lille-Flandres and Lille-Europe train stations) [?]

	Line	Direction Departure. Arrival Direction Departure Arrival
	Paris-Nord -> Lille-Flandres	Lille	8:46	9:48	Paris	10:11	11:14
	Paris-Nord -> Lille Europe	Lille	9:46	10:45	Paris	11:13	12:15
	Paris-Nord -> Lille-Flandres	Lille	10:46	11:48	Paris	11:41	12:44
	Paris-Nord -> Lille Europe	Lille	11:46	12:45	Paris	13:13	14:14
	Paris-Nord -> Lille-Flandres	Lille	12:46	13:48	Paris	15:11	16:14

  Algorithm 1 Solving the Train Timetable Generation Problem 1: Generate a random DT population ⊲ DT stands for Departure Time 2: Generate a random RT population ⊲ RT stands for Running and/or Dwelling Time 3: while Termination conditions not met do

	4:

Table 2 Case

 2 Study 1 Line Plan[START_REF] Kroon | Railway timetabling from an operations research[END_REF]

	Type	Route	Relevant stops Min. runnning time Freq.	Train
	IC	Amf -Dv	Apd	24 + 10	1	1a,1b
	IC	Amf -Zl	-	35	2	2a,2b,3a,3b
	IC	Zl -Dv	Ost	14 + 6	1	4a,4b
	stop	Amf -Zl	Hd	28 + 26	1	5a,5b

s with destination station t, is between the specified limits Le and U e previously calculated. (Lower Elapsed Time and Upper Elapsed Time, respectively). Le i,s,t ≤ A i,t -D i,s + T × Q Ai,t,Di,s ≤ U e i,s,t

Where: Le i,s,t and U e i,s,t are determined based on the details of the infrastructure between the stations s and t (including the safety system) and the running time characteristics of the involved rolling stock.

Dwell Time: This constraint will guarantee that the stopping time of train i at station s is between the Lr and U r. (Lower Rest Time and Upper Rest Time, respectively).

Where: Lr i,s and U r i,s are predefined values guaranteeing that there is enough time for passenger to board and disembark the train i. Also, if the train i have a connection at station s, this dwelling time should consider it too.

Connections: This kind of constraint will permit the passengers connections between trains, the time between the arrival of train i and the departure of train j must be within the limits of Lc and U c (Lower Connection Time and Upper connection time, respectively).

Where: Lc i,j,s and U c i,j,s are defined as the estimated time needed for a passenger from train i to disembark, find the connection platform and board the train j within station s; larger stations normally require longer connections times.

Headway: Two trains i,j departing (arriving) from (to) station s, must have a headway time between them limited by Lh and U h (Lower Headway Time and Upper Headway Time, respectively)

Where: Lh i,j,s and U h i,j,s depends on the kind of trains: Inter-City (IC) or local (stop-train), and their destination, for example: if two trains i,j departing from station s, having the same final destination, but train i stops at every station and train j only at major stations, train i cannot be overtaken by train j on the track, then the lower bound Lh i,j,s is the minimum difference in running time from station s to the next available station t (Le i,s,t -U e j,s,t ) plus the minimum headway time (predefined for train i). The upper bound U h i,j,s is the difference between the cycle T and the headway time for train j.

Headway -Single Track Opposite Trains: Similar to the Headway constraint, this one is used when there is only one track connecting two stations and there is a train i which goes and comes back, in that case the returning train is denoted by i ′ . The time between the departure of train i and the arrival of train i ′ from/to station s must be within limits Lh and U h.

Where: Lh i,i ′ ,s is two times the minimum running time from station s to the next available station t (2 × Le i,s,t ) plus the minimum headway for the train i, and the value of U h i,i ′ ,s is the difference between cycle time T and the minimum headway for the train i.

A co-evolutionary-based solving method

Our proposal to address this problem is to develop a co-evolutionary genetic algorithm inspired by Kwan et al.

[3] and Semet et al. [START_REF] Semet | An efficient memetic, permutation-based evolutionary algorithm for real-world train timetabling[END_REF] approaches. A comprehensive description about genetic algorithms can be found at [?] and [?]. Our model takes into account the objective function (1) and constraints (2) to (6) described in Section 3 and it supports the generation of periodic train timetables. There are two main parts which are detailed below, an illustrative schema is shown in Figure 2.:

(1) Data representation: In order to represent a timetable we are using two vectors:


Jonathan El

Laurent Gardes

Stéphane Girard

Nonparametric estimation of extreme risk measures from conditional heavy-tailed distributions

Keywords: Conditional tail Expectation, Heavy-tailed distributions, Kernel estimator, Asymptotic normality, Risk measures, Extreme-value statistics. AMS 2000 subject classification: 62G32, 62G30, 62E20

In this paper, we introduce a new risk measure, the so-called Conditional Tail Moment. It is defined as the moment of order a ≥ 0 of the loss distribution above the upper α-quantile where α ∈ (0, 1). Estimating the Conditional Tail Moment permits to estimate all risk measures based on conditional moments such as Conditional Tail Expectation, Conditional Value-at-Risk or Conditional Tail Variance. Here, we focus on the estimation of these risk measures in case of extreme losses (where α → 0 is no longer fixed). It is moreover assumed that the loss distribution is heavy-tailed and depends on a covariate. The estimation method thus combines nonparametric kernel methods with extreme-value statistics. The asymptotic distribution of the estimators is established and their finite sample behavior is illustrated both on simulated data and on a real data set of daily rainfalls.

Introduction

One of the most popular risk measures is the Value-at-Risk (VaR) introduced in the 1990's, see [START_REF] Jorion | Value at risk: the new benchmark for managing financial risk[END_REF] for a review. In statistical terms, the VaR at level α ∈ (0, 1) corresponds to the upper α-quantile of the loss distribution. The Value-at-Risk however suffers from several weaknesses. First, it provides us only with a pointwise information: VaR(α) does not take into consideration what the loss will be beyond this quantile. Second, random loss variables with light-tailed distributions or heavy-tailed distributions may have the same Value-at-Risk [START_REF] Tasche | Expected shortfall and beyond[END_REF]. Finally, Value-at-Risk is not a coherent risk measure [START_REF] Artzner | Coherent measures of risk[END_REF] since it is not subadditive in general 1 .

A coherent alternative risk measure is the Conditional Tail Expectation (CTE) [START_REF] Artzner | Coherent measures of risk[END_REF], also known as Tail-Value-at-Risk, Tail Conditional Expectation or Expected Shortfall in case of a continuous 2 The Regression Conditional Tail Moment: definition and estimation

A new risk measure

Let Y ∈ R be a random loss variable. For α ∈ (0, 1), the Value at Risk of level α is the quantity VaR(α) satisfying P(Y > VaR(α)) = α. The Value at Risk is the most popular risk measure [START_REF] Jorion | Value at risk: the new benchmark for managing financial risk[END_REF] but many others can be found in the literature:

-The Conditional Tail Expectation [START_REF] Artzner | Coherent measures of risk[END_REF] is defined by CTE(α) := E(Y |Y > VaR(α)).

-The Conditional Tail Variance CTV(α) := E (Y -CTE(α)) 2 |Y > VaR(α) was introduced in [START_REF] Valdez | Tail conditional variance for elliptically contoured distributions[END_REF]. It measures the conditional variability of Y given Y > VaR(α) and indicates how far away the events deviate from CTE(α).

-The Conditional Tail Skewness CTS(α) := E Y 3 |Y > q(α) (CTV(α)) 3/2 was defined in [START_REF] Hong | Conditional tail variance and conditional tail skewness[END_REF].

-The Conditional-Value-at-Risk is defined by CVaR λ (α) := λVaR(α) + (1λ)CTE(α) with 0 ≤ λ ≤ 1. It is clear that CVaR 1 (α) = VaR(α) and CVaR 0 (α) = CTE(α). This risk measure is able to quantify dangers beyond VaR(α) and is moreover coherent for λ = 1. Other fundamental properties can be found in [START_REF] Rockafellar | Optimization of conditional value-at-risk[END_REF].

-The Stop-loss Premium reinsurance risk measure with retention level equal to VaR(α) [START_REF] Cai | Optimal retention for a stop-loss reinsurance under the VaR and CTE risk measures[END_REF] is proportional to the difference between CTE(α) and VaR(α): SP(α) := E((Y -VaR(α)) + ) = α (CTE(α) -VaR(α)) , where z + = max(0, z). This measure thus permits to emphasize the dangerous cases.

The first purpose of this paper is to unify the definitions of the above risk measures. To this end, a new risk measure is introduced, the Conditional Tail Moment CTM a (α) := E(Y a |Y > VaR(α)), where a ≥ 0 is such that the moment of order a of Y exists. It is easy to check that all the above risk measures of level α can be rewritten as Φ(VaR(α), CTM 1 (α), CTM 2 (α), CTM 3 (α)), where the function Φ : R 4 → R is taken in Table 1.

Risk measure Φ(t 0 , t 1 , t 2 , t 3 ) CTE(α)

t 1 CTV(α) t 2 -t 2 1 CTS(α) t 3 /(t 2 -t 2 1 ) 3/2 CVaR(α) λt 0 + (1 -λ)t 1 , λ ∈ [0, 1] SP(α) α(t 1 -t 0 )
Table 1: Links between the new risk measure and classical risk measures More generally, the CTM can be used to define any risk measure based on conditional moments of the loss variable above the VaR of level α. For instance, one could introduce the conditional tail kurtosis thanks via the function Φ(t 0 , t 1 , t 2 , t 3 , t 4 ) = t 4 /(t 2t 2 1 ) 2 .

Extreme losses and regression case

As announced in the introduction, our second purpose is to adapt the classical risk measures to extreme losses and to the case where a covariate X ∈ R p is recorded simultaneously with the loss variable Y . To this end, the fixed level α ∈ (0, 1) is replaced by a sequence (α n ) ∈ (0, 1), such that α n → 0. Furthermore, denoting by F (.|x) the conditional survival distribution function of Y given X = x, we define the Regression Value-at Risk by RVaR(α n |x) := F ← (α n |x) = inf{t, F (t|x) ≤ α n }, and the Regression Conditional Tail Moment of order a by:

RCTM a (α n |x) := E(Y a |Y > RVaR(α n |x), X = x),
where a > 0 is such that the moment of order a of Y exists. Note that in this framework, RVaR(α n |x) is the extreme conditional quantile of level α n ∈ (0, 1), see for instance [START_REF] Beirlant | Nonparametric estimation of extreme conditional quantiles[END_REF][START_REF] Daouia | Kernel estimators of extreme level curves[END_REF][START_REF] Daouia | On kernel smoothing for extremal quantile regression[END_REF][START_REF] Gardes | Conditional extremes from heavy-tailed distributions: an application to the estimation of extreme rainfall return levels[END_REF][START_REF] Gardes | Functional kernel estimators of large conditional quantiles[END_REF][START_REF] Wang | Estimation of high conditional quantiles for heavy-tailed distributions[END_REF]. It is then quite easy to adapt the classical risk measures to extreme losses and to the presence of a covariate by applying the desired function (see Table 1) to the vector

(RVaR(α n |x), RCTM 1 (α n |x), RCTM 2 (α n |x), RCTM 3 (α n |x)).
This yields the following risk measures: RCTE(α n |x), RCTV(α n |x), RCTS(α n |x) , RCVaR(α n |x) and RSP(α n |x).

Inference

Let (X i , Y i ), i = 1, . . . , n, be independent copies of the random pair (X, Y ). To estimate the RCTM, we start from the following straightforward equality

RCTM a (α n |x) = 1 α n ϕ a (ϕ ← 0 (α n |x)|x), (1)
where for y > 0, ϕ a (y|x) = E (Y a I{Y > y}|X = x) is the conditional moment of order a ≥ 0. Estimation of the RCTM thus relies on the estimation of the conditional moment. We propose to use a classical kernel estimator (see [START_REF] Parzen | On estimation of a probability density function and mode[END_REF][START_REF] Rosenblatt | Remarks on some nonparametric estimates of a density function[END_REF]) given by φa,n (y|x) =

n i=1 K h (x -X i )Y a i I{Y i > y} n i=1 K h (x -X i ), (2)
where I{.} is the indicator function and h = h n is a non-random sequence such that h → 0 as n → ∞. We have also introduced K h (t) = K(t/h)/h p where K is a density on R p . In this context, h is called the window-width. Since φa,n (.|x) is a non increasing function, we can define an estimator of ϕ ← a (α|x) for α ∈ (0, 1) by φ← a,n (α|x) = inf{t, φa,n (t|x) < α}.

Remarking that ϕ 0 (y|x) = F (y|x), the RVaR of level α n is thus estimated by RVaR n (α n |x) = φ← 0,n (α n |x).

We thus recover the extreme conditional quantile estimator studied in [START_REF] Daouia | Kernel estimators of extreme level curves[END_REF][START_REF] Daouia | On kernel smoothing for extremal quantile regression[END_REF]. The RCTM of order a is estimated by

RCTM a,n (α n |x) = 1 α n φa,n ( φ← 0,n (α n |x)|x). (4)
An estimator of each of the above mentioned risk measures is thus given by

Φ( RVaR n (α n |x), RCTM 1,n (α n |x), RCTM 2,n (α n |x), RCTM 3,n (α n |x)), (5)
where the function Φ is chosen in

(α n |x) = φ1,n ( φ← 0,n (α n |x)|x)/α n .
The joint asymptotic distribution of the RCTM and RVaR estimators, and consequently of all the above mentioned estimators, is established in the next section.

Main results

Our main assumption is the following: (F.1) We assume that the conditional survival distribution function of Y given X = x is heavytailed and admits a probability density function.

To summarize, (F.1) amounts to assuming that the conditional distribution of Y given X = x is in the Fréchet maximum domain of attraction. Assumption (F.1) is also equivalent to stating that for all y > 0, F (y|x) = P(Y > y|X = x) is regularly varying at infinity (see [START_REF] Bingham | Regular Variation[END_REF]) with index -1/γ(x) denoted by F (.|x) ∈ RV -1/γ(x) i.e for all λ > 0,

lim y→∞ F (λy|x) F (y|x) = λ -1/γ(x) .
In this context, γ(.) is a positive function of the covariate x and is referred to as the conditional tail index since it tunes the tail heaviness of the conditional distribution of Y given X = x. It also appears that, under (F.1), a sufficient condition for the existence of RCTM a (1/.|x) is a < 1/γ(x). As established in Lemma 1, condition (F.1) also implies that, for all a

∈ [0, 1/γ(x)), ϕ a (.|x) ∈ RV a-1/γ(x) . Since, moreover, F (.|x) ∈ RV -1/γ(x) , we have RCTM a (1/.|x) ∈ RV aγ(x)
. This is equivalent to state that for a ∈ [0, 1/γ(x)) and for all y > 0,

RCTM a (1/y|x) = y aγ(x) ℓ a (y|x), (6)
with a positive index aγ(x) and, for x fixed, ℓ a (.|x) is a slowly-varying function at infinity, i.e for all λ > 0,

lim y→∞ ℓ a (λy|x) ℓ a (y|x) = 1. (7)
To establish the asymptotic normality of (4), the following additional conditions are required. First, as remarked in [START_REF] Bingham | Regular Variation[END_REF], p.15, since slowly-varying functions are of interest only asymptotically, one can assume without loosing generality that in (6)

(F.2) ℓ a (.|x) is normalized for all a ∈ [0, 1/γ(x)).
In such a case, the Karamata representation (see [START_REF] Bingham | Regular Variation[END_REF], Theorem 1.3.1) of the slowly-varying function can be written as

ℓ a (y|x) = c a (x) exp y 1 ε a (u|x) u du , (8)
where c a (.) is a positive function and ε a (y|x) → 0 as y → ∞. Thus, ℓ a (.|x) is differentiable and the auxiliary function is given by ε a (y|x) = yℓ ′ a (y|x)/ℓ a (y|x). This function plays an important role in extreme-value theory since it drives the speed of convergence in [START_REF] Chen | Extreme value behavior of aggregate dependent risks[END_REF] and more generally the bias of extreme-value estimators. Therefore, it may be of interest to specify how it converges to 0. In [START_REF] Goegebeur | Local estimation of the second order parameter in extreme value statistics and local unbiased estimation of the tail index[END_REF], the auxiliary function is supposed to be regularly varying and the estimation of the conditional regular variation index is addressed. Here, we limit ourselves to assuming that for all a ∈ (0, 1/γ(x)), (F.3) |ε a (.|x)| is continuous and ultimately non-increasing.

A Lipschitz condition on the probability density function g of X is also required. For all (x, x ′ ) ∈ R p × R p , the distance between x and x ′ is denoted by d(x, x ′ ) and the following assumption is introduced:

(L) There exists a constant c g > 0 such that |g(x) -g(x ′ )| ≤ c g d(x, x ′ ).
The next assumption is standard in the kernel estimation framework.

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

For ξ > 0, the largest oscillation at point (x, y) ∈ R p × R + * of the conditional moment of order a ∈ [0, 1/γ(x)) is given by

ω(x, y, a, ξ, h) = sup ϕ a (z|x) ϕ a (z|x ′ ) -1 , z ∈ [(1 -ξ)y, (1 + ξ)y] and x ′ ∈ B(x, h) ,
where B(x, h) denotes the ball centred at x with radius h. Finally, for all finite set E, let L(E) = {e i + e j , (e i , e j ) ∈ E × E} ∪ E. We are now in position to establish our main result.

Theorem 1 Suppose (F.1), (F.2), (L) and (K)

hold. Let us introduce 0 ≤ a 1 < a 2 < • • • < a J
where J is a positive integer. For all x ∈ R p such that g(x) > 0 and γ(x) < 1/(2a J ), let us introduce a sequence (α n ) with α n → 0 and nh p α n → ∞ as n → ∞. If there exists ξ > 0 such that

nh p α n h ∨ max a∈L({0,a1,...,a J }) ω(x, ϕ ← 0 (α n |x), a, ξ, h) 2 → 0,
then, the random vector

nh p α n    RCTM aj ,n (α n |x) RCTM aj (α n |x) -1 j∈{1,...,J} , RVaR n (α n |x) RVaR(α n |x) -1    is asymptotically Gaussian, centred, with a (J +1)×(J +1) covariance matrix K 2 2 γ 2 (x)Σ(x)/g(x)
where for (i, j) ∈ {1, . . . , J} 2 we have Σ J+1,j (x) = a j , Σ i,J+1 (x) = a i , Σ J+1,J+1 (x) = 1 and

Σ i,j (x) = a i a j (2 -(a i + a j )γ(x)) (1 -(a i + a j )γ(x)) .
Theorem 1 permits to establish the asymptotic normality for any regression estimator of a risk measure based on arbitrary moments above an extreme conditional quantile. In particular, RCTE n (α n |x), RCVaR λ,n (α n |x) and RSP n (α n |x) only involve the first order moment, their asymptotic normality can be derived under the assumption γ(x) < 1/2:

Corollary 1 Suppose (F.1), (F.2), (L) and (K) hold. For all x ∈ R p such that g(x) > 0 and γ(x) < 1/2, let us introduce a sequence (α n ) with α n → 0 and nh p α n → ∞ as n → ∞. If there exists ξ > 0 such that

nh p α n h ∨ max a∈{0,1,2} ω(x, ϕ ← 0 (α n |x), a, ξ, h) 2 → 0, then nh p α n RCTE n (α n |x) RCTE(α n |x) -1 d -→ N 0, 2γ 2 (x)(1 -γ(x)) 1 -2γ(x) K 2 2 g(x) , nh p α n RCVaR λ,n (α n |x) RCVaR λ (α n |x) -1 d -→ N 0, γ 2 (x)(λ 2 + 2 -2λ -2γ(x)) 1 -2γ(x) K 2 2 g(x) , nh p α n RSP n (α n |x) RSP(α n |x) -1 d -→ N 0, γ 2 (x) 1 -2γ(x) K 2 2 g(x)
.

The RCTV(α n |x) estimator involves the computation of a second order moment, its asymptotic normality requires the stronger condition γ(x) < 1/4.

Corollary 2 Suppose (F.1), (F.2), (L) and (K) hold. For all x ∈ R p such that g(x) > 0 and γ(x) < 1/4, let us introduce a sequence (α n ) with α n → 0 and nh p α n → ∞ as n → ∞. If there exists ξ > 0 such that

nh p α n h ∨ max a∈{0,...,4} ω(x, ϕ ← 0 (α n |x), a, ξ, h) 2 → 0, then nh p α n RCTV n (α n |x) RCTV(α n |x) -1 d -→ N 0, 8(1 -γ(x))(1 -2γ(x))(1 + 2γ(x) + 3γ 2 (x)) (1 -3γ(x))(1 -4γ(x)) K 2 2 g(x)
.

Similarly, the RCTS(α n |x) estimator involves the computation of a third order moment, its asymptotic normality requires the even stronger condition γ(x) < 1/6.

Corollary 3 Suppose (F.1), (F.2), (L) and (K) hold. For all x ∈ R p such that g(x) > 0 and γ(x) < 1/6, let us introduce a sequence (α n ) with α n → 0 and nh p α n → ∞ as n → ∞. If there exists ξ > 0 such that

nh p α n h ∨ max a∈{0,...,6} ω(x, ϕ ← 0 (α n |x), a, ξ, h) 2 → 0, then nh p α n RCTS n (α n |x) RCTS(α n |x) -1 d -→ N 0, V (x) K 2 2 g(x) ,
where

V (x) = 18(1 -13γ(x) + 50γ 2 (x) -44γ 3 (x) -23γ 4 (x) -3γ 5 (x)) (1 -3γ(x))(1 -4γ(x))(1 -5γ(x))(1 -6γ(x)) .
In Theorem 1, the condition nh p α n → 0 provides a lower bound on the level of the risk measure to estimate. This restriction is a consequence of the use of kernel estimator (2) which cannot extrapolate beyond the maximum observation in the ball B(x, h). In consequence, α n must be an order of an extreme quantile within the sample. To overcome this limitation, we propose to adapt Weissman's estimator [START_REF] Weissman | Estimation of parameters and large quantiles based on the k largest observations[END_REF], initially designed for the estimation of unconditional quantiles, to the estimation of the RCTM:

RCTM W a,n (β n |x) = RCTM a,n (α n |x) α n β n aγn(x)
, where a is a fixed value, 0 < β n < α n and γn (x) is an estimator of the conditional tail-index γ(x) (see [START_REF] Gardes | A moving window approach for nonparametric estimation of the conditional tail index[END_REF][START_REF] Gardes | Estimating the conditional tail index by integrating a kernel conditional quantile estimator[END_REF][START_REF] Gardes | Estimation of the conditional tail index using a smoothed local Hill estimator, Extremes[END_REF][START_REF] Goegebeur | Local estimation of the second order parameter in extreme value statistics and local unbiased estimation of the tail index[END_REF][START_REF] Wang | Tail index regression[END_REF]). As illustrated in the next theorem, the extrapolation factor (α n /β n )

aγn(x)
allows us to estimate RCTM of arbitrary small levels β n .

Theorem 2 Suppose the assumptions of Theorem 1 hold together with (F.3). Let us consider γn (x) an estimator of the tail index such that

nh p n α n (γ n (x) -γ(x)) d → N 0, v 2 (x) , with v(x) > 0. If, moreover, (β n ) n≥1 is a positive sequence such that β n → 0, β n /α n → 0 and √ nh p α n ε a (1/β n |x) → 0 as n → ∞, we then have nh p n α n log(α n /β n )   RCTM W a,n (β n |x) RCTM a (β n |x) -1   d → N 0, (av(x)) 2 .
Let us also note that the asymptotic normality of

RVaR W n (β n |x) = RVaR n (α n |x) (α n /β n ) γn(x)
has been established in [START_REF] Daouia | Kernel estimators of extreme level curves[END_REF]. As a consequence, replacing RVaR n by RVaR W n and RCTM a,n by RCTM W a,n in (5) provides estimators for all risk measures considered in this paper adapted to arbitrary small levels. Their asymptotic normality is a simple consequence of Theorem 2. In the next section, a procedure to select the tuning parameters h and α n is introduced and applied to the estimation of risk measures associated to extreme rainfall data.

Application: Risk measures for extreme rainfall data

The rainfall data is described in subsection 4.1. The implementation of the risk measure estimators requires the selection of two tuning parameters. An automatic procedure is proposed in subsection 4.2. Its finite sample performance is assessed on simulated data in subsection 4.3. Finally, the whole methodology is applied on the real data in subsection 4.4.

Problem and data description

The behaviour and the efficiency of our estimators are illustrated on rainfall observations in the Cévennes-Vivarais region (southern part of France). This data set is provided by the French meteorological service Météo-France and consists in daily rainfalls measured at N = 523 raingauge stations from 1958 to 2000. In this context, the variable of interest Y is the daily rainfall measured in millimeters (mm). The number of measurements at each station t ∈ {1, . . . , N } is denoted by n t , the total number of observations being n = N t=1 n t = 5, 513, 734. The covariate X is the three dimensional geographical location (longitude, latitude and altitude). A subset of the coordinates S = {x t = (x 1,t , x 2,t , x 3,t ); t = 1, . . . , N } of the raingauge stations is depicted in Figure 3. Extreme rainfall statistics are often used when a flood has occurred to assess the rarity of such an event. A typical problem is to estimate the amount that will fall on a day of exceptionally heavy rainfall which is expected to occur every T years. Usually, hydrologists are interested in the value T = 100 corresponding to a centenary event. Statistically speaking, the problem is to estimate the T -year return level which is the quantile of level β = 1/(365.25 × T ) of the daily rainfall. The goal of this study is to go further and estimate the average rainfall over the T -year return level which is the RCTE of level β = 1/(365.25 × T ).

Tuning parameters selection

Our estimators of risk measures depend on the two tuning parameters h and α n . The choice of the bandwidth h, which controls the degree of smoothing, is a recurrent issue in non-parametric statistics. Similarly, in extreme-value theory, the choice of the number of upper order statistics, or equivalently α n is of great importance since it raises a compromise between bias and variance. A high value of α n is expected to lead to a large bias (since we move out of the distribution tail) while a small value of α n leads to a large variance, see for instance Theorem 1. Here, we propose a leave-one-out cross validation type procedure to select simultaneously h and α n . To this end,

let us consider A = {α 1 ≤ • • • ≤ α R } such that α 1 > 1/ min(n j ), j = 1, . . . , N , α R < 0.1 and H = {h 1 ≤ • • • ≤ h M },
such that there is at least one observation in the ball B(x, h 1 ) for all x. The principle of the procedure is to select the empirical pair (h emp , α emp ) ∈ H × A for which two different estimations of the tail index γ(x t ) at each station t approximately coincide. The first estimator, denoted by γn,t , is the well-known Hill estimator [START_REF] Hill | A simple general approach to inference about the tail of a distribution[END_REF], it only depends on α n and is uniquely based on the rainfall measures at station t. The second estimator denoted by γn (x t ) is the conditional tail-index estimator introduced in [START_REF] Daouia | Kernel estimators of extreme level curves[END_REF]. It depends both on α n and h and is computed on all the rainfall measures in the ball B(x t , h) except the measurements at the current station t. To summarize, the main idea is to select the pair (h emp , α emp ) for which the local estimations γ(x t ) and the predicted ones γn (x t ) using the neighbour stations are coherent. To be more specific, the algorithm is the following:

1. Loop on all pairs (h i , α j ) ∈ H × A and on all stations t ∈ {1, . . . , N }.

2. Compute the Hill estimator at station t with level α j to obtain γn,t,j .

3. Compute the conditional tail-index estimator using the measures in B(x t , h i )\{x t } with level α j to obtain γn,i,j (x t ).

4. Compute the distance W hi,αj (x t ) = (γ n,t,jγn,i,j (x t )) 2 .

5. End of the loop.

6. The optimal pair is given by (h emp , α emp ) = arg min (hi,αj )∈H×A median{W hi,αj (x t ) , t ∈ {1, . . . , N }}.

Validation on simulation

The previous procedure is tested on two heavy-tailed distributions, the Fréchet distribution and the Burr distribution. The survival function of the Fréchet distribution is F (y|x) = 1-exp(-y -1/γ(x) ), for y ≥ 0, and the associated RCTE can be written

RCTE(β n |x) = 1 β n RVaR(βn|x) -1/γ(x) 0 t γ(x) exp(-t)dt with RVaR(β n |x) = (-log(1 -β n )) -γ(x) .
The survival function of the chosen Burr distribution is given by F (y|x) = 1 + y 1/γ(x) -1 for y ≥ 0, and the associated RCTE is

RCTE(β n |x) = I RVaR(β n |x) -1/γ(x) 1 + RVaR(β n |x) -1/γ(x) , 1 -γ, 1 + γ with RVaR(β n |x) = β -γ(x) n (1 -β n ) γ(x) ,
and where I(r, p, q) = B(r, p, q) B(p, q) with B(r, p, q) = r 0 w p-1 (1w) q-1 dw being the incomplete beta function and B(p, q) the beta function. In this simulation study, we choose the following conditional tail index:

γ : x ∈ (0, 1) → γ(x) = 1 2 1 10 + sin (πx) 11 10 - 1 2 exp (-64(x -1/2) 2 ) .
Note that γ(x) is close to 1/2 when x = 0.3 or x = 0.7. Let z 1 , z 2 and z 3 be respectively the latitude, longitude and altitude normalised in the unit interval. Two choices of covariates x were used for γ(x): x euc := (z 2 1 + z 2 2 )/2 and x alt := z 3 . The tuning parameters are selected in the sets A = {1/(6 × 365.25), 1/(5 × 365.25), . . . , 1/365.25, 4.10 -3 , 6.10 -3 , . . . , 10 -2 , 2.10 -2 , . . . , 10 -1 } and H = {14, 15, . . . , 30}. Recall that the conditional tail index estimator introduced in [START_REF] Daouia | Kernel estimators of extreme level curves[END_REF] is given by γn

(x) = J j=1 [log RVaR n (τ j α n |x) -log RVaR n (τ 1 α n |x)] J j=1 log(τ 1 /τ j ) ,
where (τ j ) j≥1 is a positive non-increasing sequence of weights. Two sequences are investigated:

1. the harmonic sequence defined for all j = 1, . . . , J by τ Ha j = 1/j with J = 9, 2. the geometric sequence defined for all j = 1, . . . , J by τ G j = (1/j) (j/J) with J = 15.

In both cases, the number of terms J was selected to minimize the asymptotic variance of γn (x). Finally, a bi-quadratic kernel was used:

K(x) := K(z 1 , z 2 ) = 15 16 1 -(z 2 1 + z 2 2 ) 2 I{z 2 1 + z 2 2 ≤ 1}.
To assess the performance of our procedure, it is compared to the Oracle (optimal) choice (h opt , α opt ) which is based on the knowledge of the true tail index function:

(h opt , α opt ) = arg min (hi,αj )∈H×A median{V hi,αj (x t ) , t ∈ {1, . . . , N }}, where V hi,αj (x t ) = (γ(x t )γn,i,j (x t )) 2 is the distance to the true tail index function. The selected parameters are displayed in Table 2. It appears that the cross-validation procedure approximately selects the same tuning parameters as the Oracle for all the considered choices of distribution, covariate and weights. Most importantly, one can observe on Figure 1 that the error distributions on the tail index also nearly coincide. This result indicates that the cross-validation procedure is almost as efficient as the Oracle who knows the solution.

Burr distribution

Fréchet distribution x euc and τ Ha , and computed for β = 1/(365.25 × 100) corresponding to a centenary rainfall. In this case, the quality of the estimation is assessed thanks to the relative error:

Q n (x) =   RCTE W n (β|x) RCTE(β|x) -1   2 .
The two histograms of Q n (x t ), t ∈ {1, . . . , N } obtained with (h emp , α emp ) and (h opt , α opt ) are depicted on Figure 2. Both set of parameters yield approximately the same error distribution.

Estimated risk measures on extreme rainfalls

The cross-validation procedure applied to the real data set with τ Ha j yields h emp = 24 and α emp = 1/(365.25 × 3). The estimated conditional tail index is then computed on a grid of 200 × 200 ungauged locations regularly distributed on the Cévennes-Vivarais region, see Figure 3, top panels. Using the asymptotic distribution of γn (x) established in [START_REF] Daouia | Kernel estimators of extreme level curves[END_REF], Corollary 2, pointwise confidence intervals can also be computed. It appears that, for a confidence level of 95%, one can assume that γ(x) < 1/2. At the opposite, the assumption that γ(x) does not depend on x, i.e. γ(x) is constant on the Cévennes-Vivarais region, cannot be accepted. It is then possible to estimate risk measures associated to a 100-year return period. Here, we focus on RVaR n (β n |x) and RCTE n (β n |x) with 5 Appendix: Proofs

Preliminary results

This lemma provides an equivalent of ϕ a (y|x) when y → ∞. We refer to [7, Corollary 3.2] for a similar result in the unconditional case.

Lemma 1 Under (F.1), if y → ∞, then for a ∈ [0, 1/γ(x)), ϕ a (y|x) = 1 1 -aγ(x) y a F (y|x)(1 + o(1)).
Furthermore, under the additional condition (F.2), the derivative ϕ ′ a (.|x) of the function ϕ a (.|x) exists and is a regularly varying function such that

ϕ ′ a (y|x) = aγ(x) -1 γ(x) ϕ a (y|x) y (1 + o(1)).
Proof. First, integrating by part leads to

a ∞ y z a-1 F (z|x)dz = ϕ a (y|x) -y a F (y|x). (9)
Using [28, Eq. (0.32)] together with y → y a-1 F (y|x)

∈ RV a-1/γ(x)-1 , a -1/γ(x) -1 < -1 and y → ∞ yield ∞ y z a-1 F (z|x)dz = γ(x) 1 -aγ(x) y a F (y|x)(1 + o(1)).
Replacing in ( 9) and dividing both sides by

1 1-aγ(x) y a F (y|x) lead to ϕ a (y|x) 1 1-aγ(x) y a F (y|x) -1 + aγ(x) = aγ(x)(1 + o(1)),
which concludes the first part of the proof. Next, under (F.2), derivating both sides of (9) yields

ϕ ′ a (y|x) = y a F ′ (y|x) = y a-1 F (y|x) y F ′ (y|x) F (y|x) ,
and using [28, Corollary of Theorem 0.6], it follows that

y F ′ (y|x) F (y|x) = - 1 γ(x) (1 + o(1)),
which concludes the proof.

As a consequence of Lemma 1 and (1), we obtain RCTE(α n |x)/RVaR(α n |x) ∼ 1/(1γ(x)) which is an extension of the unconditional result, see for instance [START_REF] Hua | Second order regular variation and conditional tail expectation of multiple risks[END_REF]. The second lemma is also of analytical nature. It provides a second order asymptotic expansion of the RCTM.

Lemma 2 Suppose (F.1), (F.2) and (F.3) hold and let 0 < β n < α n be two sequences such that α n → 0 as n → ∞. Then,

|log RCTM a (α n |x) -log RCTM a (β n |x) + aγ(x) log(α n /β n )| = O (log(α n /β n )ε a (1/β n )) .
Proof. Using ( 6) and (F.2), we have

log RCTM a (α n |x) = -aγ(x) log(α n ) + log(c(x)) + 1/αn 1 ε a (u|x) u du,
and consequently

∆ n := log RCTM a (α n |x) -log RCTM a (β n |x) + aγ(x) log(α n /β n ) = 1/αn 1/βn ε a (u|x) u du.
From (F.3), we obtain

|∆ n | ≤ |ε a (1/β n )| log(β n /α n )
and the conclusion follows.

Let us remark that the kernel estimator (2) of the conditional expectation can be rewritten as φa,n (y n |x) = ψa,n (y n |x)/ĝ n (x) where ψa,n (y|x

) = 1 n n i=1 K h (x -X i )Y a i I{Y i > y},
is an estimator of ψ a (y|x) = g(x)ϕ a (y|x) and ĝn (x) is the kernel estimator of the density g(x)

ĝn (x) = 1 n n i=1 K h (x -X i ). ( 10
)
Lemma 3 Suppose (F.1), (F.2), (L) and (K) hold. Let x ∈ R p such that g(x) > 0 and y n → ∞ such that nh p F (y n |x) → ∞.

(i) Let 0 ≤ a < 1/γ(x). If ω(x, y n , a, 0, h) → 0 then E( ψa,n (y n |x)) = ψ a (y n |x)(1 + O(h) + O(ω(x, y n , a, 0, h))). (ii) Let 0 ≤ a 1 < • • • < a J+1 < 1/(2γ(x))
where J is a positive integer and consider sequences (y n,j ), j = 1, . . . , J + 1 such that max j∈{1,...,J+1}

y n,j y n -1 → 0.

If there exists, ξ > 0 such that max a∈L({a1,...,a J+1 }) ω(x, y n , a, ξ, h) → 0 then, the random vector

nh p F (y n |x) ψaj,n (y n,j |x) -E( ψaj,n (y n,j |x)) ψ aj (y n,j |x) j∈{1,...,J+1}
is asymptotically Gaussian, centred, with covariance matrix K 2 2 Σ (1) (x)/g(x) where

Σ (1) i,j (x) = (1 -a i γ(x))(1 -a j γ(x)) 1 -(a i + a j )γ(x) , (i, j) ∈ {1, . . . , J + 1} 2 .
Proof. (i) Since the (X i , Y i ), i = 1, . . . , n are identically distributed, it follows that

E( ψa,n (y n |x)) = R p K h (x -t)ϕ a (y n |t)g(t)dt = S K(u)ϕ a (y n |x -hu)g(x -hu)du, under (K). Let us now consider |E( ψa,n (y n |x)) -ψ a (y n |x)| ≤ ϕ a (y n |x) S K(u)|g(x -hu) -g(x)|du (11) + ϕ a (y n |x) S K(u) ϕ a (y n |x -hu) ϕ a (y n |x) -1 g(x -hu)du. (12)
Under (L), and since g(x) > 0, we have

(11) ≤ ϕ a (y n |x)c g h S d(u, 0)K(u)du = ϕ a (y n |x)O(h). (13)
Besides, in view of ( 13),

(12) ≤ ϕ a (y n |x)ω(x, y n , a, 0, h) S K(u)g(x -hu)du = ϕ a (y n |x)g(x)ω(x, y n , a, 0, h)(1 + o(1)), ≤ ψ a (y n |x)ω(x, y n , a, 0, h)(1 + o(1)). (14)
Combining ( 13) and ( 14) concludes the first part of the proof.

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

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

var(Z 1,n ) = 1 n 2 h 2p Λ 2 n (x) var   J+1 j=1 β j K x -X h Y aj I{Y ≥ y n,j } ψ aj (y n,j |x)   = 1 n 2 h p Λ 2 n (x) β t Bβ,
where B is the (J + 1) × (J + 1) covariance matrix defined by B j,l = A j,l ψ aj (y n,j |x)ψ a l (y n,l |x) , for all (j, l) ∈ {1, . . . , J + 1} 2 and

A j,l = 1 h p cov K x -X h Y aj I{Y ≥ y n,j }, K x -X h Y a l I{Y ≥ y n,l } , = K 2 2 E 1 h p Q x -X h Y aj +a l I{Y ≥ y n,j ∨ y n,l } -h p E(K h (x -X)Y aj I{Y ≥ y n,j })E(K h (x -X)Y a l I{Y ≥ y n,l }),
with Q(.) := K 2 (.)/ K 2 2 also satisfying assumption (K). One can use part (i) of the proof to obtain

A j,l = K 2 2 ψ aj +a l (y n,j ∨ y n,l |x)(1 + O(h) + O(ω(x, y n,j ∨ y n,l , a j + a l , 0, h))) -h p ψ aj (y n,j |x)ψ a l (y n,l |x)(1 + O(h) + O(ω(x, y n,j , a j , 0, h)))(1 + O(h) + O(ω(x, y n,l , a l , 0, h))).
Let ξ > 0 such that max a∈L({a1,...,a J+1 }) ω(x, y n , a, ξ, h) → 0. Remarking that ω(x, y n,j , a j , 0, h) ≤ ω(x, y n , a j , ξ, h) for n large enough, we obtain

A j,l = K 2 2 ψ aj +a l (y n,j ∨ y n,l |x)(1 + O(h) + O(ω(x, y n , a j + a l , ξ, h))) -h p ψ aj (y n,j |x)ψ a l (y n,l |x)(1 + O(h) + O(ω(x, y n , a j , ξ, h)))(1 + O(h) + O(ω(x, y n , a l , ξ, h))).
Now, max (ω(x, y n , a j , ξ, h), ω(x, y n , a l , ξ, h), ω(x, y n , a j + a l , ξ, h)) ≤ max a∈L({a1,...,a J+1 }) ω(x, y n , a, ξ, h), leads to

B j,l = K 2 2 ψ aj +a l (y n,j ∨ y n,l |x) ψ aj (y n,j |x)ψ a l (y n,l |x) 1 + O(h) + O max a∈L({a1,...,a J+1 }) ω(x, y n , a, ξ, h) -h p 1 + O(h) + O max a∈L({a1,...,a J+1 }) ω(x, y n , a, ξ, h) .
Let us recall that, since ψ a (.|x) is regularly varying, it follows that ψ aj (y n,j |x) ∼ ψ aj (y n |x) → 0 for all j ∈ {1, . . . , J + 1}. Lemma 1 thus entails

B j,l = K 2 2 ψ 0 (y n |x) (1 -a j γ(x))(1 -a l γ(x)) 1 -(a j + a l )γ(x) (1 + o(1)) = K 2 2 ψ 0 (y n |x) Σ (1) j,l (x)(1 + o(1)).
Therefore, var(Z 1,n ) ∼ K 2 2 β t Σ (1) (x)β/n. As a preliminary conclusion, the variance of Ψ n converges to K 2 2 β t Σ (1) (x)β. Consequently, using Lyapounov theorem for the asymptotic normality of sums of triangular arrays, it remains to prove that there exists η > 0 such that:

n i=1 E |Z i,n | 2+η = nE |Z 1,n | 2+η → 0.
Straightforward calculations lead to

E |Z 1,n | 2+η = 1 nΛ n (x) 2+η E J+1 j=1 β j K h (x -X) Y aj I{Y ≥ y n,j } ψ aj (y n,j |x) -E   J+1 j=0 β j K h (x -X) Y aj I{Y ≥ y n,j } ψ aj (y n,j |x)   2+η .
Besides, for every pair of random variables (T 1 , T 2 ) with finite (2 + η)th order moments, one has

E |T 1 + T 2 | 2+η ≤ 2 2+η max i={1,2} E |T i | 2+η , leading to E |Z 1,n | 2+η ≤ 2 nΛ n (x) 2+η E J+1 j=1 β j K h (x -X) Y aj I{Y ≥ y n,j } ψ aj (y n,j |x) 2+η .
Lemma 1 and y n,j = y n (1 + o( 1)) for all j ∈ {1, . . . , J + 1} yield

E |Z 1,n | 2+η ≤ 2 nΛ n (x)ψ 0 (y n |x) 2+η × E J+1 j=1 β j K h (x -X) I{Y ≥ y n,j } Y y n,j aj (1 -a j γ(x)) 2+η (1 + o(1)).
Letting ã = max{a 1 , . . . , a J+1 } and ỹn = min{y n,1 , . . . , y n,J+1 }, it follows that for n large enough,

nE |Z 1,n | 2+η ≤ 2n 2(1 -ãγ(x)) nh p Λ n (x)ψ 0 (y n |x)ỹ ã n 2+η J+1 j=1 |β j | 2+η E K x -X h Y ãI{Y ≥ ỹn } 2+η .
Choosing η such that 0 < η < -2 + 1/(ãγ(x)), (i) implies that

E K x -X h Y ãI{Y ≥ ỹn } 2+η = h p K 2+η 2+η E N h (x -X) Y ã(2+η) I{Y ≥ ỹn } , = h p K 2+η 2+η ψ ã(2+η) (ỹ n |x)(1 + o(1)),
since N (.) := K 2+η (.)/ K 2+η 2+η also fulfils assumption (K). Using Lemma 1 and the fact that ỹn = y n (1 + o(1)), we obtain nE |Z 1,n | 2+η = O (Λ η n (x)) → 0 as n → ∞ which concludes the proof. The asymptotic behaviors of the estimators φa,n (.|x) and φ← a,n (.|x) are established in the following two propositions.

Proposition 1 Suppose (F.1), (F.2), (L) and (K) hold. Let x ∈ R p such that g(x) > 0 and

0 ≤ a 1 < • • • < a J+1 < 1/(2γ(x))
where J is a positive integer. Consider y n → ∞ such that nh p F (y n |x) → ∞ as n → ∞ and sequences (y n,j ), j ∈ {1, . . . , J + 1} such that max j∈{1,...,J+1}

y n,j y n -1 → 0.
If there exists ξ > 0 such that nh p F (y n |x) h ∨ max a∈L({a1,...,a J+1 }) ω(x, y n , a, ξ, h) 2 → 0 then, the random vector

nh p F (y n |x) φaj,n (y n,j |x) ϕ aj (y n,j |x) -1 j∈{1,...,J+1}
is asymptotically Gaussian, centred, with covariance matrix K 2 2 Σ (1) (x)/g(x).

Proof. Keeping in mind the notations of Lemma 3, the following expansion holds

Λ -1 n (x) J+1 j=1 β j φaj,n (y n,j |x) ϕ aj (y n,j |x) -1 = ∆ 1,n + ∆ 2,n -∆ 3,n ĝn (x) , (15)
where

∆ 1,n = g(x)Λ -1 n (x) J+1 j=1
β j ψaj,n (y n,j |x) -E( ψaj,n (y n,j |x))

ψ aj (y n,j |x) , ∆ 2,n = g(x)Λ -1 n (x) J+1 j=1
β j E( ψaj,n (y n,j |x))ψ aj (y n,j |x) ψ aj (y n,j |x) ,

∆ 3,n =   J+1 j=1 β j   Λ -1 n (x) (ĝ n (x) -g(x)) .
Thus, from Lemma 3(ii), the random term ∆ 1,n can be rewritten as

∆ 1,n = g(x) K 2 β t Σ (1) (x)βξ n , (16)
where ξ n converges to a standard Gaussian random variable. The non-random term ∆ 2,n is controlled with Lemma 3(i):

∆ 2,n = O hΛ -1 n (x) + O Λ -1 n (x
Finally, ∆ 3,n is a classical term in kernel density estimation, which can be bounded by [START_REF] Daouia | Kernel estimators of extreme level curves[END_REF], Lemma 4:

∆ 3,n = O(hΛ -1 n (x)) + O P (Λ -1 n (x)(nh p ) -1/2 ), = O nh p+2 F (y n |x) 1/2 + O P ( F (y n |x)) 1/2 = o P (1). (18)
Collecting ( 15)-( 18), it follows that ĝn (x)Λ - is asymptotically Gaussian, centred, with covariance matrix K 2 2 Σ (2) (x)/g(x) where

Σ (2) i,j (x) = γ 2 (x) 1 -(a i + a j )γ(x)
, (i, j) ∈ {1, . . . , J + 1} 2 .

Proof. Introduce for j ∈ {1, . . . , J + 1},

σ n,j (x) = ϕ ← aj (α n,j |x)(nh p α n ) -1/2 , v n,j (x) = α -1 n,j γ(x) 1 -a j γ(x) (nh p α n ) 1/2 , W n,j (x) = v n,j (x) φaj (ϕ ← aj (α n,j |x) + σ n,j (x)z j )|x) -ϕ aj (ϕ ← aj (α n,j |x) + σ n,j (x)z j )|x , t n,j (x) = v n,j (x) α n,j -ϕ aj (ϕ ← aj (α n,j |x) + σ n,j (x)z j )|x ,
where (z 1 , . . . , z J+1 ) ∈ R J+1 . We examine the asymptotic behavior of the cumulative distribution function defined by

Φ n (z 1 , . . . , z J+1 ) = P   J+1 j=1 σ -1 n,j (x)( φ← aj ,n (α n,j |x) -ϕ ← aj (α n,j |x)) ≤ z j   , = P   J+1 j=1 {W n,j (x) ≤ t n,j (x)}   .
Let us first focus on the non-random terms t n,j (x), j ∈ {1, . . . , J + 1}. From Lemma 1, for all a ∈ [0, 1/(2γ(x))), the function ϕ a (.|x) is differentiable and thus, for each j ∈ {1, . . . , J + 1} there exists θ n,j ∈ (0, 1) such that

ϕ aj ϕ ← aj (α n,j |x)|x -ϕ aj ϕ ← aj (α n,j |x) + σ n,j (x)z j |x = -σ n,j (x)z j ϕ ′ aj (r n,j |x), (19)
where r n,j = ϕ ← aj (α n,j |x)+θ n,j σ n,j (x)z j . It is thus clear that r n,j ∼ ϕ ← aj (α n,j |x) → ∞ and Lemma 1 yields

ϕ ′ aj (r n,j |x) = (a j γ(x) -1)α n,j γ(x)ϕ ← aj (α n,j |x) (1 + o(1)). (20)
In view of ( 19) and ( 20), we end up with

t n,j (x) = (1 -a j γ(x))v n,j (x)σ n,j (x)α n,j z j γ(x)ϕ ← aj (α n,j |x) (1 + o(1)) = z j (1 + o(1)). (21)
Let us now turn to the random terms W n,j (x), j ∈ {1, . . . , J + 1}. Clearly, sequences y n,j := ϕ ← aj (α n,j |x) + σ n,j (x)z j , j = 1, . . . , J + 1 and y n := ϕ ← 0 (α n |x) satisfy the assumptions of Proposition 1 and consequently,

W n,j (x) = γ(x) 1 -a j γ(x) ϕ aj (ϕ ← aj (α n,j |x) + σ n,j (x)z j |x) α n,j (nh p α n ) 1/2 φaj (y n,j |x) ϕ aj (y n,j |x) -1 .
Moreover, since ϕ a (.|x) is regularly varying, the following equivalences hold,

ϕ aj (ϕ ← aj (α n,j |x) + σ n,j (x)z j |x) α n,j = ϕ aj (ϕ ← aj (α n,j |x)(1 + o P (1))|x) α n,j = 1 + o P (1).
As a consequence of Slutsky's theorem, the random vector (W n,1 , . . . , W n,J+1 ) is equal to A(x)ξ n where

A(x) = diag γ(x) 1 -a 1 γ(x) , . . . , γ(x) 1 -a J+1 γ(x) ,
and ξ n is a (J + 1)-random vector converging to a centred Gaussian random variable with covariance matrix K 2 2 Σ (1) (x)/g(x). Taking account of (21), we obtain that Φ n (z 1 , . . . , z J+1 ) converges to the cumulative distribution function of a centred Gaussian distribution with covariance matrix

K 2 2 A(x)Σ (1) (x)A(x)/g(x) = K 2 2 Σ (2) (x)/g(x)
which is the desired result.

Proofs of main results

Proof of Theorem 1. Let us introduce for j ∈ {1, . . . , J},

v n,j (x) = (1 -a j γ(x))(nh p α n ) 1/2 γ(x)ϕ ← 0 (α n |x) , σ n,j (x) = ϕ aj (ϕ ← 0 (α n |x)|x)(nh p α n ) -1/2 , σ n,0 (x) = ϕ ← 0 (α n |x)(nh p α n ) -1/2 , t n,j = v n,j (x) ϕ ← 0 (α n |x) -ϕ ← aj ϕ aj (ϕ ← 0 (α n |x)|x) + σ n,j (x)z j |x , W n,j (x) = v n,j (x) φ← aj ,n ϕ aj (ϕ ← 0 (α n |x)|x) + σ n,j (x)z j |x -ϕ ← aj ϕ aj (ϕ ← 0 (α n |x)|x) + σ n,j (x)z j |x , W (0) n,j (x) = v n,j (x) φ← 0,n (α n |x) -ϕ ← 0 (α n |x) , W (0)
n,0 (x) = σ -1 n,0 (x) φ← 0,n (α n |x) -ϕ ← 0 (α n |x) ,
where (z 0 , z 1 , . . . , z J ) ∈ R J+1 . We examine the asymptotic behavior of the cumulative distribution function defined by

Φ n (z 0 , z 1 , . . . , z J ) = P      J j=1 σ -1 n,j (x)( φaj,n ( φ← 0,n (α n |x)|x) -ϕ aj (ϕ ← 0 (α n |x)|x)) ≤ z j    W (0) n,0 (x) ≤ z 0   , = P      J j=1 W n,j (x) -W (0) n,j (x) ≤ t n,j    W (0) n,0 (x) ≤ z 0   .
Let us first focus on the non-random terms t n,j (x), j = 1, . . . , J. From Lemma 1, for all a ∈ [0, 1/(2γ(x))), ϕ ← a (.|x) is a differentiable regularly varying function such that

(ϕ ← a ) ′ (y n |x) = 1 ϕ ′ a (ϕ ← a (y n |x)|x) = γ(x)ϕ ← a (y n |x) (aγ(x) -1)y n (1 + o(1)), (22)
as n → ∞. For all j ∈ {1, . . . , J}, a first order Taylor expansion leads to:

ϕ ← aj (ϕ aj (ϕ ← 0 (α n |x)|x)|x) -ϕ ← aj ϕ aj (ϕ ← 0 (α n |x)|x) + σ n,j (x)z j |x = -σ n,j (x)z j q n,j (x), where q n,j (x) = (ϕ ← aj ) ′ (ϕ aj (ϕ ← 0 (α n |x)|x) + θ n,j σ n,j (x)z j |x) with (θ n,1 , . . . , θ n,J ) ∈ (0, 1) J . Since σ n,j (x)/ϕ aj (ϕ ← 0 (α n |x)|x) = (nh p α n ) -1/2 → 0 as n → ∞, (22) entails that q n,j (x) = γ(x)ϕ ← 0 (α n |x) (a j γ(x) -1)ϕ aj (ϕ ← 0 (α n |x)|x) (1 + o(1)).
Hence,

ϕ ← aj (ϕ aj (ϕ ← 0 (α n |x)|x)|x) -ϕ ← aj ϕ aj (ϕ ← 0 (α n |x)|x) + σ n,j (x)z j |x = z j v n,j (x) (1 + o(1)), (23)
which shows that for all j ∈ {1, . . . , J}, t n,j → z j as n → ∞. Let us now turn to the random terms W n,j (x), j = 1, . . . , J. Clearly,

W n,j (x) = 1 -a j γ(x) γ(x) (nh p α n ) 1/2 φ← aj ,n ϕ aj (ϕ ← 0 (α n |x)|x) + σ n,j (x)z j |x ϕ ← aj ϕ aj (ϕ ← 0 (α n |x)|x) + σ n,j (x)z j |x -1 (1 + o(1)), since, from (23), ϕ ← aj ϕ aj (ϕ ← 0 (α n |x)|x) + σ n,j (x)z j |x ϕ ← 0 (α n |x) = 1 + z j ϕ ← 0 (α n |x)v n,j (x) (1 + o(1)) = 1 + o(1)
.

Furthermore, we have

W (0) n,j (x) = 1 -a j γ(x) γ(x) (nh p α n ) 1/2 φ← 0,n (α n |x) ϕ ← 0 (α n |x) -1 .
As a consequence, applying Proposition 2 with a J+1 = 0, α n,j = ϕ aj (ϕ ← 0 (α n |x)|x) + σ n,j (x)z j for j = 1, . . . , J and α n,J+1 = α n entails

W n,j (x) -W (0) n,j (x) j=1,...,J , W (0) n,0 (x) = M (x)ξ n ,
where M is the (J + 1) × (J + 1) matrix defined by

M (x) =      Ã(x) c 1 (x)
. . .

c J (x) 0 • • • 0 1     
with Ã(x) = diag 1a 1 γ(x) γ(x) , . . . , 1a J γ(x) γ(x) and c j = -1a j γ(x) γ(x) , j ∈ {1, . . . , J} and where ξ n is a (J + 1)-random vector asymptotically Gaussian, centred with covariance K 2 2 Σ (2) (x)/g(x). Since for each j ∈ {1, . . . , J}, t n,j → z j as n → ∞, the cumulative distribution function Φ n converges to a centred Gaussian cumulative distribution function with covariance matrix K 2 2 M (x)Σ (2) (x)M (x) t /g(x) = K 2 2 γ 2 (x)Σ(x)/g(x), which is the desired result.

Proof of Corollary 2. Clearly, from Theorem 1 one has for i = 1, 2,

RCTM i,n (α n |x) = RCTM i (α n |x) 1 + (nh p α n ) -1/2 ξ i,n ,
where the random vector (ξ 1,n , ξ 2,n ) is asymptotically Gaussian, centred with covariance matrix Σ (3) defined by Σ

i,j = ijγ 2 (x) 2 -(i + j)γ(x) 1 -(i + j)γ(x)

jx

  alt and τ G j x euc and τ G j x alt and τ Ha j h opt = 22 h opt = 24 h opt = 22 h opt = 26 h emp = 24 h emp = 22 h emp = 20 h emp = 24 α opt = 1/365.25 α opt = 1/(365.25 × 2) α opt = 1/(365.25 × 3) α opt = 1/365.25 α emp = 0.001 α emp = 1/365.25 α emp = 1/365.25 α emp = 0.004

β n = 1 /

 1 (365.25 × 100). The associated estimators RVaR W n (β n |x) and RCTE W n (β n |x) are displayed on Figure 3, bottom panels. The estimated 100-year return level RVaR W n (β n |x) is similar to the results obtained in [6] using kriging methods. More interestingly, the RCTE W n (β n |x) can be 150 millimeters higher than the RVaR W n (β n |x) on the mountains area.

K 2 2 g.Figure 1 :Figure 2 :

 212 Figure 1: Histogram of the errors V hopt,αopt (x t ) (solid line, white bars) and V hemp,αemp (x t ) (dotted lines, grey bars) computed on simulated data (Burr distribution on the top panel, Fréchet distribution on the bottom panel). Upper left: x euc and τ Ha j , upper right: x alt and τ G j , bottom left: x euc and τ G j , bottom right: x alt and τ Ha j .

Figure 3 :

 3 Figure 3: Upper left: map of the Cévennes-Vivarais region, horizontally: longitude (km), vertically: latitude (km), the color scale represents the altitude (m), the white dots represent some raingauge stations, upper right: γn (x), bottom left: RVaR W n (β n |x) for a 100-year return period, bottom right: RCTE W n (β n |x) for a 100-year return period.

Table

Table 2 :

 2 Results of the selection procedure

Similar results can be observed on the extrapolated RCTE defined by

RCTE W n (β|x) = RCTE n (α|x) α β γn(x)

  Suppose (F.1), (F.2), (L) and (K)hold. Let x ∈ R p such that g(x) > 0 and let 0 ≤ a 1 < • • • < a J+1 < 1/(2γ(x))where J is a positive integer. Consider α n → 0 such that nh p α n → ∞ as n → ∞. Let (α n,j ), j = 1, . . . , J + 1 be sequences such that

		1 n (x)	J+1 j=1	β j		φaj,n (y n,j |x) ϕ aj (y n,j |x)	-1 = g(x) K 2 β t Σ (1) (x)βξ n + o P (1).
	Finally, ĝn (x)	P -→ g(x) yields
		nh p F (y n |x)	J+1 j=1	β j	φaj,n (y n,j |x) ϕ aj (y n,j |x)	-1 = K 2	β t Σ (1) (x)β g(x)	ξ n + o P (1),
	and the result is proved.
	Proposition 2 max j∈{1,...,J+1}	ϕ ← aj (α n,j |x) ϕ ← 0 (α n |x)	-1 → 0,
	If there exists ξ > 0 such that nh p α n h ∨ max a∈L({a1,...,a J+1 }) ω(x, ϕ ← 0 (α n |x), a, ξ, h) the random vector	2 → 0 then,
						nh p α n	φ← aj ,n (α n,j |x) ϕ ← aj (α n,j |x)	-1	j∈{1,...,J+1}

Proof of Corollary 3. Clearly, from Theorem 1 one has for i = 1, 2, 3, RCTM i,n (α n |x) = RCTM i (α n |x) 1 + (nh p α n ) 1/2 ξ i,n , where the random vector (ξ 1,n , ξ 2,n , ξ 3,n ) is asymptotically Gaussian, centred with covariance matrix Σ (4) defined by

From the proof of Corollary 2, it appears that

and thus

Clearly, from Theorem 1,

is asymptotically Gaussian, centred with variance

and the result is proved.

Proof of Theorem 2. The proof is based on the following expansion:

Let us consider the three terms separately. Under the hypotheses of Theorem 2, it is clear that

As a consequence, Q n,2 P → 0 when n → ∞. Finally, Lemma 2 entails Q n,3 = O nh p n α n ε a (1/β n |x) which converges to 0 by assumption.


David Delahaye
email: david.delahaye@cnam.fr

Damien Doligez
email: damien.doligez@inria.fr

Frédéric Gilbert
email: frederic.charles.gilbert@inria.fr

Pierre Halmagrand
email: pierre.halmagrand@inria.fr

Olivier Hermant
email: olivier.hermant@mines-paristech.fr

Proof Certification in Zenon Modulo: When Achilles Uses Deduction Modulo to Outrun the Tortoise with Shorter Steps *

We present the certifying part of the Zenon Modulo automated theorem prover, which is an extension of the Zenon tableau-based first order automated theorem prover to deduction modulo. The theory of deduction modulo is an extension of predicate calculus, which allows us to rewrite terms as well as propositions, and which is well suited for proof search in axiomatic theories, as it turns axioms into rewrite rules. In addition, deduction modulo allows Zenon Modulo to compress proofs by making computations implicit in proofs. To certify these proofs, we use Dedukti, an external proof checker for the λΠ-calculus modulo, which can deal natively with proofs in deduction modulo. To assess our approach, we rely on some experimental results obtained on the benchmarks provided by the TPTP library.

Introduction

When dealing with automated theorem provers, the critical point is to ensure the soundness of the corresponding implementation: completeness is a minor concern and is often dismissed in favor of some strategies that focus efficiently on some categories of problems. To do so, two options can be considered: either formally proving the soundness of the implemented proof search method, or requiring the implementation to produce proof certificates, which can be checked to reach an appropriate degree of confidence.

The first option is clearly a difficult task, even though successful experiments have been carried out in this domain, in particular for SAT solving algorithms. This task is difficult because the tools used to formalize the soundness of algorithms provide constrained environments (frequently, functions must be total and must terminate), where it is arduous to encode the given algorithm that needs to be certified. In addition, it should be noted that the more sophisticated the algorithm is, the more difficult the soundness proof will be. As a consequence, the soundness proofs are generally made over variants of the initial algorithm, which may be far from the actual code, omitting the most error-prone parts of the algorithm.

The second option alleviates some of the difficulties inherent to the previous option, as it imposes no constraint over the proof search algorithm, which can be as clever as needed.

However, we must prove the soundness for each result produced by the algorithm, which must therefore provide enough information to build these soundness proofs. This information is generally called proof certificate, and it has become fashionable these days that certifying automated theorem provers satisfy the De Bruijn criterion (formulated by Barendregt in [START_REF] Barendregt | Autarkic Computations in Formal Proofs[END_REF]), i.e. they generate proof certificates (or even directly proofs) in a format that can be independently checked by external proof tools. In addition, proof certificates must be of good quality and some principles tend to emerge in this domain. One of them is the Poincaré principle (also formulated by Barendregt in [START_REF] Barendregt | Autarkic Computations in Formal Proofs[END_REF]), which is directly related to the size of proofs and which states that traces of computation should not be included in proof certificates. In this case, the external proof checker is expected to redo computations and must be therefore adapted.

In this paper, we present the certifying part of the Zenon Modulo automated theorem prover [START_REF] Delahaye | Zenon Modulo: When Achilles Outruns the Tortoise using Deduction Modulo[END_REF], which is an extension of the Zenon tableau-based first order automated theorem prover [START_REF] Bonichon | Zenon: An Extensible Automated Theorem Prover Producing Checkable Proofs[END_REF] to deduction modulo. The theory of deduction modulo is an extension of predicate calculus, which allows us to rewrite terms as well as propositions, and which is well suited for proof search in axiomatic theories, as it turns axioms into rewrite rules. This way, we turn proof search among the axioms into computations, avoiding unnecessary blowups, and we shrink the size of proofs by recording only their meaningful steps (see [START_REF] Burel | Efficiently Simulating Higher-Order Arithmetic by a First-Order Theory Modulo[END_REF]). Zenon Modulo, just like Zenon, adopts a certifying approach, and is able to produce proofs to be checked by an external proof tool (in this sense, Zenon Modulo satisfies the De Bruijn criterion). In addition, deduction modulo allows Zenon Modulo to compress proofs by producing proofs with no trace of computation (in this sense, Zenon Modulo verifies the Poincaré principle). However, the external proof checker must be able to deal with these proofs in deduction modulo. For that purpose, we use Dedukti [START_REF] Boespflug | The λΠ-Calculus Modulo as a Universal Proof Language[END_REF], which is a proof checker based on the λΠ-calculus modulo.

This paper is organized as follows: in Sec. 2, we present the proof generation of Zenon Modulo, and particularly focus on the proof compression aspect (compared to Zenon); in Sec. 3, we describe the proof verification of Zenon Modulo using Dedukti; for both sections, we rely on some experimental results obtained on the benchmarks provided by the TPTP library [START_REF] Sutcliffe | The TPTP Problem Library and Associated Infrastructure: The FOF and CNF Parts, v3.5.0[END_REF].

Proof Compression in Zenon Modulo

In order to validate our implementation of Zenon Modulo, we compare proofs found by both Zenon and Zenon Modulo where rewriting is used at least once in the latter (otherwise, proofs are identical). This leads us to examine a subset of 624 first order problems (FOF category) of the TPTP library (v5.5.0) out of the 1,446 problems proved by both Zenon and Zenon Modulo.

We compare the number of nodes of the resulting proof trees after pruning of the useless formulas (pruning is a technique that allows us to minimize the size of the proof search tree by determining the useful formulas to close the tree; see [START_REF] Bonichon | Zenon: An Extensible Automated Theorem Prover Producing Checkable Proofs[END_REF] for more information). The average proof size reduction obtained by Zenon Modulo over the 624 FOF problems is 6.8% (with a maximum of 91.4%). Over the 110 problems from the SET category (problems of set theory), the average proof size reduction reaches 21.6% (with a maximum of 84.6%). In Fig. 1, we present the average proof size reduction depending on the size of the initial proofs in Zenon. We cluster the proofs of Zenon in 10 deciles and compute the associated average proof size reductions. The results for the first decile (i.e. the smallest Zenon proofs, between 1 and 3 proof nodes for FOF and between 1 to 6 proof nodes for SET) are negative: about -42% for FOF and -65% for SET. These results are not shown in Fig. 1, since it is of dubious interest to turn small proofs into bigger yet still small proofs. We observe a positive correlation between the average proof size reduction and the size of Zenon proofs, i.e. the bigger Zenon proofs are, the higher the average proof size reduction is. We obtain the highest reduction for the last decile To illustrate the proof compression in Zenon Modulo, we present, in Fig. 2, the proofs found by Zenon and Zenon Modulo for the problem SET705+4, which states that every set is included in its powerset, i.e. ∀A A ∈ P(A). The axioms used to complete this proof are the following:

∀X, A (X ∈ P(A) ⇔ X ⊆ A) (powerset) ∀A, B (A ⊆ B ⇔ ∀X (X ∈ A ⇒ X ∈ B)) (subset)
Zenon uses the powerset axiom as an hypothesis and the subset axiom as a definition, called by the p-unfold ¬ rule. In contrast, Zenon Modulo transforms these two axioms into propositional rewrite rules as follows:

X ∈ P(A) -→ X ⊆ A (powerset) A ⊆ B -→ ∀X (X ∈ A ⇒ X ∈ B) (subset)
We observe a reduction of 55% of proof nodes in Zenon Modulo. The resulting proof tree is also much simpler since there is no branch creation.

Proof Verification with Dedukti

Integrating some rewriting capabilities on top of the two existing outputs of Zenon for Coq and Isabelle would lose compression of Sec. 2 by providing all the rewrite steps. We chose instead to use a proof checker with rewriting capabilities: Dedukti, based on the λΠ-calculus modulo.

Proof in Zenon: As the target logic is constructive, a translation from classical logic must be therefore implemented. Since Dedukti offers many facilities to define shallow embeddings (thanks to rewrite rules), we opted for a double-negation translation rather than for a non-computational axiomatization. Introducing double-negations at each place à la Kolmogorov [START_REF] Troelstra | Constructivism in Mathematics: An Introduction[END_REF] is a linear solution with respect to proof size, but appeared impracticable. This is why we forged a translation that introduces as few double-negation as possible, improving over previously existing work by Gödel, Gentzen, and Kuroda [START_REF] Troelstra | Constructivism in Mathematics: An Introduction[END_REF]. Moreover, we treat differently the formulas along the side of the sequent they appear, following the remark that left-rules are identical in classical and intuitionistic sequent calculi. Combined with the use of rewrite rules, this has an important impact on the size of the generated certificates, and on the time it takes to check them.

¬(∀A (A ∈ P(A))), ∀X, A (X ∈ P(A) ⇔ X ⊆ A) δ ¬∀ ¬(τ1 ∈ P(τ1)) γ ∀inst ∀A (τ1 ∈ P(A) ⇔ τ1 ⊆ A) γ ∀inst τ1 ∈ P(τ1) ⇔ τ1 ⊆ τ1 β⇔ ¬(τ1 ⊆ τ1) p-unfold ¬ ¬(∀X (X ∈ τ1 ⇒ X ∈ τ1) δ ¬∀ ¬(τ2 ∈ τ1 ⇒ τ2 ∈ τ1) α¬⇒ ¬(τ2 ∈ τ1), τ2 ∈ τ1 ⊙ ⊙ τ1 ∈ P(τ1) ⊙ ⊙ Proof in Zenon Modulo: ¬(∀A (A ∈ P(A))) ⋆ ¬(∀X (X ∈ τ1 ⇒ X ∈ τ1) δ ¬∀ ¬(τ2 ∈ τ1 ⇒ τ2 ∈ τ1) α¬⇒ ¬(τ2 ∈ τ1), τ2 ∈ τ1 ⊙ ⊙ where: τ1 = ǫ(A).¬(A ∈ P(A)) τ2 = ǫ(X).¬(X ∈ τ1 ⇒ X ∈ τ1) and: ⋆ = δ ¬∀ , powerset, subset
As shown in Tab. 1, this approach is currently able to check with success about 90% of the proofs found by Zenon Modulo, and no certificate has been explicitly rejected. As for failures, there are two main reasons for which some proofs fail to be checked properly. The first one (see the column "Dedukti Failure") is that Dedukti itself does not succeed to check its input. This is mainly due to the produced rewrite system, which is ill-formed according to Dedukti. In particular, this is the case when the rewrite system is non-terminating or contain some non-linear rules (Dedukti may also fail in presence of non-confluent rewrite systems). Since the heuristic that transforms automatically axioms into rewrite rules used by Zenon Modulo does not perform any verification over the produced rewrite system, we can imagine to reduce the proportion of these cases of Dedukti failure by carrying out a (mainly syntactical) study of the rewrite system (possibly using some dedicated external tools) before calling Dedukti.

The second reason of failures (see the column "Backend Issue") is that the backend fails to produce a result (timeout or memory overflow) due to internal inefficiencies. By making the translation depend of the side of the sequent where the formula appears, we lose modularity. In other words, the cut rule is not "directly" admissible. Unfortunately, cuts are introduced by Zenon Modulo during proof search and proof optimization, before the intervention of the backend. To treat them, we have to eliminate them "by hand" up to a point where both left and right translations become equal. This step has not been optimized yet, and may take exponential time and size. Even worse, in deduction modulo, this process may loop whatever the reduction strategy. Another increase in the complexity of the translation is introduced by the analysis of the proof tree as a whole to decide whether the double-negation introduced by the translation of

Conclusion

We have described the certifying part of the Zenon Modulo automated theorem prover, which is an extension of the Zenon tableau-based first order automated theorem prover to deduction modulo. Thanks to deduction modulo, Zenon Modulo turns axioms into rewrite rules and is therefore able to produce smaller proofs (compared to Zenon) with no trace of computation. However, the external proof checker used to certify these proofs must be adapted, and we use Dedukti, an external proof checker based on the λΠ-calculus modulo, which can deal natively with proofs in deduction modulo. To assess our approach, we have presented some experimental results obtained on the benchmarks provided by the TPTP library. These results show an increasing average proof size reduction with the size of the proof, while a large part of the proofs produced by Zenon Modulo are verified by Dedukti. As future work, we plan to improve the Dedukti backend by introducing more algorithmic optimizations, but an alternative solution could be to develop a more "syntactical" backend and unplug some optimizations, in order to avoid dependency on cut-elimination and proof inspection.

Figure 1 :

 1 Figure 1: Proof Compression in Zenon Modulo over the FOF and SET Categories

Figure 2 :

 2 Figure 2: Proofs of Problem SET705+4 in Zenon and Zenon Modulo

Table 1 :

 1 Proof Verification with Dedukti over a Subset of the FOF Category formulas is necessary: in fact, a double negation corresponds to the use of the excluded-middle principle, and many classical proofs do not rely on this.

	FOF (624 problems)	Dedukti Success	Dedukti Failure	Backend Issue
	Problems	559	5	60
	Rate	89.6%	0.8%	9.6%

* This work is supported by the BWare project [7] (ANR-12-INSE-0010) funded by the INS programme of the French National Research Agency (ANR).