Tag Archives: Pareto

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=
> 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


autrement dit


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


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

for(i in 1:length(R)){
norm.cop = normalCopula(r, dim = 2)
for(j in 1:20000){
X=cbind(qgpd(U[,1], .7 ,1),qgpd(U[,2], .5 ,1))
col="light yellow",border=NA)

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

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

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,

X=g(n);X1=g(n);X2=g(n);X3= g(n);X4=g(n)
Tp =cummax(X^p)/cumsum(X^p)
xlim=c(100,n),ylab="",col="light blue",xlab="")
lines(Tp1,col="light green")

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

for(i in 1:ns){
col="light green",ylab="",xlab="",main="")

In the case of exponentially distributed variables, we have


For variables with a lognormal distribution,


And finally, consider the case of a Pareto distribution


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

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




or even


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.




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


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




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


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

hist(U,probability=TRUE,,col="light green",

For heavy tailed distribution, or Pareto-type tails, consider Pareto samples, with distribution function https://f.hypotheses.org/wp-content/blogs.dir/253/files/2018/02/max-05.gif. Let https://f.hypotheses.org/wp-content/blogs.dir/253/files/2018/02/max-06.gif and https://f.hypotheses.org/wp-content/blogs.dir/253/files/2018/02/max-07.gif, then


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

hist(U,probability=TRUE,col="light green",

For light tailed distribution, or exponential tails, consider e.g. a sample of exponentially distribution variates, with common distribution function https://f.hypotheses.org/wp-content/blogs.dir/253/files/2018/02/max-01.gif. Let https://f.hypotheses.org/wp-content/blogs.dir/253/files/2018/02/max-02.gif and https://f.hypotheses.org/wp-content/blogs.dir/253/files/2018/02/max-03.gif, then


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

hist(U,probability=TRUE,col="light green",

Consider now a Gaussian https://f.hypotheses.org/wp-content/blogs.dir/253/files/2018/02/max-17.gif sample. We can use the following approximation of the cumulative distribution function (based on l’Hopital’s rule)


as https://f.hypotheses.org/wp-content/blogs.dir/253/files/2018/02/max-16.gif. Let https://f.hypotheses.org/wp-content/blogs.dir/253/files/2018/02/max-18.gif and https://f.hypotheses.org/wp-content/blogs.dir/253/files/2018/02/max-19.gif. Then we can get


as https://f.hypotheses.org/wp-content/blogs.dir/253/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,

hist(U,probability=TRUE,col="light green",

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://f.hypotheses.org/wp-content/blogs.dir/253/files/2018/02/max-22.gif, if




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


This gives us normalizing coefficients we should use here.

hist(U,probability=TRUE,col="light green",

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




de telle sorte que




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

  • propriétés de l’estimateur de la méthode des moments

Nous avions noté que


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



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




La variance asymptotique est alors

  • 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


i.e. pour la médiane,


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


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


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


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


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%,


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,

[1] 2.378498
[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….

Statistiques, de la théorie à la pratique

Afin d’illustrer un peu le cours, je mets en ligne ici. Il s’agit de coûts de sinistres ayant dépassé 1,000$ (qui ont été ordonnés de manière croissante). Pour importer la base en R, le code est simplement,

> X=read.table(
> X
 [1] 1.003 1.016 1.023 1.027 1.037 1.039 1.0
[13] 1.061 1.072 1.078 1.082 1.087 1.094 1.0
[25] 1.110 1.112 1.117 1.132 1.138 1.141 1.1
[37] 1.180 1.186 1.187 1.190 1.193 1.203 1.2
[49] 1.321 1.326 1.338 1.342 1.343 1.344 1.3
[61] 1.428 1.432 1.442 1.457 1.463 1.466 1.4
[73] 1.551 1.553 1.566 1.584 1.632 1.695 1.6
[85] 1.881 1.893 1.897 1.958 2.000 2.175 2.2
[97] 3.045 3.103 4.495 5.614

Pour information, l’histogramme ressemble à ça

On va supposer que l’on peut ajuster une loi de Pareto, dont la fonction de survie est donnée par


et de densité


L’espérance est


et la variance

  • Estimateur de la méthode des moments

Numériquement, on obtient

> mean(X)/(mean(X)-1)
[1] 2.946827
  • Estimateur du maximum de vraisemblance

On peut écrire la vraisemblance, et on recherche – numériquement – la maximum (ou le minimum de l’opposé de la log-vraisemblance car la fonction recherche ici les minimum),

> f=function(x,a){
+ a/x^(a+1)
+ }
> LogV = function(a,echantillon){
+ -sum(log(f(echantillon,a)))
+ }
> optim(fn=LogV,par=2,echantillon=X,method="BFGS")$par
[1] 2.845093

Mais on aurait pu faire mieux… car ici le maximum de vraisemblance peut s’écrire analytiquement

> length(X)/(sum(log(X)))
[1] 2.845093

(tout va bien, on trouve pareil).

  • Estimateur par régression linéaire

On peut noter que


ce qui veut dire qu’on peut tracer le nuage de points du log des couts, et s’ils sont alignés suivant une droite (passant par l’origine), la pente nous donnera un estimateur de notre paramètre. Le plus simple est de prendre la pente de la droite qui passe par l’origine, et par le centre de gravité du nuage,

> Z=X
> Y=log((length(Z):1)/length(Z))
> X=log(Z)
> plot(X,Y)
> -mean(Y)/mean(X)
[1] 2.753414
  • Estimation à l’aide de la médiane

Enfin, notons que la médiane est


Un estimateur naturel est alors


Numériquement, c’est alors tout simple,

> log(2)/log(quantile(X,.5))

Moralité, on a quatre estimateurs, qui nous donne quatre valeurs numériques différentes

2.9468 2.8451 2.7534 2.4137

Pour ceux qui veulent trouver d’autres estimateurs (car on peut en construire des centaines), une référence intéressante est la suivante,
Quandt, R. E. (1966). Old and new methods of estimation and the Pareto distribution, Metrika, 10, 55-82 (ici)
Le but de la statistique mathématique est de mieux connaître les propriétés (fondamentales) de ces estimateurs afin de choisir le “meilleur“. Et c’est un problème important, par exemple si on souhaite mettre en place un contrat de réassurance (et transférer les très grands sinistres). Par exemple si on veut transférer les sinistres dépassant 100,000$, la probabilité de dépasser ce seuil est 10 fois plus grande avec le dernier estimateur par rapport au premier…

> 1/100^2.4137
[1] 1.48799e-05
> 1/100^2.9468
[1] 1.277615e-06

Moralité, les résultats théoriques que l’on aborde en ce moment vont nous aider à résoudre des questions très pratiques…. à suivre.

Mais d’où sort la loi de Pareto ?

Un billet pour Jean François qui me demandait d’où sortait vraiment la loi de Pareto… Je ne reviendrais pas sur les applications de la loi de Pareto ici (je peux renvoyer à un ancien billet sur les applications en réassurance ici par exemple), si ce n’est pour noter que l’on apprend des choses surprenantes sur internet. Par exemple dans l’article sur la théorie du portefeuille de Wikipedia (ici) on apprend que dans le modèle moyenne-variance, les rendements sont supposés suivre une loi de Pareto…

Mais passons. La question de Jean François était sur l’origine de la loi de Pareto, et sur son lien avec la règle du 80-20. Pour les amateurs d’italien, on peut trouver les premières traces de la loi dans

(partiellement en ligne ici). Malheureusement, mes compétences en italien sont (très) limitées. Mais pour faire simple Vilfredo Pareto propose une loi en fonction puissance pour modéliser les revenus, en 1895,

Il fait également de l’estimation de son indice de queue (en calculant la pente dans le Pareto-plot comme on dirait maintenant)

Pour la version en français, on peut trouver la première apparition de la loi dans

(que l’on peut lire ici). Là encore, une illustration sur les revenus est proposée,

Le point de départ est une relation affine entre le logarithme de la fonction de répartition (ou plutôt de survie) et le logarithme du revenu,
et là encore une application est proposée (il me semble que la différence d’un point – environ – entre la version en italien et la version en français vient du fait qu’il travaille ici sur les fonctions de répartition, alors qu’auparavant, il travaillait davantage sur les densités),

Il propose alors une interprétation en terme de quantile (ce qui a été traduit par la loi du 80-20),

Alors il semble toutefois que Pareto n’est jamais explicitement noté que 20% des gens possédaient 80% de la richesse (qui a donné naissance au principe 80-20), ce qui pourrait sous-entendre que Pareto aurait devancé les travaux de Lorenz.Ce nom 80-20 semble avoir été attribué par Joseph Juran en 1954, dans un article “Pareto, Lorenz, Cournot, Bernoulli and others“. Il s’en est repenti (et excusé) par la suite, dans un mea culpa (ici).
Bref, tout cela a des implications évidentes sur les politiques de redistribution de la richesse, et donc sur la façon dont doit être mis en place un système d’impôts sur le revenus… Mais c’est une autre histoire…

Les modèles en réassurance

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 – 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

> sinpe = read.table("https://perso.univ-rennes1.fr/arthur.charpentier/sinpe.csv",header=TRUE,sep=";")
> head(sinpe)
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))