Digital resources in the Social Sciences and Humanities OpenEdition Our platforms OpenEdition Books OpenEdition Journals Hypotheses Calenda Libraries OpenEdition Freemium Follow us

Probabilité et contre-intuition

Je suis tombé un peu par hasard sur un vieux papier de Peter Enis, Enis (1973), qui revient sur une relation que l’on prend toujours pour acquis (moi le premier) à savoir http://freakonometrics.blog.free.fr/public/perso5/proba07.gif. Mais il faut des hypothèses complémentaires pour que cette relation soit vraie. Le fait que http://freakonometrics.blog.free.fr/public/perso5/proba09.gif existe par exemple ne suffit pas ! (je laisse les amateurs de résultats contre-intuitifs le temps de réfléchir à cette affirmation)

Supposons par exemple que http://freakonometrics.blog.free.fr/public/perso5/proba02.gif, conditionnellement à http://freakonometrics.blog.free.fr/public/perso5/proba03.gif, suive une loi normale centrée, d’écart-type http://freakonometrics.blog.free.fr/public/perso5/proba08.gif,

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

où http://freakonometrics.blog.free.fr/public/perso5/proba03.gif est une variable (strictement) positive. Presque sûrement http://freakonometrics.blog.free.fr/public/perso5/proba05.gifquelle que soit la loi de http://freakonometrics.blog.free.fr/public/perso5/proba03.gif. Et donc http://freakonometrics.blog.free.fr/public/perso5/proba04.gif. Mais si la loi de http://freakonometrics.blog.free.fr/public/perso5/proba03.gif présente des queues “trop épaisses” (sans pour autant être extrême pour autant d’ailleurs), il est fort possible que http://freakonometrics.blog.free.fr/public/perso5/proba02.gif ne soit plus intégrable ! Par exemple (c’est ce que propose Enis (1973)), on peut supposer que http://freakonometrics.blog.free.fr/public/perso5/proba03.gif suit une loi Gamma de paramètres 1/2 et 2, i.e. de densité

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

Dans ce cas, http://freakonometrics.blog.free.fr/public/perso5/proba02.gif suit une loi de Cauchy, qui n’est pas d’espérance finie. C’est pénible, hein ?

Bon, d’un autre côté, c’est très lié à la définition de l’espérance conditionnelle retenue. En l’occurrence, j’utilisais l’espérance de la loi conditionnelle. Mais je ferais peut être un billet un jour sur les diverses constructions de l’espérance conditionnelle.

Simulation d’un processus de Lévy, et discrétisation

Avec @renaudjf, on discutait l’autre jour de la simulation d’un processus de Lévy. Et on se posait la question d’un algorithme optimal pour combiner un processus de Poisson (ou un process Poisson composé) avec un processus de Wiener (avec éventuellement un drift, voire une diffusion plus générale). En fait, pour générer des processus de Poisson, j’ai toujours eu l’habitude de simuler les durées entre sauts (avec des lois exponentielles, indépendantes, comme dans des vieux billets). Jean François me suggérait d’utiliser une propriété d’uniformité des sauts sur un intervalle de temps donné, conditionnellement aux nombres de sauts.

Commençons par la première piste. On peut générer un processus de Wiener, éventuellement avec un drift, et à coté, on peut générer les lois exponentielles  (qui vont correspondre aux durées entre sauts), et éventuellement aussi des amplitudes de sauts (e.g. des pertes qui suivent des lois exponentielles). On a ici

où . On commence par générer  en notant que

où les incréments  sont Gaussiens (centrés et de variance ) et indépendants les uns des autres. Quant aux durées entre sauts, les  , ce sont des lois exponentielles indépendantes, de moyenne . Voilà le code qui permet de générer les trois composantes,

n=1000
h=1/n
lambda=5
set.seed(2)
W=c(0,cumsum(rnorm(n,sd=sqrt(h))))
W=rexp(100,lambda)
N=sum(cumsum(W)<1)
T=cumsum(W[1:N])
X=-rexp(N)

Le hic est que pour le processus de Wiener, on a du discrétiser, alors que pour le processus de Poisson composé, non. Pourtant, il va bien se ramener sur une même échelle de temps. Une première piste est de créer vraiment la fonction 

Lt=function(t){
W[trunc(n*t)+1]+sum(X[T<=t])+lambda*t
}

et pour faire un dessin ensuite, c’est un jeu d’enfant

L=Vectorize(Lt)
u=seq(0,1,length=n+1)
plot(u,L(u),type="l",col="blue")

Une autre possibilité est d’utiliser une propriété d’uniformité du processus de Poisson que j’évoquais en introduction. Car le processus de Poisson vérifie une propriété remarquable: si  est la date où survient le ième saut, , alors conditionnellement au fait que , les variables  correspondent aux statistiques d’ordres de  variables indépendantes, uniformément distribuées sur , i.e.

Cette propriété se trouve dans Wolff (1982). L’idée de la démonstration est relativement simple. Commençons avec un (unique) saut, alors pour ,

i.e. on retrouve la fonction de répartition d’une loi uniforme sur . On itère ensuite avec 2 sauts, 3 sauts, etc.

La traduction en R de cette idée est tout simplement (car on se place sur )

N=rpois(1,lambda)
T=runif(N)
X=-rexp(N)

Ensuite, une stratégie est de discrétiser le processus de Poisson, avec le même pas de temps que le processus de Wiener,

indice=trunc(T*n)+1
saut=rep(0,n+1)
saut[indice]=X
processus=W+cumsum(saut)+lambda*u

On retrouve la même trajectoire qu’auparavant

plot(u,processus,type="l",col="red")

,,

Sauf qu’on a eu de la chance. Avec cette procédure, il ne faut pas que l’on ait deux sauts dans le même intervalle de temps ! Bon, il est vrai qu’une caractérisation du processus de Poisson est que

et donc on doit avoir très peu de chance d’avoir deux sauts au même instant d’autant plus que le pas de temps est petit. Mais “peu de chance” ne veut pas dire nul, et si on génère des milliers de trajectoires, la probabilité d’avoir une fois un soucis n’est pas négligeable.

Jean-François a eu l’idée brillante de proposer de tirer non pas des lois uniformes sur , mais des lois uniformes discrètes, dans

T=c(0,sort(sample((1:(n-1)/n),size=N,replace=FALSE)))

sans remise afin d’éviter d’avoir deux sauts au même moment. L’idée était séduisante, mais il me restait à etre convaincu que les durées entre sauts… continuaient à être (presque) exponentielles.

Pour ça, on peut faire quelques tests. Par exemple, générer quelques simulations pour avoir une centaines de sauts (et donc une centaine de durées entre sauts), puis faire un test de loi exponentielle (de moyenne )

VT=0
for(ns in 1:20){
N=rpois(1,lambda)
if(N>0){
T=c(0,sort(sample((1:(n-1)/n),size=N,replace=FALSE)))
VT=c(VT,diff(T))
}}

On fait ici 20 boucles car on avait fixé

lambda=5

et j’avais dit que je voulais une centaine d’observations pour faire un test de loi (ce qui est purement arbitraire, on en conviendra). On peut ensuite faire un test d’ajustement de loi exponentielle,

ks.test(VT[-1],"pexp",lambda)$p.value

Si on répète un grand nombre de fois, en changeant le pas de temps (ou le nombre de subdivisons de l’intervalle de temps), on notera qu’effectivement, avec un grand pas de temps (à gauche ci-dessous) on rejette souvent, voire presque toujours, l’hypothèse de loi exponentielle. Mais qu’assez vite, c’est une hypothèse qui n’est pas invraisemblable,

Je ne sais pas entre ces deux fonctions ce qui serait le plus rapide, mais on a deux jolis algorithmes pour générer les processus de Lévy. Non ?

Sondages et prévisions de séries temporelles

Ce soir, @imparibus proposait sur son blog un billet passionnant sur l’élection présidentielle en France, avec des graphiques superbes, basé sur un lissage temporel des résultats des sondages pour le premier tour de l’élection présidentielle en France (qu’il a fait sous excel, on saluera la performance !). En reprenant les données du site http://www.lemonde.fr/, on peut obtenir la même chose assez simplement (j’ai repris – ou presque – les couleurs utilisées sur le site, François Hollande en rose, Nicolas Sarkozy en bleu, François Bayrou en orange, et Marine Le Pen en noir…). Pour commencer, le code en R pour importer les données et les manipuler est le suivant

sondage=read.table("http://freakonometrics.blog.free.fr/
public/data/sondage1ertour2012.csv",
header=TRUE,sep=";",dec=",")
sondage[36,1]="11/03/12"
sondeur=unique(sondage$SONDEUR)
sondage$date=as.Date(as.character(sondage$DATE),
"%d/%m/%y")

(je fais manuellement une correction de date car je me suis trompé en saisissant les chiffres à la main). A partir de là, on peut reprendre l’idée de faire une régression locale pour trouver une tendance,

CL=brewer.pal(6, "RdBu")
couleur=c(CL[1],CL[6],"orange","grey")
datefin=as.Date("2012/04/12")
k=1
plot(sondage$date,sondage[,k+2],col=couleur[k],
pch=19,cex=.7,xlab="",ylab="",ylim=c(0,40),
xlim=c(min(sondage$date),datefin))
rl=lowess(sondage$date,sondage[,k+2])
lines(rl,lwd=2,col=couleur[k])
for(k in 2:4){
points(sondage$date,sondage[,k+2],
col=couleur[k],pch=19,cex=.7)
rl=lowess(sondage$date,sondage[,k+2])
lines(rl,lwd=2,col=couleur[k])
}

Le code a l’air long mais il faut définir les couleurs, générer une fenêtre graphique, mettre les régressions locales dedans, etc. En tant que telle, la commande qui permet de faire le lissage est juste

rl=lowess(sondage$date,sondage[,k+2])

Bon, on peut d’ailleurs faire toutes sortes de lissages, avec des splines par exemple (ce que j’ai davantage tendance à utiliser),

D=data.frame(date=min(sondage$date)+0:148)
library(splines)
k=1
plot(sondage$date,sondage[,k+2],
col=couleur[k],pch=19,cex=.7,
xlab="",ylab="",ylim=c(0,40),xlim=
c(min(sondage$date),datefin))
rs=lm(sondage[,k+2]~bs(date),data=sondage)
prl=predict(rs,newdata=D)
prlse=sqrt(prl/100*(1-prl/100)/1000)*100
polygon(c(D$date,rev(D$date)),c(prl+2*prlse,
rev(prl-2*prlse)),col=CL[3],border=NA)
lines(D$date,prl,lwd=2,col=CL[1])
points(sondage$date,sondage[,k+2],
col=couleur[k],pch=19,cex=.7,)
for(k in 2:4){
points(sondage$date,sondage[,k+2],
col=couleur[k],pch=19,cex=.7)
}

avec, là encore, la commande suivante pour lisser

rs=lm(sondage[,k+2]~bs(date),data=sondage)

Cette fois le code est un peu plus long, parce que j’ai tracé un intervalle de confiance à 95% (intervalle de confiance classique, en supposant que 1000 personnes ont été interrogées, comme évoqué dans d’anciens billets),

http://freakonometrics.hypotheses.org/files/2015/12/adic10.gif

On a ici la courbe suivante

avec l’intervalle de confiance pour François Hollande, mais on peut faire la même chose pour Nicolas Sarkozy,

k=2
plot(sondage$date,sondage[,k+2],
col=couleur[k],pch=19,cex=.7,
xlab="",ylab="",ylim=c(0,40),xlim=
c(min(sondage$date),datefin))
rs=lm(sondage[,k+2]~bs(date),data=sondage)
prl=predict(rs,newdata=D)
prlse=sqrt(prl/100*(1-prl/100)/1000)*100
polygon(c(D$date,rev(D$date)),c(prl+2*prlse,
rev(prl-2*prlse)),col=CL[4],border=NA)
lines(D$date,prl,lwd=2,col=CL[6])
points(sondage$date,sondage[,k+2],col=
couleur[k],pch=19,cex=.7,)
for(k in c(1,3:4)){
points(sondage$date,sondage[,k+2],col=
couleur[k],pch=19,cex=.7)
}

Mais comme auparavant, on peut aussi visualiser les “tendances” des quatre candidats en tête,

Voilà pour le début. En fait, idéalement, on voudrait faire un peu de prévision… Le hic est que les code usuel pour faire de la prévision (avec des ARIMA, du lissage exponentiel ou tout autre modèle classique) nécessitent de travailler avec des séries temporelles, telles qu’elles sont classiquement définies, c’est à dire “régulièrement espacées dans le temps“.

Je ne vais pas commencer à réclamer plus de sondages, loin de moins cette idée, mais pour mon modèle, j’avoue que j’aurais préféré avoir plus de points. Beaucoup plus de points. Un sondage par jour aurait été idéal en fait…

Qu’à cela ne tienne, on va simuler des sondages. L’idée est simple: on a une tendance qui nous donne une probabilité jour par jour (c’est exactement ce qui a été calculé pour tracer les courbes), via

D=data.frame(date=min(sondage$date)+0:148)
prl=predict(rs,newdata=D)

On va ensuite utiliser une hypothèse de normalité multivariée sur nos 4 candidats, auquel j’en ai en rajouté un cinquième fictif (mais ça ne servait à rien) en utilisant l’analogue multivariée de l’expression suivante (évoqué dans le premier billet de l’année sur les élections)

http://freakonometrics.hypotheses.org/files/2015/12/ic-sondage01.gif

On peut alors générer, jour après jour, des sondages, en simulant des vecteurs Gaussiens. Je vais les générer indépendants les uns des autres, parce que c’est plus simple, et que cela ne me semble pas être une trop grosse hypothèse. Pour générer un ensemble de sondages, on utilise la fonction suivante

library(mnormt)
simulation=function(S){
proba=c(S,100-sum(S))/100
variance= -proba%*%t(proba)
diag(variance)=proba*(1-proba)
variance=variance/1000
return(rmnorm(1,proba[1:4],variance[1:4,1:4]))}
simulsondages=function(M){
Z=rep(NA,ncol(M))
for(i in 1:nrow(M)){
Z=rbind(Z,simulation(M[i,]))
}
return(Z[-1,])}
prediction4=matrix(NA,nrow(D),4)
for(k in 1:4){
rs=lm(sondage[,k+2]~bs(date),data=sondage)
prediction4[,k]=predict(rs,newdata=D)
}

Par exemple, si on la fait tourner une fois, on obtient le graphique suivant

set.seed(1)
S4=simulsondages(prediction4)
S100=100*S4
k=1
plot(D$date,S100[,k],col=couleur[k],
pch=19,cex=.7,
xlab="",ylab="",ylim=c(0,40),xlim=
c(min(sondage$date),datefin))
rs=lm(sondage[,k+2]~bs(date),data=sondage)
lines(D$date,predict(rs,newdata=D),
lwd=2,col=couleur[k])
for(k in 2:4){
points(D$date,S100[,k],col=couleur[k],
pch=19,cex=.7)
rs=lm(sondage[,k+2]~bs(date),data=sondage)
lines(D$date,predict(rs,newdata=D),
lwd=2,col=couleur[k])
}

(les courbes lissées sont celles obtenues sur les vrais sondages). Là on peut être heureux, parce qu’on a des vraies séries temporelles. On peut alors faire de la prévision, par exemple en faisant du lissage exponentiel (optimisé),

library(forecast)
k=1
X=S100[,k]
ETS=ets(X)
F=forecast(ETS,h=60)
fdate=max(D$date)+1:60
k=1
plot(D$date,S100[,k],col=couleur[k],
pch=19,cex=.7,
xlab="",ylab="",ylim=c(0,40),xlim=
c(min(sondage$date),datefin))
rs=lm(sondage[,k+2]~bs(date),data=sondage)
lines(D$date,predict(rs,newdata=D),
lwd=2,col=couleur[k])
for(k in 2:4){
points(D$date,S100[,k],col=couleur[k],
pch=19,cex=.7)
rs=lm(sondage[,k+2]~bs(date),data=sondage)
lines(D$date,predict(rs,newdata=D),lwd=2,col=couleur[k])
}
polygon(c(fdate,rev(fdate)),c(as.numeric(F$lower[,2]),
rev(as.numeric(F$upper[,2]))),col=CL[3],border=NA)
polygon(c(fdate,rev(fdate)),
c(as.numeric(F$lower[,1]),rev(as.numeric(F$upper[,1]))),
col=CL[2],border=NA)
lines(fdate,as.numeric(F$mean),lwd=2,col=CL[1])

Là encore le code peut paraître long, mais c’est surtout la partie associée au graphique qui prend de la place (par soucis purement esthétique, on passe du temps sur les codes graphiques depuis le début). Le code qui modélise la série, et qui la projette, est ici donné par les deux lignes suivantes

ETS=ets(X)
F=forecast(ETS,h=60)

Pour François Hollande, sur la simulation des sondages passés que l’on vient d’effectuer, on obtient la prévision suivante

On peut bien entendu faire une projection similaire pour Nicolas Sarkozy,

k=2
X=S100[,k]
ETS = ets(X)
F=forecast(ETS,h=60)
k=1
plot(D$date,S100[,k],col=couleur[k],pch=19,cex=.7,
xlab="",ylab="",ylim=c(0,40),xlim=
c(min(sondage$date),datefin))
rs=lm(sondage[,k+2]~bs(date),data=sondage)
lines(D$date,predict(rs,newdata=D),
lwd=2,col=couleur[k])
for(k in 2:4){
points(D$date,S100[,k],col=couleur[k],
pch=19,cex=.7)
rs=lm(sondage[,k+2]~bs(date),data=sondage)
lines(D$date,predict(rs,newdata=D),
lwd=2,col=couleur[k])
}
polygon(c(fdate,rev(fdate)),c(as.numeric(F$lower[,2]),
rev(as.numeric(F$upper[,2]))),col=CL[4],border=NA)
polygon(c(fdate,rev(fdate)),
c(as.numeric(F$lower[,1]),rev(as.numeric(F$upper[,1]))),
col=CL[5],border=NA)
lines(fdate,as.numeric(F$mean),lwd=2,col=CL[6])

On notera que je ne me suis pas trop fatigué sur la projection: autant sur la génération de sondages passés on a tenu compte de corrélation (i.e. si un candidat a un score élevé, ça se fera au détriment des autres, d’où la corrélation négative utilisée dans la simulation), autant ici, les projections sont faites de manière complétement indépendantes. En se fatiguant un peu plus, on pourrait se lancer dans un modèle vectoriel gaussien. Ou mieux, on pourrait travailler sur des processus de Dirichlet (surtout que R a un package dédié à ce genre de modèle) parce qu’on travaille depuis le début sur des taux, comme évoqué auparavant. Mais commençons par faire simple pour un premier billet rapide sur ce sujet.

On peut ensuite s’amuser à générer plusieurs jeux de sondages, et de lancer des prévisions dessus,

set.seed(1)
for(sim in 1:20){
Ssim=simulsondages(prediction4)
F1=forecast(ets(100*Ssim[,1]),h=60)
F2=forecast(ets(100*Ssim[,2]),h=60)
lines(fdate,as.numeric(F1$mean),col=CL[1])
lines(fdate,as.numeric(F2$mean),col=CL[6])}

Une fois qu’on a fait tout ça, on a presque fini. On peut regarder au 22 avril qui est en tête (voire, par scénario, calculer la probabilité qu’un des deux candidats soit en tête), c’est à dire au soir du 1er tour

set.seed(1)
VICTOIRE=rep(NA,1000)
for(sim in 1:1000){
Ssim=simulsondages(prediction4)
F1=forecast(ets(100*Ssim[,1]),h=60)
F2=forecast(ets(100*Ssim[,2]),h=60)
VICTOIRE[sim]=(F1$mean[30]>F2$mean[30])}
mean(VICTOIRE)

Et voilà. Ah oui, je n’ai pas laissé le résultat. Tout d’abord parce que je suis un statisticien, pas unprédicateur. Mais aussi parce que si ce genre de pronostic amuse des gens… ils n’ont qu’à se mettre à R pour faire tourner les codes !

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…

Sondages et probabilités en trois graphiques

Petite précision liminaire: je suis nul en sondage. Je n’ai jamais rien compris aux sondages. C’est la matière où j’ai eu ma pire note quand j’étais étudiant. Mais si on voit un sondage comme l’estimation d’une probabilité de succès, dans des tirages de lois binomiales (ou multinomiales), alors je me sens plus à l’aise pour commenter des choses lues ici ou là. Car en ces temps d’élections (qui approchent), on lit beaucoup de choses sur les sondages ou sur les intervalles de confiance des sondages, voire sur des comparaison de sondages… mais rares sont les justifications. Donc on va creuser un peu (histoire de voir ce qui peut être dit, dans un cadre largement simplifié).

  • “Il y a 1 chance sur 4 qu’un sondage donne une variation d’un point en + ou en – sur une proportion inchangée”

comme le disait @adelaigue (sur Twitter). Faire un sondage politique (si on résume simplement) c’est estimer la probabilité http://freakonometrics.hypotheses.org/files/2015/12/adic01.gif de voter pour A, sachant que l’on a alors une probabilité http://freakonometrics.hypotheses.org/files/2015/12/adic02.gif pour qu’une personne prise au hasard ne vote par pour A. Sur http://freakonometrics.hypotheses.org/files/2015/12/adic03.gif personnes interrogées, le nombre de personnes qui déclarera voter pour A est http://freakonometrics.hypotheses.org/files/2015/12/adic04.gif, qui suit une loi binomiale de paramètres http://freakonometrics.hypotheses.org/files/2015/12/adic03.gif et http://freakonometrics.hypotheses.org/files/2015/12/adic01.gif. Or par approximation (ou convergence) de la loi binomiale par une loi Gaussienne, on peut dire que

http://freakonometrics.hypotheses.org/files/2015/12/adic05.gif

Or l’estimateur naturel de http://freakonometrics.hypotheses.org/files/2015/12/adic01.gif est http://freakonometrics.hypotheses.org/files/2015/12/adic06.gif, i.e.

http://freakonometrics.hypotheses.org/files/2015/12/adic07.gif

Aussi, on en déduit l’intervalle de confiance classique, tel que http://freakonometrics.hypotheses.org/files/2015/12/adic01.gif appartient à cet intervalle avec une probabilité http://freakonometrics.hypotheses.org/files/2015/12/adic11.gif,

http://freakonometrics.hypotheses.org/files/2015/12/adic08.gif

http://freakonometrics.hypotheses.org/files/2015/12/adic12.gif désigne la fonction quantile de la loi normale centrée réduite; ou encore, en remplaçant les quantités inconnues par des estimateurs,

http://freakonometrics.hypotheses.org/files/2015/12/adic09.gif

Les bornes de l’intervalle de confiance de niveau http://freakonometrics.hypotheses.org/files/2015/12/adic11.gif sont alors

http://freakonometrics.hypotheses.org/files/2015/12/adic10.gif

C’est ce qu’on peut visualiser sur la figure ci-dessous,

Comme le notait @adelaigue, si on souhaite avoir

http://freakonometrics.hypotheses.org/files/2015/12/adic13.gif

(car on parle d’un point de base 100), il faudrait que

http://freakonometrics.hypotheses.org/files/2015/12/adic14.gif

Classiquement, http://freakonometrics.hypotheses.org/files/2015/12/adic03.gif est aux alentours de 1000. Et ce qu’on cherche c’est une probabilité d’appartenir à un intervalle, i.e. http://freakonometrics.hypotheses.org/files/2015/12/adic11.gif,

http://freakonometrics.hypotheses.org/files/2015/12/adi15.gif

Graphiquement, en fonction de http://freakonometrics.hypotheses.org/files/2015/12/adic01.gif on a alors

p=seq(0,1,by=.0025)
plot(p,2*(1-pnorm(.01/sqrt(p*(1-p)/1000))))

(ca sera notre premier graphique). Effectivement, entre 35% et 65%, il y a environ 50% de chances d’être autour de la fréquence observée, à un point près. Ou 50% de chances d’avoir plus d’un point de différence. Autrement dit, si on a observé avec 1000 personnes une fréquence de 40%, il y a une chance sur 4 que la probabilité soit entre 39% et 40%, une chance sur 4 qu’elle soit entre 40% et 41%, et surtout une chance sur 4 qu’elle soit inférieure à 39%, et 1 chance sur 4 qu’elle soit supérieure à 41%. Comme le notait @adelaigue, sur les bords, il semble que la formule ne soit pas valide. Mais si on prend justement les cas limites (probabilité très proche de 0, voire franchement nulle), est-ce qu’on peut dire quelque-chose ?

  • peut-on dire quelque chose (malgré tout) quand on ne peut rien dire ?

ou plus précisément, si en interrogeant http://freakonometrics.hypotheses.org/files/2015/12/adic03.gif personnes et que personne n’a déclaré vouloir voter pour A, quelle serait l’intervalle de confiance de la probabilité de voter pour A ? J’avais évoqué ce point dans un très vieux billet. La première piste est la réponse apporté par Pearson-Clopper (Pearson & Clopper (1934)), basée sur l’utilisation de la loi binomiale, et non plus d’une approximation.

library(Hmisc)
bornesup=function(N){ binom.test(x=0, n=N, p = 0,
alternative="less")$conf.int[2]*100 }
BS=Vectorize(bornesup)
vn=seq(10,500,by=10)
plot(vn,BS(vn))

Une alternative est de faire du bayésien. On utilise le fait que la loi conjuguée de la loi binomiale est la loi beta (comme évoqué dans des précédents billets). Si on suppose que la loi a priori de http://freakonometrics.hypotheses.org/files/2015/12/adic01.gifest une loi alors la loi a posteriori, sachant que sur http://freakonometrics.hypotheses.org/files/2015/12/adic03.gif lancés, on n’a observé que des 0, devient

bornesup=function(N){ qbeta(.95,1,1+N)*100 }
BS=Vectorize(bornesup)
lines(vn,BS(vn),lwd=2,col="blue")

La courbe ci-dessus correspond au calcul dans le cas du modèle de Pearson-Clopper, en rouge, et avec le quantile de la loi Beta en bleu (qui sont confondues sur ce graphique). Autrement dit, avec 300 personnes interrogées, si aucune ne dit vouloir voter pour A, l’intervalle de confiance pour http://freakonometrics.hypotheses.org/files/2015/12/adic01.gif est [0%;1%] (environ). Ca sera notre second graphique.

  • si on a trois alternatives, quelles sont les vraies probabilités de gagner ?

c’était déjà sur cette question que je concluais mon précédant billet, mais il est possible d’exploiter davantage cette idée de modèle bayésien avec un a priori suivant une loi beta (comme au dessus) ou de Dirichlet (si on a plus de deux alternatives). La première chose que l’on peut noter, c’est que s’il y a trois alternatives (A, B et C), et que la probabilité de voter pour A est 40%, la probabilité que A gagne… dépend des autres probabilités. On ne peut pas faire comme si A était seul, contre tous les autres:
– si les probabilités de voter pour B et C sont de l’ordre de 30%, on se dit que A a de fortes chances de gagner (oui, pour l’instant ce n’est pas très quantitatif).
– si les probabilités de voter pour B et C sont respectivement 55% et 5%, alors on se dit que A a de faibles chances de gagner.
Plus précisément, si on souhaite savoir si A va l’emporter ou pas, et que 200 personnes ont été interrogées, on peut essayer de faire une abaque, avec des courbes de niveau… dans le simplexe. Pour rappel, les scores obtenus par A, B et C se lisent sur le triangle ci-dessous (le point représenté correspond au triplet (40%,35%,25%))

Les probabilités sont les suivantes, avec n=200 (personnes qui se seraient exprimés)

où la région rouge correspond à une (quasi) certitude de gagner, pour A, et la région blanche, au contraire, une probabilité infinitésimale de gagner.

Ces résultats ont été obtenus par simulation,
vx=v[i]
vz=v[j]
vy=1-vz-vx
alpha=c(vx,vy,vz)
alpha=c(alpha,1-sum(alpha))*200
RD=rdirichlet(ns, alpha)
PROBA[i,j]=mean(RD[,1]>=apply(
cbind(RD[,3],RD[,2]),1,max))

Pour ceux que ne sont pas très à l’aise avec les excursions dans le simplexe, on peut se ramener dans un repère orthonormé. Ici, les scores obtenus par A et B se lisent sur le triangle ci-dessous (le point représenté correspond au triplet (40%,35%,25%), le même qu’auparavant

On a alors les probabilités de gagner suivantes, avec à droite la représentation sur le simplexe, en zoomant autour de notre triplet,

 

(ça sera le troisième et dernier graph). Si on regarde plus en détails, l’intervalle de confiance sur la probabilité de gagner dépend de la région où on se trouve dans le simplexe (un peu comme l’intervalle de confiance dépendait de la probabilité de l’emporter, dans la première partie de ce billet).

Maintenant, comme le rappelait @Christian dans un commentaire, tout cela est la version très simpliste de la pratique des sondages. Dans les vrais sondages (et pas des exercices de probabilités comme on vient d’en faire ici), on utilise des quotas, on fait du redessement, et c’est beaucoup plus complexe (enfin… c’est ce que j’essaye de me dire, pour me convaincre que la mauvaise note que j’avais obtenue lorsque j’étais étudiant était motivée).

La normalité ? avant tout une histoire de taille

Un billet rapide pour préciser le fond de ma pensée lorsque j’explique aux étudiants (qui travaillent sur leur devoir) que sur les bases de données avec beaucoup d’observations, il est très dur d’accepter l’hypothèse de normalité.

Pour ça, simulons des jeux de données: on a une seule variable explicative, un bruit gaussien, homoscédastique. Mais on va mal spécifier le modèle. En particulier, le vrai modèle ne sera pas linéaire. Voilà ce qu’on simule, et que l’on va chercher à ajuster.

n=50
X1=runif(n)
epsilon=rnorm(n)
Y=-2+8*X1-4*(X1<.2)+epsilon
base=data.frame(Y,X1)
reg=lm(Y~X1,data=base)

Autrement dit, les résidus (du modèle linéaire) ne devraient pas être gaussiens, car on n’ajuste pas le bon modèle. Maintenant, mon point est le suivant: avec peu de données, la plupart des tests de normalité ne vont pas rejeter l’hypothèse de normalité. Alors qu’avec beaucoup de données, l’hypothèse de normalité sera presque toujours rejetée.
Prenons par exemple le test Lillie, évoqué dans un vieux billet, qui est une adaptation du test de Kolmogorov-Smirnov quand on compare la distribution empirique avec la loi normale dont les paramètres ont été obtenues par maximum de vraisemblance. Sur le graphique ci-dessous, à gauche, on a en abscisse la taille de l’échantillon (qui va de 10 à 1000), et en ordonnée la p-value du test Lillie. Ou plus précisément, la courbe bleue est la région où s’est trouvée la p-value dans 90% des cas (la courbe supérieure étant le quantile à 95% de la p-value sur 10,000 échantillons simulés). La courbe rouge est la valeur moyenne de la p-value. A droite, on a – toujours en fonction de la taille de l’échantillon – la proportion des échantillons simulés pour lesquels la p-value était supérieure à 5%.

Autrement dit, plus la taille de l’échantillon est grande, plus on a de chance de rejeter (fort justement d’ailleurs) l’hypothèse de normalité. Et ce, uniquement sur un effet taille. Ou pour être plus rigoureux, sur les petits échantillons, on accepte beaucoup trop souvent (à tort cette fois) l’hypothèse de normalité.
On peut aussi regarder le test de Shapiro-Wilk

ou d’Anderson-Darling,

ou encore de Cramér-von Mises

Le code R pour faire ces tests est relativement simple,

X1=runif(n)
Y=-2+8*X1-4*(X1<.2)+rnorm(n)
base=data.frame(Y,X1)
reg=lm(Y~X1,data=base)
AD[n,i]=ad.test(residuals(reg))$p.value
SW[n,i]=shapiro.test(residuals(reg))$p.value
CV[n,i]=cvm.test(residuals(reg))$p.value
LL[n,i]=lillie.test(residuals(reg))$p.value

Pour ceux qui veulent aller plus loin, se cache derrière la notion de puissance du test…

Promesse et rationalité économique (à court terme)

Non, je vais pas faire une analyse des promesses en cette année électorale (bien que le sujet pourrait être amusant, on peut penser aux propos du candidat Newt Gingrich qui veut que la lune devienne un état américain, ici ou ). Non, en fait je voulais faire un billet sur le sujet les assureurs sont-ils vraiment des voleurs ? En effet, lorsqu’en classe on présente le fonctionnement d’un contrat d’assurance – en particulier lorsque l’on parle de tarification et de l’inversion du cycle de production, ou du provisionnement – on explique que l’assuré achète une promesse: celle de se faire indemniser les sinistres survenus pendant la période dite de couverture. Car formellement, c’est bien de cela qu’il s’agit. Sauf qu’il y a toujours quelqu’un pour faire noter que ce n’est peut-être pas aussi simple et que l’assureur, souvent, rechigne à payer.

En fait, il est parfaitement rationnel – d’un point de vue économique (on s’entend) – que l’assureur négocie. En tous les cas si l’assureur a une vision à très court terme… Et c’est ce que justifient Briegleb & Lemaire (1982), sur l’analyse du marchandage ou Lemaire (1982) (dont sera très largement inspiré ce billet).
Considérons un cas très simple: un assuré a subit une perte http://freakonometrics.blog.free.fr/public/perso5/tribunal33.gif (connue, de l’assuré comme de l’assureur). Mais l’assureur refuse de payer. Et l’assuré envisage d’aller au tribunal. Afin de voir s’il s’agit d’une décision rationnelle, on suppose que l’assuré (http://freakonometrics.blog.free.fr/public/perso5/tribunal18.gif) avait une richesse  avant que le sinistre ne survienne, qu’aller en justice lui coûte , et qu’il estime que sa probabilité de gagner le procès est http://freakonometrics.blog.free.fr/public/perso5/tribunal05.gif. De son coté, la compagnie (http://freakonometrics.blog.free.fr/public/perso5/tribunal19.gif) d’assurance a une richesse , qu’aller en justice lui coute http://freakonometrics.blog.free.fr/public/perso5/tribunal04.gif, et que la probabilité que la compagnie perde le procès est – selon la compagnie – http://freakonometrics.blog.free.fr/public/perso5/tribunal06.gif (les http://freakonometrics.blog.free.fr/public/perso5/tribunal.gif sont alors les probabilités que l’assuré gagne son procès, mais vu des deux cotés). Au tribunal, on suppose que le juge rend un verdict du genre tout ou rien (on n’autorise pas ici de compensation plus importantes que le cout réel du sinistre, et on ne parle pas non plus de remboursement des frais d’avocat en cas de victoire du procès).
L’espérance d’utilité de l’assuré est

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

alors que pour la compagnie d’assurance, son espérance d’utilité est

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

La situation est Pareto optimale si

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

Si on suppose que les agents sont risquophobe (i.e. http://freakonometrics.blog.free.fr/public/perso5/tribunal10.gifdécroissante), alors

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

alors que

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

Aussi, une condition nécessaire pour que la situation soit Pareto-optimale est que

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

i.e. http://freakonometrics.blog.free.fr/public/perso5/tribuanl16.gif. Aller au tribunal n’est pas Pareto optimal si http://freakonometrics.blog.free.fr/public/perso5/tribunal18.gif croit plus en ses raisons de gagner que http://freakonometrics.blog.free.fr/public/perso5/tribunal19.gif. Si http://freakonometrics.blog.free.fr/public/perso5/tribunal17.gif, il peut alors être optimal de ne pas aller au tribunal…. On ne dit pas ici qu’il faut que la compagnie indemnise (intégralement) l’assuré, mais que l’assureur et l’assuré peuvent avoir intérêt a trouver une solution amiable…
Si http://freakonometrics.blog.free.fr/public/perso5/tribuanl16.gif, il peut être intéressant d’aller au tribunal… mais pas forcément… Du point de vue de l’assureur, on peut supposer que sa richesse lui permet de ne pas être risquophobe, ou au moins d’être risque neutre (face a cette assuré au moins). Elle est donc indifférente entre aller au tribunal (et avoir un rendement incertain) et verser un règlement amiable si le versement correspond à son espérance de perte si elle va au procès. Autrement dit http://freakonometrics.blog.free.fr/public/perso5/tribunal20.gif. L’assuré acceptera cette somme si l’utilité qu’il en tire est plus grande que l’espérance d’utilité qu’il aura en allant au tribunal,

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

La différence entre le terme de gauche et le terme de droite est une fonction http://freakonometrics.blog.free.fr/public/perso5/tribunal22.gif. Si l’agent a une utilité exponentielle

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

alors

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

doit être une fonction positive. Un calcul rapide de dérivées montre que cette fonction est convexe, avec un minimum atteint en

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

Moralité ?

  • pour les petits et pour les gros sinistres, il est normal que l’assureur et l’assuré cherchent une solution amiable (transaction)
  • en revanche, pour les sinistres de cout intermédiaire, il peut être légitime d’aller au tribunal (procès).

Maintenant, si on a les grandes directions… peut-être peut-on essayer de mieux comprendre quelle serait le montant de la transaction dans le cas ou une solution amiable est envisagée. Nash (1950) a envisagé un jeu non-coopératif afin de décrire le marchandage. La transaction optimale est la quantité http://freakonometrics.blog.free.fr/public/perso5/tribunal28.gif qui maximise le produit des deux gains d’utilité (cf Nash (1950) ou Roth (1979))

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

La condition du premier ordre donne l’équation suivante

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

Si on fait un peu de calcul numérique, on note que la solution http://freakonometrics.blog.free.fr/public/perso5/tribunal32.gif est décroissante en http://freakonometrics.blog.free.fr/public/perso5/tribunal31.gif: plus l’assuré est risquophobe, plus faible sera la montant de la transaction…. Ce qui n’est pas surprenant, on retrouve ici un résultat noté dans Peeters (1981), Roth (1982) ou ecnore Osbourne (1985). Maintenant, sur le modèle théorique on pourrait bien entendu argumenter que le tribunal fonctionne selon une règle trop simple du tout ou rien, et qu’un procès peut couter beaucoup cher que le simple remboursement du cout du sinistre. Ou qu’il pourrait créer une jurisprudence. On pourrait aussi regarder économétriquement le lien entre le cout des sinistres et les montants des ententes amiables (c’est ce que fait Ayusoa, Bermúdezb & Santolinoc (2012)) mais c’est une autre histoire…. ou ça pourrait être le sujet pour une autre histoire… à suivre donc.

Sur le lissage exponentiel

Sur le lissage exponentiel, je pourrais renvoyer vers Gardner (1985), et la version remise au goût du jour, Gardner (2005). J’avais pris comme position, dans le cours, de présenter rapidement le lissage exponentiel, en notant qu’on les évoquerais à nouveau plus tard comme cas particulier des modèles ARIMA. Comme le note la dernière version, Gardner (2005), ce n’est pas forcément l’unique manière de voir, et le cours pourrait être orienté complètement autour des notions de lissage exponentiel (c’est d’ailleurs le point du vue de livre Hyndman, Koehler, Ord & Snyder (2008)), “when Gardner (2005) appeared, many believed that exponential smoothing should be disregarded because it was either a special case of ARIMA modeling or an ad hoc procedure with no statistical rationale. As McKenzie (1985) observed, this opinion was expressed in numerous references to my paper. Since 1985, the special case argument has been turned on its head, and today we know that exponential smoothing methods are optimal for a very general class of state-space models that is in fact broader than the ARIMA class.”

N’étant toujours pas convaincu, j’évoquerais à nouveau le lissage exponentiel quand on présentera les modèles ARIMA. En attendant, un peu de code pour mieux comprendre ce qui est fait quand on fait du lissage exponentiel.

Commençons par le lissage exponentiel simple, i.e.

http://freakonometrics.hypotheses.org/files/2015/12/lex03.gif

http://freakonometrics.hypotheses.org/files/2015/12/le04.gif désigne un poids attribué à la nouvelle observation dans la fonction de lissage http://freakonometrics.hypotheses.org/files/2015/12/lex05.gif. Le code pour lisser une série est alors assez simple,

> library(datasets)
> X=as.numeric(Nile)
> Lissage=function(a){
+  T=length(X)
+  L=rep(NA,T)
+  L[1]=X[1]
+  for(t in 2:T){L[t]=a*X[t]+(1-a)*L[t-1]}
+  return(L)
+ }
> plot(X,type="b",cex=.6)
> lines(Lissage(.2),col="red")

http://freakonometrics.hypotheses.org/files/2015/12/lissage-exp-1.gif

Sur la figure ci-dessus, non visualise l’impact de http://freakonometrics.hypotheses.org/files/2015/12/lisse16.gif sur le lissage. La prévision que l’on fait à la date http://freakonometrics.hypotheses.org/files/2015/12/lex14.gif, pour un horizon http://freakonometrics.hypotheses.org/files/2015/12/lisse12.gif est alors http://freakonometrics.hypotheses.org/files/2015/12/lisse11.gif. Il est alors possible de voir le poids comme un paramètre et on va alors chercher le poids optimal. La stratégie classique est de minimiser l’erreur de prédiction commise à un horizon de 1

http://freakonometrics.hypotheses.org/files/2015/12/lisse14.gif

> V=function(a){
+  T=length(X)
+  L=erreur=rep(NA,T)
+  erreur[1]=0
+  L[1]=X[1]
+  for(t in 2:T){
+  L[t]=a*X[t]+(1-a)*L[t-1]
+  erreur[t]=X[t]-L[t-1] }
+ return(sum(erreur^2))
+ }

Ici, on obtient comme poids optimal

> A=seq(0,1,by=.02)
> Ax=Vectorize(V)(A)
> plot(A,Ax,ylim=c(min(Ax),min(Ax)*1.05))
> optimize(V,c(0,.5))$minimum
[1] 0.246581

On notera que c’est ce que suggère la fonction de R,

> hw=HoltWinters(X,beta=FALSE,gamma=FALSE,l.start=X[1])
> hw
Holt-Winters exponential smoothing without trend an seasonal comp.

Call:
HoltWinters(x = X, beta = FALSE, gamma = FALSE, l.start = X[1])

Smoothing parameters:
alpha:  0.2465579
beta :  FALSE
gamma:  FALSE

Coefficients:
[,1]
a 805.0389

> plot(hw)
> points(2:(length(X)+1),Vectorize(Lissage)(.2465),col="blue")

Dans le cas du lissage exponentiel double, i.e.

http://freakonometrics.hypotheses.org/files/2015/12/lex50.gif

et

http://freakonometrics.hypotheses.org/files/2015/12/HW02.gif

Dans ce cas, la prédiction est http://freakonometrics.hypotheses.org/files/2015/12/lisse13.gif. Le code pour faire un lissage double est là encore assez simple,

> Lissage=function(a,b){
+  T=length(X)
+  L=B=rep(NA,T)
+  L[1]=X[1]; B[1]=0
+  for(t in 2:T){
+  L[t]=a*X[t]+(1-a)*(L[t-1]+B[t-1])
+  B[t]=b*(L[t]-L[t-1])+(1-b)*B[t-1] }
+ return(L)
+ }

Sur la figure suivante, on visualise l’évolution du lissage en fonction de http://freakonometrics.hypotheses.org/files/2015/12/lisse15.gif (http://freakonometrics.hypotheses.org/files/2015/12/lisse16.gif étant ici fixé),
http://freakonometrics.hypotheses.org/files/2015/12/lissage-exp-2.gif
(le lissage simple – avec le même poids http://freakonometrics.hypotheses.org/files/2015/12/lisse16.gif – apparaissant ici en trait clair).

Visualization in regression analysis

Visualization is a key to success in regression analysis. This is one of the (many) reasons I am also suspicious when I read an article with a quantitative (econometric) analysis without any graph. Consider for instance the following dataset, obtained from http://data.worldbank.org/, with, for each country, the GDP per capita (in some common currency) and the infant mortality rate (deaths before the age of 5),

> library(gdata)
> XLS1=read.xls("http://api.worldbank.org/datafiles
/NY.GDP.PCAP.PP.CD_Indicator_MetaData_en_EXCEL.xls", sheet = 1)
> data1=XLS1[-(1:28),c("Country.Name","Country.Code","X2010")]
> names(data1)[3]="GDP"
> XLS2=read.xls("http://api.worldbank.org/datafiles
/SH.DYN.MORT_Indicator_MetaData_en_EXCEL.xls", sheet = 1)
> data2=XLS2[-(1:28),c("Country.Code","X2010")]
> names(data2)[2]="MORTALITY"
> data=merge(data1,data2)
> head(data)
Country.Code         Country.Name       GDP MORTALITY
1          ABW                Aruba        NA        NA
2          AFG          Afghanistan  1207.278     149.2
3          AGO               Angola  6119.930     160.5
4          ALB              Albania  8817.009      18.4
5          AND              Andorra        NA       3.8
6          ARE United Arab Emirates 47215.315       7.1

If we estimate a simple linear regression – http://freakonometrics.blog.free.fr/public/perso5/logormal01.gif  – we get

> regBB=lm(MORTALITY~GDP,data=data)
> summary(regBB)

Call:
lm(formula = MORTALITY ~ GDP, data = data)

Residuals:
Min     1Q Median     3Q    Max
-45.24 -29.58 -12.12  16.19 115.83

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 67.1008781  4.1577411  16.139  < 2e-16 ***
GDP         -0.0017887  0.0002161  -8.278 3.83e-14 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 39.99 on 167 degrees of freedom
(47 observations deleted due to missingness)
Multiple R-squared: 0.2909,	Adjusted R-squared: 0.2867
F-statistic: 68.53 on 1 and 167 DF,  p-value: 3.834e-14

We can look at the scatter plot, including the linear regression line, and some confidence bounds,

> plot(data$GDP,data$MORTALITY,xlab="GDP per capita",
+ ylab="Mortality rate (under 5)",cex=.5)
> text(data$GDP,data$MORTALITY,data$Country.Name,pos=3)
> x=seq(-10000,100000,length=101)
> y=predict(regBB,newdata=data.frame(GDP=x),
+ interval="prediction",level = 0.9)
> lines(x,y[,1],col="red")
> lines(x,y[,2],col="red",lty=2)
> lines(x,y[,3],col="red",lty=2)

We should be able to do a better job here. For instance, if we look at the Box-Cox profile likelihood,

> boxcox(regBB)

it looks like taking the logarithm of the mortality rate should be better, i.e. http://freakonometrics.blog.free.fr/public/perso5/lognormal02.gif or http://freakonometrics.blog.free.fr/public/perso5/lognormal05.gif:

> regLB=lm(log(MORTALITY)~GDP,data=data)
> summary(regLB)

Call:
lm(formula = log(MORTALITY) ~ GDP, data = data)

Residuals:
Min      1Q  Median      3Q     Max
-1.3035 -0.5837 -0.1138  0.5597  3.0583

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)  3.989e+00  7.970e-02   50.05   <2e-16 ***
GDP         -6.487e-05  4.142e-06  -15.66   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.7666 on 167 degrees of freedom
(47 observations deleted due to missingness)
Multiple R-squared: 0.5949,	Adjusted R-squared: 0.5925
F-statistic: 245.3 on 1 and 167 DF,  p-value: < 2.2e-16

> plot(data$GDP,data$MORTALITY,xlab="GDP per capita",
+ ylab="Mortality rate (under 5) log scale",cex=.5,log="y")
> text(data$GDP,data$MORTALITY,data$Country.Name)
> x=seq(300,100000,length=101)
> y=exp(predict(regLB,newdata=data.frame(GDP=x)))*
+ exp(summary(regLB)$sigma^2/2)
> lines(x,y,col="red")
> y=qlnorm(.95, meanlog=predict(regLB,newdata=data.frame(GDP=x)),
+ sdlog=summary(regLB)$sigma^2)
> lines(x,y,col="red",lty=2)
> y=qlnorm(.05, meanlog=predict(regLB,newdata=data.frame(GDP=x)),
+ sdlog=summary(regLB)$sigma^2)
> lines(x,y,col="red",lty=2)

on the log scale or

> plot(data$GDP,data$MORTALITY,xlab="GDP per capita",
+ ylab="Mortality rate (under 5) log scale",cex=.5)

on the standard scale. Here we use quantiles of the log-normal distribution to derive confidence intervals.

But why shouldn’t we take also the logarithm of the GDP ? We can fit a model http://freakonometrics.blog.free.fr/public/perso5/lognormal03.gif or equivalently http://freakonometrics.blog.free.fr/public/perso5/lognormal04.gif.

> regLL=lm(log(MORTALITY)~log(GDP),data=data)
> summary(regLL)

Call:
lm(formula = log(MORTALITY) ~ log(GDP), data = data)

Residuals:
Min       1Q   Median       3Q      Max
-1.13200 -0.38326 -0.07127  0.26610  3.02212

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 10.50192    0.31556   33.28   <2e-16 ***
log(GDP)    -0.83496    0.03548  -23.54   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.5797 on 167 degrees of freedom
(47 observations deleted due to missingness)
Multiple R-squared: 0.7684,	Adjusted R-squared: 0.767
F-statistic:   554 on 1 and 167 DF,  p-value: < 2.2e-16

> plot(data$GDP,data$MORTALITY,xlab="GDP per capita ",
+ ylab="Mortality rate (under 5)",cex=.5,log="xy")
> text(data$GDP,data$MORTALITY,data$Country.Name)
> x=exp(seq(1,12,by=.1))
> y=exp(predict(regLL,newdata=data.frame(GDP=x)))*
+ exp(summary(regLL)$sigma^2/2)
> lines(x,y,col="red")
> y=qlnorm(.95, meanlog=predict(regLL,newdata=data.frame(GDP=x)),
+ sdlog=summary(regLL)$sigma^2)
> lines(x,y,col="red",lty=2)
> y=qlnorm(.05, meanlog=predict(regLL,newdata=data.frame(GDP=x)),
+ sdlog=summary(regLL)$sigma^2)
> lines(x,y,col="red",lty=2)

on the log scales or

> plot(data$GDP,data$MORTALITY,xlab="GDP per capita ",
+ ylab="Mortality rate (under 5)",cex=.5)

on the standard scale. If we compare the last two predictions, we have

with in blue is the log model, and in red is the log-log model (I did not include the first one for obvious reasons).

Examen intra, modèles de régression

L’examen intra comportera deux exercices proches de ceux faits en démonstration, et une analyse de sorties de régression. Cette dernière partie portera sur la base suivante (ou peut être une extraction de cette base)

> examen=read.table(
+ "http://freakonometrics.blog.free.fr/public/
data/basket-exam-v2.csv",
+ header=TRUE,sep=";")
> tail(examen)
EQUIPE halfdiff winner
9917                 MARYLAND       -1      0
9918             RHODE ISLAND      -15      1
9919            ROBERT MORRIS        4      1
9920 NORTH CAROLINA-ASHEVILLE       11      1
9921                    BROWN       -2      1
9922          MOUNT ST. MARYS        1      0
totalpointsdiff secondhalfpoints
9917             -10               26
9918               1               44
9919               6               43
9920              19               39
9921               6               48
9922              -9               26
nodiff_pts_1000_left hometeam halfwinner
9917                    9        1          0
9918                   18        1          0
9919                   17        1          1
9920                   20        1          1
9921                   16        1          0
9922                   13        1          1
halflooser secondhalfpointsother
9917          1                    35
9918          1                    28
9919          0                    41
9920          0                    31
9921          1                    40
9922          0                    36

qui contient des matchs de baskets universitaires.

  • EQUIPE le nom de l’équipe qui sert de référence
  • HALFDIFF le nombre de points de différence à la mi-temps (négatif signifie que l’EQUIPE perdait)
  • WINNER vaut 1 si l’équipe a gagné le match, 0 sinon,
  • TOTALPOINTSDIFF le nombre de points de différence à la fin du match
  • SECONDHALFPOINTS est le nombre de points marqués par l’EQUIPE pendant la seconde mi-temps
  • NODIFF_PTS_1000-LEFT la différence de points entre les équipes 10 minutes avant la fin du match (au milieu de la seconde période)
  • HOMETEAM vaut 0 si l’équipe joue à l’extérieur, 1 si elle joue à domicile
  • HALFWINNER vaut 1 si l’EQUIPE gagnait (supérieur ou nul) à la mi-temps
  • HALFLOOSER vaut 1 si l’EQUIPE perdait (inférieur ou nul) à la mi-temps
  • SECONDHALFPOINTSOTHER est le nombre de points marqués par l’adversaire au cours de la seconde mi-temps

Tests on tail index for extremes

Since several students got the intuition that natural catastrophes might be non-insurable (underlying distributions with infinite mean), I will post some comments on testing procedures for extreme value models.

A natural idea is to use a likelihood ratio test (for composite hypotheses). Let http://freakonometrics.blog.free.fr/public/perso5/lrtest21.gif denote the parameter (of our parametric model, e.g. the tail index), and we would like to know whether http://freakonometrics.blog.free.fr/public/perso5/lrtest21.gif is smaller or larger than http://freakonometrics.blog.free.fr/public/perso5/lrtest22.gif (where in the context of finite versus infinite mean http://freakonometrics.blog.free.fr/public/perso5/lrtest23.gif). I.e. either http://freakonometrics.blog.free.fr/public/perso5/lrtest21.gif belongs to the set http://freakonometrics.blog.free.fr/public/perso5/lrtest-10.gif or to its complementary http://freakonometrics.blog.free.fr/public/perso5/lrtest-11.gif. Consider the maximum likelihood estimator http://freakonometrics.blog.free.fr/public/perso5/lrtest24.gif, i.e.

http://freakonometrics.blog.free.fr/public/perso5/lrtest-9.gif

Let http://freakonometrics.blog.free.fr/public/perso5/lrtest25.gif and http://freakonometrics.blog.free.fr/public/perso5/lrtest-3.gif denote the constrained maximum likelihood estimators on http://freakonometrics.blog.free.fr/public/perso5/lrtest26.gif and http://freakonometrics.blog.free.fr/public/perso5/lrtest27.gif respectively,

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

http://freakonometrics.blog.free.fr/public/perso5/lrtest-2.gif

Either http://freakonometrics.blog.free.fr/public/perso5/lrtest-13.gif and http://freakonometrics.blog.free.fr/public/perso5/lrtest-6.gif (on the left), or http://freakonometrics.blog.free.fr/public/perso5/lrtest-14.gif and http://freakonometrics.blog.free.fr/public/perso5/lrtest-7.gif (on the right)

So likelihood ratios

http://freakonometrics.blog.free.fr/public/perso5/lrtest-15.gif      http://freakonometrics.blog.free.fr/public/perso5/lrtest-16.gif

 are either equal to

http://freakonometrics.blog.free.fr/public/perso5/lrtest-19.gif      http://freakonometrics.blog.free.fr/public/perso5/lrtest-18.gif

or

http://freakonometrics.blog.free.fr/public/perso5/lrtest-20.gif        http://freakonometrics.blog.free.fr/public/perso5/lrtest-17.gif

If we use the code mentioned in the post on profile likelihood, it is easy to derive that ratio. The following graph is the evolution of that ratio, based on a GPD assumption, for different thresholds,

> base1=read.table(
+ "http://freakonometrics.free.fr/danish-univariate.txt",
+ header=TRUE)
> library(evir)
> X=base1$Loss.in.DKM
> U=seq(2,10,by=.2)
> LR=P=ES=SES=rep(NA,length(U))
> for(j in 1:length(U)){
+ u=U[j]
+ Y=X[X>u]-u
+ loglikelihood=function(xi,beta){
+ sum(log(dgpd(Y,xi,mu=0,beta))) }
+ XIV=(1:300)/100;L=rep(NA,300)
+ for(i in 1:300){
+ XI=XIV[i]
+ profilelikelihood=function(beta){
+ -loglikelihood(XI,beta) }
+ L[i]=-optim(par=1,fn=profilelikelihood)$value }
+ plot(XIV,L,type="l")
+ PL=function(XI){
+ profilelikelihood=function(beta){
+ -loglikelihood(XI,beta) }
+ return(optim(par=1,fn=profilelikelihood)$value)}
+ (L0=(OPT=optimize(f=PL,interval=c(0,10)))$objective)
+ profilelikelihood=function(beta){
+ -loglikelihood(1,beta) }
+ (L1=optim(par=1,fn=profilelikelihood)$value)
+ LR[j]=L1-L0
+ P[j]=1-pchisq(L1-L0,df=1)
+ G=gpd(X,u)
+ ES[j]=G$par.ests[1]
+ SES[j]=G$par.ses[1]
+ }
>
> plot(U,LR,type="b",ylim=range(c(0,LR)))
> abline(h=qchisq(.95,1),lty=2)

with on top the values of the ratio (the dotted line is the quantile of a chi-square distribution with one degree of freedom) and below the associated p-value

> plot(U,P,type="b",ylim=range(c(0,P)))
> abline(h=.05,lty=2)

In order to compare, it is also possible to look at confidence interval for the tail index of the GPD fit,

> plot(U,ES,type="b",ylim=c(0,1))
> lines(U,ES+1.96*SES,type="h",col="red")
> abline(h=1,lty=2)

To go further, see Falk (1995), Dietrich, de Haan & Hüsler (2002), Hüsler & Li (2006) with the following table, or Neves & Fraga Alves (2008). See also here or there (for the latex based version) for an old paper I wrote on that topic.

Lissage exponentiel et séries temporelles

L’idée du lissage exponentiel est de supposer que http://freakonometrics.hypotheses.org/files/2015/12/lex01.gif, où on suppose que le bruit http://freakonometrics.hypotheses.org/files/2015/12/lex02.gif est ici un bruit indépendant (imprévisible). On va alors supposer que

http://freakonometrics.hypotheses.org/files/2015/12/lex03.gif

où http://freakonometrics.hypotheses.org/files/2015/12/le04.gif désigne un poids attribué à la nouvelle observation dans la fonction de lissage http://freakonometrics.hypotheses.org/files/2015/12/lex05.gif.

On parle de lissage pour moyenne mobile à poids exponentiel (ewma, exponentially weighted moving average, évoqué ici ou  sur ce blog). En effet

http://freakonometrics.hypotheses.org/files/2015/12/lex08.gif

i.e., en itérant

http://freakonometrics.hypotheses.org/files/2015/12/lex09.gif

(même s’il faut faire attention à cette décomposition qui, en pratique, s’arrête avec la première observation).

La méthode de Holt-Winters propose de généraliser ce lissage, en introduisant deux (voire trois) composantes.

  • http://freakonometrics.hypotheses.org/files/2015/12/lex05.gif sera interprété comme un niveau
  • http://freakonometrics.hypotheses.org/files/2015/12/lex11.gif correspondant à la variation du processus
  • http://freakonometrics.hypotheses.org/files/2015/12/lex12.gif qui correspond à un cycle (éventuel)

On suppose désormais, dans cette méthode (dite de Holt-Winters), que

http://freakonometrics.hypotheses.org/files/2015/12/lex50.gif

et

http://freakonometrics.hypotheses.org/files/2015/12/HW02.gif

D’un point de vue plus pratique, la programmation se fait comme suit,

> HW=HoltWinters(X,alpha=.2,beta=0)
> plot(HW)

Cette trois composantes seront utilisées pour faire une prévision de la série: la prédiction sera de la forme http://freakonometrics.hypotheses.org/files/2015/12/lex16.gif, à la date http://freakonometrics.hypotheses.org/files/2015/12/lex14.gif, pour un horizon http://freakonometrics.hypotheses.org/files/2015/12/lex17.gif (avec les notations du précédant billet).

> P=predict(HW,24,prediction.interval=TRUE)
> plot(HW,xlim=range(c(time(X),time(P))))
> polygon(c(time(P),rev(time(P))),c(P[,2],rev(P[,3])),
+ col="yellow",border=NA)
> lines(P[,1],col="red",lwd=3)

What is working as an actuary really like?

Just a short post to mention the page entitled “What is working as an actuary really like ?” on http://beanactuary.com/. I wanted to share that page with some students (even if it is quite incomplete) because I am always surprised when I hear that what we teach during lectures is too mathematical (or even to much based on computing skills). Practice can be much more complex than what we do in our courses (e.g. variable selection in ratemaking, integrating expert opinion in claims reserving, or hedging life insurance contracts with several complex financial options). If practitioners want to share their experience, comments are welcome… I just hope that some students will stop thinking that Barney Stinson is doing the job of an actuary… and that what we teach might – somehow – be useful someday…

the Dirichlet distribution

In the course, since we are still introducing some concepts of dependent distributions, we will talk about the Dirichlet distribution, which is a distribution over the simplex of http://freakonometrics.hypotheses.org/files/2017/07/diri11.gif. Let http://freakonometrics.hypotheses.org/files/2017/07/diri01.gif denote the Gamma distribution with density (on http://freakonometrics.hypotheses.org/files/2017/07/diri03.gif)

http://freakonometrics.hypotheses.org/files/2017/07/diri02.gif

Let http://freakonometrics.hypotheses.org/files/2017/07/diri04.gif denote independent http://freakonometrics.hypotheses.org/files/2017/07/diri05.gif random variables, with http://freakonometrics.hypotheses.org/files/2017/07/diri06.gif. Then http://freakonometrics.hypotheses.org/files/2017/07/diri07.gif where

http://freakonometrics.hypotheses.org/files/2017/07/diri08.gif

has a Dirichlet distribution with parameter

http://freakonometrics.hypotheses.org/files/2017/07/diri09.gif

Note that http://freakonometrics.hypotheses.org/files/2017/07/diri10.gif has a distribution in the simplex of http://freakonometrics.hypotheses.org/files/2017/07/diri11.gif,

http://freakonometrics.hypotheses.org/files/2017/07/diri40.gif

and has density

http://freakonometrics.hypotheses.org/files/2017/07/diri12.gif

We will write http://freakonometrics.hypotheses.org/files/2017/07/diri13.gif.

The density for different values of http://freakonometrics.hypotheses.org/files/2017/07/diri20.gif can be visualized below, e.g. http://freakonometrics.hypotheses.org/files/2017/07/diri21.gif, with some kind of symmetry,
http://freakonometrics.hypotheses.org/files/2017/07/dirichlet222.gif
or http://freakonometrics.hypotheses.org/files/2017/07/diri22.gif and http://freakonometrics.hypotheses.org/files/2017/07/diri23.gif, below
http://freakonometrics.hypotheses.org/files/2017/07/dirichlet522.gif
and finally, below, http://freakonometrics.hypotheses.org/files/2017/07/diri24.gif


Note that marginal distributions are also Dirichlet, in the sense that if

http://freakonometrics.hypotheses.org/files/2017/07/diri13.gif

then

http://freakonometrics.hypotheses.org/files/2017/07/diri14.gif

if http://freakonometrics.hypotheses.org/files/2017/07/diri15.gif, and if http://freakonometrics.hypotheses.org/files/2017/07/diri16.gif, then http://freakonometrics.hypotheses.org/files/2017/07/diri17.gif‘s have Beta distributions,

http://freakonometrics.hypotheses.org/files/2017/07/diri18.gif

See Devroye (1986) section XI.4, or Frigyik, Kapila & Gupta (2010) .This distribution might also be called multivariate Beta distribution. In R, this function can be used as follows

> library(MCMCpack)
> alpha=c(2,2,5)
> x=seq(0,1,by=.05)
> vx=rep(x,length(x))
> vy=rep(x,each=length(x))
> vz=1-x-vy
> V=cbind(vx,vy,vz)
> D=ddirichlet(V, alpha)
> persp(x,x,matrix(D,length(x),length(x))

(to plot the density, as figures above). Note that we will come back on that distribution later on so-called Liouville copulas (see also Gupta & Richards (1986)).

"sendo l'intento mio scrivere cosa utile a chi la intende…"