Modélisation statistique des valeurs extrêmes
M1 Actuariat, Université Félix Houphouët-Boigny
Nicolas Raillard ([Link]@[Link])
Année scolaire 2020 — 2021
(Dernière mise à jour: 18/03/2021)
Modélisation statistique des evènements extrêmes
Présentation
Nicolas Raillard:
Chercheur à l'IFREMER depuis 2016;
Statisticien de formation (ENSAI, Université de Rennes 1);
Domaines d'étude: modélisation des extrêmes, processus stochastiques,
apprentissage statistique, statistique environementale et géostatistique.
12h de cours: 6h de CM et 6h de TD/TP;
Objectifs:
Comprendre le contexte et les particularités du domaine;
Connaître les outils probabilistes et statistiques;
Savoir mettre en pratique avec R;
2 / 40
Introduction
Risques naturel: vagues et surcote
Tempête Xynthia, La Faute-Sur-Mer, 1er Mars 2010
3 / 40
Introduction
Risques naturel: pluies torrentielles
Inondations du 26 octobre 2012, La Garde
4 / 40
Introduction
Risques nanciers
Gauche: valeur de fermeture de l'indice S&P, du 04/01/1960 au 11/06/1993, droite: log-return.
5 / 40
Introduction
Risques nanciers: queues lourdes
6 / 40
Introduction
Besoin d'extrapoler (1/2)
Soit une grandeur X (hauteur des vagues, contrainte mécanique, montant à
rembourser par un assureur);
Etant donné un seuil (hauteur de digue, résistance de matériaux, fonds propre...),
calculer ph = IP (X > h);
Plus fréquent: étant donné p (e.g. p = 1/10000), trouver hp / p = IP (X > hp ) 7 / 40
Introduction
Besoin d'extrapoler (2/2)
Souvent, on observe X sur une période limitée (10, 20, 100 ans);
Echantillon de max annuels sur N années: X1max , … , XN
max
N
Fonction de répartition empirique des max annuels : F
ˆ
n (h) = ∑
i=1
𝟙
Xi <h
Si h > max(X
max
i
) , ou si p ,
ˆ
≤ 1/N Fn (h) = 1
Besoin d'un modèle pour extrapoler.
Par exemple X ∼ N :
Ajusté aux données ?
Calcul de la loi des maxima et estimation de hp (ou calcul de ph );
8 / 40
Introduction
Références:
Embrechts, P., Klüppelberg, C., & Mikosch, T. (1997). "Modelling Extremal Events". Springer
Berlin Heidelberg;
Coles, S. (2001). "An Introduction to Statistical Modeling of Extreme Values". Springer
Series in Statistics.;
Reiss, R.-D., Thomas, M. (2007). "Statistical Analysis of Extreme Values, with Applications
to Insurance, Finance, Hydrology and Other Fields". Birkhäuser Basel;
de Haan, L., Ferreira, A. F. (2006). "Extreme Value Theory, An Introduction". Springer-Verlag
New York.
Packages R:
extRemes
ismev
evd
POT
9 / 40
Plan
1. Introduction
2. Lois des maxima par blocs
3. Lois des dépassements de seuils
4. Cas de données dépendantes
5. Applications
10 / 40
Loi des maxima: uctuation des sommes
Dans toute la suite, on considère un échantillon i.i.d X1 , … , Xn et X la loi commune;
Usuellement, on s'intéresse à Sn = X1 + … Xn ;
Loi des grand nombre:, Sn /n → IE(X);
Sn
Théorème central-limite: √ V(X)
n
(
n
− I
E(X)) → N (0, 1)
Normalisation nécessaire: Sn → ∞:
Caractérisation de la loi limite.
Hands on R !
Comment trouver la loi limite par simulation ?
Tester différentes lois :
];
1 1
U ∼ U ([− ,
2 2
Z ∼ N (0, 1) ;
X := sign(U ) × log(1 − 2|U |) (Loi de Laplace)
T ∼ T (1.8)
...
Refaire l'expérience sur la loi des maxima.
11 / 40
Loi des maxima: uctuation des sommes
12 / 40
Loi des maxima
13 / 40
Loi des maxima
Propriété
i.i.d
Si Xi ∼ X alors:
n
I
P [max(X1 , … , Xn ) ≤ x] = ∏ I
P [Xi ≤ x] = I
P [X ≤ x]
Exemples:
U
i.i.d Mn
Un ∼ U ([0, θ]) , trouver la loi de n ( θ
U
− 1) , Mn = max(U1 , … , Un ) ;
i.i.d
X ∼ E(λ) , trouver la loi de (λMnE − log(n)), Mn
E
= max(E1 , … , En )
14 / 40
Loi des maxima
Théorème de Fischer-Tipett-Gnedenko
Soient (Xi )i≥0 des v.a. i.i.d, et Mn = maxi≤n Xi . S'il existe des suites an > 0 et bn ∈ I
R
et une v.a. Y non-dégénérée, telles que
Mn − bn d
→ Y,
an
alors ∃μ ∈ I
R, σ > 0, ξ ∈ I
R tels que
1
−
x−μ
ξ
−[(1+ξ )+ ]
σ
I
P [Y < x] = Gμ,σ,ξ (x) = e ,
avec (y)+ = max(y, 0) . Ceci dé nit une Distribution de valeurs extrêmes généralisée (GEV).
−x
Dans le cas où ξ = 0, G0,1,0 (x) = e
−e
.
15 / 40
Loi des maxima
ξ > 0 :
Support : [μ − ;
σ
, +∞[
ξ
Loi de Fréchet, heavy tailed, à queues lourdes (leptokurtique);
ξ < 0:
Support : ]−∞, μ − σξ ];
Loi de Weibull : maxima bornés supérieurement;
ξ = 0;
Support : ]−∞, +∞[;
Loi de Gumbel, queue "légère" (moins que la loi normale!). 16 / 40
Loi des maxima
Soit X ∼ GEV(μ, σ, ξ) . Alors:
σ
⎧ μ + (Γ(1 − ξ) − 1) si ξ < 1 et ξ ≠ 0
⎪
⎪ ξ
E(X) = ⎨ μ + γσ si ξ = 0
⎪
⎩
⎪
∞ si ξ ≥ 1
2
(Γ(1−2ξ)−[Γ(1−ξ)] )
⎧
⎪
⎪ 2
⎪ σ + si ξ < 1/2 et ξ ≠ 0
⎪ ξ
2
V(X) = ⎨ 2
2 π
σ si ξ = 0
⎪
⎪ 6
⎪
⎩
⎪
∞ si ξ ≥ 1/2
Fonction quantile:
⎧ μ − σ log(− log(p)) si ξ = 0 et p ∈ (0, 1)
Q(p; μ, σ, ξ) = ⎨ −ξ
σ
⎩μ + ((− log(p)) − 1) si ξ > 0 et p ∈ [0, 1)
ξ
17 / 40
Loi des maxima
Idée de la preuve:
d
Mn −bn
On a an
→ Y , i.e. F n (an x + bn ) → G(x) . Donc:
nk
∀k ≥ 1, F (an x + bn ) ⟶ G(x). (1)
n→∞
On a aussi
nk n k k
F (an x + bn ) = [F (an x + bn )] → G (x). (2)
Avec (1) et (2), on montre que G est max-stable, c.a.d :
k
∀k ≥ 1, ∃αk , β / G (αk x + βk ) = G(x)
k
En résolvant cette équation fonctionnelle, on a le théorème.
18 / 40
Domaine d'attraction
Sous les conditions du Théorème de Fisher-Tippet-Gnedenko, on dit que F est dans le
domaine d'attraction de
Fréchet si ξ > 0 :
Exemples types: Loi de Fréchet, Loi de Cauchy, Loi de Student...
Applications: précipitations, crues de rivières, pics de pollution, rendements
logarithmiques (Dow Jones...);
Weibull si ξ < 0:
Exemples types: Loi de Weibull, loi uniforme...
Applications: Résistance des matériaux, températures, durées de vie...
Gumbel si ξ = 0:
Loi de Gumbel, Loi Normale, Loi Lognormale, Loi exponentielle...
Applications: hydrologie (Vagues, débits, niveaux d'eau...), vent...
19 / 40
Estimation des paramètres
Deux méthodes principales: Maximum de vraisemblance et méthode des moments
pondérés;
Nombreux packages R ( ismev , extRemes ...) et python ( scikit extreme , [Link] ...)
20 / 40
Maximum de vraisemblance (1/2)
Soient yi=1,…,n des observations indépendantes d'une loi Fθ (y);
dFθ
θ un paramètre (vectoriel éventuellement) et fθ (y) =
dy
(y) sa densité;
n
Vraisemblance L(θ|y1 , … , yn ) = ∏
i=1
fθ (yi )
On choisi θˆ = arg max L(θ|y1 , … , yn )
θ
Sous des hypothèses générales, θ̂ est asymptotiquement normal;
Cas de la GEV:
θ = (μ, σ, ξ) ;
Méthode valide seulement pour ξ > −0.5, car le support dépend des paramètres.
On maximise généralement la log-vraisemblance:
l(θ|y1 , … , yn ) = log(L(θ|y1 , … , yn ));
n
1 yi − μ
l(θ|y1 , … , yn ) = −n log σ − (1 + ) ∑ log[1 + ξ ( )]
ξ σ
i=1
n −1/ξ
yi − μ
− ∑ [1 + ξ ( )]
σ
i=1
yi −μ
Si ∀i, 1 + ξ ( σ
) > 0 ; sinon l(θ|y1 , … , yn ) = −∞
21 / 40
Maximum de vraisemblance (2/2)
On considère toute la série du phénomène étudié;
22 / 40
Maximum de vraisemblance (2/2)
On considère toute la Série temporelle du phénomène étudié;
Les données sont aggrégées dans k blocs de même taille (p.e. une année);
23 / 40
Maximum de vraisemblance (2/2)
Les données sont aggrégées dans k blocs de même taille (p.e. une année);
On conserve le maximum Mi chaque bloc;
On suppose M1 , … , Mk ∼ GEV
iid
Perte de données : on n'utilise que une valeur par bloc ;
Compromis entre le biais (beaucup de petits blocs) et variance (peu de grands blocs).
La loi des maxima annuels peut être obtenue à partir de blocs de taille inférieure, p.e.
max annuel max mensuel 12
P(X < h) = P(X < h)
24 / 40
Exemple
library(extRemes)
fitgevpirie fevd(evd portpirie)
fitgevpirie
fevd(x = evd portpirie)
[1] "Estimation Method used: MLE"
Negative Log-Likelihood Value: -4.339058
Estimated parameters:
location scale shape
3.8747499 0.1980440 -0.0501095
Standard Error Estimates:
location scale shape
0.02793224 0.02024798 0.09825416
Estimated parameter covariance matrix. 25 / 40
Niveaux de retour
Niveau de retour (quantile) zp correspondant à une période de retour de 1
p
:
σ
^ ^
⎧
⎪
−1/ξ
^
⎪ μ
^ − [1 − yp ], si ξ ≠ 0
−1 ^
^p = G
z (1 − p) = ⎨ ξ
(μ,
^ σ, ^)
^ ξ
⎪
⎩
⎪
^
μ
^ − σ
^ log(yp ), si ξ = 0
avec yp = − log(1 − p) .
Niveau de retour centennal, z100 / IP (X max annuel
1
> z100 ) =
100
max sur 100 ans max annuel 100 1 100
I
P (X < z100 ) = I
P (X < z100 ) = (1 − ) ≈ 0.368
100
Sur 100 ans, la probabilité de dépasser le niveau de retour centennal est de 63% !
Intervalles de con ance via la delta method:
T
^ p ) ≈ ∇z V ∇z ,
Var(z
p p
avec V la matrice de covariance de (μ, ^ ξ ) et
^ σ,
^
26 / 40
Diagnostic graphique pour le domaine de Gumbel
x−μ
−
Si Fn (x) , alors Fn← (u) :
σ
−e
= e = σ[− log(− log(u))] + μ
, la i-ème plus petite observation;
← i
Fn ( ) ≈ x(i)
n+1
le graphe des (− log(− log( n+1
i
)); x(i) ) ≈ une droite
xord sort(portpirie,decreasing=F)
n length(portpirie)
inds (1:n)/(n 1)
gbquant log( log(inds))
plot(gbquant,xord,pch=20)
[Link] lm(xord~gbquant)
coeff [Link]$coefficients
abline(coeff[1],coeff[2],
col="red",lwd=2)
Le modèle de Gumbel semble adapté.
Conséquence : niveaux de retour beaucoup
plus élevés qu'avec une Weibull (comparer
avec slide 24). 27 / 40
Intervales de con ance: pro le likelihood (1/3)
On veux de "bons" intervalles de con ances sur les paramètres et les niveaux de retour;
Vraisemblance de con ance basé sur la normalité asymptotique: valide, mais pas
optimaux (vraisemblance non symétrique autour de θ^);
Vraisemblance pro lée (pro le likelihood):
θ = (μ, σ, ξ) , vraisemblance L(θ; y1 , … , yn ) = ∏ fθ (yi )
i
;
on xe θ1 (μ,
^ σ ou
^ ξ ) et on maximise la vraisemblance sur les paramètres restant
^
θ ; 2
on trace L1 (θ1 ) = max L(θ; y1 , … , yn )
2
θ
28 / 40
Intervales de con ance: pro le likelihood (2/3)
Intervalle de con ance pour θ1 : zone de forte vraisemblance pro lée;
^ 1 2 1
maxθ L(θ) = L(θ ) = max 1 {max 2 L(θ , θ )} = max 1 L(θ )
θ θ θ
Rapport de vraisemblance : 2 (log L(θ^) − log L(θ10 )) ∼ χ
2
1
⇒ région de con ance à 95%: log L(θ1 ) ^
,
> log L(θ ) − Δ Δ = 1/2q
χ
2 (0.95) .
1
29 / 40
Intervales de con ance: pro le likelihood (3/3)
Table: IC à 95% - Vraisemblance pro lée
95% lower CI location 95% upper CI
location 3.82 3.87 3.93
scale 0.16 0.20 0.24
shape -0.22 -0.05 0.17
100yr return level 4.49 4.69 5.26
Table: IC à 95% - Approximation normale
95% lower CI location 95% upper CI
location 3.82 3.87 3.93
scale 0.16 0.20 0.24
shape -0.24 -0.05 0.14
100yr return level 4.38 4.69 5.00
30 / 40
Plan
1. Introduction
2. Lois des maxima par blocs
3. Lois des dépassements de seuils
4. Cas de données dépendantes
5. Applications
31 / 40
Loi des dépassements de seuil
Dans la méthode précédente, on ne garde que le maximum par bloc: perte
d'information conséquente;
On va s'intéresser aux dépassements de seuil:
¯
I
P (X > y + u|X > u) = F u (y) ⟶ . . .
u→∞
Théorème de Pickands:
1
−
x ξ
I
P (X > x + u|X > u) ≈ (1 + ξ )
u→∞ ~
σ
~
avec σ = σ + ξ(u − μ) .
1
−
La loi dé nie par Hσ,ξ (x) est appellée Loi de Pareto Généralisée,
x ξ
= 1 − (1 + ξ )
σ +
GPD.
Méthodologie POT:
Nu
Estimation de la fonction de survie F¯ avec F¯(u) =
n
où
n
¯
Nu = ∑
i=1
1
1(−∞,u] (xi );
Pour y ¯
> u, Fu (y) ≈ H̄ σ,ξ (y) donc
1
−
Nu y−u ξ
¯ ¯ ¯
F (y) = F (u)Fu (y − u) ≈ (1 + ξ )
n σ
+
32 / 40
Lois de dépassements de seuil
Comment choisir le seuil u ?
Espérance de vie résiduelle (Redisual life):
Si Y = [X − u0 |X > u0 ] ∼ Hξ,σ(u ) , E(X − u|X
0
> u) est une fonction
linéaire en u, pour u > u0 ;
Nu
Les points {(u, doivent êtres alignés suivant une droite;
1
∑ (x(i) − u))}
Nu i=1
véri cation graphique, dif cile en pratique;
Exemple: intensité des précipitations journalières dans le Sud-Ouest de l'Angleterre
entre 1914 et 1962.
data(rain,package = "ismev")
mrlplot(rain)
33 / 40
Lois de dépassements de seuil
Alternative pour le choix du seuil stabilité des estimateurs σ(u), ξ
tcplot(rain, tlim=c(10,50),nt=30)
34 / 40
Loi des dépassements de seuil
Ajustement du modèle POT avec un seuil à 25.
fitpotrain fevd(rain,threshold=25,typ
summary(fitpotrain)
fevd(x = rain, threshold = 25, type = "GP")
[1] "Estimation Method used: MLE"
Negative Log-Likelihood Value: 900.6671
Estimated parameters:
scale shape
7.7018590 0.1077241
Standard Error Estimates:
scale shape
0.65928880 0.06222172
35 / 40
Plan
1. Introduction
2. Lois des maxima par blocs
3. Lois des dépassements de seuils
4. Cas de données dépendantes
5. Applications
36 / 40
Cas de données dépendantes
Souvent, les données présentent une dépendance temporelle;
Le théorème de Fischer-Tippet-Gnedenko est énnoncé dans le cas i. i. d ;
Extension sous certaines conditions :
Soient X1 , X2 , … une série stationnaire et X1∗ , X2∗ , … une série i.i.d de même loi
marginale que (Xi ). On note Mn = max{X1 , … , Xn } et
Mn = max{X , … , Xn }. Sous certaines conditions,
∗ ∗ ∗
1
∗
Mn − bn Mn − bn
I
P ( ≤ z) → G1 (z) ⟺ I
P ( ≤ z) → G2 (z)
an an
avec G2 (z) = G1 (z)
θ
pour 0 < θ ≤ 1 .
Pas de modi cation de la loi limite ( GEV );
Si G1 = GEV (μ, σ, ξ), alors G2 = GEV (μ − ;
σ −ξ ξ
(1 − θ ), σθ , ξ)
ξ
θ est appelé extremal index, et peut s'interpréter comme l'inverse de la taille des
clusters consécutifs d'extrêmes.
37 / 40
Cas de données dépendantes
Suivant le théorème précédent, pas de changement pour la loi des maxima par blocs: on n'a
accès qu'a l'observation de (Xi );
Pour les modèles de dépassements de seuil, deux approches:
incorporer la dépendance dans le modèle;
se ramener au cas précédent en faisant du declustering
1. Dé nir une règle pour identi er des clusters d'extrêmes;
2. identi er le maximum au sein de chaque cluster;
3. Ajuster une loi GP D à la loi des excès;
Ncluster
4. Estimer l'extremal index θ^ =
Nu
5. Le niveau de retour à m-observations est xm avec
σ ξ
= u + [(mζu θ) − 1]
ξ
ζu = I
P (X > u) .
38 / 40
Plan
1. Introduction
2. Lois des maxima par blocs
3. Lois des dépassements de seuils
4. Cas de données dépendantes
5. Applications
39 / 40
Applications
Cf TD
40 / 40
