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
We’ve seen Monday, in the MAT8595 course how to use the Generalized Pareto Distribution to estimate some downside risk measures, given a sample (assumed to be i.i.d., I will not mention here properties on extremes for stochastic processes) with distribution . The cumulative distribution function of the Pareto distribution is here
For some threshold , and
, we can write
From Pickands–Balkema–de Haan theorem, if is large enough, then
Given our sample , let
denote the number of observations over, threshold
. Then we can write
or equivalently
If we invert this function, we get the quantile of level ,
Actually, a threshold and then the implied number of observation exceeding that threshold, it is possible to consider a fixed number of observation, and then the associated threshold will be the associated order statistics.
The density of the Pareto distribution is here
which is here function of two paramters, and
.As discussed in the course, it is possible to use the Delta method to derive the asymptotic distribution of any quantile, and get then an approximated (asymptotic) confidence interval.
But since is usually not a parameter of interest, why not considering a reparametrization of our density, as a function of
and
(for some probability
that will be considered as fixed from now on). We can easily get (assuming that
) that
Tis expression is simple, and can be used to derive the likelihood (on the observations exceeding the threshold)
Numerically, let us write (and plot) that function. Consider some real data here
> X=as.numeric(danish) > Xs=sort(X,decreasing=TRUE) > n=length(X) > u=10 > nu=sum(X>u)
Consider, say, the 99.9% quantile,
> p=.999
The empirical quantile is here
> quantile(X,p) 99.9% 131.5519
The density and the loglikelihood functions are here
> gq=function(x,xi,q){ + ( (n/nu*(1-p) ) ^ (-xi)-1)/(xi*(q-u))* + (1+((n/nu*(1-p))^(-xi)-1)/(q-u)*x)^(-1/xi-1)} > loglik=function(param){ + xi=param[2];q=param[1] + lg=function(i) log(gq(Xs[i],xi,q)) + return(-sum(Vectorize(lg)(1:nu))) + }
We can try to plot this likelihood using
> h=201 > Q=seq(50,300,length=h) > XI=seq(.1,1,length=h) > XIQ=as.matrix(expand.grid(Q,XI)) > M=mapply(loglik,XIQ)
Unfortunately, it was not working, so I used the old style
> M=matrix(NA,h,h) > for(i in 1:h){for(j in 1:h){M[i,j]=loglik(c(Q[i],XI[j]))}}
The level curves of the log-likelihood are here
> hc=heat.colors(100) > image(Q,XI,-M,col=hc) > contour(Q,XI,-M,add=TRUE)
Again, since our interest is in the quantile, we can draw the profile likelihood and get the maximum of that function
> PL=function(Q){ + profilelikelihood=function(xi){ + loglik(c(Q,xi))} + return(optim(par=.8,fn=profilelikelihood)$value)} > (OPT=optimize(f=PL,interval=c(100,500))) $minimum [1] 111.1055
and the graph is
> XQ=seq(50,300,length=101) > L=Vectorize(PL)(XQ) > plot(XQ,-L,type="l") > up=OPT$objective > abline(h=-up) > abline(h=-up-qchisq(p=.95,df=1),col="red") > I=which(-L>=-up-qchisq(p=.95,df=1)) > lines(XQ[I],rep(-up-qchisq(p=.95,df=1),length(I)), + lwd=5,col="red") > abline(v=range(XQ[I]),lty=2,col="red")
which can be seen as an alternative to
> gpd.q(tailplot(gpd(X,u)),.999) Lower CI Estimate Upper CI 64.66184 94.28956 188.91752 $objective [1] 454.6481
If we want to focus on another downside risk measure, that shouldn’t be too difficult. For instance, the expected shortfall, can be estimated as
where denotes the mean excess function, which can be writen, with a Generalized Pareto Distribution
Thus, a natural estimator for the expected shortfall is
One more time, it is possible to re-parametrize the density of the Pareto distribution, using instead of
. Here, we get
The code to get the associated log-likelihood is here
> ge=function(x,xi,es){ + (xi+(n/nu*(1-p))^(-xi)-1)/(xi*(1-xi)*(es-u))*(1+(xi+(n/nu*(1-p))^(-xi) + -1)/((es-u)*(1-xi))*x)^(-1/xi-1) + } > loglik=function(param){ + xi=param[2];es=param[1] + lg=function(i) log(ge(Xs[i],xi,es)) + return(-sum(Vectorize(lg)(1:nu))) + }
and again, we can plot it
and the profile (log) likelihood is here (for the 99.9% expected shortfall)
> PL=function(ES){ + profilelikelihood=function(xi){ + loglik(c(ES,xi))} + return(optim(par=.8,fn=profilelikelihood)$value)} > (OPT=optimize(f=PL,interval=c(100,500))) $minimum [1] 143.66 $objective [1] 454.6481
which could be compared with
> gpd.sfall(tailplot(gpd(X,u)),.999) Lower CI Estimate Upper CI 96.64625 191.36972 394.87555
This week, we conclude the part on extremes with an application of extreme value theory to risk measures. We have seen last week that, if we assume that above a threshold 
, a Generalized Pareto Distribution will fit nicely, then we can use it to derive an estimator of the quantile function (for percentages such that the quantile is larger than the threshold)
It the threshold is 
, i.e. we keep the 
largest observations to fit a GPD, then this estimator can be written
The code we wrote last week was the following (here based on log-returns of the SP500 index, and we focus on large losses, i.e. large values of the opposite of log returns, plotted below)
> library(tseries)
> X=get.hist.quote("^GSPC")
> T=time(X)
> D=as.POSIXlt(T)
> Y=X$Close
> R=diff(log(Y))
> D=D[-1]
> X=-R
> plot(D,X)
> library(evir)
> GPD=gpd(X,quantile(X,.975))
> xi=GPD$par.ests[1]
> beta=GPD$par.ests[2]
> u=GPD$threshold
> QpGPD=function(p){
+ u+beta/xi*((100/2.5*(1-p))^(-xi)-1)
+ }
> QpGPD(1-1/250)
97.5%
0.04557386
> QpGPD(1-1/2500)
97.5%
0.08925095
This is similar with the following outputs, with the return period of a yearly event (one observation out of 250 trading days)
> gpd.q(tailplot(gpd(X,quantile(X,.975))), 1-1/250, ci.type = + "likelihood", ci.p = 0.95,like.num = 50) Lower CI Estimate Upper CI 0.04172534 0.04557655 0.05086785
or the decennial one
> gpd.q(tailplot(gpd(X,quantile(X,.975))), 1-1/2500, ci.type = + "likelihood", ci.p = 0.95,like.num = 50) Lower CI Estimate Upper CI 0.07165395 0.08925558 0.13636620
Note that it is also possible to derive an estimator for another population risk measure (the quantile is simply the so-called Value-at-Risk), the expected shortfall (or Tail Value-at-Risk), i.e.
The idea is to write that expression
so that we recognize the mean excess function (discussed earlier). Thus, assuming again that above (and therefore above that high quantile) a GPD will fit, we can write
or equivalently
If we substitute estimators to unknown quantities on that expression, we get
The code is here
> EpGPD=function(p){
+ u-beta/xi+beta/xi/(1-xi)*(100/2.5*(1-p))^(-xi)
+ }
> EpGPD(1-1/250)
97.5%
0.06426508
> EpGPD(1-1/2500)
97.5%
0.1215077
An alternative is to use Hill’s approach (used to derive Hill’s estimator). Assume here that 
, where 
is a slowly varying function. Then, for all 
,
Since is a slowly varying function, it seem natural to assume that this ratio is almost 1 (which is true asymptotically). Thus
i.e. if we invert that function, we derive an estimator for the quantile function
which can also be written
(which is close to the relation we derived using a GPD model). Here the code is
> k=trunc(length(X)*.025)
> Xs=rev(sort(as.numeric(X)))
> xiHill=mean(log(Xs[1:k]))-log(Xs[k+1])
> u=Xs[k+1]
> QpHill=function(p){
+ u+u*((100/2.5*(1-p))^(-xiHill)-1)
+ }
with the following Hill plot
For yearly and decennial events, we have here
> QpHill(1-1/250) [1] 0.04580548 > QpHill(1-1/2500) [1] 0.1010204
Those quantities seem consistent since they are quite close, but they are different compared with empirical quantiles,
> quantile(X,1-1/250) 99.6% 0.04743929 > quantile(X,1-1/2500) 99.96% 0.09054039
Note that it is also possible to use some functions to derive estimators of those quantities,
> riskmeasures(gpd(X,quantile(X,.975)),1-1/250) p quantile sfall [1,] 0.996 0.04557655 0.06426859 > riskmeasures(gpd(X,quantile(X,.975)),1-1/2500) p quantile sfall [1,] 0.9996 0.08925558 0.1215137
(in this application, we have assumed that log-returns were independent and identically distributed… which might be a rather strong assumption).
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…
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.
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,
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 surprenant1.
Effectivement, on peut retrouver le document originel,
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”,
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é…
1 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….
Pour reprendre le contexte général, Solvency II (l’analogue de la directive CRD pour les banques*) repose sur 3 piliers,
Cette histoire de pilliers peut s’illustrer de la manière suivante
Sur le premier pilier, assureurs et réassureurs devront mesurer les risques, et devront s’assurer qu’ils détiennent suffisamment de capital pour les couvrir. En pratique, le CEIOPS et la Commission Européenne ont retenu une probabilité de ruine de 0,5%. Les calculs de capital se font alors de deux manières, au choix,
En avril 2007, QIS3 a été lancé, afin de proposer une formule standard pour le calcul des MCR et SCR, en étudiant la problématique spécifique des groupes. En particulier, on trouve dans les documents la formule suivante (pour un calcul de basic SCR)
Cette formule sort du QIS3, mais on trouve des choses analogues dans Sandström (2004), par exemple,
Avec une contrainte forte sur la forme du SCR, il obtient alors
D’où sort cette formule ? Certains ont tenté des éléments de réponse, par exemple
Ce résultat n’est malheureusement pas très probant car il n’est jamais rien évoqué sur la dépendance entre les composantes, ce qui est troublant. Sandstôrm écrit quelque chose de similaire, même si pour lui “normalité” est ici entendu dans un cadre multivarié.
Une explication peut être trouvée dans un papier de Dietmar Pfeiffer et Doreen Straßburger (ici) paru dans le Scandinavian Actuarial Journal (téléchargeable ici). Il cherche à expliquer comment calculer le SCR,
Il note, et c’est effectivement l’intuition que l’on avait, que dans un monde Gaussien (multivarié), cette formule marche, aussi bien pour un SCR basé sur la VaR que la TVaR. En particulier, ils citent un livre de Sven Koryciorz, correspondant à sa thèse de doctorat, intitulée “Sicherheitskapitalbestimmung und –allokation in der Schadenversicherung. Eine risikotheoretische Analyse auf der Basis des Value-at-Risk und des Conditional Value-at-Risk“, publiée en 2004.
Sinon, pour aller un peu plus loin, on peut aussi noter, dans les rapports du CEIOPS des déclarations un peu troublantes, par exemple
Il est pourtant facile de montrer que ce n’est pas le cas (même si c’est effectivement ce que préconise la “formule standard“). Le graphique ci-dessous montre l’évolution de la VaR d’une somme de risques corrélés (échangeables) en fonction de la corrélation sous-jacente: sur cet exemple, les risques très très corrélés sont moins risqués que des risques moyennement corrélés.
(la loi sous-jacente est une copule de Student). En revanche pour la TVaR, sur le même exemple, la TVaR de la somme est effectivement une fonction croissante avec la corrélation,
(plus de compléments dans les slides de l’école d’été à Lyon l’été dernier, ici).
* Pour reprendre des éléments de la page de wikipedia (ici), la directive européenne CRD (Capital Requirements Directive, i.e. Fonds Propres Réglementaires) transpose dans le droit européen les recommandations des accords de Bâle II, visant à calculer les fonds propres exigés pour les établissements financiers (i.e. directives 2006/48/CEet 2006/49/CE) .
Pour le second cours à l’X, les slides sont en ligne, ainsi que les énoncés desexercices. Je ferais PC l’après midi, salle 75, bâtiment Paul Levy.
Chose promise, chose due… sur le second exercice, sur la dominance stochastique à l’ordre 2, il existe une propriété dite de single-crossing, comme je l’avais évoqué en PC (malheureusement ce n’est qu’une condition suffisante). Autrement dit, si les fonctions de répartition soient…
On parle parfois de lemme de Ohlin, ou du critère de Karlin et Novikoff (1963). Parmi les références, je peux vous mettre le résultat donné dans le livre de Michel Denuit, Jan Dhaene, Marc Goovaerts, et Rob Kaas, dans Actuarial Theory for Dependent Risks,
Sinon pour l’interprétation en terme de dominance stochastique à l’ordre 2, il faut remonter au livre de Marc Goovaerts, et Rob Kaas, Effective Actuarial Methods.
D’ailleurs, cette notion de comparaison est très liée aux notions de “stop-loss order”, “TVaR order”, “variability order” ou de “increasing convex”. Le Chapitre 3 du livre de Moshe Shaked et George Shantikumar porte d’ailleurs exclusivement sur cette notion.
En ce qui concerne le dernier point, sur l’interprétation dans le cas Gaussien, on peut utiliser ce résultat (même s’il ne donne qu’une condition suffisante).