Project_EV_OM

Installation librarie Valeur Extreme. Si ces librairies ne sont pas installés vous pouvez lancer le code suivant.

#evd Install ok
#install.packages(‘evd’) #evir des objets en conflits avec EVD #install.packages(‘evir’) #Install ok juste depend de mgcv. #install.packages(‘ismev’) #fExtremes des objets en conflits avec EVD #install.packages(‘fExtremes’) #extRemes des objets en conflits avec stats et EVD. #install.packages(‘extRemes’)

Chargement des Librairies

library(evd)
library(evir)

## 
## Attaching package: 'evir'

## The following objects are masked from 'package:evd':
## 
##     dgev, dgpd, pgev, pgpd, qgev, qgpd, rgev, rgpd

library(ismev)

## Loading required package: mgcv

## Loading required package: nlme

## This is mgcv 1.8-34. For overview type 'help("mgcv-package")'.

library(fExtremes)

## Loading required package: timeDate

## Loading required package: timeSeries

## Loading required package: fBasics

## Loading required package: fGarch

## 
## Attaching package: 'fExtremes'

## The following objects are masked from 'package:evir':
## 
##     dgev, dgpd, pgev, pgpd, qgev, qgpd, rgev, rgpd

## The following objects are masked from 'package:evd':
## 
##     dgev, dgpd, pgev, pgpd, qgev, qgpd, rgev, rgpd

library(extRemes)

## Loading required package: Lmoments

## Loading required package: distillery

## 
## Attaching package: 'extRemes'

## The following object is masked from 'package:evir':
## 
##     decluster

## The following objects are masked from 'package:evd':
## 
##     fbvpot, mrlplot

## The following objects are masked from 'package:stats':
## 
##     qqnorm, qqplot

library(fitdistrplus)

## Loading required package: MASS

## Loading required package: survival

library(TSstudio)
library(lubridate)

## 
## Attaching package: 'lubridate'

## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union

library(ggplot2)

## 
## Attaching package: 'ggplot2'

## The following object is masked from 'package:evir':
## 
##     qplot

Fin installation et chargements des librairies.

1 Chargement des données et Visualisation.

On charge les données de pluviométrie Marseille.

pluies = read.table("marseille.txt")[,1]
summary(pluies)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00    0.00    0.00   15.84    0.00 2215.00

pluies max est de 2215 ? Il y’a surement un problemes d’unite. Dans l’ennoncé les données sont en 0.1 mm Apres une recherche il semble que les données de pluie sont en mm cf https://fr.wikipedia.org/wiki/Pluviom%C3%A9trie Du coup transformons un peu le jeu de données .

pluies = pluies/10
summary(pluies)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.000   0.000   0.000   1.584   0.000 221.500

Le max est 221.5mm ce qui est plus en harmonie avec la “littérature”. Exemple https://www.infoclimat.fr/historic-details-evenement-737-orage-historique-de-marseille.html Ce chiffre de 221 correspond à 3 mois de précipitation.

pluies.ts = ts(pluies,start=c(1864,213),frequency=365) # 213: 1 aout 1864
matPluies = matrix(pluies,365,127) # decoupage en annees aout-juillet

On transforme le jeux de données.

maxPluies = apply(matPluies,2,max,na.rm = TRUE)
vecAn = 1864:1990

OK jusqu’ici tout va bien concernant le chargement des données. Des infos sur la série Temporelle.

ts_info(pluies.ts)

##  The pluies.ts series is a ts object with 1 variable and 46355 observations
##  Frequency: 365 
##  Start time: 1864 213 
##  End time: 1991 212

Et son graphe

ts_plot(pluies.ts,
        title = "Pluviométrie Marseille",
        Ytitle = "Thousands of Units",
        Xtitle = "Years",
        Xgrid = TRUE,
        Ygrid = TRUE)

On s’aperçoit qu’il y’a de nombreuses valeurs extremes notamment en 1886.Nous allons par la suite les identifier notamment en utilisant la methodologie bloxmaxima.

Mais avant étudions les jours sans pluie.

summary(pluies)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.000   0.000   0.000   1.584   0.000 221.500

#valeur max 221.5
sum(pluies>0)

## [1] 10370

#10370
sum(pluies==0)

## [1] 35985

#35985

72% des jours il ne pleut pas ,il fait plus beau à Marseille qu’a Paris .

Etude des quantiles

# On etudie le quantile 
quantile(pluies,probs=c(seq(.75,.95,.05)))

## 75% 80% 85% 90% 95% 
## 0.0 0.2 1.0 3.3 9.7

#On zoom
quantile(pluies,probs=c(seq(.90,0.99,.01)))

##  90%  91%  92%  93%  94%  95%  96%  97%  98%  99% 
##  3.3  4.0  5.0  6.2  7.7  9.7 12.1 15.8 20.9 30.8

Affichons l’empiral function survie

emplot(pluies)

## Warning in xy.coords(x, y, xlabel, ylabel, log): 35985 x values <= 0 omitted
## from logarithmic plot

Les données sont lognormés cela est du au grand nombre de jours sans pluie. Affichons la distribution

plot(sort(pluies), (1:length(pluies)/length(pluies)))

2 Methode GEV

2.1 Identifications des maximums

Dans un premier temps nous allons rechercher les maximums. Approche Block Maxima

Longueur de la chaine

str(pluies)

##  num [1:46355] 0 0 0 0 0 0 0 0 0 0 ...

length(pluies)

## [1] 46355

#On crée dataframe au bon format
Rain=data.frame(jour = seq(as.Date("1864-01-01"), by = "days",
                           length = length(pluies)),Prec=pluies)
annee=year(Rain$jour)
Rain=cbind(Rain,annee)
#Extraction des maximas
bmRain=blockmaxxer(Rain, blocks = Rain$annee, which="Prec")

Visualisation des maximas

plot(Rain$annee,Rain$Prec,xlab="Year",ylab="Prec",col = "blue",pch = 16)
points(bmRain$annee, bmRain$Prec, col="red", pch=16)

La méthode semble bien fonctionner le nuage de points rouge (maxima) #semble (en majorité) représenter les maximums de précipatations

Revenons à l’étude la méthode GEV.

2.2 Estimateur with ‘extremes’. Utilisons la fonction vu en cours

#FitRainMle=fevd(bmRain\(Prec,type="GEV",method="MLE") #FitRainMle #fevd(x = bmRain\)Prec, type = “GEV”, method = “MLE”)

On obtient un Negative Log-Likelihood Value: 1730.139 et un AIC = 3466.277 Ces valeurs me semblent élévés utilisons une fonction d’un autre package.

2.3 Estimateur with ‘ismev’. Aide fonction R gevFit(x, block = NA, type = c(“mle”, “pwm”), gumbel = FALSE, …)

Fit_GEV_MLE= gevFit(maxPluies,block=1,type = "mle")# Block à changer
par(mfcol = c(2, 2))
summary(Fit_GEV_MLE)

## 
## Title:
##  GEV Parameter Estimation 
## 
## Call:
##  gevFit(x = maxPluies, block = 1, type = "mle")
## 
## Estimation Type:
##   gev mle 
## 
## Estimated Parameters:
##         xi         mu       beta 
##  0.1003062 45.9285978 19.6778771 
## 
## Standard Deviations:
##          xi         mu       beta 
## 0.06082154 1.94555478 1.45327679 
## 
## Log-Likelihood Value:
##   585.7141 
## 
## Type of Convergence:
##   0 
## 
## Description
##   Wed Mar 09 23:39:04 2022

Voici les parametres obtenus: Estimated Parameters: xi mu beta 0.1003062 45.9285978 19.6778771 Nous obtenons de bien méilleurs résultats avec cette fonction (que celle utilisé avec Extremes) la Log-Likelihood Value: 585.7141 L’étude de la standard déviation et des graphes des résidus nous permettent de penser que le modèle est adapté.

On essaye avec PWM

Fit_GEV_MLE_Pwm= gevFit(maxPluies,block=1, type = "pwm")
par(mfcol = c(2, 2))
summary(Fit_GEV_MLE_Pwm)

## 
## Title:
##  GEV Parameter Estimation 
## 
## Call:
##  gevFit(x = maxPluies, block = 1, type = "pwm")
## 
## Estimation Type:
##   gev pwm 
## 
## Estimated Parameters:
##        xi        mu      beta 
##  0.123604 45.675999 19.269356 
## 
## Description
##   Wed Mar 09 23:39:04 2022

Voici les parametres obtenus: Estimated Parameters: xi mu beta 0.123604 45.675999 19.269356 L’etude des graphisques des résidus ne montrent pas de grandes différence significatives

On essaye avec gumbel pour la curiosité

Fit_GEV_MLE_Gum=gumbelFit(maxPluies,block=1,type = "mle")
par(mfcol = c(2, 2))
summary(Fit_GEV_MLE_Gum)

## 
## Title:
##  Gumbel Parameter Estimation 
## 
## Call:
##  gumbelFit(x = maxPluies, block = 1, type = "mle")
## 
## Estimation Type:
##   gum mle 
## 
## Estimated Parameters:
##       mu     beta 
## 47.05621 20.46544 
## 
## Standard Deviations:
##        mu     beta 
## 1.902995 1.444880 
## 
## Log-Likelihood Value:
##   587.4234 
## 
## Type of Convergence:
##   0 
## 
## Description
##   Wed Mar 09 23:39:04 2022

Pas de converge et de param xi donc pas de prévision avec Gumbel De plus l’étude des graphiques nous permet de conclure que le modèle n’est pas adapté.

Diagnostique et comparaison des 3 approches . Aide R Maximum-likelihood fitting for the generalized extreme value distribution, including generalized linear modelling of each paramete.

Fit_GEV_MLE2= gev.fit(maxPluies)

## $conv
## [1] 0
## 
## $nllh
## [1] 585.7141
## 
## $mle
## [1] 45.9274163 19.6773659  0.1003259
## 
## $se
## [1] 1.9455068 1.4532200 0.0608235

Fit_GEV_MLE_Pwm2 = gevFit(maxPluies, type = "pwm")
Fit_GEV_MLE_Gum2= gum.fit(maxPluies)

## $conv
## [1] 0
## 
## $nllh
## [1] 587.4234
## 
## $mle
## [1] 47.04623 20.47224
## 
## $se
## [1] 1.903635 1.446124

On lance un “diagnostic de la fonction”

gev.diag(Fit_GEV_MLE2)

Globalement le modele est adapté , il n’y a pas de points abérents pour les graphiques de probalilité et de densité. Concernant le graphique du Return level plot on a 2 points proches des limites et 1 points 100 ans qui franchit la limite. Pour le graphique du Quantile à partir de 120 les points dévient de plus en plus.

gum.diag(Fit_GEV_MLE_Gum2)

On voit bien que le modele avec gumbel fonctionne moins bien que la MLE.Le graphique du return level plot a un nombre plus important de points “déviants” que celui de la MLE.

Pour information on aurait pu utiliser un test de likehood Mais malheusement le test disponible avec le package “extremes” ne fonctionne pas avec les objets du package “isev”. Aprés recherche le package “isev” n’a pas d’équivalent de test.

Pour l’insant je n’ai pas réussi à adapter les fonctions. lr.test(Fit_GEV_MLE@data[[“x”]],Fit_GEV_MLE_Gumb@data[[“x”]])

A l’aide du summary on obtient les coefficient pour la prédiction.

summary(Fit_GEV_MLE2)

##       Length Class  Mode     
## trans   1    -none- logical  
## model   3    -none- list     
## link    1    -none- character
## conv    1    -none- numeric  
## nllh    1    -none- numeric  
## data  127    -none- numeric  
## mle     3    -none- numeric  
## cov     9    -none- numeric  
## se      3    -none- numeric  
## vals  381    -none- numeric

#Coeff MLE :
summary(Fit_GEV_MLE_Pwm2)

## 
## Title:
##  GEV Parameter Estimation 
## 
## Call:
##  gevFit(x = maxPluies, type = "pwm")
## 
## Estimation Type:
##   gev pwm 
## 
## Estimated Parameters:
##        xi        mu      beta 
##  0.123604 45.675999 19.269356 
## 
## Description
##   Wed Mar 09 23:39:04 2022

#Coeff PWM
summary(Fit_GEV_MLE_Gum2)

##       Length Class  Mode     
## trans   1    -none- logical  
## model   2    -none- list     
## link    2    -none- character
## conv    1    -none- numeric  
## nllh    1    -none- numeric  
## data  127    -none- numeric  
## mle     2    -none- numeric  
## cov     4    -none- numeric  
## se      2    -none- numeric  
## vals  254    -none- numeric

#Coeff Gumbel

Pour la suite nous ferons des prévions uniquement avec les méthodologies MLE et PMW.En effet on a vu que gumbel ne converge pas et que le modéle est le moins adaptés.

On fait prévison en utlisant les données auparavant Methode MLE

fExtremes::qgev( 1/100, xi = 0.1003062, mu = 45.9285978, beta = 19.6778771, lower.tail = F)

## [1] 160.9538
## attr(,"control")
##         xi      mu     beta lower.tail
##  0.1003062 45.9286 19.67788      FALSE

fExtremes::qgev( 1/1000, xi = 0.1003062, mu = 45.9285978, beta = 19.6778771, lower.tail = F)

## [1] 241.9864
## attr(,"control")
##         xi      mu     beta lower.tail
##  0.1003062 45.9286 19.67788      FALSE

Uniquement pour MLE on a une visualisation graphique

gevrlevelPlot(Fit_GEV_MLE,kBlocks=1000)

##                       min        v      max kBlocks
## GEV Return Level 187.8739 241.9864 357.1427    1000

Methode PMW

fExtremes::qgev( 1/100, xi = 0.123604, mu = 45.675999, beta = 19.269356, lower.tail = F)

## [1] 165.0591
## attr(,"control")
##        xi     mu     beta lower.tail
##  0.123604 45.676 19.26936      FALSE

fExtremes::qgev( 1/1000, xi = 0.123604, mu = 45.675999, beta = 19.269356, lower.tail = F)

## [1] 255.897
## attr(,"control")
##        xi     mu     beta lower.tail
##  0.123604 45.676 19.26936      FALSE

Le point d’attraction 100 et 1000 sont plus éléves avec PMW. Comparaison graphique Return level Plot.

RP_MLE = fExtremes::qgev( 1/seq(100,1000,50), xi = 0.1003062, mu = 45.9285978, beta = 19.6778771, lower.tail = F)
RP_PWM = fExtremes::qgev( 1/seq(100,1000,50), xi = 0.123604, mu = 45.675999, beta = 19.269356, lower.tail = F)

Visulation graphique

plot(RP_PWM~1/seq(100,1000,50), type = "l" , col = "Blue",xlab="Return Period Years",
     ylab='Return Level Rain', main="")
lines(RP_MLE~1/seq(100,1000,50), col = "Red")

Le return level plot est plus faible avec la méthode MLE, cela confime nos résultats précédents ainsi la MLE est plus adapté que la PMW.

3 La Méthode GPD

3.1 Recherche de Threshold. En utilisant le package evir

findthresh(pluies,c(1,10,100))

## [1] 169.0 100.9  56.7

#On définie U=56.7
pluies[pluies>56.7]

##   [1]  70.2  58.2  72.7  63.6  66.8  64.3  59.2  81.8  83.1 120.1  68.6 120.5
##  [13]  77.5  60.2 169.0  59.3  94.1  76.6  90.5  81.1  58.8 113.1 100.9  64.2
##  [25] 164.7  85.4 115.9 221.5  59.2 140.0  83.6  64.0 100.1  73.5  63.8  64.2
##  [37]  65.4  57.9  63.5  61.1  60.4  67.1  98.3  56.8  87.0  65.8  63.0  60.0
##  [49]  64.4 104.8  87.5  95.6  68.0  65.8  58.7  67.6 149.8  69.4  69.8  98.2
##  [61]  67.9  95.8  78.3  58.4  61.1  95.1  72.8  70.4  73.2  72.7  63.3  62.3
##  [73]  66.0  69.4  90.6  76.1  68.3  84.8  90.2  64.7  57.7  63.2  62.6  60.4
##  [85]  58.4  64.1  72.6  67.7  67.4  70.0  81.0  70.3  64.4  67.2  64.3  65.0
##  [97]  68.3  87.5  70.7  94.1

u=56.7
#On regarde  si ce choix de seuil est cohérent graphiquement.
plot(pluies,type="h",col="blue")
abline(h=56.7, col="red",lwd=2)

# Avec un calcul.
UnderU=1-length(pluies[pluies>u])/length(pluies)
#UnderU 0.997 ok .On valide notre choix du threshold.

Affichons le Mean Residual Life Plot #mrlplot(pluies) Celui du package isev plus sympa.

mrl.plot(pluies)

Graphique cohérent avec notre choix de de U.

3.2 Construction Estimateur

Fit_GPD_MLEO=fevd(pluies.ts,threshold=56.7,type="GP",method="MLE",time.units="days")
#summary(Fit_GPD_MLEO)

Probleme de Dimension @Jonathan Pourriez vous svp m’expliquer ultérieument ce message d’erreur ? Essayons une fonction autre package(“evir”) .

GPD MLE

Fit_GPD_MLE1= gpdFit(pluies.ts, u = 56.7, type = "mle")
summary(Fit_GPD_MLE1)

## 
## Title:
##  GPD Parameter Estimation 
## 
## Call:
## gpdFit(x = pluies.ts, u = 56.7, type = "mle")
## 
## Estimation Type:
##   gpd mle 
## 
## Estimated Parameters:
##         xi       beta 
##  0.1517246 18.6953210 
## 
## Standard Deviations:
##        xi      beta 
## 0.1097137 2.7602456 
## 
## Log-Likelihood Value:
##   407.9969 
## 
## Type of Convergence:
##   0

## 
## Description
##   Wed Mar 09 23:39:08 2022 by user: aidou

La Log-Likelihood Value est inferieur à celle obtenu avec la GEV. En vu des graphiques des résidus on peut penser que le modéle semble etre adapté.On a pas de valeurs aberrantes.

GPD PWM

Fit_GPD_PMW1 = gpdFit(pluies.ts, u = 56.7, type = "pwm")
summary(Fit_GPD_PMW1)

## 
## Title:
##  GPD Parameter Estimation 
## 
## Call:
## gpdFit(x = pluies.ts, u = 56.7, type = "pwm")
## 
## Estimation Type:
##   gpd pwm 
## 
## Estimated Parameters:
##         xi       beta 
##  0.1295549 19.1889613

## 
## Description
##   Wed Mar 09 23:39:08 2022 by user: aidou

Résultats globalement simpliaire avec MLE hormis quelques valeurs plus “déviantes”.

Représentation Graphique Tail Risk

Tail_GPD_MLE1 = gpdTailPlot(Fit_GPD_MLE1) 

Tail_GPD_PMW1 = gpdTailPlot(Fit_GPD_PMW1)

Pas de différence significative. Etudions les prévisions ce qui nous permettras surement de départager les 2 méthodologies.

Prevision MLE et PMW

RGPD_MLE1=gpdRiskMeasures(Fit_GPD_MLE1,p=1-1/(seq(100,1000,50)*365.25))
RGPD_PMW1=gpdRiskMeasures(Fit_GPD_PMW1,p=1-1/(seq(100,1000,50)*365.25))

Observons les quantiles

RGPD_MLE1$quantil

##  [1] 172.4943 187.6599 199.0000 208.1435 215.8474 222.5294 228.4452 233.7638
##  [9] 238.6026 243.0470 247.1609 250.9936 254.5838 257.9628 261.1557 264.1836
## [17] 267.0641 269.8118 272.4395

#PWM
RGPD_PMW1$quantile

##  [1] 169.3795 183.4453 193.8827 202.2509 209.2700 215.3353 220.6882 225.4872
##  [9] 229.8426 233.8340 237.5212 240.9500 244.1564 247.1693 250.0122 252.7043
## [17] 255.2621 257.6989 260.0266

plot(RGPD_MLE1$quantile~1/seq(100,1000,50), type = "l" , col = "blue",xlab="Return Period (années)",
     ylab='Return Level (mm de pluie)', main="")
lines(RGPD_PMW1$quantile~1/seq(100,1000,50), col = "red")

Le return level plot est plus faible avec la méthode PWM.

4 Comparaison des méthodes GEV VS GPD

Nous comparons les meilleurs méthodes celle MLE pour GEV et PWM pour la GPD

4.1 Graphique Return Level Plot GEV MLE vs GPD PWM

plot(RP_MLE~1/seq(100,1000,50), type = "l" , col = "Blue",xlab="Return Period Years",
     ylab='Return Level Rain', main="")
lines(RGPD_PMW1$quantile~1/seq(100,1000,50), col = "Red")

4.2 Previsions GEV MLE

fExtremes::qgev( 1/c(100,1000),  xi = 0.1003062, mu = 45.9285978, beta = 19.6778771, lower.tail = F)[1:2] 

## [1] 160.9538 241.9864

#160.9538 241.9864

GPD PWM

gpdRiskMeasures(Fit_GPD_PMW1,p=1-1/(c(100,1000)*365.25))

##           p quantile shortfall
## 1 0.9999726 169.3795  208.1955
## 2 0.9999973 260.0266  312.3343

#p quantile shortfall
#1 0.9999726 169.3795  208.1955
#2 0.9999973 260.0266  312.3343

Pour 100 Ans on a pour la GEV à 160.95 et la GPD à 169 la valeurs extreme sur le jeu de données pour un intevalle de 127 ans est 221 la GEV minimise le “risque”. Pour 1000 Ans on a pour la GEV à 241.96 et la GPD à 260.02 la valeurs extreme sur le jeu de données pour un intevalle de 127 ans est 221 donc la GEV semble plus approprié.

Conclusion Nous avons exploré les 2 méthodologies GEV et GPD et pour chacune conserver la meilleur approche. En comparant les 2 méthodes la GPD semble sous évaluer les valeurs extremes, Nous conservons donc celle fourni par la GPD en comparant avec un “climat” similaire La ville de Rossiglione le 4 octobre 2021 a enregistré une pluieviometrie de 740 mm(cf lien) #https://www.leprogres.fr/environnement/2021/10/06/un-metre-de-pluie-en-24h-cet-incroyable-record-qu-a-battu-l-italie Ainsi une GPD à 260 pour un horizon de 1000 ans me semble parfaitement justifié.

J’ai beaucoup apprecié travailler sur ce projet j’aurais aimé avoir plus de temps(mais d’autres projets m’appellent ..) afin de parcourir plus en profondeur Les packages R disponibles pour les calculs de valeurs extremes.