The book Recent Econometric Techniques for Macroeconomic and Financial Data, edited by Gilles Dufrénot and Takashi Matsuki has been published. With Emmanuel Flachaire, we wrote a chapter on Pareto models for risk management
The book Recent Econometric Techniques for Macroeconomic and Financial Data, edited by Gilles Dufrénot and Takashi Matsuki has been published. With Emmanuel Flachaire, we wrote a chapter on 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.
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
library(devtools) install_github("freakonometrics/TopIncomes") library(TopIncomes)
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.
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.
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.
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')
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).
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
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
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/.
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,
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.
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 intégrable, alors pour tout
et
,
. La conséquence de cette propriété (qui est juste) est que
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 et
(définies sur le même espace),
, pour tout
et
. Et ça c’est fort. En particulier, on ne dit pas qu’il faut que
et
sont des variables indépendantes. Car cette propriété est valide quelle que soit la dépendance qui pourrait exister entre
et
.
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 , et que l’on simule les couples
, 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 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 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),
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…
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
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,
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).
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")
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“.
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).
i.e.
de telle sorte que
et
Aussi l’information de Fisher est .
La variance aymptotique de est alors
Nous avions noté que
i.e. , alors
. Aussi, comme le théorème central limite limite
où
la delta method permet d’écrire, en posant ,
i.e.
La variance asymptotique est alors
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. , on en déduit que
La variance asymptotique de notre estimateur est alors (en enlevant le facteur au 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 , 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 , à l’aide de nos trois estimateurs, i.e. pour un intervalle de confiance à 90%,
où les estimations des variances asymptotiques sont obtenues en substituant à
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 ,
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 . Mais c’est une autre histoire… à suivre donc….