Freakonometrics

Suite à mon exposé d'hier, Pierre me demandait d'argumenter un point que j'avais évoqué oralement sans le justifier: " une régression médiane passe forcément par deux points du nuage ". Par exemple,

x=c(1,2,3)
y=c(3,7,8)
plot(x,y,pch=19,cex=1.5,xlim=c(0,4),ylim=c(0,10))
library(quantreg)
abline(rq(y~x,tau=.5),col="red")

Essayons de le justifier... de manière un peu heuristique. Commençons un peu au hasard, avec une droite qui passe "entre" les points (on admettra que la régression médiane sépare l'espace en deux, avec la moitié des points en dessous, la moitié au dessus, à un près pour des histoires de parité).

plot(x,y,pch=19,cex=1.5,xlim=c(0,4),ylim=c(0,10))
abline(a=-1,b=3.2,col="blue")

c'est joli, mais on peut sûrement faire mieux. En particulier, on va chercher ici à minimiser la somme des valeurs absolue des erreurs (c'est le principe de la régression médiane). Essayons de translater la courbe, vers le haut, ou vers le bas,

plot(x,y,pch=19,cex=1.5,xlim=c(0,4),ylim=c(0,10))
abline(a=-1,b=3.2,col="blue",lwd=2)
for(i in seq(-2,3,by=.5)){
abline(a=i,b=3.2,col="blue",lty=2)}

Si on regarde ce que vaut la somme des valeurs absolues des erreurs, on a

d=function(h) sum(abs(y-(h+3.2*x)))
H=seq(-4,6,by=.01)
D=Vectorize(d)(H)
plot(H,D,type="l")

> optimize(d,lower=-5,upper=5)
$minimum
[1] -0.2000091
 
$objective
[1] 2.200009

Bref, retenons cette courbe, et notons qu'elle passe par un des points,

Et c'est assez normal. On commençait avec deux points au dessus, et un en dessous. Soit la somme des valeurs absolues des erreurs . Si on translate de , on passe de à (tant que l'on a toujours deux points au dessus, car celui en dessous sera toujours en dessous). Si on translate de , (là encore tant qu'il n'y a qu'un point en dessous). Bref, on a intérêt à translater vers le haut. Si on dépasse le premier point rencontré, on se retrouve dans la situation inverse, avec un point au dessus et deux au dessous, et on a intérêt à redescendre. Bref, la translation optimale revient a s'arrêter dès qu'on croise un point. Autrement dit, la régression passe forcément par un point. Au moins.

Maintenant, essayons de faire pivoter la courbe autour de ce point, là encore afin de minimiser la somme des valeurs absolues des erreurs,

plot(x,y,pch=19,cex=1.5,xlim=c(0,4),ylim=c(0,10))
abline(a=optimize(d,lower=-5,upper=5)$minimum,b=3.2,
col="blue",lwd=2)
points(x[1],y[1],cex=1.8,lwd=2,col="blue")
for(i in seq(-1,5,by=.25)){
abline(a=(y[1]-i*x[1]),b=i,col="blue",lty=2)}

La distance est alors, en fonction de la pente de la droite

d2=function(h) sum(abs(y-((y[1]-h*x[1])+h*x)))
H=seq(-4,6,by=.01)
D=Vectorize(d2)(H)
plot(H,D,type="l")

Là encore on peut formaliser un peu

> optimize(d2,lower=-5,upper=5)
$minimum
[1] 2.500018
 
$objective
[1] 1.500018

Et si on regarde cette dernière droite, on passe par deux points,

plot(x,y,pch=19,cex=1.5,xlim=c(0,4),ylim=c(0,10))
h=optimize(d2,lower=-5,upper=5)$minimum
abline(a=(y[1]-h*x[1]),
b=h,col="purple",lwd=2)

Pour comprendre pourquoi, comparons les deux cas,

• en faisant un pivot vers le bas, pour un des points, la valeur absolue des erreurs augmente alors que pour l'autre, elle diminue

• en faisant un pivot vers le haut, là encore pour un des points, la valeur absolue des erreurs augmente alors que pour l'autre, elle diminue, mais en sens inverse par rapport au cas précédent,

Et pour faire simple, dans le premier cas, le gain (sur la baisse) compense la perte, alors que dans le second cas, c'est le contraire. Bref, on a intérêt à pivoter vers le bas, ou vers le point à droite. Jusqu'à l'atteindre. Et optimalement, on passera par ce point. Bref, on passera par deux points. Et je laisse les plus courageux regarder avec plus de trois points, mais c'est toujours pareil...

Exchangeability is an extremely concept, since (most of the time) analytical expressions can be derived. But it can also be used to observe some unexpected behaviors, that we will discuss later on with a more general setting. For instance, in a old post, I discussed connexions between correlation and risk measures (using simulations to illustrate, but in the context of exchangeable risk, calculations can be performed more accurately). Consider again the standard credit risk problem, where the quantity of interest is the number of defaults in a portfolio. Consider an homogeneous portfolio of exchangeable risk. The quantity of interest is here

or perhaps the quantile function of the sum (since the Value-at-Risk is the standard risk measure). We have seen yesterday that - given the latent factor - (either the company defaults, or not), so that

i.e. we can derive the (unconditional) distribution of the sum

so that the probability function of the sum is, assuming that

> proba=function(s,a,m,n){
+ b=a/m-a
+ choose(n,s)*integrate(function(t){t^s*(1-t)^(n-s)*
+ dbeta(t,a,b)},lower=0,upper=1,subdivisions=1000,
+ stop.on.error =  FALSE)$value
+ }
> QUANTILE=function(p=.99,a=2,m=.1,n=500){
+ V=rep(NA,n+1)
+ for(i in 0:n){
+ V[i+1]=proba(i,a,m,n)}
+ V=V/sum(V)
+ return(min(which(cumsum(V)>p))) }

> PICTURE=function(P){
+ A=seq(.01,2,by=.01)
+ VQ=matrix(NA,length(A),5)
+ for(i in 1:length(A)){
+ VQ[i,1]=QUANTILE(a=A[i],p=.9,m=P)
+ VQ[i,2]=QUANTILE(a=A[i],p=.95,m=P)
+ VQ[i,3]=QUANTILE(a=A[i],p=.975,m=P)
+ VQ[i,4]=QUANTILE(a=A[i],p=.99,m=P)
+ VQ[i,5]=QUANTILE(a=A[i],p=.995,m=P)
+ }
+ plot(A,VQ[,5],type="s",col="red",ylim=
+ c(0,max(VQ)),xlab="",ylab="")
+ lines(A,VQ[,4],type="s",col="blue")
+ lines(A,VQ[,3],type="s",col="black")
+ lines(A,VQ[,2],type="s",col="blue",lty=2)
+ lines(A,VQ[,1],type="s",col="red",lty=2)
+ lines(A,rep(500*P,length(A)),col="grey")
+ legend(3,max(VQ),c("quantile 99.5%","quantile 99%",
+ "quantile 97.5%","quantile 95%","quantile 90%","mean"),
+ col=c("red","blue","black",
+"blue","red","grey"),
+ lty=c(1,1,1,2,2,1),border=n)
+}

> PICTURE(.15)

> PICTURE(.05)

(with the modified code to visualize the quantile as a function of the underlying default correlation) or even worse,

> PICTURE(.0075)

An alternative (to profile likelihood techniques) to derive confidence intervals is to use the delta method. Consider an estimator such that

then for any differentiable transformation

The proof of that result is based on Taylor's expansion (see here or there for more details on the theory, or even on this blog - here, in French or there in English - for some codes in R). This can be used to derive an asymptotic confidence interval for a quantile. Consider the following dataset

> base1=read.table(
+ "http://freakonometrics.free.fr/danish-univariate.txt",
+ header=TRUE)
> library(evir)
> X=base1$Loss.in.DKM

It is possible to fit a Generalized Pareto distribution on observations above a given threshold,

In that case, if exceed the threshold out of a sample of size , the estimator of the quantile

i.e. . Then

> u=5
> GPD=gpd(X,u)
> theta=GPD$par.ests
> sigma=GPD$varcov
> k=GPD$n.exceed
> n=length(X)
> p=.975
> Q=u+theta[2]/theta[1]*((n*(1-p)/k)^(-theta[1])-1)
> nabla=c(-theta[2]/theta[1]^2*((1-p)^(-theta[1])-1)-
+ theta[2]/theta[1]*(1-p)^(-theta[1]*log(1-p)),
+ 1/theta[1]*((1-p)^(-theta[1])-1))
> variance=t(nabla)%*%sigma%*%nabla

Based on the assumption of normality, it is possible to derive confidence intervals, and to compare them with the one obtained in R,

> c(Q-1.96/sqrt(k)*sqrt(variance),
+   Q+1.96/sqrt(k)*sqrt(variance))
[1] 13.11562 16.82852
+ tailplot(gpd(X,5))
+ gpd.q(tailplot(gpd(X,5)), .975, ci.type =
+ "likelihood", ci.p = 0.95,like.num = 50)
Lower CI Estimate Upper CI
13.33329 14.97207 17.18094 

histoire de reprendre le titre d'un vieux billet mais dans un autre contexte, cf ici.
J'ai eu beaucoup de commentaires ou de réflexions suite à mes billets sur les accidents de la route, en France (ici, là, ou là encore). Certains me reprochaient de faire simplement de l'analyse statistique de chiffres, sans vraiment proposer de remèdes permettant de baisser le nombre de tués sur les routes. Ma première réponse, un peu simpliste consistait à dire que je ne suis que statisticien, et que faire baisser le nombre de tués sur les routes, c'est un problème politique (au sens étymologique du terme). Mais après réflexion, je me demande si c'est vraiment le cas... et si ce ne serait pas encore et toujours un problème de maths.
Par exemple, pour baisser le nombre de tués, une variable importante est le temps d'intervention des secours. Mais que vaut-il mieux ?

• que des ambulances (et des services d'urgence) attendent patiemment qu'on les appelle,
• ou que des ambulances tournent sur les routes, afin d'être au cœur de l'action...

• Cas où les accidents se produisent uniformément sur la route

• Cas où les accidents se produisent non-uniformément sur la route

• Comparaison de probabilités (et plus de valeurs moyennes)
  

soit numériquement

• Comparaison de probabilités dans le cas général
  

• si on suppose que les accident arrivent au bout de la route

a=1
f=function(x){3*x^2/(2*a^3)}
PROB1d=function(d){
1-integrate(f,-a,d)$value+integrate(f,-a,-d)$value}
PROB2d=function(d){
g1=function(x){(a-x-d)*f(x)}
g2=function(x){(x+a-d)*f(x)}
(integrate(g1,-a,a-d)$value+integrate(g2,-a+d,a)$value)/(2*a)}
P1=Vectorize(PROB1d)
P2=Vectorize(PROB2d)
dist=seq(0,1,by=.025)
plot(dist,P1(dist),lwd=3,col="red",type="l")
lines(dist,P2(dist),lwd=3,col="blue")

• si on suppose que les accident arrivent au centre la route

f=function(x){3*(a^2-x^2)/(4*a^3)}
 

A standard idea in extreme value theory (see e.g. here, in French unfortunately) is that to estimate the 99.5% quantile (say), we just need to estimate a quantile of level 95% for observations exceeding the 90% quantile.

In extreme value theory, we assume that the 90% quantile (of the initial distribution) can be obtained easily, e.g. the empirical quantile, and then, for the exceeding observations, we fit a Pareto distribution (a Generalized Pareto one to be precise), and get a parametric quantile for the 95% quantile. I.e.

library(mnormt)
library(quantreg)
library(splines)
set.seed(1)
mu=c(0,0)
r=0
Sigma <- matrix(c(1,r,r,1), 2, 2)
Z=rmnorm(2500,mu,Sigma)
X=Z[,1]
Y=Z[,2]
 
base=data.frame(X,Y)
plot(X,Y,col="blue",cex=.7)
I=(Y>qnorm(.25))&(Y<qnorm(.75))
baseI=base[I==FALSE,]
points(X[I],Y[I],col="light blue",cex=.7)
abline(h=qnorm(.25),lty=2,col="blue")
abline(h=qnorm(.75),lty=2,col="blue")
u=seq(-5,5,by=.02)
reg=rq(Y~X,data=base,tau=.05)
lines(u,predict(reg,newdata=data.frame(X=u)),lty=2)
reg=rq(Y~X,data=baseI,tau=.05*2)
lines(u,predict(reg,newdata=data.frame(X=u)))

Je poste un complément pour les étudiants qui suivent le cours de Statistique de l'Actuariat I, sur le calcul numérique d'un quantile. Comme je l'avais expliqué oralement en cours, il y a deux méthodes (au moins) pour étudier la distribution (et donc les quantiles) d'une loi composée

• des méthodes numériques (correspondant à la méthode de discrétisation de Panjer, et ses améliorations, qui sont détaillées sur ce blog, ici ou là)
• des méthodes de simulation (correspondant à ce que je demande d'implémenter pour le DM3).

converge uniformément vers la vraie fonction de répartition, c'est-à-dire que

1. correspond à l'inverse de la fonction de répartition empirique,
2. fait une moyenne des statistiques d'ordre pour corriger de la discontinuité,
3. est la définition de SAS, i.e. prendre la statistique d'ordre la plus proche,
4. est obtenue par interpolation linéaire entre la kème statistique d'ordre et p(k), où p(k) = k / n, i.e. linear interpolation of the empirical cdf.
5. est obtenue par interpolation linéaire entre la kème statistique d'ordre et p(k), où p(k) = (k - 0.5) / n, i.e.  piecewise linear function where the knots are the values midway through the steps of the empirical cdf. This is popular amongst hydrologists.
6. est obtenue par interpolation linéaire entre la kème statistique d'ordre et p(k) = k / (n + 1), i.e.  p(k) = E[F(x[k])]. Utilité par Minitab et SPSS.
7. est obtenue par interpolation linéaire entre la kème statistique d'ordre et p(k) = (k - 1) / (n - 1), i.e. p(k) = mode[F(x[k])]. Utilisé par S (et R)
8. est obtenue par interpolation linéaire entre la kème statistique d'ordre et p(k) = (k - 1/3) / (n + 1/3), i.e. p(k) =median[F(x[k])], the resulting quantile estimates are approximately median-unbiased regardless of the distribution.
9. est obtenue par interpolation linéaire entre la kème statistique d'ordre et p(k) = (k - 3/8) / (n + 1/4), i.e. the resulting quantile estimates are approximately unbiased for the expected order statistics if x is normally distributed.

Par défaut, R retient la méthode 7, comme le justifie le papier de Hyndman, R. J. et Fan, Y. (1996) Sample quantiles in statistical packages, American Statistician, 50, 361–365 (téléchargeable ici).

Il existe des dizaines d'autre méthodes pour estimer ces quantiles, en particulier dans library(quantile), ou encore hdquantile() implémentant l'estimaeur de Harrell et Davis. Pour aller plus loin, on avait fait un papier avec Abder Oulidi (ici) qui propose un comparaison d'un paquet d'estimateurs de quantiles).

Bref, pour résumer, une fois simulé le vecteur de coûts de sinistres S, il suffit d'utiliser quantile(S,.995).