Tag Archives: Pareto

Pareto models for risk management

Our paper, with Emmanuel Flachaire, “Pareto models for risk management” is now online…

The Pareto model is very popular in risk management, since simple analytical formulas can be derived for financial downside risk measures (Value-at-Risk, Expected Shortfall) or reinsurance premiums and related quantities (Large Claim Index, Return Period). Nevertheless, in practice, distributions are (strictly) Pareto only in the tails, above (possible very) large threshold. Therefore, it could be interesting to take into account second order behavior to provide a better fit. In this article, we present how to go from a strict Pareto model to Pareto-type distributions. We discuss inference, and derive formulas for various measures and indices, and finally provide applications on insurance losses and financial risks.

Pareto Models for Top Incomes

This week, The Society for the Study of Economic Inequality (ECINEQ) organised the Eighth ECINEQ Meeting 2019  in Paris, hosted by the Paris School of Economics and the World Inequality Lab. Emmanuel Flachaire was there to present our joint work on Pareto Models for Top Incomes. Slides are also available


The paper is still available on hal, and the package (TopIncomes) is also available from github,

library(devtools)
install_github("freakonometrics/TopIncomes")
library(TopIncomes)

Pareto Models for Top Incomes

With Emmanuel Flachaire, we uploaded on hal a paper on Pareto Models for Top Incomes,

Top incomes are often related to Pareto distribution. To date, economists have mostly used Pareto Type I distribution to model the upper tail of income and wealth distribution. It is a parametric distribution, with an attractive property, that can be easily linked to economic theory. In this paper, we first show that modelling top incomes with Pareto Type I distribution can lead to severe over-estimation of inequality, even with millions of observations. Then, we show that the Generalized Pareto distribution and, even more, the Extended Pareto distribution, are much less sensitive to the choice of the threshold. Thus, they provide more reliable results. We discuss different types of bias that could be encountered in empirical studies and, we provide some guidance for practice. To illustrate, two applications are investigated, on the distribution of income in South Africa in 2012 and on the distribution of wealth in the United States in 2013.

This paper was presented at and UCSB and in several workshops this spring, and this Summer, Emmanuel will present it at ECINEQ.

Note that a R package is also available on github, TopIncomes.

Extended Scale Free Networks

With Emmanuel Flachaire, we recently uploaded a short article, Extended Scale-Free Networks, on arxiv.

Recently, Broido & Clauset (2019) mentioned that (strict) Scale-Free networks were rare, in real life. This might be related to the statement of Stumpf, Wiuf & May (2005), that sub-networks of scale-free networks are not scale-free. In the later, those sub-networks are asymptotically scale-free, but one should not forget about second-order deviation (possibly also third order actually). In this article, we introduce a concept of extended scale-free network, inspired by the extended Pareto distribution, that actually is maybe more realistic to describe real network than the strict scale free property. This property is consistent with Stumpf, Wiuf & May (2005): sub-network of scale-free larger networks are not strictly scale-free, but extended scale-free.

Autres Données pour le cours de Statistiques

Encore quelques données pour faire un TP de statistique. Ce sont des données de coûts d’ouragans, aux Etats-Unis,

> library(gdata)
> db=read.xls(
+     "http://sciencepolicy.colorado.edu/publications/special/public_data_may_2007.xls",
+     sheet=1)
> stupidcomma = function(x){
+  x=as.character(x)
+  for(i in 1:10){x=sub(",","",as.character(x))}
+  return(as.numeric(x))}
> base=db[,1:4]
> base$Base.Economic.Damage=
Vectorize(stupidcomma)(db$Base.Economic.Damage)
> base$Normalized.PL05=Vectorize(stupidcomma)(db$Normalized.PL05)
> base$Normalized.CL05=Vectorize(stupidcomma)(db$Normalized.CL05)

Lire la base est un peu pénible, à cause de virgules qui traînent comme séparateur de milliers. Mais peu importe, le petit code permet de nettoyer la base.

  • Fréquence d’ouragans

La première étape pourra être de travailler sur la fréquence. On pourra ajuster des lois de Poisson, ou binomiales, ou binomiales négatives, et faire des tests de stabilité du paramètre, avec le temps (on pourra  couper en 4 ou 5 blocs).

> TB <- table(base$Year)
> years <- as.numeric(names(TB))
> counts <- as.numeric(TB)
> years0=(1900:2005)[which(!(1900:2005)%in%years)]
> db <- data.frame(years=c(years,years0),
+       counts=c(counts,rep(0,length(years0))))
> db[88:93,]
   years counts
88  2003      3
89  2004      6
90  2005      6
91  1902      0
92  1905      0
93  1907      0
> plot(years,counts,type='h')

  • Coût des ouragans

Dans un second temps, on pourra travailler sur le coût des ouragans

> plot(base$Normalized.PL05/1e9,type="h")

On pourra tenter d’ajuster des lois, log normale, gamma, Pareto, etc. Faire des tests d’ajustement. Et tester la stabilité des paramètres, dans le temps (en coupant en 2 blocs, par exemple, ou davantage).

Calcul(s) d’information de Fisher

La semaine passée, on avait fait quelques calculs pour obtenir l’information de Fisher pour des lois classiques. Je voulais juste remettre au propre les calculs pour les lois à plusieurs paramètres. Pour la loi Gamma,

la log-vraisemblance s’écrit

de telle sorte que

Ici, même pas besoin de prendre une espérance car la Hessienne est constante

Continue reading Calcul(s) d’information de Fisher

“A 99% TVaR is generally a 99.6% VaR”

Almost 6 years ago, I posted a brief comment on a sentence I found surprising, by that time, discovered in a report claiming that

the expected shortfall […] at the 99 % level corresponds quite closely to the […] value-at-risk at a 99.6% level

which was inspired by a remark in Swiss Experience report,

expected shortfall […] on a 99% confidence level […} corresponds to approximately 99.6% to 99.8% Value at Risk

Continue reading “A 99% TVaR is generally a 99.6% VaR”

Modeling Earthquake Dynamics

In 2012, with Marilou Durand, student at UQAM, we have been working on the seismic gap hypothesis, see e.g. McCann et al. (1978) or Kagan & Jackson (1991), or to be more specific, on the dynamics between earthquakes magnitude (or seismic moment) and inter-occurence durations. Our paper should appear soon in the Journal of Seismology,

In this paper, we investigate questions arising in Parsons & Geist (2012). Pseudo causal models connecting magnitudes and waiting times are consider, through generalized regression. We do use conditional model (magnitude given previous waiting time, and conversely) as an extension to joint distribution model described in Nikoloulopoulos & Karlis (2008). On the one hand, we fit a Pareto distribution for earthquake magnitudes, where the tail index is a function of waiting time following previous earthquake; on the other hand, waiting times are modeled using a Gamma or a Weibull distribution, where parameters are function of the magnitude of the previous earthquake. We use those two models, alternatively, to generate the dynamics of earthquake occurrence, and to estimate the probability of occurrence of several earthquakes within a year, or a decade.

The paper is online on https://hal.archives-ouvertes.fr/.

Income distribution and Tour de France

A few days ago, Jean-François Mignot published an interesting article entitled Tour de France 2014 : pourquoi le vainqueur gagne 100 fois plus que le 10e. In this article, we have the following graph, with the income of the cyclist, as a function of his final ranking (the data where downloaded from http://sportbuzzbusiness.fr/)

> bike=read.csv(
+ "http://freakonometrics.free.fr/tourdefrance.csv",
+ sep=";",header=TRUE,dec=" ")

> bike[1:19,"prime"]=bike[1:19,"prime"]*1000
> plot(bike,log="y",type="b",
+ xlab="(Final) rank",ylab="Bonus")

As pointed out by Jean-François, if the winner gets a lot of money, the bonus decreases fast, very fast actually. Gini index is very high here

> library(ineq)
> ineq(X,type="Gini")
[1] 0.910461

and if we look at Lorenz curve, indeed, the Tour de France is not very equalitarian,

Continue reading Income distribution and Tour de France

Earthquake dynamics

I just upload on http://hal.archives-ouvertes.fr/hal-00871883 a joint paper entitled Modeling earthquake dynamics.

In this paper, we investigate questions arising in Parsons & Geist (2012). Pseudo causal models connecting magnitudes and waiting times are consider, through generalized regression. We do use conditional model (magnitude given previous waiting time, and conversely) as an extension to joint distribution model described in Nikoloulopoulos & Karlis (2008). On the one hand, we fit a Pareto distribution for earthquake magnitudes, where the tail index is a function of waiting time following previous earthquake; on the other hand, waiting times are modeled using a Gamma or a Weibull distribution, where parameters are function of the magnitude of the previous earthquake. We use those two models, alternatively, to generate the dynamics of earthquake occurrence, and to estimate the probability of occurrence of several earthquakes within a year, or a decade.

L’espérance est un opérateur linéaire, so what ?

Oui, “l’espérance est un opérateur linéaire”. On n’arrête pas d’insister sur cette propriété dans la plupart des cours de probabilité. En fait, je pense que cette propriété a deux implications importantes. La première est que l’on interprète souvent cette phrase en disant que si http://freakonometrics.blog.free.fr/public/perso5/espl01.gif intégrable, alors pour tout http://freakonometrics.blog.free.fr/public/perso5/espl4.gif et http://freakonometrics.blog.free.fr/public/perso5/espl6.gifhttp://freakonometrics.blog.free.fr/public/perso5/espl7.gif. La conséquence de cette propriété (qui est juste) est que http://freakonometrics.blog.free.fr/public/perso5/espl09.gif dès lors que la transformation n’est pas linéaire (ou au moins, il n’y a aucune raison pour qu’il y ait égalité).

Mais cette linéarité dit davantage. Elle dit que pour toutes variables aléatoires http://freakonometrics.blog.free.fr/public/perso5/espl01.gifet http://freakonometrics.blog.free.fr/public/perso5/espl2.gif (définies sur le même espace), http://freakonometrics.blog.free.fr/public/perso5/espl3.gif, pour tout http://freakonometrics.blog.free.fr/public/perso5/espl4.gif et http://freakonometrics.blog.free.fr/public/perso5/espl6.gif. Et ça c’est fort. En particulier, on ne dit pas qu’il faut que http://freakonometrics.blog.free.fr/public/perso5/espl01.gif et http://freakonometrics.blog.free.fr/public/perso5/espl2.gif sont des variables indépendantes. Car cette propriété est valide quelle que soit la dépendance qui pourrait exister entre http://freakonometrics.blog.free.fr/public/perso5/espl01.gif et http://freakonometrics.blog.free.fr/public/perso5/espl2.gif.

Pourtant @Arnaud (qui passe souvent à l’occasion sur le blog) m’a laissé un couple de questions par courriel: “la dépendance modifie-t-elle la moyenne ?[…] sur des exemples simples en R, la moyenne est toujours perturbée”. L’idée était que les espérances étaient approchées numériquement par simulation. Or les méthodes de Monte Carlo reposent sur la loi des grands nombres et sur le théorème central limite. En particulier, si http://freakonometrics.blog.free.fr/public/perso5/espl34.gif, et que l’on simule les couples http://freakonometrics.blog.free.fr/public/perso5/espl10.gif, alors

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

autrement dit

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

i.e. la vitesse de convergence est affectée par la dépendance, pas la valeur vers laquelle on va tendre (car la moyenne tend vers l’espérance, qui est linéaire). En particulier, la convergence sera d’autant plus lente que http://freakonometrics.blog.free.fr/public/perso5/espl33.gif sera grande. Donc effectivement, l’approximation pourrait être relativement mauvaise si on est sur des couples présentant une forte dépendance.

Par exemple, si on s’amuse à simuler 100 couples http://freakonometrics.blog.free.fr/public/perso5/espl10.gif Gaussiens, avec

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

on a les valeurs suivantes pour la moyenne de ce couple (obtenu en générant 10,000 échantillons de 100 paires), avec en abscisse la corrélation entre les deuxvariables

le trait rouge est la moyenne empirique des moyennes obtenues, et les régions sont délimitées par les quantiles à 5% (et 95%), 10% (et 90%) et 25% (et 75%). En moyenne, les simulations donnent la même chose, peu importe la dépendance. Mais plus la dépendance est forte, plus la variance sera grande, et plus la convergence sera lente (et les simulations auront d’autant plus tendance à s’éloigner – ou se disperser autour – de la vraie valeur). Mais cette réponse est partielle, car on a supposé que la variance était finie. Or cette hypothèse n’est pas nécessaire pour assurer la convergence de la moyenne vers l’espérance (et donc pour utiliser des méthodes de simulations). Que se passerait-il avec des lois de variance infinie ?

Pour autre chose que les vecteurs Gaussiens (par exemple des variables de variance infinie, comme des lois de Pareto) on peut utiliser le code suivant (on continue pour l’instant à utiliser une copule Gaussienne),

library(mnormt)
library(copula)
library(evir)
set.seed(1)
ns=100
R=seq(-.95,.95,by=.05)
M=BINF1=BSUP1=rep(NA,length(R))
for(i in 1:length(R)){
r=R[i]
norm.cop = normalCopula(r, dim = 2)
S=rep(NA,20000)
for(j in 1:20000){
U=rcopula(norm.cop,ns)
X=cbind(qgpd(U[,1], .7 ,1),qgpd(U[,2], .5 ,1))
S[j]=mean(X[,1]+X[,2])}
BINF1[i]=quantile(S,.05)
BSUP1[i]=quantile(S,.95)
M[i]=mean(S)
}
plot(R,M,col="white",ylim=range(c(BINF1,BSUP1)))
polygon(c(R,rev(R)),c(BINF1,rev(BSUP1)),
col="light yellow",border=NA)
lines(R,BINF1,lty=2)
lines(R,BSUP1,lty=2)
lines(R,M,lwd=2,col="red")
esperance=1+1/(1-.7)+1+1/(1-.5)
segments(min(R),esperance,max(R),esperance,col="blue")

Sur cette somme de deux lois de Pareto (comme dans le code ci-dessus), on obtient numériquement

Autrement dit, ici, l’effet dépendance semble avoir un impact plus faible sur l’approximation de la moyenne. On peut aussi essayer de changer la copule, par exemple mettre une copule de type Clayton (avec en abscisse cette fois le tau de Kendall),

(qui présente de la dépendance dans la queue inférieure) ou la copule duale (qui présente de la dépendance dans la queue supérieure)

Bref, la structure de la dépendance ne semble pas trop changer l’estimation de la moyenne (ou plutôt la vitesse de la convergence). On peut aussi se demander s’il y a un effet dimension. Par exemple, si au lieu de simuler des couples, on simule des vecteurs de plus grande taille. Par exemple, un vecteur échangeable en dimension 10, avec des marges normales centrées réduites,

norm.cop =normalCopula(r, dim = 10,dispstr="ex")

Moralité ? effectivement, la dépendance a un impact sur la vitesse de convergence quand on travaille sur des vecteurs aléatoires. Mais de manière assez surprenante, c’est surtout le cas pour la loi normale. Et étrangement, cet impact ne semble pas lié à de la dépendance extrême, ou à la dimension…

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

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 http://freakonometrics.blog.free.fr/public/perso5/subexp-30.gif, for some http://freakonometrics.blog.free.fr/public/perso5/subexp-31.gif, 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: lethttp://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

  • http://freakonometrics.blog.free.fr/public/perso5/subexp-20.gif if and only if

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

  • http://freakonometrics.blog.free.fr/public/perso5/subexp-21.gif for some http://freakonometrics.blog.free.fr/public/perso5/subexp-23.gif if and only if the exists a non-degenerate variable http://freakonometrics.blog.free.fr/public/perso5/Z.gif such that

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

  • http://freakonometrics.blog.free.fr/public/perso5/subexp-21.gif with http://freakonometrics.blog.free.fr/public/perso5/subexp-22.gif 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 http://freakonometrics.blog.free.fr/public/perso5/subexp00.gif, define http://freakonometrics.blog.free.fr/public/perso5/suuuuuubexp.gif and http://freakonometrics.blog.free.fr/public/perso5/subexp-5.gif. Then analogous results can be derived,

  • http://freakonometrics.blog.free.fr/public/perso5/subexp-99.gif if and only if

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

  • http://freakonometrics.blog.free.fr/public/perso5/subexp-21.gif for some http://freakonometrics.blog.free.fr/public/perso5/subexp-25.gif if and only if the exists a non-degenerate variable http://freakonometrics.blog.free.fr/public/perso5/Zk.gif such that

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

  • http://freakonometrics.blog.free.fr/public/perso5/subexp-21.gif with http://freakonometrics.blog.free.fr/public/perso5/subexp-22.gif 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)

or

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)

or

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).

Fisher-Tippett theorem and limiting distribution for the maximum

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 https://freakonometrics.hypotheses.org/files/2018/02/max-00.gif. For bounded distribution, consider e.g. the uniform distribution on the unit interval, i.e. https://freakonometrics.hypotheses.org/files/2018/02/max-09.gif on the unit interval. Let https://freakonometrics.hypotheses.org/files/2018/02/max-10.gif and https://freakonometrics.hypotheses.org/files/2018/02/max-11.gif. Then, for all https://freakonometrics.hypotheses.org/files/2018/02/max-12.gif and https://freakonometrics.hypotheses.org/files/2018/02/max-13.gif,

https://freakonometrics.hypotheses.org/files/2018/02/max-14.gif

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 https://freakonometrics.hypotheses.org/files/2018/02/max-05.gif. Let https://freakonometrics.hypotheses.org/files/2018/02/max-06.gif and https://freakonometrics.hypotheses.org/files/2018/02/max-07.gif, then

https://freakonometrics.hypotheses.org/files/2018/02/max-08.gif

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 https://freakonometrics.hypotheses.org/files/2018/02/max-01.gif. Let https://freakonometrics.hypotheses.org/files/2018/02/max-02.gif and https://freakonometrics.hypotheses.org/files/2018/02/max-03.gif, then

https://freakonometrics.hypotheses.org/files/2018/02/max-04.gif

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 https://freakonometrics.hypotheses.org/files/2018/02/max-17.gif sample. We can use the following approximation of the cumulative distribution function (based on l’Hopital’s rule)

https://freakonometrics.hypotheses.org/files/2018/02/max-15.gif

as https://freakonometrics.hypotheses.org/files/2018/02/max-16.gif. Let https://freakonometrics.hypotheses.org/files/2018/02/max-18.gif and https://freakonometrics.hypotheses.org/files/2018/02/max-19.gif. Then we can get

https://freakonometrics.hypotheses.org/files/2018/02/max-20.gif

as https://freakonometrics.hypotheses.org/files/2018/02/max-21.gif. 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 https://freakonometrics.hypotheses.org/files/2018/02/max-22.gif, if

https://freakonometrics.hypotheses.org/files/2018/02/max-23.gif

then

https://freakonometrics.hypotheses.org/files/2018/02/max-24.gif

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

https://freakonometrics.hypotheses.org/files/2018/02/max-25.gif

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")

Statistiques, de la théorie à la pratique (partie 2)

Il y a une quinzaine de jours (ici) on avait vu comment construire des estimateurs pour l’indice de queue de la loi de Pareto. On avait 4 valeurs (différentes). Maintenant, on peut se demander lequel devrait être le “meilleur“.

  • propriétés de l’estimateur du maximum de vraisemblance

Maintenant que l’on a un peu avancé dans le cours, on peut continuer l’étude de la distribution des couts de sinistres (étude commencée ici).

http://freakonometrics.blog.free.fr/public/perso2/pareto-100.gif

i.e.

http://freakonometrics.blog.free.fr/public/perso2/pareto-101.gif

de telle sorte que

http://freakonometrics.blog.free.fr/public/perso2/pareto-102.gif

et

http://freakonometrics.blog.free.fr/public/perso2/pareto-204.gif

Aussi l’information de Fisher est http://freakonometrics.blog.free.fr/public/perso2/pareto-203.gif.
La variance aymptotique de http://freakonometrics.blog.free.fr/public/perso2/pareto-205.gif est alors

http://freakonometrics.blog.free.fr/public/perso2/pareto-206.gif
  • propriétés de l’estimateur de la méthode des moments

Nous avions noté que

http://freakonometrics.blog.free.fr/public/perso2/pareto-207.gif

i.e. http://freakonometrics.blog.free.fr/public/perso2/pareto-208.gif, alors http://freakonometrics.blog.free.fr/public/perso2/pareto-209.gif. Aussi, comme le théorème central limite limite

http://freakonometrics.blog.free.fr/public/perso2/pareto-210.gif

http://freakonometrics.blog.free.fr/public/perso2/pareto-211.gif

la delta method permet d’écrire, en posant http://freakonometrics.blog.free.fr/public/perso2/pareto-311.gif,

http://freakonometrics.blog.free.fr/public/perso2/correc-latex-pareto.gif

i.e.

http://freakonometrics.blog.free.fr/public/perso2/correc-latex-pareto2.gif

La variance asymptotique est alors

http://freakonometrics.blog.free.fr/public/perso2/pareto-310.gif
  • propriétés de l’estimateur de la méthode de la médiane

Pour le dernier estimateur, le théorème de Glivenko Cantelli sur les fonctions quantiles permet d’écrire

http://freakonometrics.blog.free.fr/public/perso2/pareto-309.gif

i.e. pour la médiane,

http://freakonometrics.blog.free.fr/public/perso2/pareto-308.gif

où la variance asymptotique (pour la médiane) peut se simplifier sous la forme

http://freakonometrics.blog.free.fr/public/perso2/pareto-307.gif

Là encore, en utilisant la delta method (je peux renvoyer ici pour une application numérique de cette méthode), avec

http://freakonometrics.blog.free.fr/public/perso2/pareto-306.gif

i.e. http://freakonometrics.blog.free.fr/public/perso2/pareto-305.gif, on en déduit que

http://freakonometrics.blog.free.fr/public/perso2/pareto-304.gif

La variance asymptotique de notre estimateur est alors (en enlevant le facteur http://freakonometrics.blog.free.fr/public/perso2/IC-std-6.gifau dénominateur – car il est commun à tous),

http://freakonometrics.blog.free.fr/public/perso2/pareto-303.gif

Pour trois de nos estimateurs (je n’ai pas fait le quatrième, ça serait hors programme), on a pu obtenir facilement des variances asymptotiques. Graphiquement, on a les évolutions suivantes pour les variances asymptotiques de nos trois estimateurs, en fonction de la valeur de http://freakonometrics.free.fr/blog/latex/information09.gif, avec le maximum de vraisemblance (en rouge), celui de la méthode des moments (en bleu) et celui basé sur la médiane (en mauve)

Manifestement le premier estimateur (le maximum de vraisemblance) a une variance asymptotique plus faible que les deux autres (mais on l’avait vu en cours, c’est ce que garantit l’inégalité de Cramér-Rao). Autrement dit, l’estimateur du maximum de vraisemblance est probablement le plus fiable….
En utilisant ces résultats mentionnés, on peut aussi construire des intervalles de confiance asymptotiques pour la valeur du paramètre http://freakonometrics.free.fr/blog/latex/information09.gif, à l’aide de nos trois estimateurs, i.e. pour un intervalle de confiance à 90%,

http://freakonometrics.blog.free.fr/public/perso2/pareto-302.gif

où les estimations des variances asymptotiques sont obtenues en substituant http://freakonometrics.blog.free.fr/public/perso2/pareto-205.gifà  http://freakonometrics.free.fr/blog/latex/information09.gif dans les formules au dessus.
Si on considère que l’estimateur du maximum de vraisemblance est le meilleur que l’on puisse utiliser, on obtient comme intervalle de confiance pour http://freakonometrics.free.fr/blog/latex/information09.gif,

2.845093-1.64/sqrt(100)*2.845093
[1] 2.378498
 
2.845093+1.64/sqrt(100)*2.845093
[1] 3.311688

L’étape d’après sera d’introduire des tests. On pourra ainsi se demander si l’on peut accepter l’hypothèse http://freakonometrics.blog.free.fr/public/perso2/pareto-301.gif. Mais c’est une autre histoire… à suivre donc….