Jean-François Filipot

Je souhaite tout d'abord remercier chaleureusement les membres de mon jury de thèse. Anthony Davison et Carlo Gaetan, j'ai été très honorée que vous ayez accepté de relire ce manuscrit. Je vous remercie pour vos retours pertinents et vos conseils avisés. Merci également à

L'estimation précise des probabilités d'occurrence des événements extrêmes environnementaux est une préoccupation majeure dans l'évaluation des risques. Pour l'ingénierie côtière par exemple, le dimensionnement de structures implantées sur ou à proximité des côtes doit être tel qu'elles résistent aux événements les plus sévères qu'elles puissent rencontrer au cours de leur vie. Cette thèse porte sur la simulation d'événements extrêmes multivariés, motivée par des applications aux hauteurs significatives de vagues, et sur l'évaluation de modèles de prédiction d'occurrence d'événements extrêmes.

Dans la première partie du manuscrit, nous proposons et étudions un simulateur stochastique qui génère conjointement, en fonction de certaines conditions d'état de mer au large, des extrêmes de hauteur significative de vagues (H s ) au large et à la côte. Pour cela, nous nous appuyons sur l'approche par dépassements de seuils bivariés et nous développons un algorithme de simulation non-paramétrique de lois de Pareto généralisées bivariées. À partir de ce simulateur d'événements cooccurrents, nous dérivons un modèle de simulation condition-nel. Les deux algorithmes de simulation sont mis en oeuvre sur des expériences numériques et appliqués aux extrêmes de H s près des côtes bretonnes françaises. Un autre développement est traité quant à la modélisation des lois marginales des H s . Afin de prendre en compte leur non-stationnaritée, nous adaptons une extension de la loi de Pareto généralisée, en considérant l'effet de la période et de la direction pic sur ses paramètres.

La deuxième partie de cette thèse apporte un développement plus théorique. Pour évaluer différents modèles de prédiction d'extrêmes, nous étudions le cas spécifique des classifieurs binaires, qui constituent la forme la plus simple de prévision et de processus décisionnel : un événement extrême s'est produit ou ne s'est pas produit. Des fonctions de risque adaptées à la classification binaire d'événements extrêmes sont développées, ce qui nous permet de répondre à notre deuxième question. Leurs propriétés sont établies dans le cadre de la variation régulière multivariée et de la variation régulière cachée, permettant de considérer des formes plus fines d'indépendance asymptotique. Ces développements sont ensuite appliqués aux débits de rivière extrêmes.

Title: Simulation and assessment of multivariate extreme models for environmental data Keywords: Multivariate extremes; Simulation of joint extremes; Multivariate and hidden regular variation; Evaluation of extreme classifiers; Significant wave height; River discharges Abstract: Accurate estimation of the occurrence probabilities of extreme environmental events is a major issue for risk assessment. For example, in coastal engineering, the design of structures installed at or near the coasts must be such that they can withstand the most severe events they may encounter in their lifetime. This thesis focuses on the simulation of multivariate extremes, motivated by applications to significant wave height, and on the evaluation of models predicting the occurrences of extreme events.

In the first part of the manuscript, we propose and study a stochastic simulator that, given offshore conditions, produces jointly offshore and coastal extreme significant wave heights (H s ). We rely on bivariate Peaks over Threshold and develop a non-parametric simulation scheme of bivariate generalised Pareto distributions. From such joint simulator, we derive a conditional simulation model. Both simulation algorithms are applied to numerical experiments and to extreme H s near the French Brittany coast. A further development is addressed regarding the marginal modelling of H s . To take into account non-stationarities, we adapt the extended generalised Pareto model, letting the marginal parameters vary with the peak period and the peak direction.

The second part of this thesis provides a more theoretical development. To evaluate different prediction models for extremes, we study the specific case of binary classifiers, which are the simplest type of forecasting and decision-making situation: an extreme event did or did not occur. Risk functions adapted to binary classifiers of extreme events are developed, answering our second question. Their properties are derived under the framework of multivariate regular variation and hidden regular variation, allowing to handle finer types of asymptotic independence. This framework is applied to extreme river discharges.

ticulier à Aussois. Philomène, Pauline, Romain, Gloria, Bastien, Paula, Nicolas, Jonathan (et tous ceux que j'oublie) merci pour toutes ces parties de baby-foot et ces karaoké d'anthologie. Un immense merci aux organisateurs de VALPRED Clément Dombry, Olivier Wintenberger et Philippe, de rendre ces rencontres possibles.

Si cette thèse a pu si bien se dérouler, c'est notamment grâce au soutien financier de l'ENS de Rennes que j'ai reçu, par le biais des Contrats Doctoraux Spécifiques Normaliens. Je remercie également l'IFREMER pour son soutien matériel et pour l'accès aux données d'états de mer.

Je tiens évidemment à remercier mes parents pour leur soutien mais aussi pour leur envie constante d'essayer de comprendre ce que font leurs enfants afin de pouvoir expliquer à d'autres ce qu'ils "cherchent".

Enfin, simplement merci pour tout Camille. Après deux agreg et deux thèses, je pense qu'on a eu notre quota pour les prochaines années.

Résumé en Français

Les zones côtières sont fréquemment soumises à des événements maritimes extrêmes. Ces événements sont, par définition, ceux ayant le plus d'impact sur les populations ou les activités économiques. En France, l'exemple le plus dramatique fût la tempête Xynthia qui s'est abattue en février 2010 sur la façade atlantique, engendrant un nombre élevé de victimes et de dommages matériels [START_REF] Genovese | Storm surge disaster risk management: the Xynthia case study in France[END_REF]. De part sa dimension catastrophique, cet événement est très certainement la tempête maritime qui a le plus marqué les mémoires.

Chaque année, plusieurs tempêtes touchent le littoral français, fort heureusement sans avoir de conséquences humaines aussi dramatiques que la tempête Xynthia, mais qui ont cependant de nombreux impacts que ce soit sur la fiabilité des structures installées en mer ou sur les côtes, ou encore sur l'érosion des littoraux.

Les hauteurs de vagues extrêmes sont une des caractéristiques de ces tempêtes maritimes. Afin de déterminer les potentiels impacts des phénomènes maritimes extrêmes, il est donc crucial de pouvoir caractériser et prédire les vagues les plus extrêmes pouvant survenir près des littoraux. Dans ce contexte, les travaux présentés dans ce manuscrit visent à répondre aux deux questions suivantes :

• Comment simuler de tels événements extrêmes ? (Problématique n°1)

• Comment comparer différents modèles de simulation, ou de prédiction, d'événements extrêmes ? (Problématique n°2)

Une première partie des travaux porte sur la modélisation des extrêmes de la hauteur significative des vagues, notée H s . Cette quantité permet de mesurer l'énergie des vagues et donc leur sévérité. La problématique n°1 est traitée en considérant les conditions d'états de mer au large pouvant générer des hauteurs de vagues extrêmes près des côtes. Une fois le modèle de simulation construit, une question naturelle est d'évaluer ses performances en termes de prédiction d'extrêmes. La seconde partie des travaux de cette thèse vise à répondre à cette problématique n°2, en considérant le cas particulier des classifieurs binaires, qui peuvent être vus comme des encodeurs d'événements extrêmes (un extrême a eu lieu ou n'a pas eu lieu).

Un résumé de ces différents travaux est donné ci-dessous.

Simulation jointe de H s à la côte et au large

Pour répondre à la problématique n°1, nous considérons deux points spécifiques un point au large (représenté en rouge sur la Figure 1) et un point à la côte (en vert sur la Figure 1). Nous souhaitons modéliser conjointement les extrêmes de hauteur significative de vagues en ces deux points. Afin de modéliser simultanément les H s au large et à la côte, l'outil adéquat est la théorie des extrêmes multivariés, qui permet de modéliser la structure de dépendance entre les valeurs les plus extrêmes d'un vecteur aléatoire. Dans le cadre univarié, un extrême peut être défini comme un événement dépassant un certain seuil, fixé au préalable (approche par dépassements de seuil). Dans ce contexte, pour un seuil u suffisamment grand, les dépassements de seuils, définis par [X -u | X > u] où X est une variable aléatoire réelle, peuvent être approximés par une loi de Pareto généralisée, notée GPD [START_REF] Coles | An Introduction to Statistical Modeling of Extreme Values[END_REF]. Cette caractérisation a été ensuite étendue au cadre multivarié. Nous suivons ici la définition développée par [START_REF] Rootzén | Multivariate generalized Pareto distributions[END_REF], en se restreignant au cadre bivarié (correspondant à notre application) : si l'on considère un vecteur X = (X 1 , X 2 ) dans R 2 , alors les dépassements de seuils bivariés

[X -u | X ≰ u] ,
peuvent être approximés par une GPD bivariée, où X ≰ u signifie que X 1 ≥ u 1 ou X 2 ≥ u 2 , si u = (u 1 , u 2 ) ∈ R 2 . Une illustration de ces dépassements de seuils est présentée sur la Figure 2 avec les données de H s au large et à la côte. Afin de construire de tels vecteurs GPD, Rootzén et al. (2018a) ont établi la représentation suivante A partir de la représentation (1), nous développons un générateur stochastique permettant de simuler conjointement des H s extrêmes au large et à la côte, et cela en fonction de certaines conditions au large, qui sont notamment des H s au large modérément élevées.

Z = E + T -max(T ), (1)
Pour cela, un algorithme de simulation non-paramétrique de vecteurs GPD est présenté. Cet algorithme s'appuie sur la ré-écriture suivante de l'Equation (1) :

Z 1 = E + ∆1 ∆<0 , Z 2 = E -∆1 ∆≥0 , (2)
où ∆ := Z 1 -Z 2 = T 1 -T 2 et 1 • correspond à la fonction indicatrice. Le point clef de l'algorithme est l'utilisation d'un bootstrap non-paramétrique pour simuler des valeurs de ∆. Les performances de l'algorithme sont illustrées sur plusieurs simulations numériques.

Avant d'appliquer ce générateur stochastique de vecteurs GPD aux données de H s , l'Equation (1) requiert que les variables soient standardisées à la même échelle. Pour cela, les distributions marginales des H s au large et à la côte sont modélisées par une EGPD [START_REF] Naveau | Modeling jointly low, moderate, and heavy rainfall intensities without a threshold selection[END_REF], qui permet d'éviter la sélection parfois difficile du seuil u de la GPD. L'EGPD autorise la modélisation de toute la plage de données et pas uniquement la queue de distribution. Cela nous permet notamment de modéliser comment des H s au large modérément élevées peuvent produire des H s extrêmes à la côte.

Une spécificité des H s est leur caractère non-stationnaire. En effet, [START_REF] Ewans | Statistical modelling of extreme ocean environments for marine design: A review[END_REF] ont montré que la modélisation des extrêmes de hauteur significative de vagues nécessite de prendre en compte l'effet de certaines covariables, telles que la direction des vagues ou la saisonnalité. Dans notre étude, nous considérons l'effet de la période pic (notée T p ) et de la direction pic (notée D p ) sur la relation entre les H s au large et celles à la côte. Tout comme la hauteur significative de vagues, ces deux variables contribuent à caractériser un état de mer. Elles correspondent à la période et la direction de propagation du maximum -ou "pic" -d'énergie des vagues. La dépendance en T p et D p des modèles sur les lois marginales est portée, dans les deux cas, par le paramètre d'échelle de l'EGPD à l'aide de splines de lissage [START_REF] Carrer | egpd4gamlss[END_REF].

Enfin, un second algorithme est développé, permettant de simuler des H s à la côte conditionnellement aux H s au large. Ce second générateur peut être vu comme une méthode de dowscaling statistique des H s au large (considérées comme des données provenant d'un modèle global de vagues) vers les H s à la côte (issues d'un modèle local).

Les deux algorithmes de simulation sont ensuite validés sur des expériences numériques et appliqués aux extrêmes de H s près des côtes bretonnes françaises.

Evaluation et comparaison de classifieurs d'événements extrêmes

Pour répondre à la problématique n°2, nous étudions le cas particulier des classifieurs binaires. Pour cela, nous considérons que l'on dispose d'une variable binaire Y ∈ {-1, +1} qui encode l'occurrence d'un événement extrême :

Y := Y (u) = +1, si H > u, -1, sinon,
où H ∈ [0, ∞) est une certaine variable aléatoire (par exemple un débit de rivière) et u > 0 un seuil critique élevé.

A partir d'un ensemble d'observations multivariées X ∈ [0, ∞) d , nous cherchons à déterminer, parmi différents classifieurs binaires g(X; u) ∈ {-1, +1} de ces observations, le meilleur en terme de prédictions d'occurrence d'événements extrêmes, i.e. correspondant à Y (u) = +1. Pour fixer les idées, on peut considérer que les observations X représentent des débits de rivière en différentes stations de mesure (points rouges sur la carte 3), et que Y (u) décrit l'occurrence d'un débit de rivière extrême en une station différente de X (triangle jaune sur la carte 3). Afin de comparer différents classifieurs, nous définissons la fonction de risque suivante R (u) (g) = P(g(X; u) ̸ = Y (u) ) P(Y (u) = 1 ou g(X; u) = 1)

∈ [0, 1].
Si l'on regarde la Figure 4 représentant des prédictions binaires obtenues via un arbre de décision, cette fonction de risque revient à compter le nombre moyen de points se trouvant dans les zones jaune et rouge, pondéré par la probabilité d'être dans les zones jaune, verte ou rouge. Un second développement est porté sur le comportement de dépendance asymptotique entre la variable d'intérêt H et les classifieurs binaires. Pour cela, nous considérons des classifieurs de la forme g(X; u) = +1, si g(X) > u, -1, sinon, pour une certaine fonction g : R d → (0, ∞). Dans ce cas-là, nous montrons que la mesure de risque extrême R(g) telle que définie précédemment ne permet pas de discriminer des classifieurs dans le cas où g(X) et H sont asymptotiquement indépendants.

Pour contourner ce problème, nous définissons une seconde fonction de risque, adaptée aux cas de l'indépendance asymptotique, en considérant le conditionnement suivant pour ε ∈ [0, 1) R ε (g) := lim u→∞ P(g(X; u) ̸ = Y (u) | Y (εu) = g(X; εu) = 1) P(Y (u) = 1 ou g(X; u) = 1 | Y (εu) = g(X; εu) = 1) .

Nous établissons alors un estimateur empirique de R ε (g), dont les propriétés asymptotiques sont dérivées dans le cadre de la variation régulière multivariée et de la variation régulière cachée, à l'aide principalement du modèle de Ramos and Ledford (2009), dont une adaptation est la suivante

P[g(X) > u, H > v] = L(u, v)(u -αg v -α H ) 1/2η ,
où η ∈ (0, 1] indique le vitesse de décroissance de la queue bivariée, L(•, •) est une fonction à variation lente et α H , α g > 0 sont les indices de queue respectifs de H et g(X).

Cette seconde fonction de risque est appliquée à l'analyse des débits extrêmes de rivière dans le bassin versant du Danube. L'application compare différents algorithmes de prédiction (classifieur linéaire, régression logistique, arbres de décision,. . . ) et teste leur capacité à prédire les débits de rivière extrêmes en une station donnée à partir des observations en d'autres stations (voir Figure 3).

Introduction Motivations

Coastal areas are frequently subject to extreme maritime events. These events are, by definition, those that have the greatest impact on populations or economic activities. In France, the most dramatic example was the Xynthia storm which hit the Atlantic coast in February 2010, causing a high number of fatalities and material damages [START_REF] Genovese | Storm surge disaster risk management: the Xynthia case study in France[END_REF]. Due to its catastrophic dimension, this event is undoubtedly the maritime storm that has left the greatest mark on people's memories.

Every year, numerous storms hit the French coastline, fortunately without having such dramatic human consequences as the Xynthia storm, but which nevertheless have multiple impacts, whether on the reliability of structures installed at sea or on the coasts, or on coastal erosion.

Extreme wave heights are one of the characteristics of these maritime storms. In order to determine the potential impacts of extreme maritime events, it is therefore crucial to be able to characterise and predict the most extreme waves that can occur near the coast. In this context, the work presented in this manuscript aims to answer the following questions:

• How to simulate such extreme events? (Problematic n°1)

• How to compare different simulation, or prediction, models of extreme events? (Problematic n°2)

The first part of the work concerns the modelling of extreme significant wave heights, denoted H s . This quantity measures the energy of the waves and, consequently, their severity. The problematic n°1 is addressed by considering the offshore sea state conditions that can generate extreme wave heights near the coast.

Once the simulation model is built, a natural question is to evaluate its performance in terms of prediction of extremes. The second part of the work of this thesis aims at answering this problematic n°2, by considering the particular case of binary classifiers, which can be seen as extreme event encoders (an extreme has or has not occurred).

Outline of the thesis

The outline of this manuscript is as follows. In Chapter 1, the main results of extreme value theory (EVT) are discussed. Since this theory entails many concepts and questions, a focus is made on the key elements of univariate and multivariate EVT, that are useful for the rest of this thesis.

The main goal of Chapter 2 is to answer the problematic n°1. For that, we propose and study a stochastic simulator that, given offshore conditions, produces jointly offshore and coastal extreme significant wave heights. We rely on bivariate Peaks over Threshold [START_REF] Rootzén | Multivariate generalized Pareto distributions[END_REF] and develop a non-parametric simulation scheme of bivariate GPD. To take into account non-stationarities, we also adapt the extended generalised Pareto model [START_REF] Naveau | Modeling jointly low, moderate, and heavy rainfall intensities without a threshold selection[END_REF], letting the parameters vary with specific sea state parameters.

Chapter 3 addresses the problematic n°2. We study the specific case of binary classifiers, which are the simplest type of forecasting and decision-making situation: an extreme event did or did not occur. Such classifier tailored for extremes will be called an extremal classifier and risk functions that answer our initial question will be developed. Their properties will be derived under the framework of multivariate regular variation and hidden regular variation, allowing to handle finer types of asymptotic independence.

All the codes used in this manuscript are produced with R (R Core Team, 2022).

Chapter 1

Elements of extreme value theory

Overview of Chapter 1

In this Chapter, the main results of extreme value theory (EVT) are discussed. Since this theory entails many concepts and questions, a focus is made on the key elements that are useful for the rest of this thesis. The proofs are omitted and can be found in many great reference books such as [START_REF] Coles | An Introduction to Statistical Modeling of Extreme Values[END_REF], [START_REF] Beirlant | Statistics of Extremes: Theory and Applications[END_REF] or de [START_REF] De Haan | Extreme Value Theory: An Introduction[END_REF]. The outline of the present chapter is as follows.

In Section 1.1, we give the foundations of univariate EVT. Starting from its key result that the generalised extreme value distributions are the only possible limits that can be obtained for the maximum of a random sample, under a suitable renormalisation, we introduce two modelling approaches for extremes: block maxima and peaks over thresholds. Then, a brief note is made on regular variation which is a powerful tool in EVT.

Section 1.2 is devoted to multivariate EVT. A focus is made on bivariate extremes, this enables us to avoid complex notations. Similarly to the univariate case, the bivariate extreme value distributions are introduced as the limit of component-wise maxima. Characterisations of this family are given with its spectral representation. Bivariate peaks over thresholds models are discussed, following the representation of [START_REF] Rootzén | Multivariate generalized Pareto distributions[END_REF]. Then a discussion is made on the dependence structure of bivariate extreme value distributions and we show that extreme value models fail for asymptotically independent variables. Therefore, the last section addresses this issue by looking at other types of models apart from the class of extreme value distributions.

This chapter was written with the help of the above-mentioned textbooks. Other references that have proved very valuable include [START_REF] Fougères | Multivariate extremes. Extreme values in finance[END_REF], [START_REF] Reiss | Statistical Analysis of Extreme Values: with Applications to Insurance, Finance, Hydrology and Other Fields[END_REF] and the recent review of [START_REF] Davison | Statistics of extremes[END_REF] which contains helpful animations that illustrate some of the following concepts.

In this chapter, we do not address in details the question of statistical inference methods. A comprehensive account can be found in [START_REF] Beirlant | Statistics of Extremes: Theory and Applications[END_REF] or [START_REF] Dey | Extreme Value Modeling and Risk Analysis: Methods and Applications[END_REF] among others.

Univariate extremes

As mentioned in the introduction of this dissertation, univariate extreme value theory (EVT) provides a solid theoretical basis and framework for studying the largest values of a sample. In this context different modelling methods coexist: block maxima, r largest order statistics, peaks over threshold and point process. In the following, we focus on the maxima and the threshold exceedances of random variables, see, e.g. [START_REF] Coles | An Introduction to Statistical Modeling of Extreme Values[END_REF] or [START_REF] Dey | Extreme Value Modeling and Risk Analysis: Methods and Applications[END_REF] for details on the other approaches.

Possible limits of block maxima

The strength of EVT began in 1928 with the [START_REF] Fisher | Limiting forms of the frequency distribution of the largest or smallest member of a sample[END_REF] Theorem (and later with [START_REF] Gnedenko | Sur la distribution limite du terme maximum d'une serie aleatoire[END_REF]). The idea was that if we consider a sequence of independent and identically distributed random variables X 1 , . . . , X n with common cumulative distribution function F , then there are only three possible limiting distributions for the maximum M n := max{X 1 , . . . , X n } after proper renormalisation. These three families are known as GEV (Generalised Extreme Value distributions).

Theorem 1.1 [START_REF] Fisher | Limiting forms of the frequency distribution of the largest or smallest member of a sample[END_REF]; [START_REF] Gnedenko | Sur la distribution limite du terme maximum d'une serie aleatoire[END_REF]). If there exists sequences

(a n > 0) and (b n ) such that lim n→∞ P [(M n -b n )/a n ≤ x] = G(x), (1.1)
where G is a non-degenerate distribution function. Then G belongs to the GEV family i.e. is of the form:

G ξ (x) =          exp -1 + ξ x -µ σ -1/ξ if ξ ̸ = 0, x ∈ {x; 1 + ξ(x -µ)/σ > 0} , exp -exp x -µ σ if ξ = 0, x ∈ R,
for parameters µ ∈ R (location), σ > 0 (scale) and ξ ∈ R (shape).

Univariate extremes

In EVT, an important issue concerns the shape parameter ξ, which governs the behaviour of the upper-tail distribution (see Figure 1.1):

• If ξ > 0, the upper end-point of G ξ is infinite and the distribution is heavytailed, i.e. the tail of the distribution decreases polynomially (Fréchet family).

• If ξ = 0, the upper end-point of G ξ is also infinite but the distribution is rather light-tailed, i.e. the tail decreases exponentially (Gumbel family).

• If ξ < 0, the upper end-point of G ξ is finite, the distribution tail is bounded (Weibull family). Furthermore, the upper end-point x F can be linked to the parameters of the GEV through x F = µ -σ/ξ. Example 1.2. For example, the Uniform distribution is in the Weibull domain of attraction, with ξ = -1. Indeed, let X i ∼ U (0, 1), setting a n = n -1 and b n = 1 we have:

P (M n -1)/n -1 ≤ x = F n (1 + n -1 x) = 1 + n -1 x n ---→ n→∞ exp(x).
Note that the Weibull family and the Fréchet family can be linked through

X ∈ D(G ξ<0 ) ⇔ Y := (x F -X) -1 ∈ D(G ξ>0 ), (1.2)
where G ξ<0 (resp. G ξ>0 ) is the Weibull (resp. Fréchet) extreme value distribution and x F denotes the upper end-point of the distribution function F of X (see [START_REF] Beirlant | Statistics of Extremes: Theory and Applications[END_REF]. This correspondence is one of the reason why the EVT literature often focuses on the heavy-tail case. It is also justified from an application perspective. Indeed, there is an historical tradition between hydrology and EVT, and it is widely agreed that the distributions of hydrologic variables are heavy-tailed (e.g. [START_REF] Katz | Statistics of extremes in hydrology[END_REF].

From a statistical point of view, however, Equation (1.2) is of no use since, in general, no prior information is available on the sign of ξ.

A convenient related concept to the extreme value distributions is max-stability.

A distribution function F is max-stable if, for every n ∈ N,

F n (b n + a n x) = F (x), (1.3)
for a suitable choice of constants b n and a n > 0. And it can be shown that the classes of generalised extreme value distributions and max-stable distributions coincide (see, for example, [START_REF] Beirlant | Statistics of Extremes: Theory and Applications[END_REF]).

A classification of the standard distributions according to their domain of attraction, with their respective normalising constants, can be found for example in [START_REF] Embrechts | Modelling Extremal Events: for Insurance and Finance[END_REF].

In practice, the normalising constants a n and b n are unknown and the maxima are modelled directly using the following approximation for a large enough block

size n P(M n ≤ x) ≈ G [(x -b n )/a n ] = G * (x),
where G * also belongs to the GEV family and differs from G only in the scale and location parameters.

The max-stability property (1.3) implies in particular that if block maxima with fixed block size n are well approximated by a GEV distribution G, then block maxima over a larger block size should also be well approximated by the GEV G, up to a change of scale and location. Indeed, consider m = αn ≥ n, α ∈ N, then

P(M m ≤ x) = F m (x) = [F n (x)] α = [P(M n ≤ x)] α .
And from this simple equation, if F ∈ D(G), using the max-stability property (1.3), one can write

P(M m ≤ x) ≈ {G [(x -b n )/a n ]} α = G [(x -b αn )/a αn ] .
From Theorem 1.1, and the preceding remark, one can derive directly the first approach for modelling extreme events which is the so-called block maxima approach (see Figure 1.4). For that, the data will be divided into blocks of equal sizes (typically blocks of one year or one month), we then compute the maxima for each block and we can then fit a GEV distribution to the set of block maxima. Note that the size of the block implies a bias-variance trade-off: if the block size is too small, the approximation from the limit in Theorem 1.1 will be bad inducing higher bias; if the block size is too big, we will have too few datapoints, leading to higher variance.

Peaks over thresholds approach

One limitation of the block maxima approach is if many extremes occur in the same block, only one event will be kept, which may lead to a significant loss of data. The second approach, named peaks over thresholds, is an alternative where instead of considering maxima per block, we will look at all events, or peaks, above a given -usually high -threshold. This approach relies on the following theorem where, with a similar emphasis, it states that the only possible limiting distributions of the threshold excesses is the generalised Pareto family (GP).

Theorem 1.3 [START_REF] Balkema | Residual life time at great age[END_REF]; [START_REF] Pickands | Statistical inference using extreme order statistics[END_REF]). For every ξ ∈ R, X is in the domain of attraction of a GEV distribution G ξ if, and only if, the distribution function of the excesses X -u, conditionally on X > u, converges as follows

lim u→x F sup 0<x<x F -u |P(X -u ≤ x | X > u) -H ξ,σ(u) (x)| = 0, (1.4)
for some positive function σ, where x F is the upper end-point of F and H ξ,σ(u) is the generalised Pareto distribution function and defined as

H ξ,σ(u) (x) := 1 -(1 + ξx/σ(u)) -1/ξ .
If Theorem 1.3 holds true, then the shape parameter ξ is the same for the GEV distribution G ξ and the GP distribution H ξ,σ(u) , and σ(u) = σ + ξ(u -µ).

The different behaviours of the tail distribution given the shape parameter ξ are the same as for the GEV family (see Figure 1.2), where for ξ = 0, taking the limit ξ → 0 in (1.4), we have

H 0,σ(u) (x) = 1 -exp(-x/σ(u)).
Similarly to the notion of max-stability, regarding peak over thresholds a useful property is the threshold stability (e.g. [START_REF] Falk | Laws of small numbers: extremes and rare events[END_REF]

. A distribution F is threshold stable (sometimes called POT-stable) if, for every u ∈ [0, x F ], with x F the upper-end point of F , F [u] (b u + a u x) = F (x), (1.5)
for a suitable choice of constants b u and a u > 0. F [u] denotes the distribution of exceedances above the threshold u and is defined by Note that if X has distribution function F , then the threshold exceedances X -u | X > u have distribution function F [u] :

F [u] (x) = F (x) -F (u) 1 -F (u) , x ≥ u.
P(X -u ≤ x | X > u) = F [u] (x + u).
The following proposition shows that the class of threshold stable distributions coincides with the class of GP distributions (see, e.g. [START_REF] Falk | Laws of small numbers: extremes and rare events[END_REF]).

Proposition 1.4. Let F be a distribution function with upper-end point x F .

1. If there exists functions a(u) > 0 and b(u) such that

lim u→x F F [u] (b(u) + a(u)x) = L(x), (1.6)
where L is a non-degenerate distribution function, then L is threshold stable.

Conversely, if

L is a threshold stable distribution, then L satisfies (1.6).

If, in addition to being threshold stable, L is continuous, then L is a GP distribution, up to a change of scale and location.

Example 1.5. Considering the Uniform distribution as in Example 1.2, its conditional threshold exceedances correspond to the generalised Pareto distribution with ξ = -1 and σ(u) = 1 -u

P(X -u ≤ x | X > u) = F (x + u) -F (u) 1 -F (u) = (x + u) -u 1 -u = x 1 -u .

Univariate extremes

From Theorem 1.3 we can derive the second approach to model extreme events, the so-called peaks over thresholds model (see Fig. 1.4). Instead of splitting data into blocks, we will choose a high enough threshold u and then fit a GP distribution to the thresholds exceedances. As for the block maxima approach with the choice of the block length, the choice of the threshold u implies a bias-variance trade-off. For many applications, u will typically be chosen as a high quantile of the data of interest.

Similarly to the max-stability for the block maxima approach, the threshold stability property (1.5) implies in particular that if threshold exceedances over a given threshold u are well approximated by a GP distribution H, then the threshold exceedances over a higher threshold will also be well approximated by the GP distribution H, up to a change of scale and location. To see this, consider a higher threshold v = α + u > u, α > 0, then

P(X -v > x | X > v) = P(X > x + v) P(X > u) × P(X > u) P(X > v) = P(X -u > x + α | X > u) P(X -u > α | X > u)
.

Then, if the exceedances X -u | X > u are approximated by a GP distribution H, using the threshold stability of H, one can write An example of peaks over threshold modelling is given in Figure 1.3, with the wave data used in Chapter 2. As mentioned in the Introduction, in this chapter, we will consider significant wave height, denoted H s . This quantity basically give information about the ocean surface roughness, more details are provided in Section 2.1.1. Figure 1.3 depicts the distribution of the threshold exceedances of H s given that H s > u, where u is the 0.98 quantile of H s . Two GP densities are superimposed, showing that the H s exceedances should be rather well approximated by a light-tailed distribution. This a priori analysis will be confirmed in Chapter 2, where the estimated shape parameters for H s data in two locations will be slightly negative (see Table 2.2).

P(X -v ≤ x | X > v) ≈ H(x + α) -H(α) 1 -H(α) = H((x + α -b u )/a u ).
To avoid a careful threshold selection, [START_REF] Naveau | Modeling jointly low, moderate, and heavy rainfall intensities without a threshold selection[END_REF] introduced an extension of the GP distribution, the so-called extended generalised Pareto distribution (EGPD). In a nutshell, the EGPD allows to model the entire range of a distribution with a smooth transition between the lower and the upper tails. This model is introduced in Section 2.1.3 and it has been applied to the wave data in [START_REF] Legrand | Joint stochastic simulation of extreme coastal and offshore significant wave heights[END_REF], see Chapter 2.

Figure 1.4: Illustration of the block maxima versus the peaks over thresholds approaches with one year of the wave data described in Chapter 2, Section 2.1.1.

Univariate extremes

To sum up

To model extreme events, two main strategies can be used:

1. The block maxima approach which relies on the approximation of renormalised maxima by a GEV distribution (Theorem 1.1).

⇝ Needs to choose the size of the blocks.

2. The peaks over thresholds approach which relies on the approximation of the thresholds exceedances by a GP distribution (Theorem 1.3).

⇝ Needs to choose the threshold.

Regular variation and domain of attraction

When studying EVT, one theoretical tool that quickly comes into play is regular variation. Here we briefly introduce the key concepts and show how regular variation can be used to determine domains of attraction. See, e.g. [START_REF] Bingham | Regular Variation. Encyclopedia of Mathematics and its Applications[END_REF] for more details.

Definition 1.6. A positive measurable function f is regularly varying (r.v.) with index α ∈ R if

lim t→∞ f (tx) f (t) = x α , for all x > 0.
If α = 0, f is called a slowly varying function.

Example 1.7. Examples of regularly varying functions are: all the functions converging to a positive constant, x → x α for α ∈ R, x → log β (x) for β ∈ R.

In the following, we consider regularly varying survival functions 1 -F of a random variable X. For this, an additional constraint is required by imposing α > 0. This can be seen through the following convergence, assuming that 1 -F is r.v. with index α,

lim t→∞ 1 -F (tx) 1 -F (t) = lim t→∞ P(X > tx) P(X > t) = x -α ,
(1.7) then necessarily, α must be strictly positive.

Thanks to the following result, due to [START_REF] Gnedenko | Sur la distribution limite du terme maximum d'une serie aleatoire[END_REF], the domains of attraction of the Fréchet family and the Weibull family can then be characterised in terms of regular variation.

Theorem 1.8 [START_REF] Gnedenko | Sur la distribution limite du terme maximum d'une serie aleatoire[END_REF]). In the following, G ξ>0 (resp. G ξ<0 ) denotes the Fréchet (resp. Weibull) family distribution and x F denotes the upper end-point of F .

• F ∈ D(G ξ>0 ) iif 1 -F is regularly varying (r.v.) with index 1/ξ. • F ∈ D(G ξ<0 ) iif x F < ∞ and 1 -F (x F -1/x) is r.v. with index -1/ξ.
Note that using the link between the Weibull and the Fréchet families provided in Equation (1.2), only the first statement of Theorem 1.8 is needed.

Example 1.9. Every Pareto-like distributions, i.e. such that 1 -F (x) ∼ Cx -α for large x and for some C, α > 0, are in the domain of attraction of the Fréchet family. For example the Pareto, Cauchy, and Burr distributions.

From this characterisation theorem, we see that the theory of regular variation is an appropriate framework to study heavy-tailed distributions, bringing many valuable theoretical tools (e.g. [START_REF] Resnick | Heavy-tail phenomena: probabilistic and statistical modeling[END_REF]. In Chapter 3, we will provide the multivariate extension of regular variation and this theoretical framework will be applied in particular to extreme river discharges, which are heavy-tailed data.

Multivariate extremes

The extension from univariate extremes to multivariate extremes is not straightforward, partly because it implies a point of view choice. For example, considering temperature data in different cities, one could say that an extreme (multivariate) event occurs if large values of temperature in one of the cities is observed, but it could also be that large values in all cities occur simultaneously. In the bivariate case, an illustration of these two approaches is given in Figure 1.5 with daily river discharges from the Danube river network (more details on this dataset are given in Chapter 3). The underlying question here could then be the following: Do we want to assess the risk of flooding at one (or more) location or are we more interested in assessing the overall risk of flooding -which increases if extreme river discharges occur at different locations simultaneously?

The second key aspect of multivariate extremes is the dependence structure between the different components in the asymptotic. For large values, the dependence structure between the components of a multivariate distribution can be disrupted. An example developed by [START_REF] Sibuya | Bivariate extreme statistics[END_REF] is the bivariate Gaussian distribution, for which, whenever the correlation coefficient ρ is strictly smaller than 1, its components are independent in the extremes (see Figure 1.6).

To focus on the dependence structure , the standard method is to transform the data to common margins which removes marginal effects. This is what is done in Figure 1.6, where the data are transformed to common Pareto margins (for example). Note that this transformation does not alter the dependence structure. Once the marginal effects are removed, it is then easier to seek for dependence or independence. This two step study can be elegantly illustrated through the copula theory and the following theorem: Figure 1.5: Scatter plots of daily river discharges of two stations in the Danube river network. The two possible regions of extremes are depicted in blue: when an extreme event is defined as a large value observed in at least one of the components (left) and when it is defined as a large value in both components simultaneously (right).

Theorem 1.10 [START_REF] Sklar | Fonctions de répartition à n dimensions et leurs marges[END_REF]). Let F be a joint distribution function with margins F X and F Y . Then there exists a copula C1 such that, for all x, y ∈ R,

F (x, y) = C {F X (x), F Y (y)} .
(1.8) If F X , F Y are continuous, then C is unique. Conversely, if C is a copula and F X , F Y are distribution functions, then the function F defined by (1.8) is a joint distribution function with margins F X and F Y .

From Theorem 1.10 we see that the study of the joint distribution F can be decomposed into the study of its margins and of the dependence structure, via the copula C. More details on copula theory, and its links with EVT, can be found for example in [START_REF] Nelsen | An Introduction to Copulas[END_REF].

On this basis and derived from theoretical and practical perspectives, the modelling of multivariate extreme events is in general carried out in two steps: the first step concerns the treatment of the margins and the second step the description of the dependence behaviour (e.g. [START_REF] Cai | Environmental data: multivariate Extreme Value Theory in practice[END_REF][START_REF] Tendijck | Modeling the extremes of bivariate mixture distributions with application to oceanographic data[END_REF][START_REF] Rohmer | Partitioning the uncertainty contributions of dependent offshore forcing conditions in the probabilistic assessment of future coastal flooding at a macrotidal site[END_REF].

The marginal modelling can involve classic univariate EVT (see Section 1.1). Given the application, some models will be more adapted than others. For example, in Chapter 2 the margins are modelled using a specific type of extreme value distribution: the extended generalised Pareto distribution [START_REF] Naveau | Modeling jointly low, moderate, and heavy rainfall intensities without a threshold selection[END_REF]. If there is no specific constraint on the marginal modelling, to facilitate the choice, [START_REF] Beirlant | Statistics of Extremes: Theory and Applications[END_REF] suggest that when the main interest lies in the dependence structure, a convenient choice is to transform into standard Fréchet margins as there are good theoretical properties in this case.

Regardless of the chosen standardisation distribution, once the data are transformed to common margins, different theoretical objects can be defined to model and to measure the dependence at high levels. In the remaining of this section we will briefly introduce these tools and provide the links between them. We will restrict our attention to the bivariate case, generalisation to higher dimensions can be found in the references mentioned earlier.

To lighten the notations, in the following bold face symbols will denote bivariate vectors: x = (x 1 , x 2 ) ∈ R 2 , 0 = (0, 0),. . . Operations and relations involving such vectors will be considered component-wise. For example min(x, 0) will correspond to (min(x 1 , 0), min(x 2 , 0)) ∈ R 2 .

Domain of attraction and characterisation

The definition commonly adopted in the literature for the notion of maximum in the multivariate case is the component-wise maxima. If we assume that we have a set of bivariate vectors (X 1 , Y 1 ), . . . , (X n , Y n ) that are independent and identically distributed with common bivariate distribution function F , then the sample maximum M n is defined by

M n := max 1≤i≤n X i , max 1≤i≤n Y i .

Multivariate extremes

Note that the sample maximum does not necessarily correspond to a sample point. Then, similarly to the univariate case, we can define the family of multivariate extreme value distributions as the only possible limit for the sample maxima after proper renormalisation. Recall that operations involving vectors are meant component-wise.

Definition 1.11. If there exists sequences in R 2 (a n ) and (b n ) with a n,j > 0 for j = 1, 2, such that

lim n→∞ P [(M n -b n )/a n ≤ (x, y)] = G(x, y), (1.9)
where G is a non-degenerate distribution function, then G belongs to the bivariate extreme value family and F is said to be in the (bivariate) domain of attraction of G, denoted F ∈ D(G).

This definition implies that, if X = (X, Y ) ∈ D(G), the marginal distribution F X (resp. F Y ) of X (resp. Y ) must be in the univariate domain of attraction of the GEV distribution G 1 (resp. G 2 ), where G j denotes the jth marginal distribution of G, for j = 1, 2.

As in the univariate case, there is a multivariate version of max-stability:

G is max-stable if G n (b n + a n x) = G(x),
for suitable vectors b n and a n > 0. And it can be proved that G is a bivariate extreme value distribution if, and only if, G is max-stable.

Example 1.12. Consider a bivariate vector with unit Fréchet margins. Then in this case, taking b n = 0 and a n = (n, n) for all n > 0,

P [M n /a n ≤ (x, y)] = F n (n(x, y)) = (exp(-1/(n(x, y)) n = exp(-1/(x, y)).
The standard Fréchet distribution is then in its own domain of attraction.

We now assume that both margins are transformed to standard Fréchet distributions (whose distribution function is given by x → exp(-1/x) for x > 0). Then the following theorem provides a characterisation of the bivariate extreme value family defined in Equation (1.9) (see [START_REF] Fougères | Multivariate extremes. Extreme values in finance[END_REF] or more specifically [START_REF] Resnick | Extreme Values, Regular Variation and Point Processes[END_REF] for the multivariate version).

Theorem 1.13. The following statements are equivalent:

(i) G is a bivariate extreme value distribution with unit Fréchet margins.

(ii) There exists a finite measure H on [0, 1] such that for each x, y > 0, one has that

G(x, y) = exp -2 1 0 max w x , 1 -w y dH(w)
with the mean constraints 1 0 dH(w) = 1 and 1 0 wdH(w) = 1/2.

CHAPTER 1. ELEMENTS OF EXTREME VALUE THEORY

If we are in the framework of Theorem 1.13, the following measure is also described

V (x, y) := 2 1 0 max w x ,
1 -w y dH(w), for any x, y > 0.

(1.10)

The measure V can be any positive function such that for t > 0 and x, y > 0, V (tx, ty) = t -1 V (x, y) (homogeneity property) and V (z, ∞) = V (∞, z) = 1/z for any z > 0 (marginal constraints).

Definition 1.14. The measure H is known as the spectral measure and V is sometimes called the exponent measure of the extreme value distribution function G.

Instead of V defined in (1.10), the stable tail dependence function (e.g. [START_REF] Beirlant | Statistics of Extremes: Theory and Applications[END_REF] is widely considered in the literature and is defined by

l(x, y) := V 1 x , 1 y , x, y ≥ 0.
Other alternative and equivalent characterisations can be found. For example, the Pickands dependence function [START_REF] Pickands | Multivariate extreme value distributions[END_REF] is defined as

A(t) := V 1 1 -t , 1 t , for any t ∈ [0, 1].
(1.11)

The Pickands dependence function can be linked to the spectral measure through

A(t) = 2 1 0 max {w(1 -t), (1 -w)t} dH(w).
From the marginal constraints on V , A(0) = A(1) = 1. Moreover, A lies within the set max(t, 1 -t) ≤ A(t) ≤ 1 for any t ∈ [0, 1] and is convex within this region. The measure V is completely determined by the Pickands dependence function A through

V (x, y) = 1 x + 1 y A x x + y , for x, y > 0 such that x + y > 0.
The spectral measure H can also be recovered from V , or equivalently A, see [START_REF] Beirlant | Statistics of Extremes: Theory and Applications[END_REF] for more details, so that there are one-to-one mappings between H, V and A.

Example 1.15. Historically, the first parametric family of bivariate extremes is the logistic model [START_REF] Beirlant | Statistics of Extremes: Theory and Applications[END_REF]. Its distribution G is given by G(x, y) = exp -x -1/α + y -1/α α , x, y > 0,

(1.12) 1.2. Multivariate extremes for a parameter 0 < α ≤ 1. The corresponding exponent measure is therefore

V (x, y) = x -1/α + y -1/α α ,
and it can be shown that the corresponding spectral measure H has associated spectral density h given by

h(w) = 1 -α α {w(1 -w)} 1/α-2 (1 -w) 1/α + w 1/α α-2 , w ∈ (0, 1).
The behaviour of this family given α will be discussed later in Section 1.2.3. Other parametric families can be found in Section 9.2.2 of [START_REF] Beirlant | Statistics of Extremes: Theory and Applications[END_REF]. Simulations from this model are displayed in Figure 1.7 for different values of the parameter α. Alternative ways of describing multivariate extremes to the ones presented here exist, see for instance [START_REF] Segers | Max-stable models for multivariate extremes, revstat-stat[END_REF] for further details.

To sum up

Assuming that the bivariate extreme value distribution G has unit Fréchet margins, it can be equivalently characterised by:

1. the spectral measure H G(x, y) = exp -2 1 0 max w x , 1 -w y dH(w) ,
2. the exponent measure V G(x, y) = exp(-V (x, y)),

the Pickands dependence function

A G(x, y) = exp - 1 x + 1 y A x x + y .

Bivariate threshold excess model

As for the univariate case, instead of considering maxima one can consider threshold exceedances. However, the extension to the multivariate case is not straightforward since it requires to define the conditioning event X > u in a multivariate way. For that, there is no universal choice on how to define multivariate exceedances. In this work, we follow the definition of [START_REF] Rootzén | Multivariate generalized Pareto distributions[END_REF], saying that an observation is extreme if it is extreme in at least one of its components, but other definitions exist (e.g. [START_REF] Falk | Peaks-over-threshold stability of multivariate generalized pareto distributions[END_REF].

Such event is denoted X ≰ u and means that at least one component of X is above a given (high) threshold, i.e.

{X ≰ u} := {X 1 > u 1 or X 2 > u 2 }, (1.13)
where u = (u 1 , u 2 ) and X = (X 1 , X 2 ). The subset defined in (1. In Chapter 2, a joint simulation model of significant wave heights (H s ) in two locations is developed. An application of such stochastic simulator could be on the 1.2. Multivariate extremes installation of wind farms in two relatively distant sites and the study of the extreme wave conditions that these sites could face simultaneously. In this context, the definition of [START_REF] Rootzén | Multivariate generalized Pareto distributions[END_REF] appears to be the most appropriate since it allows for the possibility that the H s are not necessarily extreme on both sites simultaneously (see Figure 1.8). [START_REF] Rootzén | Multivariate generalized Pareto distributions[END_REF] provide the following Definition 1.16 for the multivariate generalised Pareto distribution (stated here only in the bivariate case), which fulfils the two following motivating conditions for the extension from univariate to multivariate GP (Rootzén and Tajvidi, 2006, p.919

H(x) := 1 -log G(0) log G(x) G(x ∧ 0) , (1.14)
for some bivariate extreme value distribution G with non-degenerate margins and with 0 < G(0) < 1. In particular, H(x) = 0 for x < 0 and H

(x) = 1 -log G(x)/ log G(0) for x > 0.
Considering the spectral representation of Theorem 1.13, when x > 0, H(x)

can easily be linked to the spectral measure of the associated GEV distribution G, or similarly, to the Pickands dependence function A (see Propositions 5.3 and 6.1 of [START_REF] Tajvidi | Characterisation and some statistical aspects of univariate and multivariate generalised Pareto distributions[END_REF]). In the following, we will provide the explicit link between the exponent measure V and the GP distribution in a specific case. [START_REF] Rootzén | Multivariate generalized Pareto distributions[END_REF] showed that the properties (M 1) and (M 2) are fulfilled with this definition of bivariate GP distributions thanks to the two following theorems. Theorem 1.17 answers the motivation (M1), and can be viewed as a multivariate extension of Theorem 1.3. The second theorem, Theorem 1.18, is the analogue of the threshold stability property in the univariate case (see Proposition 1.4), and answers the above motivation (M2).

Hereinafter we consider a threshold function defined as an increasing continuous curve u(t) parameterised by t ∈ [1, ∞), i.e. the thresholds are defined from the set {u(t) = (u

1 (t), u 2 (t)) : t ∈ [1, ∞)}. Here "increasing" means that if t > s, with t, s ∈ [1, ∞), u(t) > u(s).
Recall that this inequality is meant componentwise, i.e. it corresponds to u 1 (t) > u 1 (s) and u 2 (t) > u 2 (s). An illustration is given in Figure 1.9. This parametrisation is needed to specify how the 2-dimensional levels increase as we look further into the tail of the bivariate distribution F of X (see Tajvidi, 1996, Chapter 4).

Figure 1.9: Illustration of the curve u(t) = (u 1 (t), u 2 (t)).

Theorem 1.17 [START_REF] Rootzén | Multivariate generalized Pareto distributions[END_REF]). Let consider X = (X 1 , X 2 ) a bivariate random vector with distribution function F .

1. If F is in the domain of attraction of a bivariate extreme value distribution G with 0 < G(0) < 1, then there exists an increasing continuous curve u(t), starting at u(1) = 0, with lim t→∞ F (u(t)) = 1 and a function σ(u(t)) > 0, such that

lim t→∞ P σ(u(t)) -1 (X -u(t)) ≤ x | X ≰ u(t) = 1 -log G(0) log G(x) G(x ∧ 0)
, for all x.

Conversely, if there exists an increasing continuous curve u with lim

t→∞ F (u(t)) = 1
and a function σ(u(t)) > 0 such that, for x > 0,

lim t→∞ P σ(u(t)) -1 (X -u(t)) ≤ x | X ≰ u(t) = H(x),
(1.15)

for some function H, where the marginals of H on R + are non-degenerate. Then the left-hand side of (1.15) converges to a limit H(x) for all x and there is a unique bivariate extreme value distribution G with G(0) = e -1 such that

H(x) = log G(x) G(x ∧ 0) , (1.16) G(x) = exp {H(x) -1} for x > 0, and F ∈ D(G).
Note that the two representations (1.14) and (1.16) of the bivariate GP distribu-

tions coincide: letting t = 1/(-log G(0)) in (1.14), H(x) = log G(x) t G(x ∧ 0) t
and by max-stability of G, G t is again an extreme value distribution with G(0) t = exp(-1).

Theorem 1.18 [START_REF] Rootzén | Multivariate generalized Pareto distributions[END_REF]). Again, let consider X = (X 1 , X 2 ) a bivariate random vector.

1. If X has a bivariate GP distribution H, then there exists an increasing continuous curve u(t) with u(1) = 0 and lim t→∞ H(u(t)) = 1, and a function σ(u(t)) > 0 such that

P σ(u(t)) -1 (X -u(t)) ≤ x | X ≰ u(t) = H(x),
(1.17)

for t ∈ [1, ∞) and all x.
2. Conversely, if there exists an increasing continuous curve u(t) with u(1) = 0 and lim t→∞ P [X < u(t)] = 1, and a function σ(u(t)) > 0 such that (1.17) holds for x > 0, and X has non-degenerate margins, then X has a bivariate GP distribution.

The marginals of a bivariate GP distribution are typically not univariate GP. However, their restriction to the positive subset are:

P(X j ≤ x | X j > 0) = 1 -(1 + ξ j x/σ j ) -1/ξ j ,
for x ≥ 0 such that σ j + ξ j x > 0 and where ξ j and σ j , for j = 1, 2, are the marginal shape and scale parameters. For any bivariate GP distribution H, we therefore consider its associated shape parameter ξ as the vector of the marginal parameters (ξ 1 , ξ 2 ), and, similarly, its scale parameter σ := (σ 1 , σ 2 ). We write GP (σ, ξ) the GP distribution H with parameters σ and ξ, and the standard GP distribution is defined for ξ = 0 and σ = 1.

Any bivariate vector X = (X 1 , X 2 ) ∼ GP (σ, ξ) can be characterised through the standard GP distribution, as follows

X = σ e ξZ -1 ξ , (1.18)
where Z ∼ GP (1, 0). From Equation (1.18), the study of bivariate GP distributions can then be reduced to the study of standard form GP distributions. Rootzén et al. (2018b,a) introduce a number of stochastic representations of such standardised GP random vectors. One of them is as follows.

Let T be a bivariate random vector and let E be a unit exponential random variable independent of T . Then

E + T -max(T ) ∼ GP (1, 0). (1.19)
Conversely, any standard GP vector can be expressed in this way. To emphasise the dependence on the bivariate vector T , it may be better to write GP (1, 0, T ) instead of GP (1, 0).

Using Definition 1.16, we can consider the exponent measure V , as defined in (1.10), of the GEV G associated to the GP distribution GP (1, 0, T ). Then Rootzén et al. (2018a) showed that V and T = (T 1 , T 2 ) can be linked through

V (x 1 , x 2 ) = E max 1 x 1 e T 1 -max(T ) E(e T 1 -max(T ) ) , 1 x 2 e T 2 -max(T )
E(e T 2 -max(T ) )

.

From (1.19), it appears that to construct bivariate GP models, a possibility is to assume a parametric distribution on T . This question was addressed by [START_REF] Kiriliouk | Peaks over thresholds modeling with multivariate generalized Pareto distributions[END_REF] where different explicit density formulas for T are derived. Two examples of such models are given in Figure 1.10 with T that follows a bivariate Gaussian distribution and, on the other hand, T 1 and T 2 being two independent Gumbel-distributed variables. In Chapter 2, we will built on model (1.19) to develop a non-parametric generator of such standard GP vectors, without assuming any specific parametric shape on T .

Finally, from a practical modelling perspective, a similar approach to the univariate case can be used as motivated by Theorem 1.17. Indeed, for a sufficiently large threshold u (i.e. large in each component), Xu | X ≰ u can be well approximated by a bivariate GP distribution. As this statistical modelling approach is still at an early stage, few papers deal with practical applications of multivariate GP models (e.g. [START_REF] Brodin | Univariate and bivariate GPD methods for predicting extreme wind storm losses[END_REF][START_REF] Kiriliouk | Peaks over thresholds modeling with multivariate generalized Pareto distributions[END_REF][START_REF] Kiriliouk | Climate extreme event attribution using multivariate peaks-over-thresholds modeling and counterfactual theory[END_REF] following the representation of [START_REF] Rootzén | Multivariate generalized Pareto distributions[END_REF].

To sum up

Bivariate GP distributions emerge as the limit distributions of bivariate threshold excesses, defined as Xu | X ≰ u.

One construction of standard forms for such bivariate GP vectors is given by

E + T -max(T ),
where E ∼ Exp(1) and T is any bivariate vector independent of E.

Asymptotic dependence and independence: Definitions and characterisations.

As already mentioned, an important feature in multivariate EVT is to describe the dependence behaviour of two variables as they become larger and larger. This question has many applications, among which we mention the following two.

• Among the components of a multivariate random vector, identifying subgroups such that large values occur simultaneously allows to reduce the dimension of the model by removing all the other components that do not contribute to the joint extremal behaviour. This is related to the notion of sparsity for extreme value, a very active research field (e.g. [START_REF] Lehtomaa | Asymptotic independence and support detection techniques for heavy-tailed multivariate data[END_REF][START_REF] Engelke | Sparse structures for multivariate extremes[END_REF][START_REF] Meyer | Sparse regular variation[END_REF]. This issue will also be addressed in a more humble way in Section 3.7.2.

• Another typical application of multivariate EVT is the probability estimation of failure regions, i.e. P [(X, Y ) ∈ D], where D is a remote set (with no observations). This can simply corresponds to the subsets {X > x, Y > y} or {X > x or Y > y} for large x, y (see Figure 1.5), or it can also corresponds to more sophisticated sets given some specific constraints such as environmental loads (e.g. de Haan and de Ronde, 1998; [START_REF] Beersma | Joint probability of precipitation and discharge deficits in the Netherlands[END_REF][START_REF] Cai | Environmental data: multivariate Extreme Value Theory in practice[END_REF][START_REF] Ewans | Evaluating environmental joint extremes for the offshore industry using the conditional extremes model[END_REF]. In this context, wrongly assuming the dependence structure between the large values of X and Y will lead to an overestimation or an underestimation of the probability of failure region (see the aforementioned references).

Given the data at hand, the dependence structure in the extremes can be very different. For example, looking at Figure 1.11, which compares the data used in Chapter 2 with those used in Chapter 3, the dependence behaviour between the largest values in the extreme region, symbolised by the upper blue hatched quadrant, appears to be quite different. Regarding the wave data, large values mainly occur simultaneously in both variables, whereas in the case of the river data, a large value in one component does not always corresponds to a large value in the other component. The limiting extreme value distribution G provides this type of information, i.e. the dependence structure of the variables X and Y . In the following, we keep the previous notations, i.e. the bivariate distribution function F of a vector (X, Y ) is in the domain of attraction of a bivariate extreme value distribution G characterised by a spectral measure H, an exponent measure V or, equivalently, a Pickands dependence function A.

The random variables X and Y are said to be asymptotically independent if, equivalently,

• H puts mass 1/2 on the boundaries w = 0 and w = 1,

• V (x, y) = 1/x + 1/y for any x, y > 0, • A(t) = 1 for any t ∈ [0, 1].
These properties correspond to be in the domain of attraction of the independence, i.e. the limiting extreme value distribution decomposes into

G(x, y) = G 1 (x) × G 2 (y),
where G 1 , G 2 are the marginal distributions of G.

Conversely, X and Y are said to be completely dependent if, equivalently,

• H puts mass 1 at the barycenter w = 1/2,

• V (x, y) = max(1/x, 1/y) for any x, y > 0,

• A(t) = max(t, 1 -t) for any t ∈ [0, 1].
Example 1.19. Continuing with Example 1.15, Figure 1.12 depicts the spectral measure H and the Pickands dependence function A for the bivariate logistic model with different values of α. Figure 1.12 shows that α measures the strength of dependence between the two components. In particular, when α → 0 in (1.12),

G(x, y) → exp -max x -1 , y -1 ,
which corresponds to complete dependence, and when α → 1,

G(x, y) → exp -x -1 + y -1 ,
and this corresponds to asymptotic independence. A for the logistic model defined in (1.12) [START_REF] Coles | An Introduction to Statistical Modeling of Extreme Values[END_REF], with α = 0.2 (red solid line), α = 0.5 (blue dashed line) and α = 0.99 (green dotted line). The spectral measures on the left panel have been normalised between 0 and 1 for visual purposes. On the right panel, the black dotted line represents the lower bound of A. Estimations are performed using the R packages evd [START_REF] Stephenson | evd: Extreme value distributions[END_REF] for H and copula [START_REF] Hofert | copula: Multivariate Dependence with Copulas[END_REF] for A.

In practice, asymptotic independence means that the largest values between the components of a random vector are unlikely to occur simultaneously. If this is not the case, i.e. large values occur simultaneously with a high probability, we say that X and Y are asymptotically dependent. The goal is then to determine the strength of this dependence.

From Theorem 1.13 and the spectral representation, all the dependence structure of the limiting distribution is contained in the spectral measure H. From this, it can be very useful to have a summary of this information. Different coefficients have been proposed to measure the strength of asymptotic dependence between two variables. The best-known is probably the dependence measure χ [START_REF] Coles | Dependence measures for extreme value analyses[END_REF], defined as

χ := lim x→x F P(X > x | Y > x)
(1.20)

A generalisation in the case where the marginal distributions of X and Y are non-identical is given by

χ := lim u→1 P(F X (X) > u | F Y (Y ) > u).
The coefficient χ defined in (1.20) can be written as the limit of a function

χ(u) such that χ = lim u→1 χ(u) with χ(u) := 2 - log {P(F X (X) < u, F Y (Y ) < u)} log(u) .
The different properties for χ and χ(u) are the following:

• 0 ≤ χ ≤ 1.

• χ = 0 if, and only if, X and Y are asymptotically independent. This characterisation corresponds to the definition given by [START_REF] Sibuya | Bivariate extreme statistics[END_REF] for asymptotically independent variables.

• If χ > 0, X and Y are asymptotically dependent, and the value of χ increases with the strength of asymptotic dependence.

• The sign of χ(u) determines whether X and Y are positively or negatively associated at the quantile level u.

Example 1.20. Looking at the wave data used in Chapter 2, we can consider the strength of asymptotic dependence between the H s at the two locations of interest.

From Figure 1.13, it seems that χ = lim u→1 χ(u) > 0, indicating that the two variables might be asymptotically dependent.

Remark 1.21. Looking at the right-hand side of Figure 1.13, as u gets closer to 1, the point-wise confidence interval increases, up to contain 0, due to the smaller number of points in the extreme. It is therefore not so evident to be in favour of asymptotic dependence.

In this context, one could use other graphical tools to check for extremal dependence. For example, and as already mentioned, it is possible to standardised data in common margins in order to focus on the extremal dependence structure. This is shown in Figure 1.14, where we compare the asymptotic dependence behaviour of the H s data and the simulated Gaussian vector of Figure 1.6. Unlike the Gaussian vector, for the H s data some points are lying in the interior of the positive quadrant, supporting asymptotic dependence.

Note that more precise diagnostics can be achieved with the help of statistical tests.

There is an extensive literature on this subject, a review can be found in [START_REF] Bacro | Measuring and modelling multivariate and spatial dependence of extremes[END_REF]. For example [START_REF] Draisma | Bivariate tail estimation: dependence in asymptotic independence[END_REF] proposed a test based on an estimator of the coefficient of tail dependence η (see Equation (1.21)). With an extension to the spatial context, [START_REF] Bacro | Testing the independence of maxima: from bivariate vectors to spatial extreme fields[END_REF] developed a test linked to the F-madogram, a quantity often considered in spatial extremes and defined by 1

2 E|F X (X) -F Y (Y )|
(see the references therein).

If we assume that X and Y have common standard Fréchet margins, χ can also be written as χ = lim x→+∞ χ(x) where

χ(x) := 2 + x log {P(X < x, Y < x)} .
From this, one can easily see that for the bivariate extreme value distribution given by G(x, y) = exp(-V (x, y)), χ(x) ≡ χ = 2 -V (1, 1). More generally, if (X, Y ) is in the domain of attraction of a bivariate extreme value distribution G with exponent measure V , then χ = 2 -V (1, 1). In terms of the Pickands dependence function this can be written as χ = 2(1 -A(1/2)). From this, we retrieve that asymptotic independence is equivalent to χ = 0.

Example 1.22. Following Examples 1.15 and 1.19, the dependence measure for the logistic model is given by

χ = 2 -2 α .
The role of α is once again highlighted: when α gets closer to 0, χ tends to 1 (i.e. stronger asymptotic dependence) and when α = 1, χ = 0 (i.e. asymptotic independence).

Yet, in the case of asymptotic independence (i.e. χ = 0), no information is available on the strength of dependence at large but finite levels. In this case, an alternative measure of dependence, denoted χ, is often introduced (e.g. [START_REF] Coles | Dependence measures for extreme value analyses[END_REF] and provides information on the strength of dependence in case of asymptotic independence. This coefficient is defined through the following limit, provided that it exists, χ := lim u→1 χ(u), where

χ(u) := 2 log(1 -u) log P(F X (X) > u, F Y (Y ) > u) -1, u ∈ (0, 1).
Its properties are the following:

• -1 ≤ χ ≤ 1.

• χ = 1 if, and only if, X and Y are asymptotically dependent, and in this case χ informs us on the strength of asymptotic dependence between the variables.

• If χ < 1, X and Y are asymptotically independent, and the value of χ provides a measure that increases with dependence strength.

All in all, the pair (χ, χ) provides complete information on the extreme dependence behaviour between two variables. This has been applied for example in [START_REF] Towe | The extremal dependence of storm severity, wind speed and surface level pressure in the northern North Sea[END_REF] to study the dependence structure between extreme significant wave height and extreme wind speed (such data will be considered in Chapter 2).

The disadvantage of χ is that, unlike χ, it is not easily interpreted. A third dependence coefficient was introduced by [START_REF] Ledford | Statistics for near independence in multivariate extreme values[END_REF], which is perhaps easier to interpret as it is directly related to the way the joint tail decreases.

Assuming that X and Y have unit Fréchet margins, they suggested the following tail model

P(X > x, Y > x) ∼ l(x) x 1/η , for large x, (1.21)
where η ∈ (0, 1] and l is a slowly varying function (see Definition 1.6). The coefficient η in (1.21) is called the coefficient of tail dependence and constitute another measure of extremal dependence.

Using (1.20), one can show that, provided that the limit exists,

χ = lim x→∞ x 1-1/η l(x).
(1.22)

From (1.22), if η = 1 and lim x→∞ l(x) = c ∈ (0, 1], X and Y are asymptotically dependent with strength of dependence given by χ = c. On the other hand, if η < 1 or if η = 1 and lim x→∞ l(x) = 0, X and Y are asymptotically independent.

Within the class of asymptotic independence, three types of dependence can be distinguished:

To sum up

The coefficient of extremal dependence χ defined in (1.20) measures the strength of asymptotic dependence between the components of a random vector. Unlike the measures introduced in Section 1.2.1, χ is a finite number, easier to handle.

If X ∈ D(G), χ can be linked to the previous measures through

χ = 2 -V (1, 1) = 2(1 -A(1/2)).
In case of asymptotic independence (i.e. χ = 0), other measures have been developed to quantify the dependence structure at high but finite levels. Namely the coefficient of extremal dependence χ and the coefficient of tail dependence η. The relation between the two is given by χ = 2η -1.

Limitations of bivariate EVT

From Theorem 1.17, the bivariate threshold excess model presented in Section 1.2.2, and used in Chapter 2, requires that the underlying distribution F belongs to the domain of attraction of a bivariate EV distribution. For asymptotically dependent variables, such models are fully adequate. Therefore, looking at Figure 1.13, a bivariate threshold excess model seems well adapted to the wave data discussed in Chapter 2.

For some applications, asymptotically independent data can also be encountered (e.g. [START_REF] Heffernan | Hidden regular variation and the rank transform[END_REF][START_REF] Wadsworth | Dependence modelling for spatial extremes[END_REF][START_REF] Towe | Statistical downscaling for future extreme wave heights in the North Sea[END_REF]. For example, in Chapter 3, we will consider asymptotically independent measurements of river discharges from distant stations.

However, if asymptotic independence arises, the class of bivariate extreme value distributions is no longer appropriate. This can be illustrated in several ways. For example, assume that we are in the framework of Section 1.2.1, i.e. we have a set of bivariate vectors {(X i , Y i )} 1≤i≤n , independent and identically distributed, with unit Fréchet margins and common distribution F . Moreover, assume that F ∈ D(G), i.e.

lim n→∞ P [M n /n ≤ (x, y)] = G(x, y). Then, using that log x ∼ x -1 for x → 1, log P [M n /n ≤ (x, y)] = n log F (nx, ny) ∼ n (F (nx, ny) -1) .
Therefore, the above convergence can then be rewritten as

lim n→∞ n P(X 1 /n > x or Y 1 /n > y) = -log G(x, y). 50 1.2. Multivariate extremes This leads to lim n→∞ n P(X 1 /n > x and Y 1 /n > y) = log G(x, y) -log G 1 (x) -log G 2 (y),
where G 1 and G 2 are the marginals of G (see the remark following Definition 1.11). Now going back to our question, assume that we are in the domain of attraction of the asymptotic independence. Then in this case

lim n→∞ n P(X 1 > nx and Y 1 > ny) = 0.
This shows that in the case of asymptotic independence, we can not extrapolate joint probabilities at high but finite levels. Thus, if we only rely on bivariate EVT, asymptotic independence is similar to perfect independence and this is certainly an excessive assumption.

To fill this gap, models adapted to the case of asymptotic independence have been developed. We present two of them in the following section.

Models coping with asymptotic independence

To overcome the above mentioned issue, models that accommodate both asymptotic dependence and asymptotic independence have been introduced. In the following, X and Y are assumed to have unit Fréchet margins.

Note that compared to the bivariate GP model introduced in Section 1.2.2, the following models focus on regions where the two components are simultaneously large. This corresponds for example to the situation in the right-hand side of Building on their first study [START_REF] Ledford | Statistics for near independence in multivariate extreme values[END_REF], [START_REF] Ledford | Modelling dependence within joint tail regions[END_REF] introduced an extension of model (1.21) with a more flexible model given by

P(X > x, Y > y) ∼ L(x, y) x c 1 y c 2 , for large x and y, (1.23)
where c 1 , c 2 > 0 are such that c 1 + c 2 ≥ 1 and L is a bivariate slowly varying function, that is, there exists a function g such that, for all x, y > 0 and c > 0, g(x, y) = lim r→∞ L(rx, ry) L(r, r) and g(cx, cy) = g(x, y).

(1.24)

The homogeneity property in (1.24) implies that g is constant along any ray y = ax for a > 0 and one can introduce a ray dependence function g * defined as g(x, y) = g * (x/(x + y)) = g * (w) for w = x/(x + y) ∈ (0, 1).

Setting η = 1/(c 1 + c 2 )
shows that this latter model (1.23) contains the previous model (1.21).

Example 1.23. Assume that X and Y are two independent Fréchet-distributed variables. Then

P(X > x, Y > y) = [1 -exp(-1/x)] [1 -exp(-1/y)] ∼ L(x, y)(xy) -1 ,
for large x and y, where L(x, y) ≡ 1 is a bivariate slowly varying function. In particular c 1 = c 2 = 1 and η = 1/2 (exact independence).

Finally, building on a simplified version of model (1.23) by considering [START_REF] Ramos | A new class of models for bivariate joint tails[END_REF] proposed the following joint tail model

c 1 = c 2 ,
P(X > x, Y > y) = L(x, y) (xy) 1/(2η) ,
(1.25)

with the same assumptions on L and η as before, i.e. η ∈ (0, 1] and L is a bivariate slowly varying function as in (1.24). This new model aimed to overcome some issues of model 1.23 for which, in some special cases, the joint density obtained was valid only when x/(x + y) ∈ (0, 1).

The main result of [START_REF] Ramos | A new class of models for bivariate joint tails[END_REF] is an alternative spectral representation as in Theorem 1.13 which provides a unified framework covering both asymptotic dependence and asymptotic independence, unlike the classical multivariate extreme value case.

Theorem 1.24 [START_REF] Ramos | A new class of models for bivariate joint tails[END_REF]). If (1.25) holds and if

lim u→∞ P (X > us, Y > ut) P (X > u, Y > u) = g * (s/(s + t)) (st) 1/(2η) ,
for all s, t ≥ 1 where g * is the limit function that is defined following equation (1.24) and η is the coefficient of tail dependence. Then g * (w) satisfies

η -1 g * (w) = 1 -w w 1/(2η) w 0 z 1/η dH η (z) + w 1 -w 1/(2η) 1 w (1 -z) 1/η dH η (z),
where H η is a (non-negative) measure on [0, 1] satisfying the normalisation condition

η -1 = 1/2 0 w 1/η dH η (w) + 1 1/2 (1 -w) 1/η dH η (w).
The measure H η introduced in Theorem 1.24 is the analogue of the spectral measure H defined in Theorem 1.13 and a relationship between the two measures can be derived [START_REF] Ramos | Multivariate joint tail modelling and score tests of independence[END_REF][START_REF] Ramos | A new class of models for bivariate joint tails[END_REF]. For example, if (X, Y ) follows a bivariate extreme value distribution with standard Fréchet margins and spectral measure H, assuming that X and Y are not exactly independent, so that η = 1 (or equivalently χ > 0), then one can show that H = {2 -V (1, 1)} H 1 (see [START_REF] Ramos | A new class of models for bivariate joint tails[END_REF].

In Chapter 3,model (1.25) of [START_REF] Ramos | A new class of models for bivariate joint tails[END_REF] is used to develop a risk measure adapted to a wide range of dependence models, including both 1.3. Tools considered in the following asymptotic dependence and asymptotic independence. In particular, links between model (1.25) and hidden regular variation are discussed in Section 3.1.2.

In this section, only joint tail models based on multivariate regular variation have been introduced. But many other models for asymptotically independent data have been developed in the literature such as the conditional tail model of [START_REF] Heffernan | A conditional approach for multivariate extreme values (with discussion)[END_REF], which avoids the restriction to regions where both components are simultaneously large. de Carvalho and Ramos (2012) provide a review on the existing methods for the statistical modelling of asymptotically independent data.

We give here some possible limitations of the aforementioned appealing models that deal with both asymptotic dependence and asymptotic independence. Firstly, and as already mentioned, [START_REF] Ramos | A new class of models for bivariate joint tails[END_REF] model focuses on regions where both variables are large, this can be a significant limitation for some applications, such as the application considered in Chapter 2 with the wave data. On the other hand, for the flexible model of [START_REF] Heffernan | A conditional approach for multivariate extreme values (with discussion)[END_REF], which removes this constraint, some authors have raised some theoretical issues (e.g. [START_REF] Drees | Conditional extreme value models: fallacies and pitfalls[END_REF].

Finding models that deal with both asymptotic dependence and asymptotic independence is a current intensive research area, and recent advances in this direction have been made (e.g [START_REF] Wadsworth | Modelling across extremal dependence classes[END_REF].

A final remark is that, in our case, the data at hand used in Chapter 2, are asymptotically dependent. This is one of the reasons why we choose to focus and develop bivariate threshold excess models, which are theoretically valid models, relying on EVT.

Tools considered in the following

To sum up Among the different tools presented in this chapter, the ones that we will use in the continuation of this manuscript are the following.

• In Chapter 2, a non-parametric simulator of bivariate GP vectors (see [START_REF] Rootzén | Multivariate generalized Pareto distributions[END_REF] and Section 1.2.2) will be developed to simulate extreme wave events.

• In Chapter 3, a risk measure adapted to binary classification of extreme events in case of both asymptotic dependence and independence will be considered, relying on [START_REF] Ramos | A new class of models for bivariate joint tails[END_REF] model (see also Section 1.2.5).

Chapter 2

Joint stochastic simulation of extreme coastal and offshore significant wave heights

Overview of Chapter 2

The main goal of the following chapter is to propose and study a stochastic simulator that, given offshore conditions, produces jointly offshore and coastal extreme significant wave heights. For that, we rely on bivariate Peaks over Threshold (see Section 1.2.2 and [START_REF] Rootzén | Multivariate generalized Pareto distributions[END_REF]) and develop a non-parametric simulation scheme of bivariate GPD.

To take into account non-stationarities, we also adapt the extended generalised Pareto model [START_REF] Naveau | Modeling jointly low, moderate, and heavy rainfall intensities without a threshold selection[END_REF], letting the parameters vary with specific sea state parameters. Section 2.2 gives some background motivations for the study of extreme coastal wave heights. More specifically, precise definitions of the parameters used to describe a sea state are detailed in Section 2.1.1, along with the databases used in Paper I [START_REF] Legrand | Joint stochastic simulation of extreme coastal and offshore significant wave heights[END_REF]. Section 2.1.2 briefly describe the evolution of the conducted study, motivating the modelling developed in Paper I. This preliminaries are ended in Section 2.1.3, with a brief introduction on the extended generalised Pareto distribution (EGPD), which is adapted in Paper I for a non-stationary modelling.

From Section 2.2 to Appendix 2.B, Paper I is reproduced as it is. In Section 2.3, the sea state data considered in this study are presented. We then develop in Section 2.4 the marginal regression models, incorporating the effects of the peak period and the peak direction on the EGPD parameters. In Section 2.5, our non-parametric method to simulate MGP vectors is presented, along with some numerical experiments. Two algorithms are outlined, one for bivariate simulations and a second one for conditional simulations. Finally, both algorithms are applied in Section 2.6 to the sea state data. The results show that these algorithms successfully simulate new realistic extreme H s events.

Preamble to Paper I

The characterisation of future extreme wave events is crucial because of their multiple impacts, covering a broad range of topics such as coastal flood hazard, coastal erosion, reliability of offshore and coastal structures . . . For instance, for coastal flood risk assessment, Idier et al. ( 2020) combined historical data, statistical and physical models to reconstruct past coastal flood events in a specific site on the French coastline. They showed that the return period associated to historical flood events tends to decrease due to sea-level rise. On a global scale (i.e. along the world coastlines), [START_REF] Marcos | Increased extreme coastal water levels due to the combined action of storm surges and wind waves[END_REF] considered the compound effects of extreme wave heights and storm surges. As mentioned in Section 1. 2.3, Marcos et al. (2019) found that return periods of coastal extreme water levels were significantly overestimated if the dependency structure between wave heights and storm surges was not considered.

Characterisation of wave climatology is also of great concern regarding the implementation of Marine Renewable Energy structures (often shortened to MRE). For example, [START_REF] Stopa | Wave energy resources along the Hawaiian Island chain[END_REF] used a numerical wave model (see Section 2.1.1) to characterise the wave climate in Hawaii for potential implementation sites of wave energy converters.

Concerning the variability of the intensity of the extremes, authors have shown an increase in extreme significant wave heights in the last decades over the North Atlantic Ocean [START_REF] Bertin | A significant increase in wave height in the North Atlantic Ocean over the 20th century[END_REF][START_REF] Young | Multiplatform evaluation of global trends in wind speed and wave height[END_REF], and which could be attributed to climate change [START_REF] Rohmer | Attribution of Extreme Wave Height Records along the North Atlantic Coasts using Hindcast Data: Feasibility and Limitations[END_REF]. Regarding future wave climate projections, under different emission scenarios1 (RCP4.5 and RCP8.5), their is a consensus on a decrease of the mean significant wave heights in the North Atlantic Ocean [START_REF] Aarnes | Projected changes in significant wave height toward the end of the 21st century: Northeast atlantic[END_REF][START_REF] Bricheno | Future wave conditions of Europe, in response to high-end climate change scenarios[END_REF][START_REF] Lemos | On the decreases in North Atlantic significant wave heights from climate projections[END_REF], but it has also been shown a slight increase in the most extreme wave events [START_REF] Bricheno | Future wave conditions of Europe, in response to high-end climate change scenarios[END_REF]. Caution should be exercised here as there is no general consensus on this subject (e.g. [START_REF] Lobeto | Future behavior of wind wave extremes due to climate change[END_REF]. Note that in order to address such issues, initiatives such as COWCLIP (Coordinated Ocean Wave Climate Project) have been carried out with the aim, among others, of producing a unified database of global wave climate projections [START_REF] Morim | A global ensemble of ocean wave climate projections from CMIP5-driven models[END_REF].

All these studies show the importance of predicting extreme wave heights accurately, and more particularly in coastal regions. In Paper I [START_REF] Legrand | Joint stochastic simulation of extreme coastal and offshore significant wave heights[END_REF], we propose a first attempt in this direction by looking at how moderately high offshore wave heights can produce large coastal significant wave heights. For this, we will propose a joint stochastic simulator which, given offshore sea states conditions, produces extreme offshore and coastal significant wave heights.

Wave data

Spectral definition of sea state parameters

In this section, we provide a precise definition of significant wave height, denoted H s . More details on ocean waves can be found, for example, in [START_REF] Holthuijsen | Waves in oceanic and coastal waters[END_REF]. Recall that in Chapter 1, we said that H s gave a measure of the ocean surface roughness. More precisely, this quantity can be defined from the energy of a sea state.

A sea state is the characterisation, in a given area and for a limited period, of the sea surface elevation, which is assumed to be a stationary process. Because waves are generated from various weather systems, mainly wind sea (waves generated locally from the wind) and swell (that have travelled from a remote generating area), the description of this sea surface can be complex.

To overcome this issue, a sea state can be decomposed into several simple wave trains (or harmonic waves), each having a particular direction of propagation and wave period, the so called spectral decomposition (see Figure 2.1). From this spectral decomposition, one can consider the spectral energy density function of each simple wave train involved, denoted E(f, θ), where θ is the direction of propagation and f is the wave frequency (i.e. the inverse of the wave period). Note that it is a convention to consider the frequency rather than the period. A representation of the total spectral energy of a sea state is then given by the so-called wave spectrum, which gives the spectral energy of each wave train given its direction θ and frequency f . An example is given in Figure 2.2. From this representation, one can identify the different wave regimes:

• swells that correspond to a relatively narrow directional dispersion and to low frequencies: they are visually rather regular and long-crested,

• wind sea corresponding to a much broader directional dispersion and to higher frequencies: corresponding to irregular and short-crested waves. To synthesise the information contained in this 2-dimensional spectrum, several statistics are defined. These parameters are then used to describe a sea state. We give here the formal definitions of the spectral parameters used in Paper I. More details are given, for example, in Chapter 4 of [START_REF] Holthuijsen | Waves in oceanic and coastal waters[END_REF].

• The significant wave height H s (in meters) is estimated from the total variance of the moving sea surface elevation, denoted ⟨η 2 ⟩, which is shown to be the integrated spectral density in both frequency and direction:

H s = 4 ⟨η 2 ⟩ = 4 ∞ 0 2π 0 E(f, θ)df dθ.
The significant wave height is usually interpreted as a measure of the wave energy of the sea state. Historically, it is also defined as the average of the highest one-third of wave heights.

• The peak period T p (in seconds) is the period associated to the maximum of the omnidirectional spectrum, it corresponds to the period of the dominant waves in the sea state and is defined as

T p = arg max f E(f ) -1 = arg max f 2π 0 E(f, θ)dθ -1
.

• The peak direction D p (in degrees) is the direction associated to the maximum of energy of the wave spectrum, defined as

D p = arg max θ ∞ 0 E(f, θ)df.
As an example, for the wave spectrum in Figure 2.2, H s = 0.66m, T p = 3.57s and D p = 117 • , indicating that the corresponding sea state is rather dominated by a wind sea regime and with dominant waves coming from the North-West.

These statistics are limited as they give a limited picture of a sea state, as explained by [START_REF] Holthuijsen | Waves in oceanic and coastal waters[END_REF]

(p.25):

For instance, wave conditions may well be similar in the sense that the significant wave height and period are equal, but they may still be very different in detail: a mixed sea state of wind sea (short, irregular, locally generated waves) and swell (long, smooth waves, generated in a distant storm) may have the same significant wave height and period as a slightly higher wind sea without swell.

However, for a long-term description of the wave climate, these statistics are sufficient and are conventionally used, especially for coastal studies.

Wave data sources

Significant wave height is obtained either from direct measurements of the sea surface elevation, using in situ measurements (such as buoys) or remote-sensing techniques (e.g. satellites), or from numerical wave models. Numerical wave models produce simulated wave spectrum based on the physical equation of wave energy balance [START_REF] Filipot | A unified spectral parameterization for wave breaking: From the deep ocean to the surf zone[END_REF], taking into account different forcing (such as winds, currents, sea ice, bathymetry,. . . ).

Note that in situ measurements generally provide estimates of the omnidirectional spectrum E(f ), whereas numerical models can estimate the full wave spectrum E(f, θ), but they are considered to be closer to the truth than satellite or numerical data and therefore used often as reference (e.g. [START_REF] Ailliot | Space-time models for moving fields with an application to significant wave height fields[END_REF]. Historical outputs from numerical wave models are often called hindcast data (in contrast to forecast data), since it consist of numerical predictions in the past. In this thesis, two hindcast databases based on the numerical wave model WAVEWATCH-III [START_REF] Tolman | User manual and system documentation of WAVEWATCH III version 4[END_REF] are considered (more details are given later, Section 2.3). Both have been developed by IFREMER (French Research Institute for Exploitation of the Sea) and provide sea state statistics (as H s , T p and D p ). For the two databases, the numerical wave model is forced by the CFSR wind reanalysis2 [START_REF] Saha | The NCEP Climate Forecast System Reanalysis[END_REF]. Their main characteristics are the following:

• IOWAGA [START_REF] Ardhuin | IOWAGA[END_REF]: covers the whole globe with a regular 0.5 • resolution grid in latitude and longitude (see Figure 2.

On the choice of appropriate covariates for modelling extreme coastal sea states

The starting point of this study was to characterise the conditional distribution of extreme significant wave height, considered as exceedances over a high threshold, close to the French coast.

We chose to focus on one specific point from the HOMERE database, represented in yellow in Figure 2.4, which corresponds to the SEM-REV sea test site [START_REF] Mouslim | Development of the french wave energy test site SEM-REV[END_REF]. The main advantage of this specific location is that it is well-documented and is used for the development and optimisation of renewable energy structures, it is therefore of great interest to describe extreme sea states at this specific site.

Since waves are generated by winds, a first part of the study was devoted to find a conditional model for extreme H s given wind fields. Indeed, waves observed at the point of interest depends on wind conditions over the North Atlantic Ocean.

The large spatial extent of wind fields implied the need to identify the characteristics of the spatio-temporal wind processes associated with extreme H s and to derive a representation in a finite dimensional space, in order to facilitate the calculation of the above mentioned model. For that, we considered the 10m wind components U 10 and V 10 from the ERA5 global reanalysis database [START_REF] Hersbach | The ERA5 global reanalysis[END_REF]. This dataset was used in a preliminary approach, it was chosen mainly for accessibility reasons, although in a second approach it would have been more relevant to consider the CFRS reanalysis dataset since it corresponds to the wind forcing of the wave data. Figure 2.4 illustrates the spatial footprint of the wind field corresponding to the maximum of H s in the dataset. This preliminary work highlighted, for extreme H s , dominant sectors of wind direction. But it also raised some limitations due to the large spatial extent of wind fields and the complex processes involved in the generation of extreme sea states, which imply the superposition of very different wave regimes (swells and wind sea). [START_REF] Obakrim | Statistical modeling of the space-time relation between wind and significant wave height[END_REF] addressed this issue using a downscaling approach based on weather-types. This approach was also investigated by [START_REF] Michel | Deep learning for statistical downscaling of sea states[END_REF], relying on deep learning methods. One drawback of these studies is that it requires to take into account the whole wind field over the North Atlantic Ocean.

To avoid this, we decided to exploit the somewhat space-time Markovian behaviour of extreme H s . Indeed, there is a strong spatial and temporal dependency, notably in the extremes, among the nearest data points due to the wave propagation (see Figure 2.5). Furthermore, the benefits of considering a wave field instead of a wind field is that the former contains most of the information resulting from the complex wave generation processes. The advantage of the global IOWAGA database is that it can be run rapidly. Whereas HOMERE is slower to run, due to its strong refinement. Thus, finding a relationship between the large-scale model (represented by IOWAGA) and the small-scale model (represented by HOMERE) would allow rapid calculations of the probabilities of extreme H s at the location of interest.

In Paper I [START_REF] Legrand | Joint stochastic simulation of extreme coastal and offshore significant wave heights[END_REF], we will built on the relation between a largescale grid point and a small-scale grid point of interest, near the coast (see Figure 2.7). The proposed model will allow for a statistical downscaling of extreme H s from the global wave model to the coastal model.

Extended generalised Pareto distribution

Among extreme value methodologies, GP models (see Section 1.1.2) are powerful tools for modelling extreme climate observations. However a major drawback of such methods is the threshold choice. To answer this problem, threshold se-lection methods have been proposed (see, e.g. Raymond-Belzile (2019), Chapter 1).

Another drawback is that by definition of GP models, only values that are above the threshold are considered. But for some applications, one needs to model the entire range of the data. For instance, this is the case in Paper I, where we want to model H s exceedances given moderately high values of H s , i.e. we look at quantities of the following type

H s -u | H s > u,
but for moderate threshold u.

In this context, different models have been proposed to describe the range of observations. For example, studies have considered piece-wise models, for which one consider a GP model above the threshold and a parametric model for the bulk below the threshold, for wave applications people considered for instance a truncated Weibull [START_REF] Randell | Bayesian inference for nonstationary marginal extremes[END_REF] or a truncated Gamma distribution [START_REF] Ross | Efficient estimation of return value distributions from non-stationary marginal extreme value models using Bayesian inference[END_REF].

However such models are discontinuous at the threshold. To alleviate this, [START_REF] Carreau | A hybrid pareto model for asymmetric fat-tailed data: the univariate case[END_REF] developed a mixture GP model, with a smooth transition at the threshold, from an hybrid GP distribution that combines a Gaussian and a GP tail. Yet this model is rather complicated to use in practice due to the strong continuity constraints.

Recently, [START_REF] Naveau | Modeling jointly low, moderate, and heavy rainfall intensities without a threshold selection[END_REF] (see also the references therein) proposed a model that bypasses the threshold selection and models the entire range of data, the so called extended generalised Pareto distribution (EGPD).

The EGPD is a distribution defined on (0, ∞) whose distribution function is given by

F (x) = G [H ξ,σ (x)] ,
where G is a continuous distribution function on [0, 1] and H ξ,σ is a GP distribution function as defined in Section 1.1.2 (Chapter 1). In [START_REF] Naveau | Modeling jointly low, moderate, and heavy rainfall intensities without a threshold selection[END_REF], several constraints are required for G in order to ensure the lower and upper tails behaviours. Four different parametric models for G are derived.

In this study we will only use the first one, which in practice appears to be flexible enough [START_REF] Tencaliec | Flexible semiparametric generalized Pareto modeling of the entire range of rainfall amount[END_REF][START_REF] De Carvalho | An extreme value bayesian lasso for the conditional left and right tails[END_REF][START_REF] Rivoire | A comparison of moderate and extreme ERA-5 daily precipitation with two observational data sets[END_REF][START_REF] Gall | Improved regional frequency analysis of rainfall data[END_REF], and given by G(v) = v κ . This model has therefore three parameters: a scale σ > 0, a shape ξ ∈ R (controlling the upper-tail behaviour) and a second shape parameter κ > 0 (which controls the shape of the lower tail). Figure 2.6 depicts the density of the EGPD for different values of κ.

In Paper I bellow, the EGPD will be used for marginal modelling, as described in Equation (2.2). The review of [START_REF] Ewans | Statistical modelling of extreme ocean environments for marine design: A review[END_REF] on the statistical modelling of extreme significant wave heights highlights the necessity of taking into account covariates effects. For this, Generalised Additive Models (GAMs) provide a nice flexible framework [START_REF] Hastie | Generalized Additive Models[END_REF]. Recall that the simplest regression model, the linear model, can be written in compact form as

Y ∼ N (µ, σ 2 ), µ = Xβ,
where X is the design matrix and β the vector of parameters. A GAM is a generalisation of the linear model, that allows for response distributions other than Gaussian, for a certain degree of non-linearity in the model and with linear predictors specified as smooth function of the covariates. It can be written as

       Y ∼ some exponential family distribution, g(E[Y ]) = Xβ + J j=1 s j (x j ),
where g is a smooth monotonic link function and s j are smooth functions of some covariates (x 1 , . . . , x J ).

The smooth functions s j are typically represented as linear combinations of basis functions b j,k

s j (x) = K k=1 α j,k b j,k (x), 1 ≤ j ≤ J,
where α j,k are some coefficients to be estimated. The basis functions b j,k are completely known, and different basis (also called splines) can be considered for a same model, such as natural cubic splines, thin plate splines, tensor product splines (see [START_REF] Perperoglou | A review of spline function procedures in r[END_REF] for a comprehensive review, with applications in R). A complete introduction on GAMs can be found in [START_REF] Wood | Generalized Additive Models: An Introduction with R[END_REF].

The above GAM considers only single-parameter distributions with mean parameter that vary as a smooth function of some covariates. Such models have then been extended to distribution families with more than one parameter. For example, regarding extreme value models, Chavez-Demoulin and [START_REF] Chavez-Demoulin | Generalized additive modelling of sample extremes[END_REF] proposed to model the behaviour of extremes with GAMs by allowing the GPD parameters to be represented as smooth functions of covariates. [START_REF] Yee | Vector generalized linear and additive extreme value models[END_REF] linked an extended version of GAMs (named Vector Generalised Additive Models, VGAM3 ) with EVT. All these models are built for extreme value models. However, as discussed previously, we wish to model values that are not necessarily extreme, relying on the EGPD.

To model non-stationarity with EGPD models, Le Carrer ( 2022) developed an add-on to the gamlss package [START_REF] Stasinopoulos | Instructions on how to use the gamlss package in R Second Edition[END_REF]. GAMLSS (Generalised Additive Models for Location, Scale and Shape) allows for GAM forms of (up to) four-parameter distributions. In the case of the EGPD, it can be written as

         Y ∼ EGP D(ξ, κ, σ) g 1 (ξ) = X 1 β 1 + s 1,1 (x 1,1 ) + • • • + s 1,J 1 (x 1,J 1 ), g 2 (κ) = X 2 β 2 + s 2,1 (x 2,1 ) + • • • + s 2,J 2 (x 2,J 2 ), g 3 (σ) = X 3 β 3 + s 3,1 (x 3,1 ) + • • • + s 3,J 3 (x 3,J 3 ),
where, for i ∈ {1, 2, 3}, g i is a link function, X i is the design matrix containing the linear additive terms, β i is the vector of linear parameters and (s i,j ) 1≤j≤J i are the smoothing functions of the explanatory variables (x i,1 , . . . , x i,J i ). The implementation of Le Carrer (2022) allows to consider negative shape parameters ξ by setting g 1 as an identity link, allowing an unconstrained search domain for the inference of ξ. Note that the model fitting is performed using maximum penalised likelihood estimation. This is used in Paper I to let the scale parameter σ of the EGPD vary as a smooth function of the peak period and peak direction (see Equation (2.2)).

A final remark is that besides GAM related models, other recent approaches have been considered to let the parameters of the EGPD vary given some covariates. For example, de Carvalho et al. ( 2021) developed a Bayesian approach to learn the effect of covariates on an EGPD-based model. * * * Below, Paper I [START_REF] Legrand | Joint stochastic simulation of extreme coastal and offshore significant wave heights[END_REF] is reproduced as it will be submitted. It was written and developed by Juliette L. under the supervision of Pierre Ailliot, Philippe Naveau and Nicolas Raillard.

Introduction

French coastlines have been particularly affected by extreme maritime events in the past (Nicolae [START_REF] Lerma | Spatial variability of extreme wave height along the Atlantic and channel French coast[END_REF]. Co-occurrence of high tidal coefficients, atmospheric surge conditions and specific sea states can lead to extreme maritime events. These events are particularly crucial for assessing flooding risks and their consequences [START_REF] Genovese | Storm surge disaster risk management: the Xynthia case study in France[END_REF][START_REF] Bertin | Importance of wave age and resonance in storm surges: The case Xynthia, Bay of Biscay[END_REF]. According to the special IPCC report (see [START_REF] Collins | Extremes, abrupt changes and managing risk[END_REF]), extreme wave heights, which contribute to these extreme maritime events, have increased over the past few years. The recent IPCC report [START_REF] Seneviratne | Weather and climate extreme events in a changing climate[END_REF] indicates, with high confidence, an increase in the occurrence and magnitude of such coastal events in the future. More specifically, [START_REF] Caires | 100-year return value estimates for ocean wind speed and significant wave height from the ERA-40 data[END_REF] showed that the most extreme wave conditions were expected to occur in the North Atlantic, which include the Bay of Biscay, our study area.

Sea surface elevation over a geographical area results from the superposition of waves generated by local winds and by remote swell (generated in distant regions). The characterisation of this complex surface is called a sea state and to describe it, various parameters are available. In this work, we will focus on three variables: the significant wave height denoted H s [m], the peak period, T p [s], and the peak direction, D p [°] (see e.g. [START_REF] Holthuijsen | Waves in oceanic and coastal waters[END_REF] for more details).

From a coastal risk point of view, a fundamental question is to determine how moderately high offshore significant wave heights can produce large coastal H s .

Peak direction and peak period influence the relationship between coastal and offshore H s . In this context, our main goal is to propose and study a stochastic simulator that, given offshore conditions (T p , D p , H o moderately high), produces jointly offshore and coastal extreme significant wave heights. The left-hand side map of Figure 2.7 shows the two locations of interest.

From such a stochastic generator (1st goal), many products can be derived. In particular our second objective is to propose a conditional simulation model (2nd goal). The framework for each step of this study is summarised in Table 2.1. In order to make the two simulation models as flexible as possible, non-parametric algorithms are derived using resampling techniques (or non-parametric bootstrap [START_REF] Efron | Bootstrap methods: Another look at the jackknife[END_REF]). While this study focuses on simulation of extreme H s , and as illustrated by the numerical simulations in Section 2.5, the two non-parametric algorithms developed could be applied to a broad range of data.

In multivariate extreme value analysis, one is often interested in the joint behaviour of the variables as they become large. As illustrated by the right-hand side of Figure 2.7, which depicts a scatter plot between H o and H c , large values tend to occur simultaneously. For this specific type of dependence, called asymptotic dependence [START_REF] Coles | An Introduction to Statistical Modeling of Extreme Values[END_REF], models from the class of multivariate Extreme Value Theory (EVT) can be used. To achieve the two objectives (1st Goal and 2nd Goal), we will therefore propose two simulation algorithms based on multivariate Peaks over Thresholds (Sec. 8.3.1 Beirlant et al., 2004). Note that to assess whether data fall within the class of asymptotic dependence or not, summary statistics have been developed such as the dependence measures χ and χ [START_REF] Coles | Dependence measures for extreme value analyses[END_REF] and plots of these measures (not shown here) confirm strong dependence between large values of H o and H c .

For weakly dependent extremes, conditional models based on [START_REF] Heffernan | A conditional approach for multivariate extreme values (with discussion)[END_REF] should be favoured to deal with our 2nd goal, see for example [START_REF] Towe | Statistical downscaling for future extreme wave heights in the North Sea[END_REF]; [START_REF] Shooter | On spatial conditional extremes for ocean storm severity[END_REF]; [START_REF] Tendijck | Modeling the extremes of bivariate mixture distributions with application to oceanographic data[END_REF].

Before modelling the joint behaviour of large values, it is necessary to model margins [START_REF] Beirlant | Statistics of Extremes: Theory and Applications[END_REF]. To visualise this task with respect to our data, the empirical histograms displayed in Figure 2.8 indicate that, given {H o > v o } (moderately high offshore significant wave heights), a traditional univariate extreme value approach based on fitting a generalised Pareto distribution (GPD) to the exceedances [START_REF] Coles | An Introduction to Statistical Modeling of Extreme Values[END_REF] is not appropriate. To tackle this issue, we use the extended generalised Pareto distribution (EGPD) introduced by [START_REF] Naveau | Modeling jointly low, moderate, and heavy rainfall intensities without a threshold selection[END_REF] which handles this type of setting (see also [START_REF] Papastathopoulos | Extended generalised Pareto models for tail estimation[END_REF]), more details are given in Section 2.4. Furthermore, like many other environmental data, extreme H s are non--stationary with respect to covariates [START_REF] Ewans | Statistical modelling of extreme ocean environments for marine design: A review[END_REF] and marginal models need to take into account this non-stationarity (e.g. [START_REF] Ewans | The effect of directionality on northern north sea extreme wave design criteria[END_REF][START_REF] Méndez | Seasonality and duration in extreme value distributions of significant wave height[END_REF][START_REF] Casas-Prat | A physical-based statistical method for modeling ocean wave heights[END_REF]. To incorporate non-stationarities, Chavez-Demoulin and [START_REF] Chavez-Demoulin | Generalized additive modelling of sample extremes[END_REF] proposed to let the parameters of an extreme value model vary as smooth functions of covariates. This has been intensively applied to oceanographic data (e.g. [START_REF] Feld | Estimation of storm peak and intra-storm directional-seasonal design conditions in the north sea[END_REF][START_REF] Jonathan | Non-stationary conditional extremes of northern North Sea storm characteristics[END_REF][START_REF] Ross | Efficient estimation of return value distributions from non-stationary marginal extreme value models using Bayesian inference[END_REF]. However, there are only a few papers dealing with nonstationary EGPD [START_REF] De Carvalho | An extreme value bayesian lasso for the conditional left and right tails[END_REF][START_REF] Haruna | Performance-based comparison of regionalization methods to improve the at-site estimates of daily precipitation[END_REF]. In this study, and as illustrated in Figure 2.9, the marginal EGPD models parameters will be described as smooth functions of the covariates T p and D p .

The key steps of our study are the following: (1) marginal regression modelling within the class of EGPD, (2) transformation of the data to common exponential margins, (3) modelling extremal dependence between the variables using multivariate generalised Pareto model (hereafter MGP models) [START_REF] Rootzén | Multivariate generalized Pareto distributions[END_REF], (4) non-parametric simulation of bivariate extreme H s within the class of MGP distributions. In our modelling scheme, different steps are novelties.To our knowledge, little attention has been paid in the literature to the modelling of mul-2.3. Sea state data Figure 2.9: Dependence of H s with respect to T p and D p (left) Offshore significant wave heights H o given peak period T p . (right) Estimated quantiles of H o given the peak direction D p for different quantile levels q ∈ {0.25, 0.5, 0.75, 0.9, 0.99}, estimation is performed using smoothed quantile regression [START_REF] Koenker | Quantile smoothing splines[END_REF]. tivariate non-stationary extremes using EGPD and to the non-parametric simulation within the MGP class.

Our paper is organised as follows. In Section 2.3, the sea state data are presented and the marginal inference incorporating covariates in the EGPD modelling is described in Section 2.4. In Section 2.5, the non-parametric method to simulate MGP vectors is presented and some numerical experiments are shown. Two algorithms are outlined, one for bivariate simulations and a second one for conditional simulations. Both algorithms are applied in Section 2.6 to the sea state data.

Sea state data

Our study is carried out in the northern part of the Bay of Biscay in France. The specificity of this region is that it is exposed to the Atlantic Ocean and therefore subject to complex superpositions of wind generated waves and swell. The data are extracted from two different wave hindcasts provided by IFREMER and consist of simulations of sea states by a numerical model. First, the IOWAGA database [START_REF] Ardhuin | IOWAGA[END_REF] corresponds to sea states parameters that are generated by the wave model WAVEWATCH-III and forced by CFSR winds on a Global grid (0.5 • resolution grid in latitude and longitude). HOMERE is the second database [START_REF] Accensi | HOMERE[END_REF], also based on WAVEWATCH-III model and forced by IOWAGA on the wet boundaries, but on an unstructured grid covering only the English Channel and the Bay of Biscay, more refined close to the coast and with the inclusion of currents and water levels.

As mentioned in the introduction, we restrict our attention to two specific locations: an offshore grid point (47 2.7, highlighting a strong dependence structure between the variables. Data are then split into two sets:

• Set 1 contains the first 70% of the data and is used for the inference of the marginal regression models and the preliminary steps for the simulation of H s (see Section 2.6);

• Set 2 contains the remaining 30% of the data and is used for the simulation of extreme H s . 2013)), which appears to be flexible enough, and whose cumulative distribution function is of the form

Marginal regression analysis

Marginal regression

F (x) = 1 -1 + ξx σ -1/ξ κ .
(2.1)

Marginal regression analysis

The model has three parameters: scale σ > 0, shape ξ ∈ R and an additional parameter κ > 0 which controls the shape of the lower tail. In [START_REF] Naveau | Modeling jointly low, moderate, and heavy rainfall intensities without a threshold selection[END_REF], the authors developed the EGPD for non-negative shape parameter ξ. Indeed, the main applications were the modelling of daily rainfall (e.g. [START_REF] Naveau | Modeling jointly low, moderate, and heavy rainfall intensities without a threshold selection[END_REF][START_REF] Tencaliec | Flexible semiparametric generalized Pareto modeling of the entire range of rainfall amount[END_REF][START_REF] De Carvalho | An extreme value bayesian lasso for the conditional left and right tails[END_REF][START_REF] Rivoire | A comparison of moderate and extreme ERA-5 daily precipitation with two observational data sets[END_REF], which are heavytailed (ξ > 0). In our case, as reported in [START_REF] Ewans | Statistical modelling of extreme ocean environments for marine design: A review[END_REF], extreme H s data are generally described by upper-bounded tail distributions. Still, the case ξ < 0 can be handled by model (2.2).

As mentioned in the Introduction, H s data are non-stationary (see Figure 2.9), see also [START_REF] Jonathan | Non-stationary conditional extremes of northern North Sea storm characteristics[END_REF]; [START_REF] Feld | Estimation of storm peak and intra-storm directional-seasonal design conditions in the north sea[END_REF]; De [START_REF] De Leo | Non-stationary extreme value analysis of sea states based on linear trends. Analysis of annual maxima series of significant wave height and peak period in the Mediterranean Sea[END_REF]. Therefore, we regress H s on the peak direction D p and the peak period T p . We choose here to put the dependency on the scale parameter. The regression marginal models can then be written as follows

           P(H o -v o ≤ x|H o > v o , T p , D p ) = 1 -1 + ξ o x σ o (T p , D p ) -1/ξo κo , P(H c -v c ≤ x|H o > v o , T p , D p ) = 1 -1 + ξ c x σ c (T p , D p ) -1/ξc κc . (2.2)
The coastal marginal threshold v c introduced in (2.2) is defined by v c := min(H c ; H o > v o ) so that the minimum of H c -v c is equal to zero. The regression marginal models (2.2) are estimated using the R package gamlss [START_REF] Stasinopoulos | Instructions on how to use the gamlss package in R Second Edition[END_REF] with the EGPD family (Le Carrer, 2022). The inference is performed using maximum penalised likelihood estimation (note that the model fitting is achieved with the CG algorithm [START_REF] Cole | Smoothing reference centile curves: The LMS method and penalized likelihood[END_REF]). In our model (2.2), we assume that the parameters σ o and σ c vary smoothly with T p and D p . This is achieved using tensor product of cubic regression splines. The parameter estimates for the regression marginal models are reported in Table 2.2. Asymptotic 95% confidence intervals for each parameter are given in brackets and derived from the asymptotic variance-covariance matrix of the fitted models. Table 2.2: Estimated parameters for the regression marginal models. 95% asymptotic confidence intervals are given in brackets.

Both estimated shape parameters are negative but close to zero, which suggests light-tailed or bounded distributions. This is in accordance with previous studies and the physical behaviour of wave heights in shallow waters [START_REF] Castillo | Engineering analysis of extreme value data: Selection of models[END_REF][START_REF] Vanem | A truncated, translated Weibull distribution for shallow water sea states[END_REF]. The goodness of fit of the model (2.2) are shown in Appendix 2.A.

Covariates effects

To visualise the effect of the covariates on the scale parameters, we choose to consider the theoretical expectation of the fitted EGPD models. As can be seen from Equation (2.3) (and similarly for H o ), the theoretical expectation of model (2.2) is directly proportional to the scale parameter (see [START_REF] Naveau | Modeling jointly low, moderate, and heavy rainfall intensities without a threshold selection[END_REF]):

E(H c |H o > v o , T p , D p ) = σ c (T p , D p ) 1 ξ c [κ c B(κ c , 1 -ξ c ) -1] , (2.3)
where B denotes the Beta function defined by The first panel in Figure 2.10 represents the empirical ratio H c /H o given T p and D p values. Local polynomial interpolation (LOESS, [START_REF] Cleveland | Locally weighted regression: An approach to regression analysis by local fitting[END_REF]) is performed between the observed data points to get values on a regular grid of D p and T p , ranging from 230 to 300 degrees for the peak direction and from 10 to 20 seconds for the peak period. This first panel is then compared to the estimated ratio of the two conditional expectations (second panel in Figure 2.10), which, from Equation (2.3), is proportional to

B(a, b) = 1 0 t a-1 (1 -t) b-1 dt.
σc (T p , D p ) /σ o (T p , D p ) .
This ratio can give us an idea of the propagation of the wave energy from the offshore to the coast as a function of the covariates. These results can be physically interpreted: the loss of wave energy between the offshore and the coast is lower for small periods but also for waves coming from the SW rather than the NW due to the bathymetry (see the map on Figure 2.7). In the third panel of Figure 2.10, the dark blue region, corresponding to T p > 16s and D p < 255 • , has no data and consequently, produces large differences due to extrapolation issues.

Using the estimated σc (T p , D p ) and σo (T p , D p ), the H s data are then transformed to common exponential scale using the probability integral transform

H E o := -log 1 -Fo [(H o -v o )/σ o (T p , D p )] , H E c := -log 1 -Fc [(H c -v c )/σ c (T p , D p )] , (2.4)
where Fo (resp. Fc ) is the fitted EGP D( ξo , κo , 1) cdf's (resp. EGP D( ξc , κc , 1)) from Equation (2.2).

Multivariate Pareto model

In this section, the threshold exceedances of H s transformed to common exponential margins are modelled. This vector is denoted

H E := (H E c , H E o ) in the fol- lowing.
For that, we adapt the definition of [START_REF] Rootzén | Multivariate generalized Pareto distributions[END_REF] of bivariate threshold exceedances given as

H E -u|H E ≰ u (2.
Z = E + T -max(T ), (2.6)
with T a random vector and E a unit exponential random variable independent of T .

Note that the standard MGP distribution is supported by the set

L := {x ∈ R d ; x ≰ 0}.
In our study, Equation (2.6) is adapted to Z = (Z 1 , Z 2 ) defined as

Z 1 := H E o -u o |H E o > u o or H E c > u c , Z 2 := H E c -u c |H E o > u o or H E c > u c .
(2.7) [START_REF] Kiriliouk | Peaks over thresholds modeling with multivariate generalized Pareto distributions[END_REF] established several parametric MGP models by setting explicit densities for T in a multivariate setting. In the following, we consider only vectors of dimension 2, i.e. Z = (Z 1 , Z 2 ) and T = (T 1 , T 2 ). To bypass the choice of the underlying distribution for T , we start from the following rewriting of Equation (2.6),

Simulation of bivariate standard generalised Pareto distributed vectors

Z 1 = E + ∆1 ∆<0 , Z 2 = E -∆1 ∆≥0 , (2.8)
where ∆ := Z 1 -Z 2 = T 1 -T 2 and 1 A denotes the indicator function, equals to 1 if A is true and 0 otherwise. Equation (2.8) is the basis for our simulation algorithms. From this equation, we see that we need to simulate values of ∆ and E independently, instead of (T 1 , T 2 ). Generating independent, and identically distributed, unit exponential is trivial, so the main difficulty is to simulate ∆. This can be achieved by bootstrapping (see [START_REF] Efron | Bootstrap methods: Another look at the jackknife[END_REF]). Our approach is then described in Algorithm 1 and a theoretical proof can be found in Appendix 2.B.

Define ∆ i := Z 1,i -Z 2,i for 1 ≤ i ≤ n 5: Generate m realisations E (m) k ∼ Exp(1), independently of (∆ i ) 1≤i≤n , for 1 ≤ k ≤ m 6: Bootstrap m realisations ∆ (m) k , 1 ≤ k ≤ m, from (∆ 1 , . . . , ∆ n ) 7: end procedure 8: return Z (m) 1,k := E (m) k + ∆ (m) k 1 ∆ (m) k <0 and Z (m) 2,k := E (m) k -∆ (m) k 1 ∆ (m) k >0 , for 1 ≤ k ≤ m

Numerical experiments

In the following, we simulate MGP vectors Z = (Z 1 , Z 2 ) from the representation (2.6) with different parametric models on (T 1 , T 2 ) and we compare with our simulation algorithm. The different experiments are reported in Table 2.3 and some graphical results are shown in Figure 2.11 which displays for each model a scatter plot of the data, the measure of extremal dependence χ(u) for increasing values of u, and the marginal quantile-quantile plots. We use the measure χ(u) which gives a measure of asymptotic dependence between two variables X and Y (for more details see e.g. [START_REF] Coles | Dependence measures for extreme value analyses[END_REF]) and which is defined by

χ(u) := P Y > F -1 Y (u) | X > F -1 X (u) , u ∈ (0, 1).

Bivariate model

Joint distribution of T Parameters (a) Gaussian symmetric

N (µ 1 , µ 2 ), 1 ρ ρ 1 µ 1 = 0 µ 2 = 0 ρ = 0.4 (b) Gaussian asymmetric N (µ 1 , µ 2 ), 1 ρ ρ 1 µ 1 = 0 µ 2 = 2 ρ = 0.4 (c) Logistic F (x 1 , x 2 ) = 1 + e -x 1 /σ 1 + e -x 2 /σ 2 -1 σ 1 = 1 σ 2 = 5 (d) Gumbel F (x 1 , x 2 ) = exp [-exp{-x 1 /σ 1 }] exp [-exp{-x 2 /σ 2 }] σ 1 = 1 σ 2 = 4 (e) Exponential S(x 1 , x 2 ) = exp {-λ 1 x 1 -λ 2 x 2 -λ 3 max(x 1 , x 2 )} λ 1 = 2 λ 2 = 10 λ 3 = 1 Table 2.
3: Overview of the different experiments carried-out. For each, we give the joint distribution F (x 1 , x 2 ) when it writes easily or the survival function S(x 1 , x 2 ).

In the third column, we give the different parameters values used in the numerical experiments.

In Table 2.3, (a) and (b) are two bivariate Gaussian models with same correlation coefficient ρ < 1 but with µ 1 ̸ = µ 2 for (b), leading to asymmetry. Model (c) corresponds to the Type I bivariate logistic distribution proposed by [START_REF] Gumbel | Bivariate logistic distributions[END_REF]. For model (d) we consider two independent Gumbel distributed variables with different scale parameters. And lastly, (e) corresponds to a bivariate exponential distribution as defined in [START_REF] Marshall | A multivariate exponential distribution[END_REF].

As seen in Figure 2.11, Algorithm 1 successfully simulates draws from the parametric simulations in terms of the marginals Z 1 and Z 2 , but also recovers well the dependence structure when looking at the measure of dependence χ(u).

Conditional simulation within the MGP class

From an application perspective, we also want to be able to simulate conditionally on one of the two variables. In this section we describe the conditional simulation algorithm for the MGP model. From Equation (2.6) we can derive

Z 2 = Z 1 + T 2 -T 1 = Z 1 -∆.
(2.9)

From Equation ( 2.9), we can design a simulation strategy but caution is required because ∆ and Z 1 are not necessarily independents. But for some values of Z 1 , this will be the case. To see this, one can compute the conditional distribution of ∆ given Z 1 = z 1 starting from the joint distribution function of (Z 1 , Z 2 )

which is given by

f (Z 1 ,Z 2 ) (z 1 , z 2 ) = e -max(z 1 ,z 2 ) f ∆ (z 1 -z 2 )1 (z 1 ,z 2 )∈L , for z 1 , z 2 ∈ R,
where f ∆ denotes the distribution function of ∆. 1st case: If z 1 > 0. In this case, noting that if z 1 > 0 then 1 (z 1 ,z 2 )∈L = 1, the marginal distribution of Z 1 is given as follows

f Z 1 (z 1 ) = e -z 1 K, where K := 0 -∞ e u f ∆ (u)du + ∞ 0 f ∆ (u)du and does not depend on z 1 . Therefore, the conditional distribution of Z 2 given Z 1 = z 1 when z 1 > 0 is given by f Z 2 |Z 1 (z 2 | z 1 ) = 1 K f ∆ (z 1 -z 2 )1 z 1 ≥z 2 + e z 1 -z 2 f ∆ (z 1 -z 2 )1 z 1 <z 2
From this, the conditional distribution of ∆ given Z 1 = z 1 > 0 can then be deduced:

f ∆|Z 1 (δ | z 1 ) = 1 K f ∆ (δ)1 δ≥0 + e δ f ∆ (δ)1 δ<0 .
(2.10)

This shows that, conditionally on Z 1 > 0, ∆ does not depend on Z 1 . 2nd case: If z 1 < 0. Then, noting that

1 (z 1 ,z 2 )∈L = 1 z 2 >0 if z 1 < 0, we get f Z 1 (z 1 ) = e -z 1 K(z 1 ),
where K(z 1 ) := z 1

-∞ e u f ∆ (u)du. And we can derive the conditional distribution of

Z 2 given Z 1 = z 1 < 0 as follows f Z 2 |Z 1 (z 2 | z 1 ) = 1 K(z 1 ) e z 1 -z 2 f ∆ (z 1 -z 2 )1 z 2 >0 .
The conditional distribution of ∆ given Z 1 = z 1 < 0 is then given by

f ∆|Z 1 (δ | z 1 ) = 1 K(z 1 ) e δ f ∆ (δ)1 δ<z 1 .
(2.11)

From Equations (2.10) and ( 2.11) we derive the conditional simulation algorithm described in Algorithm 2, where the simulation procedure is split into two cases:

1. If z 1 > 0, we can sample values of ∆ independently of Z 1 , 2. otherwise, if z 1 < 0, we use a rejection sampling approach to approximate the targeted conditional density in Equation (2.11).

Similarly, we could also derive a simulation scheme of Z 1 given Z 2 = z 2 .

Algorithm 2 Non-parametric conditional MGP simulation

1: input A sample (∆ i ) 1≤i≤n ; a realisation z 1 of Z 1 2: output A simulated sample (Z (m) 2,k ) 1≤k≤m conditionally on Z 1 = z 1 , potentially with m ̸ = n 3: procedure 4: if z 1 > 0 then 5: Define ∆ |Z + 1 the subset of (∆ i ) 1≤i≤n such that Z 1 > 0 6: Bootstrap m realisations ∆ (m) k , 1 ≤ k ≤ m, from ∆ |Z + 1 independently of Z 1 7: else 8: for 1 ≤ k ≤ m do 9: Sample one realisation ∆ (m) k from (∆ i ) 1≤i≤n independently of Z 1 10: Generate a random number u ∈ [0, 1] 11: while u > exp(∆ (m) k )1 ∆ (m) k <z 1 do 12:
Repeat steps 9 and 10 13: end for 14: end procedure

15: return Z (m) 2,k := z 1 -∆ (m) k for 1 ≤ k ≤ m

Numerical experiment continued

As for the bivariate simulations, we can illustrate Algorithm 2 with numerical experiments. We choose here to show the results only for Model (a) (Symmetric Gaussian) since for this specific model we have an explicit form for the theoretical distribution of ∆. The results are presented in Figure 2.12 where we simulated the conditional distribution of Z 2 for eight different conditioning values. The sampled and theoretical conditional distributions appear to be in close conformity.

Application to extreme significant wave height

The methodology presented in Section 2.5 is applied to the joint and the conditional simulations of extreme significant wave heights. For that, the sample of bivariate threshold exceedances (Z 1 , Z 2 ) defined in Equation (2.7) is used as input data for Algorithm 1 or Algorithm 2. The thresholds u o and u c in (2.7) are defined as the 0.8 quantile of H E o , or equivalently of H E c . Recall that in both cases, simulations are performed on the exponential scale.

A final step of back transformation is then necessary to get simulations of H s on the original scale. This final step corresponds to part 3 (resp. 5) in the following procedure for the joint (resp. conditional) simulation of H s . For the sake of clarity we now divide the joint and the conditional simulation scheme of H s in two

Joint simulation of significant wave heights

The joint simulation scheme for extreme H s is described hereafter. In the following we fix the pair value (t p , d p ) ∈ R 2 which may be taken from Set 2.

1. Compute σo (t p , d p ) and σc (t p , d p ) from the marginal EGPD models fitted on Set 1 (see Section 2.4).

2. Simulate m pairs of (z 1 , z 2 ) applying Algorithm 1 with input data (Z 1 , Z 2 ) as defined in (2.7).

We therefore obtain m simulated pairs ((z 1,1 , z 2,1 ), . . . , (z 1,m , z 2,m )) for a fixed value (t p , d p ).

Transform the simulated values to the original scale

h o,i := σo (t p , d p ) F -1 o (1 -e -(z 1,i +uo) ) + v o ∈ R m , h c,i := σc (t p , d p ) F -1 c (1 -e -(z 2,i +uc) ) + v c ∈ R m ,
where

F -1 o (resp. F -1 c ) is the inverse cdf of the EGP D( ξo , κo , 1) (resp. EGP D( ξc , κc , 1)) estimated in Section 2.4.
This procedure is then applied to four selected pairs (t p , d p ) from the Set 2 corresponding to the four largest H c of this dataset. Figure 2.13 depicts simulated pairs of offshore and coastal H s with simulation sample size m = 1000.

Date-Time

H c [m] H o [m] T p [s] D p [°]
1/m m i=1 1{h o,i > h o , h c,i > h c },
applying the above procedure with simulation sample size m = 1.10 6 . This estimated probability quantify the observations made with Figure 2.13: for the most extreme events, the associated joint probabilities are expected to be lower.

Conditional simulation of coastal significant wave heights

The conditional simulation scheme for extreme coastal H s is described bellow. Note that the procedure is very symmetrical for simulating offshore H s . In the following we fix the triplet value (h o , t p , d p ) ∈ R 3 which may be taken from Set 2.

1. Compute σo (t p , d p ) and σc (t p , d p ) from the marginal EGPD models fitted on Set 1 (see Section 2.4).

2. Transform h o to the standardised space using the probability integral transform:

h E o = -log 1 -Fo [(h o -v o )/σ o (t p , d p )]
where Fo is the EGP ( ξo , κo , 1) cdf's.

Set z

1 := h E o -u o ,
:= (Z 1 -Z 2 )1{Z 1 > 0}, with Z 1 , Z 2
as defined in (2.7), and bootstrapped m times. We therefore obtain m simulations of z 2 = (z 2,1 , . . . , z 2,m ) for a fixed triplet (h o , t p , d p ), given z 1 > 0.

Transform the predicted values to the original scale

h c,i := σc (t p , d p ) F -1 c (1 -e -(z 2,i +uc) ) + v c where F -1 c
is the inverse cdf of the EGP D( ξc , κc , 1) (see Section 2.4).

Note that in

Step 4 above, the simulation is restricted to the case when z 1 > 0 for convenience, since our focus is on the simulation of extreme H s . The pseudo-algorithm described above is applied with all the triplet values (h o , t p , d p ) from Set 2, with simulation sample size m = 1000. These conditional simulations of coastal significant wave heights are then compared to the true values of H c from Set 2. The overall coverage probability (i.e. the number of times the actual value of coastal H s is within the 95% range of the predicted distribution) is equal to 95% and the simulations are shown in Figure 2.14. The simulations and the true H c values (red dots) are most of the time in good agreement. Since no declustering approach has been adopted, consecutive observations, which belong to the same storm event, are kept. They are depicted with identical colour, highlighting a temporal dependence structure between each storm cluster. The effect of the covariates T p and D p in the conditional simulations is depicted in Figure 2.15, showing that the simulation model is able to simulate both the most intense and the more moderate H c . This plot also highlight for which sea state conditions the simulations are far from the observed values. It appears that the two predictions such that the observed H c value does not fall within the 95% simulation range correspond to small H c and H o .

Discussion and conclusions

Simulation of extreme events in a multivariate setting is of great interest to capture not only the statistical behaviour of the extremes, but also the dependence between large values of complex processes. Based on the multivariate EVT, this work presents two non-parametric simulation algorithms of bivariate generalised Pareto distributed variables, without assuming any specific parametric shape for the MGP model. Thanks to Algorithm 1, one can simulate joint extremes. As for Algorithm 2, it allows the simulation of conditional extremes. Both methods have been validated with numerical simulations.

We would like to point out that in the context of bivariate extremes, other simulation algorithms have been developed. For example [START_REF] Marcon | A semi-parametric stochastic generator for bivariate extreme events[END_REF] proposed a simulation method with a semi-parametric structure for the extremal dependence function, but it was not based on a MGP model and did not cover non-stationarities. [START_REF] Michel | Simulation and estimation in multivariate generalized Pareto models[END_REF] derived a non-parametric simulation framework of bivariate generalised Pareto variables using a different representation of MGP vectors than the one used in this paper.

For application purposes and as a byproduct, a non-stationary marginal modelling with the EGPD was also developed, adding covariate effects on the scale parameter of the EGPD using smoothing splines.

We applied this work to the simulation of extreme significant wave heights near the Brittany coast given specific offshore sea state conditions (T p , D p ) with compelling results. In both joint and conditional settings, thanks to the presented algorithms, we are able to simulate realistic extreme H s events.

Note that in possible extensions of this work to climate projections, it is assumed that data will not be available at the coastal location but only on a coarse grid, similar to the IOWAGA Global hindcast. This argument is illustrated in Table 2.1, and is favour of the first pre-selection of the H s data through {H o > v o } for the marginal regression analysis (see Section 2.4).

Extensions to the multivariate case will be the subject of future works. Considering more than two locations raises different modelling issues. It would also be interesting to apply this methodology to other locations in order to ensure the proper generalisation of the methodology.

2.B. Proof of Algorithm 1

We show here the goodness of fit for the marginal regression models defined in Equation (2.2). As our models depend on some covariates, the diagnostic plots presented here are built for the standardised H s exceedances which are defined as (H c -v c ) /σ c (T p , D p ) (and similarly for H o ). From Figure 2.16, one can see that the fits seem to be fairly good, a slight discrepancy in the lower values can be noticed for the coastal model but this is not a major issue as the interest lies mainly in the larger values.

2.B Proof of Algorithm 1

With the same notations as in Algorithm 1, let F be the common distribution function of ∆ 1 , . . . , ∆ n and F (m) nm be the empirical distribution function of the bootstrap sample

∆ (m) 1 , . . . , ∆ (m) m . Lemma 2.1. If F (m)
nm converges in distribution to F , as n and m tend to infinity, then (Z

(m) 1,k , Z (m) 2,k ) 1≤k≤m converge in distribution to a bivariate GPD G where G is the com- mon distribution function of the sample (Z 1,i , Z 2,i ) 1≤i≤n . Proof. As P (E ≤ u) = 1 -min(1, exp(-u)) for any u ∈ R if E ∼ Exp(1), the bivariate distribution function of (Z (m) 1 , Z (m) 2 ) is equal to P (Z (m) 1 , Z (m) 2 ) ≤ (z 1 , z 2 ) = 1 -E min 1, e -min(z 1 -∆ (m) 1 ∆ (m) <0 ,z 2 +∆ (m) 1 ∆ (m) >0 ) = 1 -E min 1, e -min(z 1 -∆ (m) ,z 2 )-max(∆ (m) ,0) , for any (z 1 , z 2 ) ∈ {x ∈ R 2 ; x ≰ 0}.
Then, one can show that the function x → min(1, x), defined for x ≥ 0, is Lipschitz and bounded by 1. And applying the Portmanteau theorem, we have, letting min(n, m) → ∞,

P (Z (m) 1 , Z (m) 2 ) ≤ (z 1 , z 2 ) | ∆ 1 , . . . , ∆ n → 1 -E min 1, e -min(z 1 -∆,z 2 )-max(∆,0) = 1 -E min 1, e max(T 1 -z 1 ,T 2 -z 2 )-max(T 1 ,T 2 ) .
Which is the cumulative distribution function of the MGP vector (Z 1 , Z 2 ) as defined in Rootzén et al. (2018a) (Prop. 8).

The assumption in Lemma 2.1 is linked to the bootstrap asymptotic theory (e.g. [START_REF] Bickel | Some asymptotic theory for the bootstrap[END_REF]).

CHAPTER 2. JOINT STOCHASTIC SIMULATION OF EXTREME COASTAL AND OFFSHORE SIGNIFICANT WAVE HEIGHTS

Summary of Paper I

To sum up

• Paper I addressed the specific issue of joint simulation of extreme environmental variables.

• In the context of asymptotic dependence, multivariate EVT provides a solid mathematical framework.

→ Marginal treatment: non-stationary EGPD modelling, taking into account covariate effects on the scale parameter through smoothing splines.

→ Dependence treatment: development of non-parametric simulation of bivariate GP vectors, based on a simple rewriting.

→ Two algorithms for the simulation of extreme H s given offshore conditions (T p , D p , H o moderately high), are outlined:

• A first stochastic simulator that produces jointly offshore and coastal extreme significant wave heights.

• A second stochastic simulator that produces coastal extreme significant wave heights given offshore significant wave heights.

̸ → Underlying hypothesis of asymptotic dependence between H o and H c , otherwise other models should be considered.

̸ → Focus on a specific site: to ensure the generalisation of the method, other locations should be investigated.

Chapter 3

Evaluation of binary classifiers for asymptotically dependent and independent extremes

Overview of Chapter 3

This chapter addresses another issue, that can be seen as a subsequent step to the work of Chapter 2. Instead of building one prediction model, assume that we have several prediction models available -or rather several forecasts -and we want to compare the performance of each model. Since this thesis focuses on EVT, this comparison should be made given their capacity to capture extreme occurrences.

In the following, we study the specific case of binary classifiers, which are the simplest type of forecasting and decision-making situation: an extreme event did or did not occur. Such classifier tailored for extremes will be called an extremal classifier and risk functions that answer our initial question will be developed. Their properties will be derived under the framework of multivariate regular variation and hidden regular variation, allowing to handle finer -or more specific -types of asymptotic independence.

Section 3.1 gives some preliminary notions on multivariate regular variations and the refinement to hidden regular variation. Their definition and some examples are provided.

Section 3.2 up to Section 3.6 reproduces Paper II [START_REF] Legrand | Evaluation of binary classifiers for asymptotically dependent and independent extremes[END_REF]. First, the construction of our risk functions is developed by investigating some specific cases: beginning with two naive classifiers and then looking at a particular type of asymptotic independence (corresponding to hidden regular variation). We then propose an empirical estimator of our risk functions, whose inferential properties are derived under the framework of multivariate regular variation and hidden regular variation. A simulation study compares different classifiers and indicates their performance with respect to our risk functions. To conclude, we apply our framework to the analysis of extreme river discharges in the Danube river basin. The application compares different predictive algorithms and tests their capacity at forecasting river discharges from other river stations.

Finally, Section 3.7 gives some further developments. We study in detail the special class of linear classifiers and show that the optimisation of our risk functions leads to a consistent solution. A tool to identify the explanatory variables that contribute the most to extremal behaviour is developed and applied to the river network data. To conclude, ideas for extending this pre-selection tool to the case of hidden regular variation are presented.

Preamble to Paper II

In the previous chapter, we proposed a stochastic simulator for bivariate extremes in the case of asymptotic dependence. Specifically, the conditional Algorithm 2 proposed in Paper I could be viewed as a statistical method to predict extreme values of a given variable (e.g. extreme coastal significant wave heights). Now imagine that we have a prediction model for the entire range of data, not only the largest values. Then, a natural question in EVT is whether the model performs well in predicting the most extreme values, and can we discriminate between different models.

This question is in some ways similar to the notion of scoring rules used in the forecasting literature where the aim is to provide summary measures of how well a probabilistic forecast performs compared to the reality, but also in comparison with other competing forecasts. A good introduction can be found, for example, in [START_REF] Jolliffe | Forecast Verification: a Practitioner's Guide in Atmospheric Science[END_REF]. Evaluation of probabilistic forecasts is beyond the scope of this manuscript, but we can mention some recent advances in forecast of extreme events, reflecting various challenges that it raises. For example, considering rare events as binary outputs (which is actually analogous to our approach as explained below), [START_REF] Stephenson | The extreme dependency score: a non-vanishing measure for forecasts of rare events[END_REF] discussed the advantages and drawbacks of various risk functions. To obtain accurate extreme wind gusts predictions, [START_REF] Friederichs | Forecast verification for extreme value distributions with an application to probabilistic peak wind prediction[END_REF] derived closed-form expressions of the continuous ranked probability score (CRPS) (e.g. [START_REF] Gneiting | Probabilistic forecasts, calibration and sharpness[END_REF] for extreme value-distributed data. A weighted CRPS [START_REF] Gneiting | Comparing density forecasts using thresholdand quantile-weighted scoring rules[END_REF] has been applied to compare spatio-temporal predictions of extreme sea surface temperatures during the EVA 2019 data competition [START_REF] Huser | EVA 2019 data competition on spatio-temporal prediction of Red Sea surface temperature extremes[END_REF]. Finally, [START_REF] Brehmer | Why scoring functions cannot assess tail properties[END_REF] and [START_REF] Taillardat | Extreme events evaluation using CRPS distributions[END_REF] showed the limitations of scoring functions to assess tail properties.

Binary classifiers and extremal risk

In Paper II [START_REF] Legrand | Evaluation of binary classifiers for asymptotically dependent and independent extremes[END_REF]) hereunder, we will not consider forecast models directly but rather binary classifiers. In this context, the previous question now becomes how to assess and compare binary classifiers, based upon a set of multivariate observations X and denoted g(X; u), with respect to their capacity to capture 'extreme' occurrences of a binary response Y (u) ∈ {-1, +1}. To set the scene, one can imagine that the labels of Y (u) encode the fact that a given random variable H is above (i.e. extreme) or below (i.e. non-extreme) a given threshold u (hence the above notations which depend on u):

Y (u) = +1, if H > u, -1, otherwise.
As is common practice in EVT, we will consider high thresholds u and call g(X; u)

an extremal classifier if it verifies some convergence conditions as u becomes large. But then, and due to the inherent rarities of extreme events, there will be many more -1 labels than +1. This leads somehow to a binary classification problem with imbalanced classes and a wealth of models have been developed in the machine learning community to cope with imbalanced data classification (see, e.g. [START_REF] Haixiang | Learning from class-imbalanced data: Review of methods and applications[END_REF] for a recent review).

In this context, our objective is not to construct such classifiers but rather to design an appropriate risk function allowing for their comparison. This risk function will be defined as

R (u) (g) = P(g(X; u) ̸ = Y (u) ) P(Y (u) = 1 or g(X; u) = 1)
and will depend on the threshold u considered. Subsequently, and in the same vein as the tail dependence coefficient χ(u) [START_REF] Coles | Dependence measures for extreme value analyses[END_REF], one may want to look at R (u) for high levels of u > 0 to focus on the performance of g(X; u) in the extremes.

This risk function R (u) can then be used to compare different predictive models and test their capacity at forecasting extreme events, as for instance extreme river discharges (see Section 3.5). In Section 3.3, the statistical inference properties of our risk function R (u) are derived under the framework of multivariate regular variation, which is an alternative way to describe the tail behaviour of multivariate extremes. The following section gives a brief overview of this concept, for more information see, e.g. [START_REF] Resnick | Heavy-tail phenomena: probabilistic and statistical modeling[END_REF].

Refinement of asymptotic independence

As illustrated in Chapter 2, when working in the framework of asymptotic dependence, models based on bivariate extreme value theory usually suit well. However, and as discussed in the following, in the case of binary classifiers, it can be complex to find such dependent models.To overcome this problem, a finer risk measure adapted to the case of asymptotic independence is developed. For this purpose, we will rely on the framework of [START_REF] Ramos | A new class of models for bivariate joint tails[END_REF] with joint tail models adapted to the case of asymptotic independence (see Section 1.2.5). This specific framework handles the case of hidden regular variations which can be viewed as a second-order regular variation when asymptotic independence arises [START_REF] Resnick | Hidden regular variation, second order regular variation and asymptotic independence[END_REF]. Hereinafter, we give more details on hidden regular variation, starting from the definition of multivariate regular variation, and its connection with the model of [START_REF] Ramos | A new class of models for bivariate joint tails[END_REF].

Multivariate regular variation

Formally, a d-dimensional non-negative random vector X, with identical margins (for simplicity), is multivariate regularly varying with limit measure ν if there exists a function b(t) → ∞ as t → ∞ and a non-negative Radon measure ν ̸ = 0 such that It can be shown that condition (3.1) implies that there exists a constant α > 0, the tail index, such that, for all relatively compact sets A ⊂ [0, ∞] d \ {0} and t > 0, ν(tA) = t -α ν(A). This property is fundamental in order to obtain the polar decomposition of Equation (3.3). Moreover, b is a regularly varying function with index 1/α (see [START_REF] Resnick | Heavy-tail phenomena: probabilistic and statistical modeling[END_REF].

t P X b(t) ∈ • v - → ν, (3.1)
The link between the univariate case as discussed in Section 1.1.3 and the above definition of multivariate regular variation may not be directly evident. Yet, equivalent definitions of univariate and multivariate regular variation in terms of convergence of measures and tail probabilities can be derived (see Theorems 3.6 and 6.1 in Resnick (2007)). Theorem 3.1. Let X be a non-negative d-dimensional random vector. The following statements are equivalent.

1. X, with distribution function F , is multivariate regularly varying with measure ν and tail index α > 0. for all points x > 0 which are continuity points of the function x → ν ([0, x] c ), where [0, x] c denotes the complement set taken with respect to [0, ∞] d \ {0} (see Figure 3.1 for an illustration in dimension 2).

The following convergence holds

lim t→∞ 1 -F (tx) 1 -F (t1) = lim t→∞ P(X/t ∈ [0, x] c ) P(X/t ∈ [0, 1] c ) = ν ([0, x] c ) , (3.2)
3. There exists an α-Pareto random variable P and, independently of P , a random vector Θ ∈ [0, ∞) d on the unit sphere {x ∈ R d ; ∥x∥ = 1}, such that

P ∥X∥ t , X ∥X∥ ∈ • | ∥X∥ > t v - → P [(P, Θ) ∈ •] , as t → ∞, (3.3
)

for an arbitrary norm ∥•∥ on R d . Equation (3.
2) is the analogue version of univariate regularly varying random variables as defined in (1.7). It implies that the survival function 1-F is multivariate regularly varying (see [START_REF] Resnick | Heavy-tail phenomena: probabilistic and statistical modeling[END_REF], Chapter 6).

In the polar decomposition (3.3), Θ is referred to as the spectral tail vector and its distribution as the spectral measure. A nice interpretation is that the spectral measure, or equivalently ν, places mass in directions where large events occur [START_REF] Meyer | High-dimensional Learning for Extremes[END_REF][START_REF] Meyer | Sparse regular variation[END_REF]. This decomposition will be used in Section 3.7 to develop a special case of a linear classifier. One further comment is that, as for the univariate case, multivariate regular variation and multivariate EVT are two strongly connected notions. Indeed, if X is multivariate regularly varying as in (3.1) with limit measure ν and tail index α > 0, then X is in the domain of attraction of an extreme value distribution G with [START_REF] Beirlant | Statistics of Extremes: Theory and Applications[END_REF]). Thus, the limit measure ν and the exponent measure V (defined in (1.10) for the bivariate case) are in the end the same quantity, i.e. V (x) = ν ([0, x] c ).

G(x) = exp {-ν ([0, x] c )} , (see, for example

Hidden regular variation

In view of the foregoing, hidden regular variation is defined as a refinement of regular variation on the following subset

E 0 := s ∈ [0, ∞] d \ {0} : For some 1 ≤ i < j ≤ d, min(s i , s j ) > 0 .
This set E 0 corresponds to the points of [0, ∞] d \ {0} such that at most d -2 coordinates are 0 (see Figure 3.1). Then, we say that X is hidden regularly varying if in addition to (3.1), there exists a non-decreasing function b 0 (t) → ∞ such that b(t)/b 0 (t) → ∞, as t → ∞, and a Radon measure ν 0 ̸ = 0 on E 0 and such that on the set E 0 . Similarly to multivariate regular variation, if (3.4) holds, there exists α 0 ≥ α such that b 0 is regularly varying with index 1/α 0 and ν 0 (t•) = t -α 0 ν(•), for all t > 0. Now, if X is hidden regularly varying, then necessarily the components of X are asymptotically independent, this is shown in [START_REF] Resnick | Hidden regular variation, second order regular variation and asymptotic independence[END_REF] and with Chapter 5 of [START_REF] Resnick | Extreme Values, Regular Variation and Point Processes[END_REF]. For simplicity, we illustrate this for d = 2, considering X = (X 1 , X 2 ). In this case,

t P X b 0 (t) ∈ • v - → ν 0 (3.4)
E 0 = (0, ∞] 2 . Assuming (3.1) and (3.4) hold, then, since b(t)/b 0 (t) → ∞ as t → ∞, lim t→∞ t P X 1 b(t) > δ, X 2 b(t) > δ = lim t→∞ t P min(X 1 , X 2 ) b 0 (t) > b(t) b 0 (t) δ = 0
for any δ > 0. That is, ν (E 0 ) = 0 and ν concentrates on the axes, i.e. on the lines {(x 1 , 0); x 1 > 0} and {(0, x 2 ); x 2 > 0}.

Therefore, the terminology hidden can be understood by the fact that the dependence structure is hidden by the mass put on the axes, hence a normalisation of smaller order, with b 0 , is needed to capture the finer structure that may be present away from the axes.

Example 3.2. If X := (X 1 , X 2 ) with X 1 and X 2 two i.i.d. Pareto variables with param- eter 1/2, that is P(X i > x) = x -1/2 . Then, setting b 0 (t) = 1/t, X is hidden regularly varying with ν 0 ((x 1 , ∞] × (x 2 , ∞]) = (x 1 x 2 ) -1/2 . This is illustrated in Figure 3.2.
Example 3.3. Consider X := BY + (1 -B)U , with B a Bernoulli variable, Y having regularly varying marginals with index -1, and U multivariate regularly varying with tail index 1 < α < 2, B, Y and U being independents. Then X is hidden regularly varying [START_REF] Resnick | Hidden regular variation, second order regular variation and asymptotic independence[END_REF].

Hidden regular variation implies multivariate regular variation and asymptotic independence, but the converse is in general not true (see [START_REF] Resnick | Heavy-tail phenomena: probabilistic and statistical modeling[END_REF]. For inference or simulation purposes, [START_REF] Das | Models with hidden regular variation: Generation and detection[END_REF] review different generative models exhibiting hidden regular variation, namely the mixture model developed by [START_REF] Maulik | Characterizations and examples of hidden regular variation[END_REF] and the additive model of [START_REF] Weller | A sum characterization of hidden regular variation with likelihood inference via expectation-maximization[END_REF]. [START_REF] Das | Models with hidden regular variation: Generation and detection[END_REF] also suggest diagnostic tools in order to detect regular variation and hidden regular variation in multivariate data.

Link with Ramos and Ledford (2009) model

In dimension d = 2, it appears that the model of [START_REF] Ramos | A new class of models for bivariate joint tails[END_REF] discussed in Section 1.2.5 is a special case of hidden regular variation, where the coefficient of tail dependence η is given by η = 1/α 0 (see, e.g. [START_REF] Heffernan | Hidden regular variation and the rank transform[END_REF]).

Indeed, assume that

X = (X 1 , X 2 ) satisfies (1.25). Consider b 0 (t) := 2U ← (t) 1 where U (t) := (2t) 1/η L(2t, 2t) . Then for any x 1 , x 2 > 0, lim t→∞ t P X 1 b 0 (t) > x 1 , X 2 b 0 (t) > x 2 = g(x 1 , x 2 )(x 1 x 2 ) -1/(2η) ∼ (x 1 x 2 ) -1/(2η) .
In the above, we assumed that the components of X were identically distributed. However, with real data analysis, such behaviour is rarely observed. To overcome this issue, [START_REF] Heffernan | Hidden regular variation and the rank transform[END_REF] showed that the use of the rank transform method, resulting in equal marginals, preserved hidden regular variation. Yet, in the following Paper II [START_REF] Legrand | Evaluation of binary classifiers for asymptotically dependent and independent extremes[END_REF], equal marginal transformation is not desired as explained later. Instead an extension of [START_REF] Ramos | A new class of models for bivariate joint tails[END_REF] is developed, allowing for different tail indexes between the components. * * * Below, Paper II [START_REF] Legrand | Evaluation of binary classifiers for asymptotically dependent and independent extremes[END_REF] is reproduced as submitted. It was written in collaboration with Philippe Naveau and Marco Oesting. In this paper, Juliette L. wrote the Sections 3.4 and 3.5. Juliette L. wrote the R code used and produced the figures in the article. The theoretical proofs were mainly derived by Marco Oesting. All three authors contributed equally to the writing of the remaining text.

The proofs of all lemmas and propositions are placed in Appendix A, as originally done in Paper II.

Introduction

In binary classification, one typically considers data of the form (X, Y ) ⊤ where Y ∈ {-1, 1} represents a binary response to the input X ∈ [0, ∞) d . In this paper, we focus on the case that Y = Y (u) represents the occurrence of an extreme event, Y (u) = 1 indicating that a random quantity H crosses a level u, called the threshold, and Y (u) = -1 otherwise, that is

Y (u) = +1, if H > u, -1, otherwise.
(3.5)

In the following, for simplicity, we focus on the case that H is a non-negative random variable such that P(H > u) > 0 for all u > 0 and its upper end point is infinite.

In extreme value analysis, one is interested in the behaviour of Y (u) for high levels, that is for u → ∞, and, therefore, also any classifier g needs to be adapted to the threshold u. Thus, for every u > 0, let g(•; u) be a measurable function from R d to {-1, 1}. In order to evaluate the quality of the classification at a certain level u, we consider a loss function l u that assigns a cost to a classifier g(•; u) and a realisation (x, y (u) ). Here, it is important to note that, by definition of rare events, P(Y (u) = 1) is very small and P(Y (u) = -1) is close to one as u gets large. This imbalance can lead to atypical and/or undesirable comparisons of classifiers. For example, the "always optimistic" classifier that never forecasts an extreme can be defined as g(X; u) ≡ -1, almost surely. To see how to handle this naive classifier, the classical risk function defined as the expectation of the indicator 1{g(X; u) ̸ = Y (u) } can be written as P g(X; u) ̸ = Y (u) = P(H > u, g(X; u) = -1) + P(H ≤ u, g(X; u) = 1). (3.6) critical success index. In particular, these authors linked forecast scoring rules with two dependence indices used in EVT (see [START_REF] Coles | Dependence measures for extreme value analyses[END_REF]

χ = lim u→1 P(U > u | V > u) and χ = lim u→1 2 log(P(U > u) log(P(U > u, V > u)) -1 ,
where the two random variables U and V follow the same continuous uniform distribution on [0, 1]. The choice of uniform marginals can be made whenever the forecast can be assumed to be calibrated, i.e. observations and forecasts follow the same marginal distributions and can be transformed into uniforms. Concerning the extremal dependence strength between U and V , if χ > 0, then the variables U and V are said to be asymptotically dependent and χ = 1. If χ = 0, then the variables U and V are said to be asymptotically independent and χ < 1 captures some second order extremal dependence information. [START_REF] Stephenson | The extreme dependency score: a non-vanishing measure for forecasts of rare events[END_REF] advocated the use of χ and called it the extreme dependency score. Later on, [START_REF] Ferro | Extremal dependence indices: Improved verification measures for deterministic forecasts of rare binary events[END_REF] proposed two different scores and studied their properties. But the link with the concept of asymptotic independence was not clear and the convergence results of their estimators were not fully developed. In contrast to χ, one drawback of χ is that its formula is not easy to explain to practitioners. In comparison, R(g) as a type of the critical success index can be interpreted with ease. Hence, it is of interest to extend this definition to the asymptotic independent case. In the machine learning literature, [START_REF] Jalalzai | Climate Change 2014: Synthesis Report. Contribution of Working Groups I, II and III to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change[END_REF] also worked on binary classifiers for extremes. But they did not focus on R (u) (g). Instead, they studied a different setting where the object of interest was

P(g(X) ̸ = Y | ∥X∥ > u)
where ∥X∥ represents a norm with u large. Hence, their conditioning event was {∥X∥ > u}, while our conditioning depends on Y with the set {Y (u) = 1 or g(X; u) = 1}, see Equation (3.7). So, their interest was centred on the classifier performance when the norm of the explanatory vector X was large. Our focus is on large values of H in the production of extreme events of the type Y (u) = 1 when H > u, see Equation (3.5). [START_REF] Jalalzai | Climate Change 2014: Synthesis Report. Contribution of Working Groups I, II and III to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change[END_REF] provided various theoretical results based on the main assumption that the conditional distribution of X given Y = ±1 was regularly varying with an angular measure that depends on Y = ±1.

In this study, one part of our results is based on the concept of hidden regular variation (see, e.g. [START_REF] Ledford | Statistics for near independence in multivariate extreme values[END_REF][START_REF] Heffernan | Hidden regular variation and the rank transform[END_REF][START_REF] Ferro | A probability model for verifying deterministic forecasts of extreme events[END_REF]. In particular, we take advantage of the model of [START_REF] Ramos | A new class of models for bivariate joint tails[END_REF] to derive the asymptotic properties of our estimators.

Our paper is organised as follows. In Section 3.3, we propose and study a risk function that can handle both the asymptotic dependent and independent cases. Estimators are also constructed and their asymptotic properties derived. Section 3.4 focuses on a simulation example that highlights the difficulty to compare common classifiers in the case of asymptotic independence. In Section 3.5, we revisit the well studied example of the Danube river application and see how the choice of the metric can change the ranking of classifiers. Note that, besides the proofs of all propositions, the appendix addresses the questions of how to optimise the linear classifier for extremes and how to choose the relevant features, see Section 3.7.

Risk, upper tail equivalence and extremal dependence

The following lemma provides a flexible blueprint to link risk functions with probabilities based on general sets. We will apply it under different setups linked to extreme events.

Lemma 3.4. Let A ε be a sequence of measurable sets of increasing sizes with decreasing ε ∈ [0, 1], in particular A 1 ⊆ A ε ⊆ A 0 . Let B ε be the same type of set sequence such that P(A 1 ∩ B 1 ) > 0. The following ratio R(A ε , B ε ) can be written as

R(A ε , B ε ) := P(A 1 △B 1 | A ε ∩ B ε ) P(A 1 ∪ B 1 | A ε ∩ B ε ) = 1- 1 P(B 1 |A 1 ∩ B ε ) + 1 P(A 1 | A ε ∩ B 1 ) -1 -1 , (3.8)
where A 1 △B 1 denotes the difference set. In addition, we have the three following properties for R ε :

1. R(A ε , B ε ) is non-increasing in ε with R(A 1 , B 1 ) = 0.

Let A ′

ε be another sequence of measurable sets of increasing sizes with decreasing ε.

If A 1 = A ′ 1 and A ε ⊆ A ′ ε for some ε ∈ [0, 1) then R(A ε , B ε ) ≤ R(A ′ ε , B ε ).
3. If for any ε ∈ [0, 1] and ε ′ ∈ [0, 1], there exists some positive constants a and b and some positive function c ε,ε ′ such that

P(A ε ∩ B ε ′ ) = c ε,ε ′ P a (A ε ) P b (B ε ′ ), (3.9) then R(A ε , B ε ) = 1 - c ε,1 c 1,1 (P(A 1 | A ε )) -a + c 1,ε c 1,1 (P(B 1 | B ε )) -b -1 -1 . (3.10)
We deduce from Equation (3.8) that R(A ε , B ε ) = 0 if and only if

P(B 1 | A 1 ∩ B ε ) = P(A 1 | A ε ∩ B 1 ) = 1.
The second property of this lemma indicates that, for a given ε, the risk function becomes smaller if the set A ε is as small as possible.

Equation (3.9) can be viewed as a mixing condition that leads to a simple expression of R(A ε , B ε ′ ) based on disjoint events. We want to use the ratio (3.8) in order to generalise the risk R (u) (g) defined in (3.7). In the latter definition, adapted to extremes, the set {H > u} was considered and one could set B 1 = {H > u}, for instance, and consequently, B ε = {H > εu} would be a natural sequence if H is regularly varying. The choice of A ε is open and it can play an important role2 . In the next paragraphs, we set ε = 0 and discuss the choice of A 1 and A 0 . In Section 3.3.1, we will work with ε > 0. In that case, we will call the corresponding risk the conditional risk.

In the rest of this paper, we restrict our attention on a particular form of classifiers g(X; u) = +1, if g(X) > u, -1, otherwise, (3.11) for some function g : R d → (0, ∞). The function g(.) does not have to be a norm.

It can be understood as any projection/summary of the explanatory variables X onto the positive real line (see, e.g. [START_REF] Aghbalou | Tail inverse regression for dimension reduction with extreme response[END_REF], for projection techniques for extremes). This corresponds to the set A 1 = {g(X) > u} in Lemma 3.4. A special case of this lemma is to set ε = 0 and when A 0 and B 0 are equal to the full set, i.e. P(A 0 ) = P(B 0 ) = 1, and A 1 = {g(X) > u} and B 1 = {H > u}. In this case, Equation (3.8) tells us that R (u) (g) defined by (3.7) satisfies

R (u) (g) = 1 - 1 P(H > u | g(X) > u) + 1 P(g(X) > u | H > u) -1 -1
.

This leads to the following expression of R (u) (g)

R (u) (g) = 1 - P(g(X) > u | H > u) 1 -P(g(X) > u | H > u) + P(g(X) > u)/ P(H > u)
.

(3.12)

Within the class defined by Equation (3.11), the effect of the marginal distributions of H and g(X) on R(g) 3 can be explained. When g(X) has a lighter tail than H, i.e. P(g(X) > u)/ P(H > u) → 0, then we also have P(g(X) > u | H > u) → 0, and consequently

R (u) (g) → 1, as u → ∞.
In the case where g(X) possesses a heavier tail than H, i.e. P(g(X) > u)/ P(H > u) → ∞, we can also show that

R (u) (g) → 1, as u → ∞.
This lemma indicates that c(g) = 1, i.e. ḡ(X) and H are "asymptotically calibrated", is a necessary condition to have R(ḡ) = 0. In the above example with δ = 0, we have c(g 2 ) = 1 and R(g 2 ) = 1/2 ̸ = 0. This means that χ * (g 2 ) < 1. Note that c(g) = 1 implies that the constant χ * (g) simply corresponds to the aforementioned tail dependence coefficient χ(g), and so, the case c(g) = 1 simplifies the expression of the risk

R(g) = 1 - χ(g) 2 -χ(g) .
This equality tells us that any asymptotically calibrated classifier with χ(g) = 0 always produces a risk function R(g) = 1. Consequently, any asymptotically independent classifier is as uninformative as the two naive classifiers. A reasonable strategy will be to dismiss all asymptotically independent classifiers and find/construct new asymptotically dependent classifiers with positive χ(g). But, finding asymptotically dependent classifiers can be complex in practice, and in addition, in some not so exotic setups, this is not always possible. To see this, we consider the simple non-linear regression model in the following lemma.

Hence, no classifier can outperform naive classifiers for this regression model.

Note that even if the forecaster knows exactly the function f (.) and has drawn from the explanatory X, the "ideal" classifier g(x) =f (x) will perform badly, i.e. R(f ) = 1. In addition, the classical trick of using ranks to avoid the problem of marginals discrepancy cannot be applied here. For example, suppose that H is unit Fréchet distributed, then transforming the marginals of X into unit Fréchet random variables, say into X, does not remove the issue as the unobserved noise N has still have heavier tails than f ( X) = f (X) for some function f (.). So, a finer risk measure is needed that is able to distinguish different classifiers in case of asymptotic independence.

Conditional risk and hidden regular variation

The choice of the conditioning set in Lemma 3.4 brings new possibilities to construct finer risk measures for extremes than R (u) . To do so, we opted for the sets: A 1 = {g(X) > u} and B 1 = {H > u}, (when ε = 1) and A 0 = {g(X) > 0} and B 0 = {H > 0} (when ε = 0). But, changing the size of the sets A ε and B ε to make them closer to A 1 and B 1 , will increase the conditional probabilities P(B 1 | A 1 ∩B ε ) and P(A 1 | A ε ∩B 1 ). A simple choice when ε ̸ = 0 or ε ̸ = 1 is to set A ε = {g(X) > εu} and B ε = {H > εu} with ε > 0. This modelling strategy is at the core of hidden regular variation and asymptotic independent models. More precisely, we first need to fix marginal features. We assume that both g(X) and H possess regularly varying tails with indices α g > 0 and α H > 0, respectively. This means that for any ε ∈ (0, 1),

lim u→∞ P(g(X) > u | g(X) > εu) = ε αg and lim u→∞ P(H > u | H > εu) = ε α H .
These limits have to be understood with respect to Equation (3.10), i.e. the terms P(A 1 | A ε ) and P(B 1 | B ε ). To apply (3.10), the mixing condition (3.9) needs to be satisfied. To do so, we opt for an extended version of the framework of [START_REF] Ramos | A new class of models for bivariate joint tails[END_REF], i.e. (3.13) where η ∈ (0, 1] indicates the rate of decay of the joint survival function and L(•, •) is bivariate slowly varying function, i.e. there exists a limit function ℓ : (0, ∞) × (0, ∞) → (0, ∞) defined as ℓ(s, t) = lim u→∞ L(us, ut) L(u, u) , s, t > 0 and satisfying ℓ(cs, ct) = ℓ(s, t) for all c, s, t > 0. The parameter η measures the dependence strength. The case η = 1 corresponds to the asymptotic dependence while η < 1 to the asymptotic independence case. In particular, if η = .5, then independence appears in the extremes. If .5 < η < 1 (0 < η < .5) the extremal dependence is said to be positively (negatively) associated. Now, noticing that (3.13) corresponds to the mixing condition (3.9), we can apply Lemma 3.4, see Appendix A.1 for a proof. Proposition 3.7. Under the Ramos and Ledford model defined by (3.13), the following risk function

P[g(X) > u, H > v] = L(u, v)(u -αg v -α H ) 1/2η ,
R ε (g) := lim u→∞ P(g(X; u) ̸ = Y (u) | Y (εu) = g(X; εu) = 1) P(Y (u) = 1 or g(X; u) = 1 | Y (εu) = g(X; εu) = 1)
, (3.14) which will henceforth also be called conditional risk, can be expressed as

R ε (g) = 1 - 1 ℓ(ε, 1)ε -αg/2η + ℓ(1, ε)ε -α H /2η -1
, for any ε ∈ [0, 1).

Note that η ∈ (0, 1] takes a similar role as χ in the case of the unconditional risk R.

For fixed ε ∈ [0, 1), the risk function R ε (g) decreases with increasing η. So, given all parameters are fixed but η, the forecaster should aim at maximising η. In practice, two forecasters, say g 1 and g 2 , may produce different ℓ(., .) and α g . Consequently, the minimisation of R ε (g) can also depend, besides η, on other parameters.

Risk function inference

Concerning the estimation of R ε (g) defined by (3.14), the empirical estimator can be easily computed from the sample (X i , H i ) i=1,...,n . The following proposition describes the asymptotic property of such an estimator. Proposition 3.8. Assume that the risk function R ε (g) defined by (3.14)

exists for a sequence of u n → ∞ such that np g,ε (u n ) → ∞ with p g,ε (u n ) := P(max{g(X; u n ), Y (un) } = 1, H > εu n , g(X; εu n ) = 1). If lim n→∞ np g,ε (u n ) P(g(X; u n ) ̸ = Y (un) , H > εu n , g(X; εu n ) = 1) p g,ε (u n ) -R ε (g) = 0,
then the empirical estimator based on a sample (g(X i ; u), {H i > u}) i=1,...,n and defined by .15) converges in distribution in the following way

R n,ε (g) = n i=1 1{g(X i ; u n ) ̸ = Y (un) i , H i > εu n , g(X i ; εu n ) = 1} n i=1 1{max{g(X i ; u n ), Y (un) i } = 1, H i > εu n , g(X i ; εu n ) = 1} . ( 3
np g,ε (u n ) R n,ε (g) -R ε (g) n→∞ -→ N (0, R ε (g)(1 -R ε (g))) .
3.4. Simulations

Simulations

A simple linear setup

Our main simulated example is based on a simple linear regression model but with the feature that the explanatory variables X do not have the same tail behaviour and the noise is regularly varying, see Lemma 3.6. More precisely, the multivariate vector X is defined as follows

         X 1 ∼ Pareto(3) X 2 ∼ Pareto(2) X 3 ∼ Exp(1) X 4 ∼ Exp(2) , (3.16)
where all X i are independent with X 1 and X 2 Pareto distributed with respective tail index 2 and 3, and X 3 and X 4 exponentially distributed with respective scale parameters 1 and 2. The variable of interest H is simply a linear transform of X 1 tainted by an additive noise

H d = X 1 + N, (3.17)
where N ∼ Pareto( 2) represents an independent noise with heavier tail than X 1 . So, given a sample ({X j,i } 1≤j≤4 , H i ) i=1,...,n with n = 10000, our goal is to compare different classifiers in terms of predicting extreme occurrences, here defined as {H > u} with u equal to the 97th percentile of H. In this simulation setup, it is clear from (3.17) that all variables but X 1 are useless to explain H. In addition, Lemma 3.6 tells us that the relevant information contained in the variable X 1 is hidden by the heavier noise N , i.e. we are in the case of asymptotic independence. An example of such simulation is given in Figure 3.3. The left panel displays a scatter plot between H (left axis) and g(X) = X 1 (right axis). As expected, no sign of asymptotic dependence can be found in the upper corner. In the right panel, we remove the mass along the axis (grey points) by conditioning on the joint event

A ε ∩ B ε = {g(X) > εu} ∩ {H > εu} with ε = .
7, see all dark points. The right panel zooms on these black points and highlights a clear dependence between H and X 1 that was hidden by the heavier noise N in H = X 1 + N .

In practice, we do not know the optimal choice for g(.) and we need to introduce different classifiers and compare them.

Classifiers descriptions

Table 3.1 below provides the list of classifiers that we compare with our metric (3.15). This list contains some of the most standard classifiers found in the literature (see, e.g. [START_REF] Hastie | The Elements of Statistical Learning: Data Mining, Inference, and Prediction[END_REF]: logistic regression (Logistic), decision tree (Tree), random forest (RF) and support vector machine (SVM).

Method Main features

Linear classifier Simple binary classifier, parameters estimation based on the minimisation of the risk function over the set of contributing variables, theoretical value of R(g θ ) inferred from spectral decomposition.

Logistic regression

Parametric linear model with a lasso penalty, coefficients of less contributing variables are set to zero.

Decision trees

Easy to interpret, gives relative importance of each variables, learns simple decision rules inferred from the input.

Random forests

Builds multiple decision trees combined by majority vote, better predictive power than decision trees.

Support vector machines Finds the best hyperplane to separate two overlapping classes, generally performs better than the other classifiers.

Table 3.1: Summary and key features of the different classifiers studied. See for example [START_REF] Hastie | The Elements of Statistical Learning: Data Mining, Inference, and Prediction[END_REF] for a comprehensive review of the last four classification methods.

Except for the linear classifier, we apply them with their built-in cost function that is not necessarily fine-tuned to forecast extremes. This is not an issue because our main goal is to compare existing forecasters, and not to create new ones (see, e.g. Jalalzai et al., 2018, for such developments). Still, to fix a baseline in terms of performance, the linear classifier defined as

g θ (X; u) = +1, θ ⊤ X > u, -1, θ ⊤ X ≤ u, θ ∈ [0, ∞) d .
should be optimal for the linear model (3.17), especially if the regression parameters are estimated by minimising our cost function (3.14). In such a context, we expect the linear classifier to be the best. In Appendix 3.7, Proposition 3.9 provides the condition of the consistency of the estimator θn,un based on minimising Rn (g θ (•; u n )) under a regularly varying framework.

The binary outputs from the decision tree classifier are explained in Figure 3.4. The light blue and light green regions represent the set of points that are well predicted by the classifier. On the contrary, wrongly classified points belong to the light yellow and light red regions: either an extreme is predicted when it is not (light red region), or an extreme event is missed (light yellow region). The difference between the left and right panels corresponds to the training set based The results shown are the predicted values made upon the testing set whose sample size is equal to 3000 (30% of the data). on either {H > u} or on {H > ε u }, i.e. mass removed in the latter case, see also (3.15).

Implementation and results

We split our simulated data set in two: 70% for a training part, over which we train our different classifiers to get good predictive power; 30% for a testing part, which we use to estimate the risks R 0 (g) and R ε (g). Note that each algorithm has the same inputs, in particular the same binary sequence describing the events {H > u} with u set to be equal to the 97th percentile of H. This cross-validation procedure has been repeated 50 times. The sample used to compute our risk function is based on the bivariate binary vector (g({X j,i } 1≤j≤4 ; u), {H i > u}) i=1,...,n where the output of the classifier g is binary. In addition, the binary outputs of the classifier g are obtained under the threshold u and the threshold εu, so the training part has to be performed twice (once for each threshold). Then, the empirical risk estimator defined by (3.15) can be computed. Figure 3.5 shows the sensitivity of the classifier ranking with respect to the value of ε.

As expected from Lemma 3.6, the top-left panel, that corresponds to the case ε = 0, clearly indicates that our five classifiers cannot outperform naive classifiers as all classifiers have a risk near to one, the worst possible value. To start dis- criminating classifiers, we need to remove the masses along the axes by setting a value for ε. As we increase ε, the size of sets needed to compute (3.15) becomes smaller, see the values of n ε in the legend of each panel. Hence, the blue box plots becomes wider as ε increases: a classical bias-variance trade-off. As we know from (3.17) that the true generative process is linear, ε = 0.4 appears as a reasonable value to balance the bias-variance trade-off. More importantly, the overall ranking is not sensitive to the values of ε > 0. In all cases, our linear classifier tailored to handle linear asymptotic independence cases outperforms all the other classifiers. Among the other classifiers, decision tree appears to be the best, but it is still far from the optimal linear solution. Other simulations concerning the regular variation case are available upon request. We now apply our assessment approach to summer daily river discharges (measured in m 3 /s) at 31 stations spread over the upper Danube basin, see Figure 3.6, and recorded over the time period 1960-2010 in June, July and August. These observations have been studied by the EVT community (see, e.g. [START_REF] Asadi | Extremes on river networks[END_REF][START_REF] Mhalla | Causal mechanism of extreme river discharges in the upper danube basin network[END_REF][START_REF] Gnecco | Causal discovery in heavy-tailed models[END_REF]. This dataset was made available by the Bavarian Environmental Agency (http://www.gkd.bayern.de). To remove temporal clustering in extreme river discharges, [START_REF] Mhalla | Causal mechanism of extreme river discharges in the upper danube basin network[END_REF] in their Section 5 implemented a declustering step. Each station then contains n = 428 observations that we will consider temporally independent. In order to reduce the large discrepancy in terms of discharges magnitude among stations, we force the starting value of all 31 time series to equal zero by subtracting to each station its minimum.

Danube river discharges

Then, we re-normalise each time series by its range (i.e. the difference between the maximum and the minimum of each time series). These post processing treatments are useful to display and interpret the data at hand and do not impact the classifiers performance. Although all 31 station recordings are available, we can artificially remove one station and try to predict its values from a given subset of other weather stations.

In this section, we remove station 1 (downstream) and try to predict its values from only stations 23 and 24, which are indirect tributaries to the main river flow. So, this setup is complex4 for two reasons. First, station 1, as a downstream point that accumulates all discharges, has a much heavier tail than the two tributaries. Second, it is difficult to determine if we are in the asymptotic dependent or independent case, see Figure 3.7 that displays the scatter plot between the hidden station (y-axis with station 1) and the two tributaries (x-axis, stations 23 and 24). In this graph, the threshold u is taken to be equal to the 85th quantile of X 1 and we choose ε = 0.6. Figure 3.8 summarises our findings. Removing the mass on the axes when thresholding by ε u implies that only 190 points remain from the original length of 428 data points per station. This can explain why, looking at Figure 3.8, the uncertainty in the risk estimate increases when considering R ε (.) (blue boxes) instead of R 0 (.) (red boxes). Figure 3.8: Estimation of R ε (g) for five different classifiers (cross-validation with 50 repetitions, 70% train and 30% test), threshold is the 0.85 quantile of H. In red (left) are the estimates when ε = 0 and in blue (right) for ε = 0.6. The length of the testing set is equal to 129, this leads to around 60 points such that H > ε u and nearly 20 points such that H > u.

Unlike the simulation example in Section 3.4.1, it is not clear to assess whether our river discharges analysis of our three selected weather stations belongs to the framework of asymptotic independence or not. Still, it is reassuring that the ranking of the classifiers in Figure 3.8 appears to be insensitive to the values of R 0 or R ε . i.e. whether the data are asymptotically dependent or not. This hints that, among all the classifiers, the logistic regression with lasso penalty seems to perform better than the four other classification methods. This ranking of classifiers is specific to this particular example. No general conclusions about lasso techniques for extremes should be drawn.

Besides this river example, we advocate practitioners to compute risk functions that can both handle the asymptotic dependent and independence cases. This also complements the recent tools used to discriminate between the two cases (see, e.g. [START_REF] Ahmed | Recognizing a spatial extreme dependence structure: A deep learning approach[END_REF]. In addition, the linear classifier could provide a simple benchmark with well understood properties with respect to R 0 , see Proposition 3.9.

Supplementary Materials

A R package is available on GitHub that implements the empirical estimation of the risk function developed in this paper (https://github.com/jlegrand35/Extr emesBinaryClassifier) and can be used either to reproduce the results of the conducted classifier comparisons or to perform new comparisons using other binary classifiers. The data used in the application are available in the R package graphicalExtremes [START_REF] Engelke | graphicalExtremes: Statistical Methodology for Graphical Extreme Value Models[END_REF].

Acknowledgement

Supplements to Paper II

This section was originally part of the appendix of [START_REF] Legrand | Evaluation of binary classifiers for asymptotically dependent and independent extremes[END_REF]. We reproduce it below and give further details.

In the following, the focus is on linear classifiers. We first give its definition and basic properties, then provide a way to find such optimal classifiers with Proposition 3.9. It turns out that it is sufficient to consider only the components contributing to the extremes. This leads to the notion of sparsity and from this we Here, we additionally assume that P(∥Γ∥ ∞ > 0) > 0 and P(Ω > 0) > 0. In this setup, ∥X∥ ∞ and H are tail equivalent in the sense that

lim u→∞ P(∥X∥ ∞ > u) P(H > u) = lim u→∞ P(∥X∥ ∞ > u | ∥(X, H)∥ ∞ > u) P(H > u | ∥(X, H)∥ ∞ > u) = P(P • ∥Γ∥ ∞ > 1) P(P • Ω > 1) = E(∥Γ∥ α ∞ ) E(Ω α ) ∈ (0, ∞),
This is the minimal requirement on the link between the covariates X and the unobserved extremes of H essentially saying that at least one component of X is tail-equivalent to H. It is important to highlight that we do not exclude the case that Γ i = 0 a.s. for some i ∈ {1, . . . , d} which means that X i possesses a lighter tail than H. This property can be read off from the quantity

c i = lim u→∞ P(X i > u) P(H > u) = E(Γ α i ) E(Ω α ) . (3.18)
Thus, Γ i = 0 a.s. if and only if c i = 0. By including this case, we therefore admit that most of the components of X may not contribute to the extremes of the vector H. This feature is essential when the question of sparsity will be addressed.

Under these conditions, we obtain that, for all θ ∈ [0, ∞) d , the classifier g θ is an extremal classifier as .19) Equation (3.19) implies that the function θ → R(g θ ) is well-defined and continuous on [0, ∞) d . Its value does not depend on those components θ i for which Γ i = 0 a.s., which is equivalent to c i = 0 as discussed above. Thus, in the following, we will consider this function only on the parameter set

R(g θ ) = lim u→∞ P(H > u or θ ⊤ X > u) -1 P[max{θ ⊤ X, H} > u] -P[min{θ ⊤ X, H} > u] = E max{θ ⊤ Γ, Ω} α -E min{θ ⊤ Γ, Ω} α E (max{θ ⊤ Γ, Ω} α ) = 1 - E min{θ ⊤ Γ, Ω} α E (max{θ ⊤ Γ, Ω} α ) ∈ [0, 1]. ( 3
C = {θ ∈ [0, ∞) d : θ i = 0 for all i s.t. c i = 0},
containing all the relevant information -here, note that, in practice, identifying the components i ∈ {1, . . . , d} such that c i = 0 a.s., from a given data set is a necessary step for the correct specification of the set C.

From the consideration in the introduction, it can be easily seen that R(g 0 ) = 1 -the case θ = 0 corresponds to the trivial always optimistic classifier. Furthermore, denoting the set of indices j with c j > 0 by J, we can see that

R(g θ ) ≥ 1 - E(Ω α ) ∥θ∥ α ∞ • min j∈J E(Γ α j ) → 1 (3.20)
as ∥θ∥ ∞ → ∞. By the continuity of θ → R(g θ ), we obtain that the function attains a global minimum on the domain C. Proposition 3.9. Additionally to the assumptions above on the joint distribution of (X, H) with α > 1, assume that there exists a function a(u) with a(u) → 0 as u → ∞ such that is consistent, i.e. θn,un → p θ * .

P(u -1 X ∈ A | ∥(X, H)∥ ∞ > u) ≤ (1 + a(u)) P(P Γ ∈ A) (3.21) for all A ⊂ [0, ∞). Furthermore, let u n → ∞ and n P(H > u n ) → ∞ such that, for every compact subset K ⊂ C, sup θ∈K n P(H > u n ) P(g θ (X; u n ) ̸ = Y un ) P(H > u n ) - E max{θ ⊤ Γ, Ω} α -E min{θ ⊤ Γ, Ω} α E(Ω α ) = 0 (3.22) and sup θ∈K n P(H > u n ) P(max{g θ (X; u n ), Y un } = 1) P(H > u n ) - E max{θ ⊤ Γ, Ω} α E(Ω α ) = 0. ( 3
Given the set C, this result provides a strategy to find the optimal θ, i.e., the best linear classifier. Determining the set C requires the identification of the relevant features, i.e. the index set J such that c j > 0 if and only if j ∈ J. This is discussed in more detail in the following subsection.

Feature Selection

The notion of sparsity quickly comes into play when doing classification. This is all the more true when one is only interested in the extremes. Among the whole data set, only a small proportion will truly contribute to the extremal behaviour of the variable of interest. Here, we develop a method to identify the informative signals in terms of extremes among a large data set, assuming that (X, H) is jointly regularly varying. For a comprehensive review of existing methods on sparsity and multivariate extremes we highly recommend the work of [START_REF] Engelke | Sparse structures for multivariate extremes[END_REF].

As we have seen above, for linear classifiers, all the relevant features X i necessarily satisfy c i > 0. Thus, feature selection can be based on estimation of the c i which can be done according to the following proposition. Proposition 3.10. Assume that

c i = lim u→∞ P(X i > u) P(H > u) exists. If u n → ∞ and n P(H > u n ) → ∞, then n j=1 1{X j,i > u n } n j=1 1{H j > u n } → p c i . If, additionally, n P(H > u n ) P(X i > u n ) P(H > u n ) -c i → 0 and χ * i = lim u→∞ P(X i > u | H > u) ∈ [0, 1]
exists, then, we have

n P(H > u n ) n j=1 1{X j,i > u n } n j=1 1{H j > u n } -c i → N (0, c i • [1 -2χ * i + c i ]) .
The above Proposition 3.10 helps us to discriminate between informative coordinates and non-informative ones in terms of extremes.

River network

We return to the Danube river network. Contrary to what was done in Section 3.5, here the data are pre-processed differently: we only subtract the minimum to each station, no re-normalisation is performed. Thus some stations have a heavier tail than others.

Before applying the methodology outlined above, an a priori analysis is performed to check if the data are well regularly varying. From Table 3.2, we can not reject the hypothesis that the river discharges at the stations considered are heavy-tailed. Note that we show the results only for Station 1 (the station of interest in this study), Stations 2, 13 and 30 (the stations that appear to contribute the most to the extremes of Station 1, see below), and Stations 23 and 24 (the stations considered in Section 3.5).

This section deals with a simpler case than the application in Section 3.5. Here an application could be the following: we want to know which stations should continue to be maintained to prevent extreme floods and maybe some stations are not necessary.

As before, the goal is to predict the extreme events at Station 1, denoted X 1 , where an extreme event is defined as an event exceeding the 85th quantile of 3.2: Estimates of the shape parameter when fitting a GPD. The threshold was selected using mean residual life plot and parameter stability plots (see, e.g. [START_REF] Coles | An Introduction to Statistical Modeling of Extreme Values[END_REF]). Estimation is performed using the R package ismev [START_REF] Heffernan | ismev: An Introduction to Statistical Modeling of Extreme Values[END_REF].

X 1 . Unlike the study of Section 3.5, we assume that the whole set of remaining stations is available. In this case, strong dependencies among Station 1 and other stations can be observed. Therefore, the main issue is to identify and select these stations following the procedure presented above.

The stations that may not contribute to the extremes of X 1 are identified through the estimation of the coefficients c i . The estimation of the set C on all the data is presented in Table 3.3. Among the 30 stations, only three stations are relevant: Stations 2, 13 and 30. Looking at Figure 3.6, these stations correspond to the stations closest to X 1 . Figure 3.9 shows the scatter plots between these stations and Station 1, reflecting strong dependencies between the variables. c i 0 0 0 0 0 0 0 0 0 0 0 0 0.02 0 Table 3.3: Empirical estimates of ĉi (as defined in Proposition 3.10) for each sta- tion. The values different from zero are highlighted in red.

Once the contributing variables have been identified, we compare the performance of several classifiers, on the one hand keeping all the data and on the other hand keeping only the informative stations. Since there is a strong dependence between the data, we assume that it is sufficient here to look at the risk R 0 . Comparison results are shown in Figure 3.10. The comparison is performed among the same classes of classifier as in Section 3.5, see Table 3.1 for a summary.

By definition of the linear classifier, the estimation is already done by keeping only the informative variables, which is why the estimates are identical for this specific classifier. As for the other classifiers, we see some improvements when keeping only the informative variables: the risk estimates are slightly smaller. This means that even if we remove a lot of information by going from 30 explanatory variables to 3, these 3 remaining stations contain all the information in terms of extremes of Station 1. However, caution must be exercised here, as no theoretical Figure 3.9: Correlation plot between X 1 and the stations considered to contribute to the extremes of X 1 according to Table 3.3 (ie for which c i ̸ = 0). The blue dotted lines represent the threshold u defined by the 85th percentile of X 1 .

Predictive power of the linear classifier

We provide some additional remarks on the special case of the linear classifier.

Of the foregoing, only Stations 2, 13 and 30 are relevant for the estimation of the optimal θ (see Proposition 3.9). The left-hand side of Figure 3.11 depicts the distribution estimates of θ2 , θ13 and θ30 . The estimation is performed on 70% (training set) of the data drawn 50 times at random. From Figure 3.11, X 30 has less effect than X 2 and X 13 . Looking at the river map (Fig. 3.6), Station 30 is the furthest station from Station 1 out of the three.

Then one of the 50 estimates is retained (at random, but sensitivity test showed that it does not impact our results). Let's denote it θ = ( θ2 , θ13 , θ30 ). And the predicted output of our linear classifier is computed through θ⊤ X test . The vector X test corresponds to the river discharges at Station 2, 13 and 30 from the remaining 30% data (that has not been used for the estimation of θ). In the right-hand side of Figure 3.11, the predicted output θ⊤ X test is plotted against the observed river discharges at Station 1. From this, we find that the linear classifier predicts the extremes quite well (upper right quadrant), but also, it predicts well the river discharges that are not extreme (i.e. outside the extreme region). This is most likely due to the specificity of our data, which are highly correlated. Figure 3.11: (left) Estimates of the optimal θ i of the three contributing variables X 2 , X 13 and X 30 (see Table 3.3), with 50 replications of the optimisation procedure. (right) Predicted versus observed river discharges at Station 1. The predicted values are defined as θ⊤ X, for one estimated triplet θ among the 50 estimations.

The blue delimited upper region represents the region where both coordinates are above the 0.85 quantile of X 1 . The red dotted line depicts the first diagonal.

Possible extension to hidden regular variation, a simulation study

Up to now, Section 3.7 only addressed the asymptotically dependent case. However, some similar results could be derived in the case of asymptotic independence.

For that, we consider a modified version of the c i defined in (3.18). Let ε ∈ [0, 1)

and define the following:

c i,ε = lim u→∞ P(X i > u | X i > εu) P(H > u | H > εu)
.

Then different cases arise according to the value of ε.

Case ε = 0: We get c i,0 ∈ (0, ∞) if, and only if, X i and H are tail equivalent. We get c i,0 = 0 if the tail of X i is less heavy than the tail of H. This corresponds to the previous coefficient defined in (3.18).

Case ε > 0:

If the tail function of H is regularly varying with index α H and the tail function of X i is regularly varying with index α i , we obtain

c i,ε = ε α i -α H .
Thus, c i,ε = 0 only if X i possesses light tails (assuming that H is heavy-tailed).

The following Lemma and the inference properties of c i,ε from Proposition 3.12

are still an open question, but from the simulation results shown hereinafter, we are confident that this holds true.

Lemma 3.11. Assume that the joint distribution of (g(X), H) satisfies the Ramos and Ledford model

P[g(X) > u, H > v] = L(u, v)(F g (u)v -α H ) -1/2η , u, v > 0,
where F g denotes the tail distribution of ḡ(X). Then, the following holds true:

1. If

lim u→∞ P(g(X) > u | g(X) > εu) P(H > u | H > εu) = 0,
we have that R ε (g) = 1 for all ε ∈ (0, 1].

2. If R ε (g) < 1 and c i,ε = 0, then, for all c > 0, we have that

R ε (g) = R ε (h i (c)),
where h i (c) = g(X) + cX i .

True value n = 1000 n x = 300 q = 70 n = 5000 n x = 300 q = 94 n = 10000 n x = 300 q = 97 n = 20000 n x = 300 q = 98.5 n = 1e5 n x = 300 q = 99.7 n and 95% confidence interval coverage percentages computed among the 1000 replicates. The number of exceedances is denoted by n x and the percentile chosen for the threshold is given by q. The value ε is equal to 0.7.

c
The simulation results are presented in Table 3.4. It appears that the sample size does not impact the estimates. We also see that X 1 and X 2 are detected as contributing to the extremes of H and that we could exclude X 3 and X 4 from our classification task but, in that case, confidence intervals appear to be difficult to find.

Summary of Paper II

To sum up

• Paper II addressed the specific issue of comparing binary classifiers in the case of rare event prediction.

• Asymptotic dependent models: considering our first risk function R(•) (which counts the number of mistakes), we have a benchmark with the naive classifiers and we can rank different classifiers.

• Asymptotic independent models: for regression models of the form H = f (X)+N , with tail of N heavier than tail of f (X), any classifier will be as good as the naive classifiers.

→ To tackle this issues, development of risk functions adapted to extremal classifiers.

Summary of Paper II

→ Flexible framework that allows to address both asymptotic dependence and asymptotic independence.

→ Tool to identify the explanatory variables that contribute the most to extremal behaviour, with promising results in the case of hidden regular variation.

̸ → Assumption of multivariate regular variation implies to work with heavy-tailed data.

̸ → Transforming the data into identical margins would, in particular, defeat the optimisation scheme of the linear classifier.

Chapter 4

Some directions for future work

In three years, many avenues are explored, offering various topics for future research. Hereinafter we list some potential extensions of the work presented in this thesis.

Extension of the simulation model to the multivariate case

Recall that in Chapter 2, we considered the stochastic representation of a standard GP vector through

Z = E + T -max(T ), (4.1)
with E a unit exponential variable and T a multivariate vector. For a bivariate vector Z = (Z 1 , Z 2 ), defining ∆ = Z 1 -Z 2 allowed to rewrite the representation (4.1) in a simple form, which was then easy to use for nonparametric simulations:

Z i = E ± ∆,
given the sign of ∆.

Now consider for example that max(T ) = T 1 . Then the above can be written as follows

Z 1 = E and Z 2 = E -∆.
This gives hints for an extension to the multivariate case. Indeed, if we consider Z = (Z 1 , . . . , Z d ), d > 2, a first step would be to determine the index j such that max(T ) = T j , this could be achieved through a random sampling between {1, . . . , d}. Then defining ∆ i = T j -T i for i ̸ = j, this leads to Z j = E and Z i = E -∆ i , for i ̸ = j.

(4.2)

Simulating independent and identically distributed (iid) unit exponential variables E is costless. As for the bivariate case, the main point is the simulation of the ∆ i . An avenue to explore is to bootstrap observations from the (d-1)-dimensional 123 CHAPTER 4. SOME DIRECTIONS FOR FUTURE WORK set of points ∆ -j := (∆ 1 , . . . , ∆ j-1 , ∆ j+1 , . . . , ∆ d ). This could be achieved by resampling on the indexes: assuming that we have n iid replicates of ∆ -j , random draws with replacements in {1, . . . , n} leads to m bootstrap samples From this, we could then perform non-parametric bootstrap MGP simulations in dimension greater than 2.

Perhaps the most challenging point for our application to extreme wave heights is the multivariate extension of the marginal regression models. Indeed, if we assume that instead of having only one offshore site, two or three sites are considered (represented in red in left-hand side of Figure 4.1). From an application perspective, this could certainly improve the coastal predictions as we would include more information in the covariates. But up to now, it is unclear how to define the conditioning part as in Equation (2.2). For the coastal marginal model, a possibility could be of the following )). This certainly leaves rooms for future developments.

H c -v c | H 1 o , T 1 p , D 1 p , H
In terms of applications of such multivariate simulation model, in addition to the above-mentioned one, which would allow for more refined simulations, another interesting application would be for the recovery of historical data. A possible approach is schematised in the right-hand side of Figure 4.1 and would be to consider two well documented sites such as buoys (represented by the two red dots in Figure 4.1), and to simulate data between this two measurement sites (represented by the blue dot). This idea was suggested by Jérémy Rohmer during a visit at BRGM, and would deeply increase the understanding of past extreme wave events, for which in-situ measurements are rather sparse.

Comparison of coastal extreme wave warning models

For coastal risk assessment, it is crucial to have efficient extreme sea state warning systems. From a simulation methods of extreme H s as developed in this manuscript, one could then derive a warning (or forecast) model by looking at the probability of the event {Y u = +1}, for a given high threshold u, where Y u is defined as follows In order to have a "good" warning model, one needs to ensure that the forecast does not produce to many false alarms or miss. In this context, the risk function developed in Paper II appears to be entirely relevant. The performances of several prediction models could be assess on historical/hindcast data in order to select the most efficient in terms of extreme wave events forecast. Such future developments are definitively something worth trying.

Y u = +1, if H s > u, -1, otherwise.

Application to climate projections

Finally, regarding adaptation strategies for future climate conditions, it is of fundamental importance to be able to derive information on the changes in the occurrence of extreme wave events between current and future climates, but also on the variability between different climate projection models. As discussed in the introduction of Chapter 2, such questions have been addressed on a global scale (e.g. [START_REF] Aarnes | Projected changes in significant wave height toward the end of the 21st century: Northeast atlantic[END_REF], but local variability has been less explored.

Assuming that the relation between the offshore point and the local point remains identical, some future extension of the presented work could then be the simulation of extreme coastal H s given large scale wave climate projections. This is somewhat linked to the idea of statistical downscaling, where it is assumed that there exists a transformation between the large scale and the local scale, and that this transformation will remain valid in the future climate. Such methods have been applied for instance in [START_REF] Bechler | A spatial hybrid approach for downscaling of extreme precipitation fields[END_REF] in spatial context to downscale extreme rainfalls, or in [START_REF] Towe | Statistical downscaling for future extreme wave heights in the North Sea[END_REF] to produce future extreme H s from downscaled wind fields. We can write that

P(A 1 ∪ B 1 | A ε ∩ B ε ) = P((A 1 ∪ B 1 ) ∩ (A ε ∩ B ε )) P(A ε ∩ B ε ) , = P(A 1 ∩ B ε ) + P(A ε ∩ B 1 ) -P(A 1 ∩ B 1 ) P(A ε ∩ B ε )
.

In the same way, we have

P(A 1 △B 1 |A ε ∩ B ε ) = P(A 1 ∩ B ε ) + P(A ε ∩ B 1 ) -2 P(A 1 ∩ B 1 ) P(A ε ∩ B ε )
.

Hence, we deduce that

R(A ε , B ε ) = P(A 1 △B 1 |A ε ∩ B ε ) P(A 1 ∪ B 1 |A ε ∩ B ε ) = 1 - P(A 1 ∩ B ε ) P(A 1 ∩ B 1 ) + P(A ε ∩ B 1 ) P(A 1 ∩ B 1 ) -1 -1
.

The expression given by (3.8) follows.

Item (a) of the lemma is based on the following inequality

P(U | V ) ≥ P(U | W ), if the sets U , V and W satisfy U ⊂ V ⊂ W .
For item (b), note that

R(A ε , B ε ) ≤ R(A ′ ε , B ε ) ⇐⇒ (1 -R(A ε , B ε )) -1 ≤ (1 -R(A ′ ε , B ε )) -1 ⇐⇒ P(A 1 ∩ B ε ) P(A 1 ∩ B 1 ) + P(A ε ∩ B 1 ) P(A 1 ∩ B 1 ) ≤ P(A ′ 1 ∩ B ε ) P(A ′ 1 ∩ B 1 ) + P(A ′ ε ∩ B 1 ) P(A ′ 1 ∩ B 1 )
, ⇐⇒ P(A 1 ∩ B ε ) + P(A ε ∩ B 1 ) ≤ P(A 1 ∩ B 1 )

P(A ′ 1 ∩ B 1 ) [P(A ′ 1 ∩ B ε ) + P(A ′ ε ∩ B 1 )].
A.1. Proofs

The Ramos and Ledford model corresponds to the special of case of (3.9) Then, from (3.10)

P[g(X) > u, H > v] = L(u, v)(u -αg v -α H ) 1/2η = P(A ε ∩ B ε ′ ) = c ε,
R(A ε , B ε ) = 1 - c ε,1 c 1,1 (P(A 1 |A ε )) -a + c 1,ε c 1,1 (P(B 1 |B ε )) -b -1 -1
.

Letting u gets large provides the required result. □.

A.1.4 Proof of Proposition 3.8:

As we assume that R ε (g) exists (for some ε ≥ 0), for u n → ∞ such that np g,ε (u n ) → ∞, we obtain that }=1,H i >εun,g(X i ;εun)=1} npg,ε(un)   -P(g(X,un)̸ =Y (un) ,H>εun,g(X;εun)=1) pg,ε(un)

1 n→∞ -→ N 0 0 , R ε (g) R ε (g) R ε (g) 1
Provided that the bias is negligible, i.e. lim n→∞ np g,ε (u n ) P(g(X, u n ) ̸ = Y (un) , H > εu n , g(X; εu n ) = 1) p g,ε (u n ) -R ε (g) = 0, the Delta method yields

np g,ε (u n ) R n,ε (g) -R ε (g) n→∞ -→ N (0, R ε (g)(1 -R ε (g)))
.

□

Titre:

  Simulation et évaluation de modèles d'extrêmes multivariés pour des données environnementales Mots clés: Extrêmes multivariés; Simulation d'extrêmes cooccurrents; Variation régulière multivariée et cachée; Evaluation de classifieurs extrêmes; Hauteur significative de vagues; Débits de rivière Résumé:

Figure 1 :

 1 Figure 1: Carte représentant le point au large (en rouge) et le point près des côtes (point vert).

Figure 2 :

 2 Figure 2: H s à la côte en fonction des H s au large, la zone hachurée rouge correspond aux dépassements bivariés des hauteurs significatives de vagues.

Figure 3 :

 3 Figure 3: Carte du bassin versant du Danube, montrant les 31 stations de jaugeage le long du Danube et de ses affluents.

Figure 4 :

 4 Figure 4: Exemple de prédictions binaires obtenues avec un arbre de décision, contre les observations H. La ligne horizontale bleue représente le seuil critique u, la ligne verticale sépare les deux classes prédites -1 et +1.

Figure 1 . 1 :

 11 Figure 1.1: Densities of the three generalised extreme value distributions. For the Fréchet distribution ξ = 0.5 and for the Weibull distribution ξ = -0.5. The shape parameter σ is equal to 1 and the location parameter µ is equal to 0.

Figure 1 . 2 :

 12 Figure 1.2: Densities of the GP distribution for different values of the shape parameter ξ. For the Fréchet family ξ = 0.5 and for the Weibull family ξ = -0.5. The scale parameter σ is equal to 1.

Figure 1 . 3 :

 13 Figure 1.3: Upper tail distribution of significant wave heights (H s ), corresponding to the threshold exceedances above the 0.98 quantile of H s . GP density with shape parameter ξ = 0 (resp. ξ > 0) is superimposed as a red solid line (resp. blue dashed line).

Figure 1

 1 Figure 1.6: Illustration of asymptotic independence with a bivariate Gaussian vector on the original scale (left) and transformed to common Pareto scale (right).

Figure 1 . 7 :

 17 Figure 1.7: Simulations of n = 1000 points from the bivariate logistic extreme value model (1.12) for different values of α ∈ {0.2, 0.5, 0.99} (left to right).

  13) corresponds to the blue L-shaped region in the left-hand side of Figure 1.5 or the red L-shaped region in Figure 1.8.

Figure 1

 1 Figure 1.8: (left) Map illustrating the locations considered in Chapter 2. The red dot corresponds to the "offshore" point and the green dot corresponds to the "coastal" point. (right) Scatter plot of the significant wave heights at the two locations considered, more details on the data are provided in Chapter 2. The black dots belong to the region where the multivariate model will be fitted.

Figure 1 .

 1 Figure 1.10: Simulations of bivariate GP models from the representation (1.19) with (left) T following a Gaussian distribution with zero-mean and covariance matrix Σ = [(1, 0.4), (0.4, 1)] and (right) T = (T 1 , T 2 ) where T 1 and T 2 are independent variables with T 1 ∼ Gumbel(0, 1) and T 2 ∼ Gumbel(0, 4).

Figure 1 .

 1 Figure 1.11: Asymptotic dependence behaviour for the different application data sets considered in this thesis. (left) Scatter plot between the significant wave heights in the two locations of interest (see Chapter 2). (right) Scatter plot between the river discharges in two distant stations (see Chapter 3). The blue hatched region is just a visual indication of where both components are large simultaneously.

Figure 1 .

 1 Figure1.12: (left) Spectral measure H and (right) Pickands dependence function A for the logistic model defined in (1.12)[START_REF] Coles | An Introduction to Statistical Modeling of Extreme Values[END_REF], with α = 0.2 (red solid line), α = 0.5 (blue dashed line) and α = 0.99 (green dotted line). The spectral measures on the left panel have been normalised between 0 and 1 for visual purposes. On the right panel, the black dotted line represents the lower bound of A. Estimations are performed using the R packages evd[START_REF] Stephenson | evd: Extreme value distributions[END_REF] for H and copula[START_REF] Hofert | copula: Multivariate Dependence with Copulas[END_REF] for A.

Figure 1 .

 1 Figure 1.13: (left) Scatter plot of significant wave heights in two nearby locations off the coast of Brittany, more details are given in Chapter 2. (right) Plot of the corresponding estimated dependence measure χ for increasing values of quantile levels u with 95% point-wise confidence intervals (estimation is performed using the R package evd (Stephenson, 2002)).

Figure 1 .

 1 Figure 1.14: Asymptotic dependence comparison on Pareto scale between H s data (left) and a Gaussian vector (right).

  Figure 1.5.

Figure 2 .

 2 Figure 2.1: Spectral decomposition of a sea state as a sum of harmonic waves (taken from Holthuijsen (2007)).

Figure 2 . 2 :

 22 Figure 2.2: Example of wave spectrum at the coastal location considered in Paper I (see Figure 2.7). The corresponding date is 21/12/2021 04:00. The wave spectrum is provided by https://marc.ifremer.fr/.

  3) and 3-hourly temporal resolution from 01/01/1990 to 31/03/2020. • HOMERE (Accensi and Maisondieu, 2015): covers the Channel and the Bay of Biscay on an unstructured grid, refined close to the coast (see Figure 2.3). It takes into account a high resolution bathymetry and is forced by IOWAGA on the wet boundaries. The temporal resolution is one hour time step, data are available from 01/01/1994 to 30/04/2021.

Figure 2 . 3 :

 23 Figure 2.3: Grids of the hindcast databases considered. (left) Portion of the IOWAGA grid. (right) Total HOMERE grid.

Figure 2 . 4 :

 24 Figure 2.4: Wind field corresponding to the highest H s at the location of interest, depicted by the yellow dot. Arrows represent the wind direction and colours its intensity. Wind data are extracted from the ERA5 database.

Figure 2 . 5 :

 25 Figure 2.5: Measure of asymptotic dependence χ between wave data points from IOWAGA and the location of interest extracted from HOMERE (green triangle point).

Figure 2 .

 2 Figure 2.6: Densities of the extended GP distribution for different values of the shape parameter κ, with ξ = 0 and σ = 1.

  Summary of available data for each step of this study. A tick v (resp. a cross x) indicates the availability (resp. non-availability) of the data.

Figure 2 .

 2 Figure 2.7: (left) Portion of IOWAGA hindcast database grid, the red dot corresponds to the "offshore" point (data extracted from the IOWAGA database) and the green dot corresponds to the "coastal" point (data extracted from the HOMERE database). (right) Scatter plot of the coastal significant wave heights versus the offshore significant wave heights with the different thresholds considered. The black dots belong to the region where the multivariate model is fitted.

Figure 2 .

 2 Figure 2.8: (left) Empirical histogram for the coastal significant wave height threshold exceedances and (right) similarly for the offshore significant wave heights threshold exceedances, illustrating that fitting a generalised Pareto distribution (Coles, 2001) is not suitable for the coastal data. The coastal marginal threshold is defined by v c := min (H c ; H o > v o ).

Figure 2 .

 2 Figure 2.10: (left) Interpolated ratio of empirical extreme coastal significant wave heights H c and extreme offshore significant wave heights H o . The interpolated surface between the data points is performed using local polynomial interpolation of degree 2. (middle) Ratio of the predicted conditional expectations of extreme coastal significant wave heights H c and extreme offshore significant wave heights H o , conditionally to the offshore peak direction D p and the offshore peak period T p . (right) Absolute difference between the interpolated empirical ratio and the predicted ratio. On the three plots, observed data points are superimposed.

Figure 2 .

 2 Figure 2.11: Each panel line corresponds to one of the parametric model (a) to (e), and shows from left to right: (1) Scatter plot of simulated data with the parametric model with sample size n = 10000 (black dots) and sampled data from one simulation using Algorithm 1 with sample size m = 10000 (blue dots); (2) Empirical estimates of the measure of asymptotic dependence χ(u) for the simulated data with the parametric model (black line) and for the sampled data from Algorithm 1 (red line), with associated 95% pointwise confidence intervals based on 1000 bootstrap replications; (3) and (4) Quantile-quantile plots for Z 1 and Z 2 with associated 95% pointwise confidence intervals based on 1000 bootstrap replications.

Figure 2 .

 2 Figure 2.12: Sampled conditional distribution of Z 2 given Z 1 = z 1 using Algorithm 2 for the asymmetric Gaussian model with sample size n = 10000. Eight experiments are presented for different quantiles of Z 1 whose values are reported in each panel title. The sample size for each simulation is m = 10000. The theoretical conditional density is superimposed in red.

Figure 2 .

 2 Figure 2.13 shows that for these specific conditions, Algorithm 1 successfully generates extreme H c and H o . Note that since the four points considered are the largest observations, they are expected to be among the extremes of the simulated distributions.

Figure 2 .

 2 Figure 2.13: Sampled values of coastal versus offshore significant wave heights from Algorithm 1, m = 1000 pairs of points are simulated. Each panel corresponds to a fixed value of (t p , d p ) corresponding to the four largest H c from Set 2 (in decreasing order). Next to each scatter plot, marginal distributions of simulated H o and H c are depicted with boxplots. On both the scatter plots and the boxplots, the red diamond-shaped points represent the true values of the coastal and offshore significant wave heights for the four pairs.

Figure 2 .

 2 Figure 2.14: Boxplot of predicted H c conditionally on (H o , T p , D p ) using Algorithm 2. The simulation sample size is equal to m = 1000 for each observation. Red dots represents the observed H c values from Set 2. The alternating colours depict different storms: consecutive boxes with same colour correspond to observations that belong to the same storm (i.e. separated by less than 3 hours).

2. 7 .

 7 Figure 2.15: Scatter plot from Set 2 with the peak direction D p on the y-axis and the peak period T p on the x-axis. The dots' colour corresponds to the value of the coastal significant wave heights H c . The size of the dots corresponds to the value of the offshore significant wave heights H o . The shape indicates if the observed H c falls within the 95% range of the predicted distribution from the conditional model where the simulation sample size is equal to m = 1000 for each observation.

on [ 0 ,

 0 ∞] d \ {0}. The notation v -→ stands for vague convergence, a definition and its properties can be found in Chapter 3 of Resnick (2007) for example.

Figure 3 . 1 :

 31 Figure 3.1: Illustration for d = 2 of the different spaces considered. The red lines represent the excluded parts. (left) [0, ∞] 2 \{0, 0} grey shaded area, (middle) [0, x] c with x = (x 1 , x 2 ) blue hatched area, (right) E 0 = (0, ∞] 2 grey shaded area.

Figure 3 . 2 :

 32 Figure 3.2: Example of hidden regular variation with X := (X 1 , X 2 ) as defined in Example 3.2. On the original scale (left) the variables clearly appear as asymptotically independent, but on the log-scale (right) second-order regular variation emerges.

Lemma 3 . 6 .=

 36 Assume that the variable H in(3.5) is generated by the non-linear rerepresents the equality in distribution, N corresponds to a random noise and X corresponds to the explanatory variables, independent of N . If P(f (X) > u) = o (P(N > u)) , then for any classifier of the type defined by(3.11), we always have R(g) = 1.

Figure 3 . 3 :

 33 Figure 3.3: Simulated example for Lemma 3.6, for more details see equations (3.16) and (3.17).

Figure 3 .

 3 Figure 3.4: Example of predicted binary output for the decision tree classifier versus the true values of H (on a logarithmic scale). The classifier has been trained twice as explained below. On the left-hand side the classifier is trained with the subset {H > u} and the horizontal blue dashed line represents the threshold value u. On the right-hand side the classifier is trained with the subset {H > ε u } and the horizontal red dashed line represents the threshold value ε u = 0.4 × u.

Figure 3 . 5 :

 35 Figure 3.5: Estimation of R ε (g) for different classifiers (cross-validation with 50 repetitions). In red (top-left) are the estimates when ε = 0 and in blue for different values of ε > 0 (ε ∈ {0.4, 0.6, 0.8}). At the top of each plot, the value of ε and the number of points such that H > ε u from the testing set are given.

Figure 3 .

 3 Figure 3.6: River map of the upper Danube basin, showing sites of the 31 gauging stations along the Danube and its tributaries. Water flows toward gauging station 1. The stations represented by a green triangle shaped dot are the three stations of interest as described in Section 3.5.

Figure 3 . 7 :

 37 Figure 3.7: Summer daily measurements of river discharges from station 1 (y-axis) against station 23 (left x-axis) and station 24 (right x-axis) of Figure 3.6. The blue dotted lines are the threshold value u (the 85th percentile of station 1). The black dots on the graphs in the top row are the values such that min(X j , X 1 ) > ε u where ε u = 0.6u. The graphs in the bottom row corresponds to the same data plots but on unit Frechet scale highlighting potential asymptotic independence in the data.

  Part of this work was supported by the DAMOCLES-COST-ACTION on compound events, the French national program (FRAISE-LEFE/INSU and 80 PRIME CNRS-INSU), and the European H2020 XAIDA (Grant agreement ID: 101003469) and funded by Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany's Excellence Strategy -EXC 2075 -390740016. The authors also acknowledge the support of the French Agence Nationale de la Recherche (ANR) under reference ANR-20-CE40-0025-01 (T-REX project), the ANR-Melody and the Stuttgart Center for Simulation Science (SimTech). * * *

. 23 )

 23 If the function θ → R(g θ ) has a unique minimiser θ * in C, then the estimator θn,un = argmin θ∈C Rn (g θ (•; u n )).

Figure 3 .

 3 Figure 3.10: Estimation of R 0 (g) for different classifiers g (cross-validation -70% train, 30% test -with 50 repetitions). The red distributions come from the estimation with all the stations, and the green distributions represent the estimations with only the variables having ĉi ̸ = 0.

  ), 1 ≤ k ≤ m.

Figure 4

 4 Figure 4.1: Two possible applications of a multivariate simulation model: (left) Refinement of coastal simulations; (right) Gap-filling in historical data.

  ε ′ P a (A ε ) P b (B ε ′ ) with P(A ε ) = P(g(X) > εu) = L g (εu)ε -αg u -αg , P(B ε ) = P(H > εu) = L H (εu)ε -α H u -α H and c ε,ε ′ = L(εu, εu)/L(u, u), a = α g /(2η) and b = α H /(2η).

Table of contents

 of A.1.2 Proof of Lemma 3.6: . . . . . . . . . . . . . . . . . . . . . . . . A.1.3 Proof of Proposition 3.7: . . . . . . . . . . . . . . . . . . . . . A.1.4 Proof of Proposition 3.8: . . . . . . . . . . . . . . . . . . . . . A.2 Proofs of the appendix . . . . . . . . . . . . . . . . . . . . . . . . . . A.2.1 Proof of Proposition 3.9 . . . . . . . . . . . . . . . . . . . . . . A.2.2 Proof of Lemma A.1 . . . . . . . . . . . . . . . . . . . . . . . . Carte des données de vagues . . . . . . . . . . . . . . . . . . . . . . 2 H s à la côte en fonction des H s au large . . . . . . . . . . . . . . . . . 3 Carte du bassin versant du Danube . . . . . . . . . . . . . . . . . . . 4 Prédictions issues d'un arbre de décision . . . . . . . . . . . . . . . . Grid of the wave data and scatter plot of significant wave heights . . 2.8 Empirical histograms of H s threshold exceedances . . . . . . . . . . 2.9 Dependence of H s with respect to T p and D p . . . . . . . . . . . . . . OF FIGURES 2.13 Sampled values of coastal versus offshore H s from Algorithm 1 . . . 2.14 Predicted H c conditionally on (H o , T p , D p ) using Algorithm 2. . . . . 2.15 Non-stationary illustration of the predicted H c using Algorithm 2 . . 2.16 Conditional marginal EGPD fits . . . . . . . . . . . . . . . . . . . . . . 3.1 Space considered for multivariate regular variation . . . . . . . . . . 3.2 Example of hidden regular variation . . . . . . . . . . . . . . . . . . . 3.3 Simulated example for Lemma 3.6 . . . . . . . . . . . . . . . . . . . . 3.4 Predicted binary output for the decision tree classifier . . . . . . . . 3.5 Estimation of R ε (g) for different classifiers . . . . . . . . . . . . . . . 3.6 River map of the upper Danube basin . . . . . . . . . . . . . . . . . . 3.7 Summer daily measurements of river discharges from station 1 (yaxis) against station 23 (left x-axis) and station 24 (right x-axis) . . . . 3.8 Estimation of R ε (g) for five different classifiers with the Danube river discharges . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .3.9 Correlation plot between Station 1 and the stations considered to contribute to the extremes of Station 1 . . . . . . . . . . . . . . . . . 3.10 Estimation of R 0 (g) for five different classifiers . . . . . . . . . . . . 3.11 Predictive power of the linear classifier . . . . . . . . . . . . . . . . .

	List of Figures
	Bibliography
	1

1.1 Densities of the three generalised extreme value distributions . . . 1.2 Densities of the GP distribution for different values of the shape parameter ξ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3 Upper tail distribution of significant wave heights . . . . . . . . . . . 1.4 Block maxima versus peaks over thresholds . . . . . . . . . . . . . . 1.5 Two possible definitions of bivariate extremes . . . . . . . . . . . . . 1.6 Asymptotic independence of a bivariate Gaussian vector . . . . . . . 1.7 Bivariate logistic extreme value model . . . . . . . . . . . . . . . . . 1.8 Scatter plot of significant wave heights . . . . . . . . . . . . . . . . . 1.9 Illustration of u(t) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.10 Bivariate generalised Pareto distributions . . . . . . . . . . . . . . . 1.11 Asymptotic dependence behaviour . . . . . . . . . . . . . . . . . . . 1.12 Spectral measure and Pickands dependence function . . . . . . . . . 1.13 Dependence measure χ for wave data . . . . . . . . . . . . . . . . . 1.14 Asymptotic dependence comparison on Pareto scale . . . . . . . . . 2.1 Spectral decomposition of a sea state . . . . . . . . . . . . . . . . . . 2.2 Wave spectrum . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3 Hindcast databases grid . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4 Illustration of a wind field . . . . . . . . . . . . . . . . . . . . . . . . . 2.5 Measure of asymptotic dependence between wave data . . . . . . . 2.6 Densities of the extended GP distribution for different values of the shape parameter κ . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.7 2.10 Maps illustrating the covariates effects . . . . . . . . . . . . . . . . . 2.11 Simulation results for Algorithm 1 . . . . . . . . . . . . . . . . . . . . 2.12 Simulation results for Algorithm 2 . . . . . . . . . . . . . . . . . . . .

LIST 4.1 Possible applications of a multivariate model . . . . . . . . . . . . .

  responding to the SEM-REV sea test site[START_REF] Mouslim | Development of the french wave energy test site SEM-REV[END_REF] (see Figure2.7). Among all the sea states parameters, the significant wave heights H s [m] from the two locations, the peak period T p [s] and the peak direction D p [ • ] only from the offshore location are used. Data are available at 3-hour intervals spanning from 1994 to 2016. More precisely, the original HOMERE database has a 1-hour resolution time step but the IOWAGA database is sampled every 3 hours, so to obtain data at the same time scale a sub-sampling of the HOMERE database every 3 hours is performed. In the following, H c (resp. H o ) denotes the coastal (resp. offshore) significant wave heights. A scatter plot between H c and H o can be found in Figure

• N, 3 • W) from the IOWAGA database and a coastal point from the HOMERE database, near the French coast (47 • 24N, 2 • 78W) cor-

Table 2 .

 2 4 describes the eight largest h c of Set 2, giving the time event and the corresponding h o , t p and d p values. The last column of Table 2.4 gives the empirical estimate of the joint survival probability defined as

Table 3 .

 3 4: Empirical mean estimations of the limits c * i,ε for different sample sizes

	1,ε	0.34 0.34 94.8% 0.34 93.7%	0.34 93.7%	0.34 93.6% 0.35 94.2%
	c * 2,ε	0.49 0.49 95.8% 0.49 95.0%	0.49 93.8%	0.49 94.3% 0.49 95.0%
	c * 3,ε	0	0.36 0%	0.16 0%	0.09 41.7%	0.05 99%	0 100%
	c * 4,ε	0	0.12 92.8% 0.02 98.6%	0 99.7%	0 100%	0 100%

*

  2 o , T 2 p , D 2 p ∼ EGP D, and to consider that the scale parameter of the EGPD vary as a smooth function of both (T 1 p , D 1 p ) and (T 2 p , D 2 p ) (assuming that we have two offshore grid point, site 1 with parameters (H 1 o , T 1 p , D 1 p ) and site 2 with (H 2 o , T 2 p , D 2 p

1.1. Univariate extremes

A copula (in dimension

2) is a bivariate distribution function on [0, 1] × [0, 1] with uniform margins in [0, 1][START_REF] Nelsen | An Introduction to Copulas[END_REF].

1.2. Multivariate extremes

If 1/2 < η ≤ 1, X and Y are positively associated: large observations for both X and Y occur more frequently than under exact independence.

If η = 1/2, extremes of X and Y are near independent. If additionally l(x) = 1, the variables are exactly independent.

If 0 < η < 1/2, X and Y are negatively associated: large observations for both X and Y occur less frequently than if X and Y were independent.In[START_REF] Heffernan | A directory of coefficients of tail dependence[END_REF] values of the coefficient of tail dependence η and slowly varying function l can be found for a wide range of bivariate distributions, including both asymptotic dependence and asymptotic independence. Finally, χ and η can be linked through χ = 2η -1.

IPCC (2014) defined different scenarios given the evolution of greenhouse gas concentrations in the atmosphere, named RCP, according to the climate policies adopted.

2.1. Preamble to Paper I

Reanalysis data are obtained by combining historical observations and numerical model outputs through data assimilation (e.g.[START_REF] Evensen | Data assimilation: the ensemble Kalman filter[END_REF], providing more consistent databases.

VGAMs are more general models than GAMs, that allow for multiple linear predictors and encompass models outside the exponential family distribution[START_REF] Yee | Vector generalized linear and additive models: with an implementation in R[END_REF]. These models encompass the cases where Y is a vector.

2.6. Application to extreme significant wave height

U ← (t) := inf{s; U (s) ≥ t} denotes the left-continuous inverse of U , and U (U ← (t)) ∼ t.

Although we will apply Lemma

3.4 to sets A ε that are rare events, this is a not a necessity.3 the shortcut notation R(g) corresponds to the case where g(X; u) = +1 is built from the event g(X) > u in the associated R(g).

Section 3.7.3 treats a simpler case where station 1 is predicted from the whole set of remaining stations. In this case, strong dependencies among station 1 and other stations can be observed. So, the main issue is to select these stations, a problem discussed in Section 3.7.

3.7. Supplements to Paper II

Remerciements Acknowledgements

The authors would like to thank Noémie Le Carrer for helpful discussion on the implementation of the EGPD family within the gamlss package. The authors are indebted to the Laboratoire d'Océanographie Physique et Spatiale, IFREMER for providing the data used in this study. Part of this work was supported by the DAMOCLES-COST-ACTION on compound events, the French national program (FRAISE-LEFE/INSU and 80 PRIME CNRS-INSU), and the European H2020 XAIDA (Grant agreement ID: 101003469. The authors also acknowledge the support of the French Agence Nationale de la Recherche (ANR) under reference ANR-20-CE40-0025-01 (T-REX project), the ANR-Melody (ANR-19-CE46-0011).

2.A Marginal regression modelling

If g(X; u) ≡ -1, then P(g(X; u) ̸ = Y (u) ) = P(H > u) goes to zero as u gets large. Hence, the classical risk function E(1{g(X; u) ̸ = Y (u) }) will systematically favour the always optimistic classifier for extremes. To avoid this undesirable feature, the loss function has to be modified. One natural idea is to re-scale by P(H > u) and introduce the loss function 1{g(X; u) ̸ = Y (u) }/ P(H > u). In this case, the risk, i.e., the expected loss E(1{g(X; u) ̸ = Y (u) }/ P(H > u)), goes towards one as u gets large.

Another trivial but also interesting case is the "crying wolf" forecaster who always issues g(X; u) ≡ +1, see also the forecaster's dilemma (e.g. [START_REF] Lerch | Forecaster's Dilemma: Extreme Events and Forecast Evaluation[END_REF]. In this case, Equation (3.6) implies that P g(X; u) ̸ = Y (u) = P(H ≤ u) and, consequently, the risk E(1{g(X; u) ̸ = Y (u) }/ P(H > u)) goes towards infinity as u gets large. This limiting cost indicates that the "crying wolf" forecaster is much worse than the overly optimistic one. Both of them are unreasonable in practice and there is no reason to strongly favour one over the other one. For this reason, we propose to use a following weighted loss function

and the associated risk

By construction, the event {g(X; u) ̸ = Y (u) } implies that {Y (u) = 1} or {g(X; u) = 1} and therefore, necessarily, R (u) (g) ∈ [0, 1]. In particular, the naive classifier g(X; u) ≡ -1 possesses unit risk with R (u) (g) = 1 at each level u > 0.

Similarly, the risk of the "crying wolf" classifier g(X; u) ≡ +1 is then equal to R (u) (g) = P(H ≤ u) and converges to one as u → ∞.

Hence, the value of one is reached by the two worst cases scenarios in terms of classifiers. This unit value provides a clear benchmark that can be compared to any other classifier satisfying the existence of the limit

We call such classifiers extremal.

In the weather forecast literature (e.g. [START_REF] Schaefer | The critical success index as an indicator of warning skill[END_REF], the definition of R (u) (g) can be linked to the critical success index, also called the threat score. The critical success index computes the total number of correct event forecasts (hits) divided by the total number of forecasts plus the number of misses (hits + false alarms + misses). Hence, 1 -R (u) (g) can be understood as a critical success index for extremes. In the context of rare events forecasts, [START_REF] Stephenson | The extreme dependency score: a non-vanishing measure for forecasts of rare events[END_REF] highlighted some advantages and drawbacks of various risk functions, including the This indicates that, whenever g(X) and H are not tail equivalent, the classifier g(X) cannot outperform naive classifiers with respect to R (u) (g) for large u. Hence, the marginal behaviour of g(X) has a direct impact of its predicting capacity in terms of R(g). This is a widely known fact in forecast verification. In particular, a paradigm promoted by Gneiting and his co-authors (see, e.g. [START_REF] Gneiting | Probabilistic forecasts, calibration and sharpness[END_REF] is to work only with calibrated forecasts. In our case, calibration means that g(X) and H have the same distributions, and consequently P(g(X) > u)/ P(H > u) = 1 for all u. In practice, it may be difficult to ensure that this constraint holds for extremes (see, e.g. [START_REF] Lerch | Forecaster's Dilemma: Extreme Events and Forecast Evaluation[END_REF][START_REF] Taillardat | Extreme events evaluation using CRPS distributions[END_REF]. To illustrate this, suppose that P(X > u) = u -1 for all u ≥ 1, and the variable H, independently of the value of X, is either equal to δX or (2 -δ)X with probability .5 and the constant δ ∈ (0, 1). Then, we have

Hence, X and H are more than tail equivalent, they have identical tail behaviour for large u. Concerning classifiers, linear ones of the type g a (X) = aX with a > 0 belong to the class defined by (3.11). They are tail equivalent to H, and R(g a ) < 1.

Although X and H have identical tail behaviours, the choice of a = 1 is not optimal with respect to R(g a ). In particular, one can show that

This is not surprising. By construction, the largest values of H are more likely to be produced by (2 -δ)X than X, especially if δ is small. In this context, the following lemma (that is a rewriting of Lemma 3.4) explains that the risk function R(g) both depend on the upper tail dependence between g(X) and H and their marginal behaviours.

then the limiting risk based on (3.7) has the following expression

where χ * (g) denotes the limit of P(g(X) > u | H > u). In particular,

build a simple tool to detect relevant features. In Paper II, this methodology was then applied to the Danube river network, as reproduced below.

Hereinafter, we return to the bold notation for the vectors of R d , i.e. x = (x 1 , . . . , x d ) ∈ R d .

Linear Classifiers: Definition, Basic Properties and Inference

In this section, we consider a specific type of classifiers which in Paper II is referred to as linear classifiers, i.e. classifiers of the form

This is a simple binary classifier that is based on the idea that a large value along an appropriate linear combination could produce extremes in the hidden variable

To obtain an optimal linear classifier of g θ (X; u), i.e. some weight vector θ * such that the classification risk R(g θ * ) gets minimal, we need to impose some joint extremal dependence structure on X and H from (3.6).

Even though some of the results can also be obtained in a similar manner in a more general framework for the conditional risk R ε (see Section 3.7.4), henceforth, we will focus on the asymptotically dependent case where we might find some optimal classifier with unconditional risk R(g θ * ) < 1. As discussed before, in this case, at least one component of X needs to have a similar tail behaviour as H. A natural assumption is therefore that (X, H) is jointly regularly varying on [0, ∞) d+1 with index α > 0, i.e. there exists an α-Pareto random variable P and, independently of P , a random vector (Γ, Ω) ∈ [0, ∞) d × [0, ∞), the so-called spectral tail vector, on the unit sphere

Note that assuming joint regular variation, implies in particular that H and all the relevant components of X are heavy-tailed. For example, this can be the case for rainfall data (e.g. Le [START_REF] Gall | Improved regional frequency analysis of rainfall data[END_REF] or river discharges (see Table 3.2). But other environmental data sometimes do not show such behaviour, as in the case of wind gusts [START_REF] Friederichs | A probabilistic analysis of wind gusts using extreme value statistics[END_REF], temperatures [START_REF] Toulemonde | Applications of extreme value theory to environmental data analysis[END_REF] or significant wave heights [START_REF] Legrand | Joint stochastic simulation of extreme coastal and offshore significant wave heights[END_REF], which are more light-tailed or even upper-bounded. Some precautions should therefore be taken before considering a dataset, by applying for instance an univariate extreme value analysis (e.g. Embrechts et al., 2013, Chapter 6).

Supplements to Paper II

The second result in Lemma 3.11 emphasises that for partially linear classifiers, the components X i for which c i,ε = 0, do not have any effect, and can thus be removed as in the asymptotic dependent case.

Note that if we assume that the tail function of H is regularly varying with index α H , then, lim u→∞ P(H > u | H > εu) = ε α H . Thus, instead of estimating c i,ε , we can equivalently estimate

Therefore, a good candidate is the following empirical estimator ĉi,ε :=

, for large u n .

The inference results for the estimation of the c * i,ε , when ε > 0, are obtained through the following proposition, which is the analogue of Proposition 3.10 for ε > 0.

Proposition 3.12. Assume that

exists. If u n → ∞ and n P(X i > εu n ) → ∞ and, additionally,

.

The same simulation framework as in Section 3.4 is carried out, i.e. 2) and H = X 1 + N , where X 1 , . . . , X 4 and N are independent.

The theoretical values of the limits c * i,ε are then given by

Now to compare with the theoretical values, we generate 1000 samples with sample size n and threshold u defined as a given quantile of H. Different sample sizes are considered n ∈ {5000, 1e4, 2e4, 1e5} and u is chosen such that we get 300 exceedances each time. The 95% confidence intervals are obtained using the convergence result of Proposition 3.12.

As P(A 1 ∩B 1 ) = P(A ′ 1 ∩B 1 ) and P(A

. This provides the second statement (b) since we assume that

Item (c) is a direct consequence of (3.8).

□

A.1.2 Proof of Lemma 3.6:

From Eq. ( 3.12), we know that we can only get R(g) < 1, only if H and g(X) are tail equivalent. Thus, this will be assumed in the following. Let w such that δ := lim u→∞ w(u)/u ∈ (0, 1), then for any positive u we can write that

.

Since g(X) and N are independent, the second term reduces to P(N > uw(u)) which converges to 0 since u -w(u) ∼ (1 -δ)u as u gets large.

For the first term, we rewrite the ratio as follows

.

Since w(u) → ∞ and P(f (X) > u) = o (P(N > u)), the ratio P(f (X)>w(u)) P(N >w(u)) goes to 0 as u gets large. From the assumption w(u) ∼ δu, P(N >w(u)) P(N >u) behaves as a constant when u → ∞.

The only remaining term is P(N >u) P(g(X)>u) which converges to a constant due to tail equivalence. So, lim P(g(X) > u | H > u) = 0. □

A.1.3 Proof of Proposition 3.7:

In Lemma 3.4 we fix A ε = {g(X) > εu} and B ε = {H > εu} and A 1 = {g(X) > u} and B 1 = {H > u}.

A.2 Proofs of the appendix

A.2.1 Proof of Proposition 3.9

The proof is based on the following lemma which is proven in Subsection A.2.2.

Lemma A.1. Under the assumptions from Proposition 3.9, for every compact subset K ⊂ C, the sequences of processes {A n (θ), θ ∈ K} and {B n (θ), θ ∈ K} defined by

converge to centered Gaussian processes {A(θ), θ ∈ K} and {B(θ), θ ∈ K}, respectively, weakly in ℓ ∞ (K).

If the function θ → R(g θ ) has a unique minimizer θ * , then, necessarily, R(g θ * ) < 1. Now, similarly to the notation above, let J denote the set of indices j with c j > 0, and let us consider θ ∈ C such that ∥θ∥ ∞ > k 0 for some constant k 0 > 0. Then,

where the right-hand side goes to 1 as k 0 → ∞. Thus, as R n (g θ * ) → p R(g θ * ) < 1, we obtain that, for sufficiently large k 0 ≫ ∥θ * ∥, with probability going to one,

and, consequently,

Now, we note that, by Lemma A.1, the bias conditions (3.22) and (3.23) and the functional delta method, R n (g θ ) converges in probability to R(g θ ) uniformly on every compact subset of C. In particular,

Thus,

Proof of Lemma A.1

We will proof the lemma by applying the Central Limit Theorem 2.11.9 in [START_REF] Van Der Vaart | Weak Convergence and Empirical Processes[END_REF]. To this end, we define the function spaces A = {a θ , θ ∈ K} and B = {b θ , θ ∈ K} where

Then, with

for l = 1, . . . , n, we have that

for all l = 1, . . . , n and n ∈ N. Consequently, we check the Lindeberg condition: For k ∈ N, we have

, we obtain a Lindeberg type condition that ensures convergence of A n and B n to A and B, respectively, in terms of finite-dimensional distributions. For k = 1, we obtain the Lindeberg type condition of Theorem 2.11.9 in [START_REF] Van Der Vaart | Weak Convergence and Empirical Processes[END_REF].

It remains to check the equi-continuity condition. In the following, to simplify notation, we assume that C

and the probability that any of those two expressions is equal to one is bounded by the probability

where we use that b ⊤ X > u n implies that ∥X∥ > u n /(∥K∥d) with ∥K∥ = sup x∈K ∥x∥ ∞ . Making use of the fact that P (∥(X, H)∥ > u n /(∥K∥d)) ≤ C 0 (∥K∥d) α P(H > u n ) for some constant C 0 > 0 and the bound given by Equation (3.21), we obtain that E (Z nl (f ) -Z nl (g)) 2 = 2C 0 δ, which tends to 0 as δ → 0. From this inequality, it can also be seen that any partition of K into hypercubes with length ε 2 /(2C 0 ) leads to a valid ε-bracketing, i.e. the number N ε ∝ 1/ε 2d grows with a power rate and, so, log(N ε ) is integrable.

Thus, by Theorem 2.11.9, the processes A n and B n converge to Gaussian processes A and B, weakly in ℓ ∞ (K).