Category Archives: Extreme Value

MAT8886 Extremes and sums (of i.i.d. random variables)

20/01/2012 Arthur Charpentier Leave a comment

Yesterday, we have discussed briefly sums and maximas of i.i.d. random variables using the concept of subexponential distributions. Today, we will introduce the concept of regular variation: a positive function is said to be regularly varying (at infinity), denoted , for some , if

http://freakonometrics.blog.free.fr/public/perso5/subexp-33.gif

for all

http://freakonometrics.blog.free.fr/public/perso5/subexo_34.gif

. An this concept can be related to sums and maxima (see section 6.2.6 in Embrechts et al. (1997)). Consider i.i.d. positive random variables

http://freakonometrics.blog.free.fr/public/perso5/subsexp-01.gif

: let

http://freakonometrics.blog.free.fr/public/perso5/subexp-2.gif

and

http://freakonometrics.blog.free.fr/public/perso5/subexp-3.gif

. Then it can be shown easily that

if and only if

http://freakonometrics.blog.free.fr/public/perso5/subexp-10.gif

for some if and only if the exists a non-degenerate variable such that

http://freakonometrics.blog.free.fr/public/perso5/subexp-13.gif

with if and only if

http://freakonometrics.blog.free.fr/public/perso5/subexp-14.gif

If is not that simple to check for such convergences, it is still possible to use graphs to study the behavior of the empirical version of those quantities. Consider the following function to visualize convergence of empirical ratios,

CONVERGENCE=function(g,p=1,n=500000){
set.seed(1)
X=g(n);X1=g(n);X2=g(n);X3= g(n);X4=g(n)
Tp =cummax(X^p)/cumsum(X^p)
Tp1=cummax(X1^p)/cumsum(X1^p)
Tp2=cummax(X2^p)/cumsum(X2^p)
Tp3=cummax(X3^p)/cumsum(X3^p)
Tp4=cummax(X4^p)/cumsum(X4^p)
plot(Tp4,type="l",ylim=c(0,1),log="x",
xlim=c(100,n),ylab="",col="light blue",xlab="")
lines(Tp1,col="light green")
lines(Tp2,col="yellow")
lines(Tp3,col="pink")
lines(Tp,lwd=2)
abline(h=0:1,col="red",lty=2)
}

or the following to study the “asymptotic” distribution of the ratio on simulated samples

LIMITDIST=function(g,p=1,n=500000,ns=1000){
set.seed(1)
T=rep(NA,ns)
for(i in 1:ns){
X=g(n)
T[i]=max(X^p)/sum(X^p)
}
hist(T,breaks=seq(0,1,by=.05),probability=TRUE,
col="light green",ylab="",xlab="",main="")
}

In the case of exponentially distributed variables, we have

CONVERGENCE(rexp)

For variables with a lognormal distribution,

CONVERGENCE(rlnorm)

And finally, consider the case of a Pareto distribution

rpareto=function(n){runif(n)^(-1/1.5)-1}
CONVERGENCE(rpareto)

Here, it looks like those three distributions have finite variance (and actually, they do). To go one step further, for , define and . Then analogous results can be derived,

if and only if

http://freakonometrics.blog.free.fr/public/perso5/subexp-11.gif

for some if and only if the exists a non-degenerate variable such that

http://freakonometrics.blog.free.fr/public/perso5/subexp-12.gif

with if and only if

http://freakonometrics.blog.free.fr/public/perso5/subexp-15.gif

Again, it is possible to use the function defined above,

CONVERGENCE(rexp,p=2)

CONVERGENCE(rexp,p=3)

or even

CONVERGENCE(rexp,p=10)

If the power is not too high, it looks like the ratio goes to zero. But when it becomes larger, it looks like more simulations might be necessary to say something relevant.

CONVERGENCE(rlnorm,p=2)

CONVERGENCE(rlnorm,p=3)

Here also, it looks like we have a light tailed distribution (and actually, it is the case). And finally, if we consider the case of a Pareto distribution

CONVERGENCE(rpareto,p=2)

Then it looks like it is an heavy tailed distribution. In order to get a better understanding, plot the distribution of the ratio obtained from 1,000 simulated samples (of size 500,000),

LIMITDIST(rpareto,p=1)

versus

LIMITDIST(rpareto,p=2)

So obviously, something is going on between 1 and 2 (recall that the power parameter of the Pareto distribution is 1.5).

Extreme Value, MAT8886

Fisher-Tippett theorem with an historical perspective

18/01/2012 Arthur Charpentier 1 Comment

A couple of weeks ago, Rafael asked me if I had something on the history of extreme value theory. Since I will get back to fundamental results about extremes in my course, I promised I will write down a short post on all that issue.

To start from the beginning, in 1928, Ronald Fisher and Leonard Tippett formulated the three types of limiting distributions for the maximum term of a random sample (Fisher & Tippett (1928)). The problem was to characterize function such that

http://freakonometrics.hypotheses.org/files/2015/12/ext-2.gif

where where ‘s are i.i.d. with cumulative distribution function . They had supporting arguments, but no (rigorous) proof. Nevertheless, the obtained that the only possible types for G were

http://freakonometrics.hypotheses.org/files/2015/12/ext-6.gif

i.e. Fréchet type (Pareto-type tails), or

http://freakonometrics.hypotheses.org/files/2015/12/ext-7.gif

i.e. Weibull type (bounded distribution type), or

http://freakonometrics.hypotheses.org/files/2015/12/ext-8.gif

i.e. Gumbel type (exponential-type tails). Emil Gumbel has been intensively using the so-called Gumbel distribution on river flows, since (as he explained in 1958), “it seems that the rivers know the theory. It only remains to convince the engineers of the validity of this analysis“.
Independently of that work (published in 1928), Maurice Fréchet considered in 1927 (in Sur la loi de probabilité de l’écart maximum) possible limits of

http://freakonometrics.hypotheses.org/files/2015/12/ext-9.gif

and obtained only as possible limit. Richard von Mises gave in 1936 sufficient, but not necessary conditions for their (max) domain of attraction, i.e. characterization of function such that the maxima converges to some specific function (von Mises (1936)). E.g. he noticed that a sufficient condition on to be in the (max) domain of attraction of the Gumbel distribution is that

http://freakonometrics.hypotheses.org/files/2015/12/ext-13.gif

Then in 1943, Boris Gnedenko gave a complete characterization of those three types, with a complete characterization for two of them (heavy tails, i.e. Fréchet type and bounded support, i.e. Weibull) but his necessary and sufficient condition was based on a function that was not explicitly defined (see Gnedenko (1943)). Laurens de Haan in the 70’s derived checkable condition for Gumbel’s type.
Boris Gnedenko proved (in Section 4 of his paper) that F is the (max) domain of attraction of if and only if is regularly varying at infinity, with index (even if the term “regular variation” was not mentioned in the paper). Similar results were derived to characterize functions in the (max) domain of attraction of Weibull. For the (max) domain of attraction of , Boris Gnedenko obtained that a necessary and sufficient condition was that there exists a function such goes to 0 at infinity and

http://freakonometrics.hypotheses.org/files/2015/12/ext-20.gif

Several papers have discussed what function could be e.g. David Mejzler in 1949 (in Russian, but see also his 1965 paper), and Laurens de Hann in 1970 and 1971 (following the dramatic flood in the Netherlands in 1953, researchers in the Netherlands have focuses on dikes, and extreme value applications).

Mejzler’s idea was to work on quantiles, and not on the cumulative distribution function. I.e. define

http://freakonometrics.hypotheses.org/files/2015/12/ext-21.gif

Then a necessary and sufficient condition for F to be in the (max) domain of attraction of is that

http://freakonometrics.hypotheses.org/files/2015/12/ext-23.gif

Laurens de Haan proved in 1971 that function can be – in general – given by

http://freakonometrics.hypotheses.org/files/2015/12/ext-25.gif

And in 1976, Laurens de Haan obtained a three-type convergence working on quantile function (with a much shorter proof).
There have been many many papers extending Fisher-Tippett’s theorem, e.g. on non-independent sequences, like exchangeable ones (in a paper by Simeon Berman in 1962, or on stationary Gaussian sequences in 1964).

Extreme Value, MAT8886

Fisher-Tippett theorem and limiting distribution for the maximum

12/01/2012 Arthur Charpentier 1 Comment

Tomorrow, we will discuss Fisher-Tippett theorem. The idea is that there are only three possible limiting distributions for normalized versions of the maxima of i.i.d. samples . For bounded distribution, consider e.g. the uniform distribution on the unit interval, i.e. on the unit interval. Let and . Then, for all and ,

i.e. the limiting distribution of the maximum is Weibull’s.

set.seed(1)
s=1000000
n=100
M=matrix(runif(s),n,s/n)
V=apply(M,2,max)
bn=1
an=1/n
U=(V-bn)/an
hist(U,probability=TRUE,,col="light green",
xlim=c(-7,1),main="",breaks=seq(-20,10,by=.25))
u=seq(-10,0,by=.1)
v=exp(u)
lines(u,v,lwd=3,col="red")

For heavy tailed distribution, or Pareto-type tails, consider Pareto samples, with distribution function . Let and , then

which means that the limiting distribution is Fréchet’s.

set.seed(1)
s=1000000
n=100
M=matrix((runif(s))^(-1/2),n,s/n)
V=apply(M,2,max)
bn=0
an=n^(1/2)
U=(V-bn)/an
hist(U,probability=TRUE,col="light green",
xlim=c(0,7),main="",breaks=seq(0,max(U)+1,by=.25))
u=seq(0,10,by=.1)
v=dfrechet(u,shape=2)
lines(u,v,lwd=3,col="red")

For light tailed distribution, or exponential tails, consider e.g. a sample of exponentially distribution variates, with common distribution function . Let and , then

i.e. the limiting distribution for the maximum is Gumbel’s distribution.

library(evd)
set.seed(1)
s=1000000
n=100
M=matrix(rexp(s,1),n,s/n)
V=apply(M,2,max)
(bn=qexp(1-1/n))
log(n)
an=1
U=(V-bn)/an
hist(U,probability=TRUE,col="light green",
xlim=c(-2,7),ylim=c(0,.39),main="",breaks=seq(-5,15,by=.25))
u=seq(-5,15,by=.1)
v=dgumbel(u)
lines(u,v,lwd=3,col="red")

Consider now a Gaussian sample. We can use the following approximation of the cumulative distribution function (based on l’Hopital’s rule)

as . Let and . Then we can get

as . I.e. the limiting distribution of the maximum of a Gaussian sample is Gumbel’s. But what we do not see here is that for a Gaussian sample, the convergence is extremely slow, i.e., with 100 observations, we are still far away from Gumbel distribution,

and it is only slightly better with 1,000 observations,

set.seed(1)
s=10000000
n=1000
M=matrix(rnorm(s,0,1),n,s/n)
V=apply(M,2,max)
(bn=qnorm(1-1/n,0,1))
an=1/bn
U=(V-bn)/an
hist(U,probability=TRUE,col="light green",
xlim=c(-2,7),ylim=c(0,.39),main="",breaks=seq(-5,15,by=.25))
u=seq(-5,15,by=.1)
v=dgumbel(u)
lines(u,v,lwd=3,col="red")

Even worst, consider lognormal observations. In that case, recall that if we consider (increasing) transformation of variates, we are in the same domain of attraction. Hence, since , if

then

i.e. using Taylor’s approximation on the right term,

This gives us normalizing coefficients we should use here.

set.seed(1)
s=10000000
n=1000
M=matrix(rlnorm(s,0,1),n,s/n)
V=apply(M,2,max)
bn=exp(qnorm(1-1/n,0,1))
an=exp(qnorm(1-1/n,0,1))/(qnorm(1-1/n,0,1))
U=(V-bn)/an
hist(U,probability=TRUE,col="light green",
xlim=c(-2,7),ylim=c(0,.39),main="",breaks=seq(-5,40,by=.25))
u=seq(-5,15,by=.1)
v=dgumbel(u)
lines(u,v,lwd=3,col="red")

Computer, Extreme Value, Statistics

C’est agaçant cette manie de battre des records !

27/09/2011 Arthur Charpentier 4 Comments

En janvier dernier, j’avais fait un billet pour noter que 2010 était une année record pour les catastrophes naturelles, mais que 2009 aussi, et 2008 aussi. Bref, tous les ans, on bat des records. Mais il n’y a pas que sur le climat que l’on bat sans cesse des records. En sport aussi.
Par exemple, dimanche, le record du monde du marathon a (encore) été battu, avec 42,195 kilomètres parcourus (ou courus, tout simplement) en ~~2:03:58~~ 2:03:38 par Patrick Makau ~~(bon, certes, avec une amélioration d’une seconde par rapport au précédant record)~~. On devrait être content (en tous les cas Xi’an semble l’être). Sauf qu’à l’hiver 2012, je vais donner un cours sur copules et valeurs extrêmes à l’UQÀM, alias MAT8886. Et je pensais faire travailler les étudiants sur les articles de John Einmahl, sur les records au 100 mètres, ou plus généralement les records en athlétisme. Ces papiers proposent une très jolie application de la théorie des valeurs extrêmes, en particulier si on est dans le domaine d’attraction de la loi de Weibull (on peut alors chercher la borne supérieur – si on travaille sur le max – du support de la distribution, appelé ici endpoint). On apprend ainsi que selon la théorie des valeurs extrêmes, le temps minimal pour parcourir les 42,195 km du marathon devrait être de 2:04:06.

I.e. le temps pour courir un marathon devrait être borné par cette valeur. Sauf que ce record a été battu. Et à plusieurs reprises me semble-t-il.
Pour aller un peu plus loin, et comprendre ce qui se passait, j’ai été chercher les données sur wikipedia, pour le marathon de Boston, Chicago, Paris, Berlin, NYC, Stockholm, Fukuoka, Rotterdam, Amsterdam ou encore Londres. J’ai ainsi le temps du vainqueur pour les dernières années. Les données ont été copiées puis collées et donc un peu de remise en forme est ici nécessaire,

> base=read.table("http://freakonometrics.blog.free.fr/public
/data/topmarathon.csv",sep=";",header=TRUE)
> base=base[is.na(base$TIME)==FALSE,]
> base=base[(base$TIME=="")==FALSE,]
> n=nchar(as.character(base$YEAR))
> base$Y=as.numeric(substr(as.character(base$YEAR),n-3,n))
> base$Y[is.na(base$Y)]=1987
> h1=as.numeric(substr(as.character(base$TIME),1,1))
> h2=as.numeric(substr(as.character(base$TIME),1,2))
> i =which(is.na(t2))
> h=h2; h[i]=h1[i]
> m1=as.numeric(substr(as.character(base$TIME),3,4))
> m2=as.numeric(substr(as.character(base$TIME),4,5))
> m=m2; m[i]=m1[i]
> s1=as.numeric(substr(as.character(base$TIME),6,7))
> s2=as.numeric(substr(as.character(base$TIME),7,8))
> s=s2; s[i]=s1[i]
> base$T=h+m/60+s/60/60
> base=base[base$T>0,]
> base0=base[base$Y>=1982,]
> plot(base0$Y,base0$T,xlab="",ylab="")
> library(splines)
> reg=lm(T~bs(Y,6),data=base0)
> lines(1983:2010,predict(reg,newdata=
+ data.frame(Y=1983:2010)),col="red",lwd=2)

Si on regarde attentivement on observe plusieurs choses intéressantes. La première est qu’on ne court pas à la même vitesse dans toutes les villes. Par exemple à Stockholm, le temps du plus rapide est souvent bien éloigné des temps records. Mais surtout, le temps moyen du vainqueur ne cesse de baisser avec le temps,

http://freakonometrics.hypotheses.org/files/2015/12/marathon-cities_m.gif

On a alors une forme de non-stationnarité de la série qui laisse à penser qu’il conviendrait d’étudier davantage cette non-stationnarité avant d’utiliser la théorie des valeurs extrêmes. Et cette non-stationnarité et l’analyse des records n’est pas sans rappeler la discussion que l’on avait sur les catastrophes naturelles. Il serait peut-être temps de creuser davantage…

Computer, Extreme Value, Maths, Statistics

Tennis and risk management

16/12/2010 Arthur Charpentier 3 Comments

As mentioned already here, while we were going to Québec City for the workshop, we had interesting discussions in the car, and Maciej mentioned an article recently published in The Actuary,

Hence, I wanted to discuss (extremely) rare event probabilities in tennis. The story is simple: in June 2010, at Wimbledon, Nicolas Mahut and John Isner have played the longest match ever. 980 points, 11 But first of all, we need a dataset. Thanks to Duncan Murdoch, I have been able to run a short code to build up a dataset:

CITIES=c("berlin","madrid","paris","rolandgarros","wimbledon","sydney",
"beijing","shanghai","singapore","tokyo","melbourne","melbourne-indoor")
YEARS=1970:2009
BASE0=data.frame(YEAR=NA,TRNMT=NA,LENGTH=NA,SETS=NA)
for(i in 1:length(CITIES)){
for(j in 1:length(YEARS)){
city=CITIES[i]
year=YEARS[j]
localization = paste("http://www.resultsfromtennis.com/",
year,"/atp/",city,".html",sep="")
essai = try(readLines(localization), silent=TRUE)
ERROR404=FALSE
if(inherits(essai, "try-error")){ERROR404=TRUE}
if(ERROR404==FALSE){
B=scan(localization,"character")
SETS=NA
LENGTH=NA
if(length(B)>270){
I=(substr(B,1,10)=="class=rez>")
sum(I)
X0=B[I]
X3=as.numeric(substr(X0,11,13))
X2=as.numeric(substr(X0,11,12))
X1=as.numeric(substr(X0,11,11))
X0=X3
X0[is.na(X3)==TRUE]=X2[is.na(X3)==TRUE]
X0[is.na(X2)==TRUE]=X1[is.na(X2)==TRUE]
JL=c(which(substr(B,1,9)=="class=nl>"),length(B))
IL=which(substr(B,1,10)=="class=rez>")
IC=cut(IL,JL)
base=data.frame(IC,X0)
LENGTH=as.numeric(tapply(X0,IC,sum))
SETS=as.numeric(tapply(X0,IC,length))/2}
BASE=data.frame(YEAR=year,TRNMT=city,LENGTH,SETS)
BASE0=rbind(BASE0,BASE)}}}
write.table(BASE0,"BASE-TENNIS-TOTAL.txt")

Here I consider only tournaments where players have to win 3 sets (and actually more tournaments than those in the code above), and I have something like a bit more than 72,000 matches,

> I=is.na(TENNIS$LENGTH)==FALSE
> BT=TENNIS[I,]
> nrow(BT)
[1] 72754
> maxr=function(x){max(x,na.rm=TRUE)}
> T=paste(BT$TRNMT,BT$YEAR)
> DUREE=tapply(BT$SETS,T,maxr)
> LISTE=names(DUREE[DUREE>3])
> BT=BT[T%in%LISTE,]

so, if we look briefly at matches over 35 years, we have the following boxplot (one boxplot per year),

The red line being the epic Isner-Mahut match in June 2010 (4-6, 6-3, 7-6, 6-7, 70-68, i.e. 183 games, here for the score card).

If we study theory (e.g. from Paul Newton and Kamran Aslam), a lot of results can be obtained for the expected value of the number of games, but if we want to study extremely rare events, we should generate Markov chains (with a lot of generation since the probability should be extremely small). But how many ? Consider below matches with more than 50 games,

The tail plot (over 50), i.e. the log-log Pareto plot indicates that it will be difficult to study tails,

and similarly with the Hill plot (assuming that tails are Pareto type….)

Anyway, if we want to study tails, we should consider a threshold high enough. For instance, with a threshold at 68 (we keep only 24 match), we have

> seuil=68+0.25
> GPD1=gpd(X,seuil,method = "ml")
> GPD2=gpd(X,seuil,method = "pwm")
>
> xi=GPD1$par.ests[1]
> mu=seuil
> beta=GPD1$par.ests[2]
> x=180
> P=exp((-1/xi)*log(1 + (xi * (x - mu))/beta))
> as.numeric((1-GPD1$p.less.thresh)*P)
[1] 5.621281e-09
>
> xi=GPD2$par.ests[1]
> mu=seuil
> beta=GPD2$par.ests[2]
> x=180
> P=exp((-1/xi)*log(1 + (xi * (x - mu))/beta))
> as.numeric((1-GPD2$p.less.thresh)*P)
[1] 3.027095e-09

I.e. the probability that one match last more than 183 games is 1 chance over a billion… With, say, 2500 match per year, that gives us a return period of 400 years. So yes, we might say that this way a rare event… So perhaps, generating several billions of chains, it should be possible to get a more precise estimation of the probability to play 183 games in a single match…

Actuarial science, CIRM, Conferences, Copulas, Extreme Value, Publications

Lecture notes on risk and insurance

02/11/2010 Arthur Charpentier Leave a comment

I just finished some lectures notes on risk and insurance. The notes, that can be downloaded [pdf], are in French, and will be used at the JES (Journées d’Etudes Statistiques), organised at the CIRM (mentioned here). Previous notes on risk measures [pdf] and copulas [pdf]. Again, all comments are welcome…

Computer, Extreme Value, Tests

Detecting distributions with infinite mean

29/09/2010 Arthur Charpentier Leave a comment

In a post I published a few month ago (in French, here, based on some old paper,there), I mentioned a statistical procedure to test if the underlying distribution $https://perso.univ-rennes1.fr/arthur.charpentier/latex/paretoRO02.png$ of an i.i.d. sample $https://perso.univ-rennes1.fr/arthur.charpentier/latex/paretoRO01.png$ had a finite mean (based on extreme value results). Since I just used it on a small dataset (yes, with real data), I decided to post the R code, since it is rather simple to use. But instead of working on that dataset, let us see what happens on simulated samples. Consider $https://perso.univ-rennes1.fr/arthur.charpentier/latex/paretoRO03.png$ =200 observations generated from a Pareto distribution

$https://perso.univ-rennes1.fr/arthur.charpentier/latex/paretoRO04.png$

(here $https://perso.univ-rennes1.fr/arthur.charpentier/latex/paretoRO05.png$ =2, as a start)

> a=2
> X=runif(200)^(-1/a)

Here, we will use the package developped by Mathieu Ribatet,

> library(RFA)

Using Generalized Pareto Distribution (and LR test)

A first idea is to fit a GPD distribution on observations that exceed a threshold $https://perso.univ-rennes1.fr/arthur.charpentier/latex/paretoRO11.png$ >1.
Since we would like to test $https://perso.univ-rennes1.fr/arthur.charpentier/latex/paretoRO06.png$ (against the assumption that the expected value is finite, i.e. $https://perso.univ-rennes1.fr/arthur.charpentier/latex/paretoRO07.png$ ), it is natural to consider the likelihood ratio

$https://perso.univ-rennes1.fr/arthur.charpentier/latex/paretoRO08.png$

Under the null hypothesis, the distribution of $https://perso.univ-rennes1.fr/arthur.charpentier/latex/paretoRO09.png$ should be chi square distribution with one degree of freedom. As mentioned here, the significance level is attained with a higher accuracy by employing Bartlett correction (there). But let us make it as simple as possible for the blog, and use the chi-square distribution to derive the p-value.
Since it is rather difficult to select an appropriate threshold, it can be natural (as in Hill’s estimator) to consider $https://perso.univ-rennes1.fr/arthur.charpentier/latex/paretoRO10.png$ , and thus, to fit a GPD on the $https://perso.univ-rennes1.fr/arthur.charpentier/latex/paretoRO13.png$ largest values. And then to plot all that on a graph (like the Hill plot)

> Xs=rev(sort(X))
> s=0;G=rep(NA,length(Xs)-14);Gsd=G;LR=G;pLR=G
> for(i in length(X):15){
+ s=s+1
+ FG=fitgpd(X,Xs[i],method="mle")
+ FGc=fitgpd(X,Xs[i],method="mle",shape=1)
+ G[s]=FG$estimate[2]
+ Gsd[s]=FG$std.err[2]
+ FGc$fixed
+ LR[s]=FGc$deviance-FG$deviance
+ pLR[s]=1-pchisq(LR[s],df=1)
+ }

Here we keep the estimated value of the tail index, and the associated standard deviation of the estimator, to draw some confidence interval (assuming that the maximum likelihood estimate is Gaussian, which is correct only when n is extremely large). We also calculate the deviance of the model, the deviance of the constrained model ( $https://perso.univ-rennes1.fr/arthur.charpentier/latex/paretoRO06.png$ ), and the difference, which is the log likelihood ratio. Then we calculate the p-value (since under $https://perso.univ-rennes1.fr/arthur.charpentier/latex/paretoRO12.png$ the likelihood ratio statistics has a chi-square distribution).
If $https://perso.univ-rennes1.fr/arthur.charpentier/latex/paretoRO05.png$ =2, we have the following graph, with on top the p-value (which is almost null here), the estimation of the tail index he $https://perso.univ-rennes1.fr/arthur.charpentier/latex/paretoRO13.png$ largest values (and a confidence interval for the estimator),

If $https://perso.univ-rennes1.fr/arthur.charpentier/latex/paretoRO05.png$ =1.5 (finite mean, but infinite variance), we have

i.e. for those two models, we clearly reject the assumption of infinite mean (even if $https://perso.univ-rennes1.fr/arthur.charpentier/latex/paretoRO05.png$ gets too close from 1, we should consider thresholds large enough). On the other hand, if $https://perso.univ-rennes1.fr/arthur.charpentier/latex/paretoRO05.png$ =1 (i.e. infinite mean) we clearly accept the assumption of infinite mean (whatever the threshold),

Using Hill’s estimator

An alternative could be to use Hill’s estimator (with Alexander McNeil’s package). See here for more details on that estimator. The test is simply based on the confidence interval derived from the (asymptotic) normal distribution of Hill’s estimator,

> library(evir)
> Xs=rev(sort(X))
> HILL=hill(X)
> G=rev(HILL$y)
> Gsd=rev(G/sqrt(HILL$x))
> pLR=1-pnorm(rep(1,length(G)),mean=G,sd=Gsd)

Again, if $https://perso.univ-rennes1.fr/arthur.charpentier/latex/paretoRO05.png$ =2, we clearly rejct the assumption of infinite mean,

and similarly, if $https://perso.univ-rennes1.fr/arthur.charpentier/latex/paretoRO05.png$ =1.5 (finite mean, but infinite variance)

Here the test is more robust than the one based on the GPD. And if $https://perso.univ-rennes1.fr/arthur.charpentier/latex/paretoRO05.png$ =1 (i.e. infinite mean), again we accept $https://perso.univ-rennes1.fr/arthur.charpentier/latex/paretoRO12.png$ ,

Note that if $https://perso.univ-rennes1.fr/arthur.charpentier/latex/paretoRO05.png$ =0.7, it is still possible to run the procedure, and hopefully, it suggests that the underlying distribution has infinite mean,

(here without any doubt). Now you need to wait a few days to see some practical applications of the idea (there was on in the paper mentioned above actually, on business interruption insurance losses).

Computer, Copulas, Extreme Value, Statistics

Tails of copulas, une lecture graphique

08/07/2010 Arthur Charpentier 3 Comments

Suite à une formation que je faisais en fin de semaine à Brest (les slides sont ici et là), je voulais revenir sur les histoires de tails of copulas, pour reprendre le titre de l’article (ici) de Gary Venter (et qui correspond à des choses que j’avais pu présenter il y a quelques années à Berlin, les slides étant en ligne ici).

Slides astin from Arthur Charpentier

Quantifier la dépendance de queue

L’idée est de noter qu’il est noter qu’il existe deux manières de quantifier la dépendance de queue. La première est liée à l’approche de Joe (1990, ici, ou 1997 pour le livre), qui a introduit un (strong) tail dependence index. Par exemple pour la queue inférieure,

https://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso3/toc3latex2png.2.php.png

soit

https://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso3/toc3latex2png.3.php.png

La seconde est liée à une idée que l’on retrouve dans les travaux de Janet Heffernan, Stuart Coles ou Jonathan Tawn. L’intuition est la suivante (on peut la retrouver en ligne ici). Si $https://perso.univ-rennes1.fr/arthur.charpentier/latex/toclatex2png-2.2.php.png$ et $https://perso.univ-rennes1.fr/arthur.charpentier/latex/toclatex2png-3.2.php.png$ ont la même loi et que l’on suppose les variables indépendantes, alors

$https://perso.univ-rennes1.fr/arthur.charpentier/latex/toclatex2png-1.2.php.png$

En revanche, si les variables sont comonotones (c’est à dire égales comme on suppose les lois identiques),

$https://perso.univ-rennes1.fr/arthur.charpentier/latex/toclatex2png-4.2.php.png$

Aussi, on peut supposer qu’il existe un indice $https://perso.univ-rennes1.fr/arthur.charpentier/latex/toclatex2png-6.2.php.png$ tel que

$https://perso.univ-rennes1.fr/arthur.charpentier/latex/toclatex2png-5.2.php.png$

Le soucis est que le cas d’indépendance correspond à $https://perso.univ-rennes1.fr/arthur.charpentier/latex/toclatex2png-6.2.php.png$ =2, alors que le cas de dépendance forte correspond au cas $https://perso.univ-rennes1.fr/arthur.charpentier/latex/toclatex2png-6.2.php.png$ =1. Il est alors usuel de faire une transformation affine pour se ramener sur [0,1], et que la force de la dépendance soit croissante avec $https://perso.univ-rennes1.fr/arthur.charpentier/latex/toclatex2png-6.2.php.png$ , e.g.

$https://perso.univ-rennes1.fr/arthur.charpentier/latex/toclatex2png-8.2.php.png$

Posons alors

$https://perso.univ-rennes1.fr/arthur.charpentier/latex/toc2latex2png.2.php.png$

qui pourra être interprété comme un (weak) tail dependence index.
Bref, ces deux mesures donnent de l’information sur le comportement dans les queues de distribution.

Les fonctions de concentration dans les queues

L’idée est de noter qu’il est possible d’étudier ces fonctions afin de mieux comprendre le comportement dans les queues. En s’inspirant de Gary Venter, on peut définir

https://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso3/Llatex2png.2.php.png

pour étudier le comportement dans la queue inférieure, et

pour la queue supérieure,où $https://perso.univ-rennes1.fr/arthur.charpentier/latex/toclatex2png-12.2.php.png$ est la copule de survie associée à $https://perso.univ-rennes1.fr/arthur.charpentier/latex/toclatex2png-13.2.php.png$ , au sens où

$https://perso.univ-rennes1.fr/arthur.charpentier/latex/toclatex2png-14.2.php.png$

$https://perso.univ-rennes1.fr/arthur.charpentier/latex/toclatex2png-15.2.php.png$

Cet outil permettra de modéliser la dépendance forte. On peut également poser, afin d’étudier la dépendance faible,

$https://perso.univ-rennes1.fr/arthur.charpentier/latex/toc2latex2png.3.php.png$

$https://perso.univ-rennes1.fr/arthur.charpentier/latex/toc2latex2png.4.php.png$

Application statistique

L’idée est de noter qu’il est facile d’estimer ces fonctions. Ces outils peuvent être utiles pour mieux comprendre le comportement dans les queues.
Par exemple pour une copule Gaussienne de corrélation 0,5, on a la forme théorique suivante pour les fonctions de concentration (au sens fort)

Statistiquement, il est possible d’estimer ces quantités en comptant simplement le nombre d’observations dans le coin inférieur gauche, ou le coin supérieur droit. Si on dispose d’un échantillon, on peut alors regarder ce que donnent les versions

$http://perso.univ-re<br /><br /> nnes1.fr/arthur.charpentier/latex/toclatex2png-18.2.php.png$

$https://perso.univ-rennes1.fr/arthur.charpentier/latex/toclatex2png-19.2.php.png$

Pour un échantillon de taille n=500, on obtient les intervalles de confiance à 90% de la forme suivante,

Le code R ressemble à ça

> library(evd); data(lossalae)
> cor(lossalae,method="spearman")
         Loss     ALAE
Loss 1.000000 0.451872
ALAE 0.451872 1.000000

avec le code suivant pour la version empirique,

> z=seq(0,.5,by=.001)
> U=rank(v[,1])/(nrow(v)+1)
> V=rank(v[,2])/(nrow(v)+1)
> Lemp=rep(NA,length(z))
> Remp=rep(NA,length(z))
> for(i in 1:length(z)){
+  Lemp[i]=sum((U<=z[i])&(V<=z[i]))/sum(U<=z[i])
+  Remp[i]=sum((U>=1-z[i])&(V>=1-z[i]))/sum(U<=z[i])
+ }

et pour la version théorique,

> Lg=(pcopula(copclayton,cbind(z,z)))/(z)
> Rg=((1-2*(1-z)+pcopula(copclayton,cbind(1-z,1-z))))/(z)
> plot(c(1-z,z),c(Lg,Rg))

De plus, on a des fonctions similaires pour la dépendance au sens faible, avec le code suivant pour la version théorique,

> Lg=log(pcopula(cop,cbind(z,z)))/log(z)
> Rg=log((1-2*(1-z)+pcopula(cop,cbind(1-z,1-z))))/log(z)
> Lg=1/Lg*2-1
> Rg=1/Rg*2-1

et celui là pour la version empirique

> z=seq(0,.5,by=.001)
> v <- lossalae
> U=rank(v[,1])/(nrow(v)+1)
> V=rank(v[,2])/(nrow(v)+1)
> Lemp=rep(NA,length(z))
> Remp=rep(NA,length(z))
> for(i in 1:length(z)){
+  Lemp[i]=log(mean((U<=z[i])&(V<=z[i])))/log(mean(U<=z[i]))
+  Remp[i]=log(mean((U>=1-z[i])&(V>=1-z[i])))/log(mean(U<=z[i]))
+ }
> Lemp=1/Lemp*2-1
> Remp=1/Remp*2-1

Bref, on peut utiliser ces fonctions sur des vrais échantillons. Considérons l’exemple classique loss-alae (où l’on couple les frais dans des sinistres assurés, et les frais payés par l’assureur). On souhaite ajuster une copule, sans trop savoir laquelle. On peut commencer par étudier la dépendance forte, et comparer avec une copule Gaussienne. La copule Gaussienne de référence possède ici le même rho de Spearman que l’échantillon dont on dispose,

> cor(lossalae,method="spearman")
         Loss     ALAE
Loss 1.000000 0.451872
ALAE 0.451872 1.000000
> library(copula)
> paramgauss=.47
> paramclayton=.9
> paramgumbel=1.45
> copgauss=normalCopula(paramgauss)
> copclayton=claytonCopula(paramclayton, dim = 2)
> copgumbel=gumbelCopula(paramgumbel, dim = 2)

On obtient ici

La courbe verte est l’intervalle de confiance (ponctuel) à 95% pour une copule Gaussienne et un échantillon de même taille. On voit qu’on modélise mal la structure de dépendance. Avec une copule duale de Clayton, on obtient

et enfin pour une copule de Gumbel,

Bref, la copule de Gumbel semble réellement bien adaptée… Si on creuse en étudiant la dépendance au sens faible, on peut valider là aussi ce modèle. En effet, si la référence est la copule Gaussienne,

ou pour une copule de Clayton,

alors qu’une copule de Gumbel donnerait

Extreme Value, Statistics

Some historical remarks on extreme values

12/02/2010 Arthur Charpentier 1 Comment

I will start here a short post on extreme values, with some historical perspective. In a recent paper (in French), I mentioned the use of the Pareto distribution as a standard model for extremes, but if reinsurers have been using the Pareto distribution for a long time (see here e.g.), the oldest mathematical models when dealing with extreme value should be related to work on maximum values in finite samples.

The work of Ronald Fisher and Leonard Tippett

Leonard Henry Tippett, a former student of Karl Pearson published in Biometrika a note on extremes, in 1925. The goal was “the determination of the distribution of the range and the extremes for a large number of samples“. In 1925, everyone was looking for the Gaussian distribution everywhere, and Leonard Tippett observed that the distribution of the largest value did not have a Gaussian distribution.
A few years after, a joint work with Ronald Fisher was presented to the Cambridge Philosophical Society. The starting point was the idea of “stability” (even if the term did not appear explicitely in their work): the limiting distribution the maximum should be of the “same type” as the underlying distribution. Thus, if $https://perso.univ-rennes1.fr/arthur.charpentier/latex/ext-01.png$ stands for the cumulative distribution function, it should satisfy functional equation

$https://perso.univ-rennes1.fr/arthur.charpentier/latex/ext-02.png$

Solutions of that functional equation will give all possible limiting distributions. Thus, Fisher and Tippett obtained three possible limits,

solutions of $https://perso.univ-rennes1.fr/arthur.charpentier/latex/ext-03.png$ , i.e. $https://perso.univ-rennes1.fr/arthur.charpentier/latex/ext-04.png$
solutions of $https://perso.univ-rennes1.fr/arthur.charpentier/latex/ext-05.png$ , i.e. $https://perso.univ-rennes1.fr/arthur.charpentier/latex/ext-06.png$ with $https://perso.univ-rennes1.fr/arthur.charpentier/latex/ext-07.png$ (i.e. finite lower bound for the support), i.e. $https://perso.univ-rennes1.fr/arthur.charpentier/latex/ext-08.png$
solutions of $https://perso.univ-rennes1.fr/arthur.charpentier/latex/ext-05.png$ , i.e. $https://perso.univ-rennes1.fr/arthur.charpentier/latex/ext-10.png$ if $https://perso.univ-rennes1.fr/arthur.charpentier/latex/ext-11.png$ (i.e. finite upper bound for the support), i.e. $https://perso.univ-rennes1.fr/arthur.charpentier/latex/ext-09.png$

Based on those possible limiting distributions, Fisher and Tippett wanted to derive what has been called later on the “domain of attraction” of those distributions.

The work of Maurice Fréchet, at the same time

In 1926, Maurice Fréchet wrote a paper on “la loi de probabilité de l’écart maximum“. That paper, as well as the one by Fisher and Tippett (wrote at the same time), investigated asymptotic limits. Both obtained functional equations, but only Maurice Fréchet understood the importance of the stability concept, pointed out by Paul Levy in the context of sums. Thus, Maurice Fréchet introduced the concept of what is called now “max-stability“. But Fréchet solve only functional equation $https://perso.univ-rennes1.fr/arthur.charpentier/latex/ext-05.png$ . The point is that Fréchet studied absolute values of errors, i.e. strictly positive random variables. Thus, Maurice Fréchet considered distribution

$https://perso.univ-rennes1.fr/arthur.charpentier/latex/ext-12.png$

where $https://perso.univ-rennes1.fr/arthur.charpentier/latex/ext-92.png$ is an arbitrary positive constant. The “2” comes from the fact that Fréchet considered errors with respect to the median. But he did not introduced that new distribution function, he also proved that the distribution appears as a limit when the underlying distribution of the $https://perso.univ-rennes1.fr/arthur.charpentier/latex/ext-13.png$ ‘s has an algebraic behavior at infinity, i.e. equivalent to $https://perso.univ-rennes1.fr/arthur.charpentier/latex/ext-90.png$ , for some $https://perso.univ-rennes1.fr/arthur.charpentier/latex/ext-91.png$ . I.e. he proved that Pareto-type tailed distibutions where in the domain of attraction of the Fréchet distribution.

Later on, the work of Emil Gumbel

In 1932, Emil Gumbel gave a talk in France on the “âge limite“. But as he wrote it “on peut donc supposer que la distribution de l’âge limite – c’est à dire la probabilité que la probabilité de cet âge ait une valeur donnée – soit Gaussienne“. But a few years after, he read about Fisher’s work, and observed also that “la distribution d’une valeur extrêmes peut être représentée pour un nombre suffisant d’observations par la formule doublement exponentielle, pourvu que la distribution initiale se comporte asymptotiquement comme une exponentielle. La formule devient rigoureuse si la distribution initiale est exponentielle“, as he wrote in 1935. Thus, as Fréchet proved that Pareto type distribution were in the max-domain of attraction of Fréchet’s distribution, Gumbel obtained that exponential type distributions were in the max-domain of attraction of Gumbel’s distribution. He also introduced the term “distribution de type exponentiel”
For Emil Gumbel, it was natural to study the logarithmic derivative of the distribution, since it is the mortality rate in demography (area that Emil Gumbel studied previously). As he mentioned “d’un point de vue théorique, il est intéressant de noter que M. Fréchet a construit une distribution initiale d”une variable aléatoire pour laquelle la valeur absolue de la dérivée logarithmique diminue sans limite“. But since it was not a valuable property for practical applications, he decided that “nous nous bornerons au traitement des données de type exponentiel“. Emil Gumbel always tried to relate his work on extremes and what he did on demograpy.
For instance in 1937, he wrote a paper on “les centennaires” that can also be related to the work of Bortkiewicz on rare events. He also applied his work on radioactivity, and hydrology.
In the 30’s, hydrographs as Hazen or Graszberger introduced the concept of “yearly maximum” of
a river level. They actually proposed to look for actuarial models to study decennial or centennial floods. But they only used the lognormal distribution to model yearly maxima. In 1936, French hydrologist Aimé Coutagne met Emil Gumbel (who was teaching at the ISFA, in Lyon). At that time, Emil Gumbel was looking for possible applications (outside demography) for his doubly exponential distribution. As as pointed out by Aimé, “sa formule devait être applicable au cas des crues; c’est à dire des plus grands débits, problème analogue à celui des plus grands âges“. Not only Gumbel’s distribution gave better empirical results, but also it came with a theoritical justification.

Gumbel’s distribution properties

Consider the Gumbel distribution, with location and scale parameters alpha and beta respectively, i.e.

$https://perso.univ-rennes1.fr/arthur.charpentier/latex/ext-40.png$

Note that the associated quantile function is

$https://perso.univ-rennes1.fr/arthur.charpentier/latex/ext-41.png$

with mean

$https://perso.univ-rennes1.fr/arthur.charpentier/latex/ext-43.png$

and variance

$https://perso.univ-rennes1.fr/arthur.charpentier/latex/ext-44.png$

The work of Waloddi Weibull

Waloddi Weibull, a Swedish physict proposed a distribution in 1939, to represent the distribution of breaking strength of materials. He used it in the 50’s in reliability concept. Actually, Weibull appeared late in the story of extremes, since Fréchet, Fisher and Tippett mentioned it already in the mid-20’s.

From the central limit theorem (on the average) to Fisher-Tippett theorem (on the maxima)

In order to visualize those two theorem, consider the following animation, where samples of 20 exponential variables are generated. From those 20 values, we plot the maximum in blue, and the average in red, on top. Just below, be rescale those points by considering $https://perso.univ-rennes1.fr/arthur.charpentier/latex/ext-16.png$ , and below again, $https://perso.univ-rennes1.fr/arthur.charpentier/latex/ext-15.png$ }. When then look at the position of $https://perso.univ-rennes1.fr/arthur.charpentier/latex/ext-14.png$ and the one of the mean of $https://perso.univ-rennes1.fr/arthur.charpentier/latex/ext-15.png$ . We then build an histogram to visualize the distribution of the rescaled maximum (in blue) and the rescale average (in red).

For those who might be busy, after 1000 generations of samples, we obtain the following histograms (below), including the Gaussian distribution below (i.e. the average of exponential variables looks Gaussian, even with only 20 observations, actually the Gaussian distribution is only asymptotic, i.e. we should consider samples of size 2000), and the maximum over 20 observations of exponential variables (on top) looks like a Gumbel distribution (actually, here it is the exact distribution, and it is the asymptotic distribution for exponential type variables).

The GEV distribution

The unified expression of those three distributions is call the GEV distribution. The generalized extreme value distribution has cumulative distribution function

$https://perso.univ-rennes1.fr/arthur.charpentier/latex/ext-20.png$

for $https://perso.univ-rennes1.fr/arthur.charpentier/latex/ext-21.png$ , where $https://perso.univ-rennes1.fr/arthur.charpentier/latex/ext-22.png$ is the location parameter, $https://perso.univ-rennes1.fr/arthur.charpentier/latex/ext-23.png$ the scale parameter and $https://perso.univ-rennes1.fr/arthur.charpentier/latex/ext-24.png$ the shape parameter. Note that the expected value is
$https://perso.univ-rennes1.fr/arthur.charpentier/latex/ext-30.png$

Actuarial science, Computer, Extreme Value, Publications, Reinsurance, Research, Statistics

Les modèles en réassurance

02/02/2010 Arthur Charpentier 2 Comments

Publication d’un papier sur les modèles en réassurance dans la revue Risques (ici), suite à un questionnement sur la pertinence des modèles classiques utilisés par les réassureurs. La question initiale était partie de la constatation – que l’on retrouve ici ou là – sur l’utilisation (ou la mauvaise utilisation) de modèles sophistiques en finance de marché, en essayant d’expliquer que – d’un point de vue épistémologique au moins – les modèles des réassureurs étaient plus robustes. En particulier, on notera que les modèles les plus anciens utilisés par les réassureurs (en particulier la loi de Pareto comme je l’avais évoqué ici) ont eu une légitimité pratique durant plusieurs décennies avant d’être justifiés par la théorie des valeurs extrêmes. Je ferais d’ailleurs bientôt un billet sur l’histoire des valeurs extrêmes en statistiques, en revenant en particulier sur les travaux de Gumbel ou de Fréchet.

Sinon le code utilisé dans le papier est en ligne ici, et la base là

> sinpe = read.table("https://perso.univ-rennes1.fr/arthur.charpentier/sinpe.csv",header=TRUE,sep=";")
> head(sinpe)
      DSUR  MNTPE
1 19850206 240439
2 19851228 125674
3 19850504 488331
4 19851118 457347
5 19850220 990919
6 19851214 182939
> annee=as.numeric(substr(as.character(sinpe$DSUR),1,4))
> sinistres=sinpe$MNTPE[annee>1992]
> XS=sinistres/100000

On se limite ici aux sinistres survenus après 1992. Si l’on visualise ces montants de sinistres, on obtient

> datesur=as.Date(as.character(sinpe$DSUR),"%Y%m%d")
> jour=datesur[annee>1992]
> plot(jour,sinistres/100000,xlab="",ylab="Coût individuel",cex=.5,ylim=c(0,600))
> ded=50
> abline(h=ded)

On peut aussi faire le graphique Pareto-log-log,

> library(evir)
> n=length(X)
> plot(log(sort(X)),log((n:1)/(n+1)),
+ xlab="Coûts des sinistres (logarithme)",ylab="Fonction de survie (logarithme)",cex=.8)
> out <- gpd(XS, 15)
> XI=as.numeric(out$par.ests[1]); BETA=as.numeric(out$par.ests[2]) 
> x0=seq(2,8,.01)
> lines(x,-1/XI*(x-log(15)),col="red")

(l’ajustement de la loi de Pareto permettant de tracer la droite rouge) ou encore visualiser l’estimateur de Hill de l’indice de queue,

> hill(X)

Pour finir, le code suivant permet de calculer la prime pure (ou plus généralement une prime de Wang) soit de manière non-paramétrique (burning cost) ou bien en utilisant l’ajustement d’une loi de Pareto. Le papier étant un papier de vulgarisation, j’ai pris le seuil de manière arbitraire, sans aucune recherche de valeur “optimale”

> DEDUC = seq(10,50,by=5)
> lambda=0;  seuil=5
> WG1=WG2=rep(NA,length(DEDUC))
> for(k in 1:length(DEDUC)){
+ deductible=DEDUC[k]
+ out <- gpd(XS, seuil)
+ XI=as.numeric(out$par.ests[1]); BETA=as.numeric(out$par.ests[2]) 
+ G0=function(x){1-pgpd(x+seuil, xi = XI, mu = seuil, beta = BETA)}
+ G=function(x){(G0(x+deductible-seuil))/(G0(deductible-seuil))}
+ F=function(x){pnorm(qnorm(G(x))+lambda)}
+ (wang1=integrate(F, 0, Inf))
+ X=XS[XS>deductible]
+ n=length(X)
+ FS= function(z){
+ m=rep(NA,length(z))
+ for(i in 1:length(m)){
+ m[i]=sum(X>z[i]+deductible)/n}
+ return(m)
+ }
+ G=function(x){pnorm(qnorm(FS(x))+lambda)}
+ (wang2=sum(G(seq(0,800,.01))*.01))
+ WG1[k]=as.numeric(wang1$value)
+ WG2[k]=wang2
+ }
> plot(DEDUC,WG2-DEDUC,type='b',xlab="Niveau de la priorité ('00 000 euros)",ylab="Prime pure par sinistres réassuré",ylim=c(0,50))
> lines(DEDUC,WG1-DEDUC,type='b',col="red",pch=4)
> legend(10,50,c("Aujustement d'une loi de Pareto","'Burning cost'"),
+ col=c("red","black"),lty=1,cex=.8,pch=c(4,1))

Actuarial science, Extreme Value, Publications, Reinsurance, Risk, Statistics, Tests

Peut on faire l’économie du formalisme quand on parle d’extrêmes ?

19/11/2009 Arthur Charpentier Leave a comment

Tous les blogs économiques saluent la parution en poche du joli petit livre de Daniel Zajdenweber, Economie des Extrêmes. En particulier, beaucoup de monde salue ce livre qui explique simplement des choses complexes…. Par exemple Alexandre dès 2001 “Passé le premier chapitre, un peu ardu, et qui nécessite du lecteur des connaissances de base en statistique et probabilités (notion de lois de probabilité, d’espérance, de variance…) qui décrit en termes littéraires les caractéristiques de ces lois, l’auteur applique ces résultats à un grand nombre de phénomènes concrets, et en tire les conséquences“. Mais peut-on parler d’économie des extrêmes sans être technique ?

Histoire que mon message ne soit pas déformé, je trouve passionnant ce petit livre introductif à la problématique des risques extrêmes (qui est un de mes dadas depuis quelques années) mais j’espère qu’il servira d’encouragement à une lecture d’ouvrages plus détaillés sur le sujet. Car la vulgarisation a des limites que l’on atteint vite quand on parle de sujets aussi complexes.

L’exemple que j’ai le plus étudié est celui des sinistres de perte d’exploitation (longuement évoqué par Daniel Zajdenweber dans son livre). Il y a quelques années j’avais utilisé cette partie du livre comme base pour faire un sujet d’examen pour le cours de “réassurance et grands risques” que je donnais alors à l’ENSAE¹. Et malheureusement, mes compétences littéraires sont très limitées, donc je vais faire des maths. Dans le livre, la figure suivante est présentée,

qui correspond effectivement à la fonction tracée dès 1925 par Karl Gustav Hagstroem (j’avais souligné (ici) ses travaux précurseurs où l’intérêt de la loi de Pareto pour modéliser les très grands riques apparaissait pour la première fois). C’est en effet assez naturel: si on a une loi de Pareto, i.e.

$https://perso.univ-rennes1.fr/arthur.charpentier/latex/z01.png$

alors on pourrait écrire, en passant au logarithme

$https://perso.univ-rennes1.fr/arthur.charpentier/latex/z02.png$

Si on représente la version empirique, c’est à dire le nuage de points

$https://perso.univ-rennes1.fr/arthur.charpentier/latex/z03.png$

alors pour une loi de Pareto, les points devraient être alignés suivant une droite, et la pente doit correspondre au paramètre de la fonction puissance. C’est visiblement l’idée exploitée ici.

Autrement dit, les pointillés ne sont un intervalle de confiance, mais juste un outils graphique pour se demander si la pente vaut 1, ou pas. Daniel Zajdenweber affirme que la pente doit ici être -1.

Le fait que la valeur soit unitaire ou pas a en effet un impact très important en terme d’assurabilité du risque de perte d’exploitation. Rappelons que pour une variable positive (et c’est le cas ici). Et si on a une telle loi de Pareto (de puissance unitaire), alors la prime pure d’un traité de réassurance, couvrant entre m et M s’écrit

$https://perso.univ-rennes1.fr/arthur.charpentier/latex/z06.png$

soit

$https://perso.univ-rennes1.fr/arthur.charpentier/latex/z07.png$

ce qui correspond aux calculs de Daniel Zajdenweber… Mais encore une fois “l’absence d’espérance mathématique de la distribution des sinistres” est une conclusion très forte sur laquelle on peut essayer de revenir.

$https://perso.univ-rennes1.fr/arthur.charpentier/latex/z04.png$

aussi ici

$https://perso.univ-rennes1.fr/arthur.charpentier/latex/z05.png$

autrement dit, l’espérance est finie si la pente est strictement plus grande (en valeur absolue) que 1. Si la pente est inférieure (ou égale) à 1, le risque n’est pas assurable ! Ce qui est une conclusion très très forte pour les assureurs.
J’ai donc demandé à la FFSA la base de données utilisée ici, et pour éviter des problèmes d’inflation des coûts de sinistres entre 1992 et 2000. Si je prends tous les sinsitres, on obtient l’ajustement de Pareto suivant

soit une pente (en valeur absolue de 1.47). Mais encore une fois, l’ajustement de Pareto se fait sur les grands sinistres. Hill a proposé un estimateur très populaire pour estimer ce coefficient, où on ne prend en compte que les k observations les plus grandes, et on regarde l’estimation de la pente du graphique de Pareto pour ces quelques valeurs. On représente alors l’estimation en fonction du nombre de grands sinistres, ou du seuil définissant les graphs sinistres. Numériquement, en posant

$https://perso.univ-rennes1.fr/arthur.charpentier/latex/z08.png$

on peut écrire comme estimateur de la pente

$https://perso.univ-rennes1.fr/arthur.charpentier/latex/z09.png$

soit, en simplifiant le numérateur,

$https://perso.univ-rennes1.fr/arthur.charpentier/latex/z10.png$

tel que l’a construit Hill en 1975. Graphiquement, on a ici

Bref, la question est ce savoir si on atteint la valeur 1 pour les grands sinistres. Graphiquement, on a malgré tout envie de rejeter cette hypothèse.
Une solution peut être de faire un test statistique, basé sur de ratio de vraisemblance, comme le suggèrent Reiss & Thomas (2001) or Coles (2001). En fait, on peut même utiliser d’autres estimateurs que celui de Hill, comme celui obtenu en faisant un ajustement de loi GPD (Pareto généralisée
) sur la loi des Excès, ou une loi GEV sur des maximas par blocs (Generalized Extreme Value). On introduit alors la statistique de test suivante

$https://perso.univ-rennes1.fr/arthur.charpentier/latex/z12.png$

et on regarde les p-value (ainsi que la correction de Bartlett à droite),

On peut aussi, plus simplement, estimer plusieurs coefficients de pentes pour des seuils différents, et regarder la borne supérieure de l’intervalle de confiance,

Bref, même si avec un des ajustements de loi GPD on hésite à retenir une pente unitaire, la plupart des tests rejettent cette hypothèse, et donc le risque de perte d’exploitation semble être assurable, d’espérance mathématique finie. Bref, les dessins c’est très bien pour faire passer une idée, mais ne retenir que ça pour en tirer des conclusions aussi fortes me laisser sceptique….

Actuarial science, Extreme Value, Risk

Les généralités sont généralement fausses

16/09/2009 Arthur Charpentier Leave a comment

Suite à un commentaire sur un vieux billet (ici), concomitante avec une demande d’intervention pour la formation ERM proposée par l’Institut des Actuaires, je me suis replongé dans les calculs de SCR (solvency capital requierement). Plus je lis la documentation technique, plus je suis surpris par l’importance de l’hypothèse de normalité ! J’ai repris mon Loss Model, qui est le livre de base pour l’examen d’actuariat non-vie dans les pays anglo-saxons, et j’ai vu des dizaines de lois bizarres… mais rien sur la loi normale…

de l’hypothèse de normalité dans la formule d’agrégation

Comme je le disais dans un vieux billet (ici), on ne peut montrer la validité de la formule de calcul du SCR que dans le cas Gaussien (ou plus précisément dans le cas elliptique, voire ici). En effet, la preuve est assez simple. Je renvois vers un papier tout chaud de Laurent Devineau et Stéphane Loisel sur le sujet, ici.

où l’on continue à entrevoir l’hypothèse de normalité

En fait l’hypothèse de normalité est vraiment partout dans les calculs de SCR. Mais bien sûr, sans le dire. Par exemple dans un papier de la SCOR,

https://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/.doc14-1_m.jpg

Comme je le disais l’autre jour (ici), le problème des documents “professionels” (par opposition à “académiques“), c’est que comme on veut que le CEO le lise, on évite de parler d’hypothèses, ou pire, de mettre un peu de maths…. Bref, on apprend que “generally” (terme d’une grande rigueur scientifique, que je recaserai un jour dans mes démonstrations), “the expected shortfall[…] at the 99 % level corresponds quite closely to the[…] value-at-risk at a 99.6% level “. Mais on retrouve la même phrase dans un rapport de l’edhec

https://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/.SCR-edhec_m.jpg

Là aussi, aucune source… Damned ! En creusant un peu, Actuaris m’a mis sur la piste (ici), car on apprend que c’est l'”autorité suisse” qui a énoncé ce résultat que je qualifierais de surprenant¹.

https://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/.actuaris-var_m.jpg

Effectivement, on peut retrouver le document originel,

https://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/.doc13-1_m.jpg

J’avoue avoir essayé – sans grand succès – de lire l’annexe D. Je n’ai pas tout compris, mais j’y ai lu un “assuming that the distribution functions follows a normal law”,

https://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/.doc13-2_m.jpg

Donc si j’ai bien tout compris, “generally”, “the distribution functions follows a normal law” ! Bon, je crois que je vais pouvoir envoyer un mail à Stuart, Harry et Gordon pour leur expliquer qu’ils peuvent virer des pages à leur livre ! Ou peut-être Jean-Luc Besson (pour lequel j’ai une immense estime !) qui présente beaucoup de lois utiles en assurance dans l’ouvrage qu’il a coécrit avec Christian Partrat…
Mais peut-être a-t-il raison…? On peut se dire qu’en actuariat, on avait l’habitude de présenter l’approximation “normal power“, et peut être que ce résultat est vrai pour d’autres lois….
Regardons vite-fait sur quelques lois….. Par exemple, pour la loi lognormale, si on regarde en fonction de l’écart-type de la loi normale sous-jacente,

Effectivement, on colle au quantile à 99.6%. L’ordonnée est en échelle logarithmique pour que le dessin soit présentable, avec en pointillé la VaR (formule exacte, à 99.6% en rouge, et à 99.8% en bleu), et les ronds sont des estimations de la TVaR à 99%. Par contre, par exemple pour la loi lognormale, si on regarde en fonction de l’écart-type de la loi normale sous-jacente, plus la volatilité est grande, plus l’approximation est mauvaise.

Bon, essayons une loi Gamma (classique en assurance, car c’est une loi de la famille exponentielle, et souvent comparée à la loi lognormale).

Visiblement, là aussi l’approximation semble marcher…. mais rappelons que la loi Gamma est une loi à queue fine… Allons-y franchement, comparons avec le cas d’une loi de Pareto…

cette fois, l’abscisse se lit dans l’autre sens: à gauche, c’est le cas le plus risqué (avec une variance infinie entre 1 et 2). On notera que pour les risques les plus importants, là aussi on sous-estime la TVaR en considérant la VaR (à un seuil plus important). On retrouve d’ailleurs exactement le même type de graphique avec une loi de Student,

Bref, on me dira que je pinaille, n’empêche qu’à mon avis, les compangies sous-estiment “généralement” la TVaR (à supposer que ça soit la mesure de risque qui les intéresse) en prenant une hypothèse de normalité…

¹en tant qu’enseignant, ce genre de conclusion me gêne toujours, car je me fatigue à faire des cours sur la mesure de risque, j’essaye de présenter des notions aussi proprement que possible, en expliquant que les hypothèses c’est important, qu’il faut définir différentes notions de mesures de risques (entre les convexes, les cohérentes, les dynamiques, etc), et là tout d’un coup, on nous apprend que “generally” c’est kif kif….

Actuarial science, Copulas, Extreme Value, Insurance, Multivariate Analysis, Risk

Ecole d’été en actuariat à Brest

30/06/2009 Arthur Charpentier Leave a comment

Exposé sur mesures de risques et dépendance, à Brest, dans le cadre de l’école d’été, la première partie est en ligne

ainsi que la seconde

Computer, Copulas, Extreme Value, Seminar, Spatial

Séminaire Probabilité et Statistique, UBO, Brest

24/04/2009 Arthur Charpentier Leave a comment

Talk at the statistical seminar at the Université de Bretagne Occidentale, in Brest, ~~Wednesday May 6th~~ Tuesday May 5th, 14h (in 10 days), on “multivariate extremes“. Slides can be found here.

The talk will give a detailed introduction on multivariate extremes and related concepts. Then the case of Archimedean copula will be fully described (following the paper with Johan Segers).

[04/05/2009]: some applications in risk management will be shown at the end of talk, as well as some news things on spatial correlation.

and in order to illustrate tail convergence of Archimedean copulas, I have uploaded two animations, with tail independence below,

with tail dependence (or asymptotic dependence),

Extreme Value

The return period concept

05/04/2009 Arthur Charpentier Leave a comment

In a paper with David Sibaï (that can be found here), we had a discussion about the concept of “return period” in hydrology. Actually, it looks like this concept has been introduced by Emil Gumbel in his book on Statistics of Extremes.

Graphs are also proposed,

The link between a probability and time can be established clearly using the geometric distribution. The time of the first success has the following distribution

$\Pr(X = k) = (1 - p)^{k-1}\,p\,$

where $0< p \leq 1$ denotes the success probability. Then
$\mathrm{E}(X) = \frac{1}{p}, \qquad\mathrm{var}(X) = \frac{1-p}{p^2}.$

This means that the time we have to wait – on average – before the first success is then simply the inverse of the success probability.
Note that this distribution is simply the discrete version of the exponential distribution, satisfying the memoryless property (in the context of continuous variates). And actually, this geometric distribution is the only discrete distribution satisfying such a property.
For instance, to illustrate this idea, the 1910 flood in Paris was suppose to be a centenial event. It does not mean that we must have a similar event next year, it means that we have to wait (still) 100 year – on average – before having a similar event. Or similarly, such an event occurs every year with a 1% probability, assuming temporal independence (this was actually the point we discussed in our paper).

Freakonometrics