Arthur Charpentier

Aller au contenu | Aller au menu | Aller à la recherche

lundi 8 février 2010

Le billet de la Saint Valentin

Avec une petite semaine d’avance, un billet rapide pour préparer la Saint Valentin, i.e. récupérer des conseils d’experts. Il existe un site http://www.okcupid.com/ de Free Online Dating pour reprendre le titre. Mais son grand intérêt, c’est surtout le blog associé, qui propose des analyses statistiques sur les matching, http://blog.okcupid.com/.

Par exemple, les fondateurs (qui se sont amusés à faire beaucoup de statistiques) regardent les photos, et notent que pour espérer obtenir beaucoup de messages (ou tout du moins plus que la moyenne), les femmes doivent regarder l’objectif tout en souriant, en adoptant une moue un peu dragueuse ("flirty face") (cf Myth 1 ici). Les hommes, eux, doivent plutôt tirer la tronche et regarder ailleurs.... Mais le plus amusant, c’est - je m’adresse aux hommes qui lisent mon blog - il y a une étude sur ce qui doit apparaître sur la photo,


Les résultats sont sans ambiguité: il faut être avec un animal (je ne sais pas si c’est un animal de compagnie, genre un gentil chiot du calendrier de La Poste, ou bien avec un boa autour du cou...). Et il peut être intéressant de montrer ses muscles ! Effectivement, on peut montrer ses abdos, mais avec l’âge, il commence à être plus sage de les cacher sous un gros pull....

alors que la tendance est inversée pour les décolletés des femmes,

qui ont d’autant plus intérêt à se découvrir qu’elles sont âgées... Ils ont aussi fait des études sur les correspondances raciales, en étudiant le pourcentage de matching,

qui est assez uniforme, et en le comparant au pourcentage de réponses au messages de dragues...

On observe que la "race" a un impact sur le fait qu’on se laisse draguer ou pas, mais ensuite, on y devient indifférent. On retrouve aussi l’étude sur la religion,

De manière plus amusante, l’étude est faite aussi sur les correspondances zodiacales (ici)

Bref, enfin des statistiques amusantes et instrucives sur un sujet qui intéressent tout le monde (pour citer le mémoire de nos étudiants qui avaient fait leur projet d’économétrie sur les relations amoureuses...)

samedi 6 février 2010

Provisionnement et surdispersion

Un billet rapide suite à une question que j’ai reçu par mail, au sujet des triangles de liquidation, "pourquoi tout le monde fait du ODP alors que le résultat est le même que Poisson normale ?". Le plus simple est de regarder effectivement, sur un triangle, ce qui se passe si l’on prend des modèles intégrant de la surdispersion, au lieu d’un modèle de Poisson plus classique. Pour commencer, ajustons un modèle de Poisson sur le triangle des incréments,
>  library(ChainLadder)
>  an <- 10; ligne = rep(1:an, each=an); colonne = rep(1:an, an)
>  passe = (ligne + colonne - 1)<=an; n = sum(passe)
>  PAID=GenIns; INC=PAID
>  INC[,2:an]=PAID[,2:an]-PAID[,1:(an-1)]
>  Y = as.vector(INC)
>  lig = as.factor(ligne)
>  col = as.factor(colonne)
>  base = data.frame(Y,col,lig)
>  reg1=glm(Y~col+lig,data=base,family="poisson")
>  sum(exp(predict(reg1,newdata=base))[passe!=TRUE])
[1] 18680856
Si on utilise un modèle de Poisson surdispersé, la prédiction sera la même,
>  reg2=glm(Y~col+lig,data=base,family="quasipoisson")
>  sum(exp(predict(reg2,newdata=base))[passe!=TRUE])
[1] 18680856
En effet, seule la variance change, i.e. l’incertitude sera plus grande si l’on utilise un modèle de Poisson surdispersion. Notons que les tests permettent de valider l’hypothèse de surdispersion,
> dispersiontest(reg1)
        Overdispersion test
data:  reg1
z = 4.3942, p-value = 5.558e-06
alternative hypothesis: true dispersion is greater than 1

Une autre idée peut être de faire une régression binomiale négative (qui permet de prendre en compte la surdispersion),
>  library(MASS)
>  reg3=glm.nb(Y~col+lig,data=base)
> summary(reg3)
(Dispersion parameter for Negative Binomial(13.8349) family taken to be 1)
              Theta:  13.83
          Std. Err.:  2.61
 2 x log-likelihood:  -1460.766

On peut alors calculer à nouveau le montant de provision total,
>  sum(exp(predict(reg3,newdata=base))[passe!=TRUE])
[1] 18085795
qui est légèrement différent du cas Poissonnien.

On peut d’ailleurs aller un peu plus loin, et regarder les mse de prédiction, pour chacun des modèles (que j’avais évoqué ici). Mais auparavant, pour des raisons techniques, il est plus simple de se débarrasser des valeurs nulles dans le triangle,
>  an <- 10; ligne = rep(1:an, each=an); colonne = rep(1:an, an)
>  passe = (ligne + colonne - 1)<=an; np = sum(passe)
>  futur = (ligne + colonne - 1)> an; nf = sum(passe)
>  base$Y2=base$Y; base$Y2[is.na(Y)]=.001
>  reg1b=glm(Y2~lig+col, family=poisson,data=base)
>  reg2b=glm(Y2~lig+col, family=quasipoisson,data=base)
>  reg3b=glm.nb(Y2~lig+col, data=base)

La fonction suivante calcule l’estimateur du montant total de réserve, ainsi que la prediction error (que l’on exprime également en fonction du montant de provisions)
> predCL=function(reg=reg1,regb=reg1b){
+  p = 2*6-1;
+  phi.P = sum(residuals(reg,"pearson")^2)/(np-p)
+  Sig = vcov(reg)
+  X = model.matrix(regb)
+  Cov.eta = X%*%Sig%*%t(X)
+  mu.hat = exp(predict(reg,newdata=data.frame(lig,col)))*futur
+  pe2 = phi.P * sum(mu.hat) + t(mu.hat) %*% Cov.eta %*% mu.hat
+  cat("Total reserve =", sum(mu.hat), "prediction error =", sqrt(pe2),sqrt(pe2)/sum(mu.hat),"\n")
+ }

Avec nos trois modèles, Poisson, ODP et binomiale négative, on obtient,
> predCL(reg1,reg1b)
Total reserve = 18680856 prediction error = 896876.9 0.04801048
> predCL(reg2,reg2b)
Total reserve = 18680856 prediction error = 4736425 0.2535443
> predCL(reg3,reg3b)
Total reserve = 18085795 prediction error = 2058134 0.1137984
Moralité, oui la surdispersion a un impact, ainsi que la méthode que l’on retient pour ma modéliser....

vendredi 5 février 2010

from two to three...

A short post to give more details about the final remark in the course of Financial Econometrics, and more precisely the formula that can be found in the book of Philip Jorion,

Note that this formula can be found (perhaps written with slight changes) in several papers, e.g. the following sentence (on the http://www.bis.org/ website),

or the following formula, on documents from the Bank of England website,

I recently pulished (in French, here) a paper on the Value-at-Risk, including the following graph,

Usually, three times the average over 60 trading days is the larger component, but during the financial crisis, it turned out that the daily component was almost three times higher than the average value over the past the months (this fact was mention by Paul Embrechts in some conference in Paris on risk measures).
The interpreation of the multiplicative k coefficient (which is from 2 to 3 in some publications, or which exceeds 3 in others) has been proposed in a paper of Gerhard Stahl, entitled three cheers. The idea is to use the Bienaymé-Tchebychev inequality. For random variables with finite variance, then
http://perso.univ-rennes1.fr/arthur.charpentier/latex/BTC-01.png
Recall that this inequality is simply a corrolary of Markov’s inequality
http://perso.univ-rennes1.fr/arthur.charpentier/latex/BTC-02.png
or for any increasing function http://perso.univ-rennes1.fr/arthur.charpentier/latex/BTC-99.png
http://perso.univ-rennes1.fr/arthur.charpentier/latex/BTC-03.png
(taking function http://perso.univ-rennes1.fr/arthur.charpentier/latex/BTC-04.png, applied to http://perso.univ-rennes1.fr/arthur.charpentier/latex/BTC-05.png). This upper bound can be far away from the true probability, see e.g. the gaussian case below, i.e. if  http://perso.univ-rennes1.fr/arthur.charpentier/latex/BTC-06.png,
http://perso.univ-rennes1.fr/arthur.charpentier/latex/BTC-07.png
> z = seq(0,3,by=.01)
> P = 2*dnorm(k)
> U = 1/z^2
> plot(z,P,type="l",lwd=2,col="red",xlab="",ylab="")

The ratio between the two is given below,
> plot(z,U/P,type="l",lwd=2,col="purple",xlab="",ylab="",ylim=c(0,10))

Note that it is possible to interprete the axis values as probabilities values, taking quantiles of the gaussian distribution
> plot(pnorm(z),U/P,type="l",lwd=2,col="purple",xlab="",ylab="",ylim=c(0,10),xlim=c(.9,1))
> abline(h=3,lty=2)

The interpretation is that the upper bound is 3 times higher than the true probability in the Gaussian case when z is the quantile of the http://perso.univ-rennes1.fr/arthur.charpentier/latex/BTC-06.png distribution associated with probability level 99%.
Note that
  • if z is the 95% quantile of the \mathcal{N}(0,1) distribution, the ratio is 2 (1.92)
  • if z is the 99% quantile of the \mathcal{N}(0,1) distribution, the ratio is 3 (3.04)
  • if z is the 99.55% quantile of the \mathcal{N}(0,1) distribution, the ratio is almost 4 (3.88)
  • if z is the 99.75% quantile of the \mathcal{N}(0,1) distribution, the ratio is 5 (5.04)
A more formal explaination is to assume that X is symmetric, and then
http://perso.univ-rennes1.fr/arthur.charpentier/latex/BTC-09.png
Thus, if http://perso.univ-rennes1.fr/arthur.charpentier/latex/BTC-10.png, i.e. http://perso.univ-rennes1.fr/arthur.charpentier/latex/BTC-11.png, we have an upper bound for the  http://perso.univ-rennes1.fr/arthur.charpentier/latex/BTC-12.png Value-at-Risk,
http://perso.univ-rennes1.fr/arthur.charpentier/latex/BTC-20.png
where the upper bound is the upper bound for the http://perso.univ-rennes1.fr/arthur.charpentier/latex/BTC-12.png Value-at-Risk for any distribution with finite variance and centred.
If  http://perso.univ-rennes1.fr/arthur.charpentier/latex/BTC-31.png, then http://perso.univ-rennes1.fr/arthur.charpentier/latex/BTC-32.png, i.e. http://perso.univ-rennes1.fr/arthur.charpentier/latex/BTC-33.png. But since, http://perso.univ-rennes1.fr/arthur.charpentier/latex/BTC-33.png for a http://perso.univ-rennes1.fr/arthur.charpentier/latex/BTC-36.png distribution, then
http://perso.univ-rennes1.fr/arthur.charpentier/latex/BTC-21.png
and further
http://perso.univ-rennes1.fr/arthur.charpentier/latex/BTC-22.png

Quelques éléments de référence en actuariat...

Afin d’apporter des compléments au cours, je mets en ligne quelques liens (non exhaustifs, loin de là). Il existe des séminaires de la CAS dédiés à la tarification, comme ici ou . En particulier, on peut trouver des exemples illustrés d’utilisation de techniques d’économétrie nonlinéaire (ici ou ) ou de présentations des modèles GLM (ici ou voire encore ici). Sinon je peux aussi renvoyer sur le Practitioner’s Guide to Generalized Linear Models de Watson Wyatt (ici), ou la CAS qui a mis en ligne () un très joli ouvrage Basic Ratemaking. Sur le provisionnement, la CAS là aussi vient de publier un très joli ouvrage (ici), Estimating Unpaid Claims Using Basic Techniques.

mardi 2 février 2010

Les modèles en réassurance

Publication d’un papier sur les modèles en réassurance dans la revue Risques (ici), suite à un questionnement sur la pertinance des modèles classiques utilisés par les réassureurs. La question initiale était partie de la constatation - que l’on retrouve ici ou - sur l’utilisation (ou la mauvaise utilisation) de modèles sophistiques en finance de marché, en essayant d’expliquer que - d’un point de vue épistémiologique au moins - les modèles des réassureurs étaient plus robustes. En particulier, on notera que les modèles les plus anciens utilisés par les réassureurs (en particulier la loi de Pareto comme je l’avais évoqué ici) ont eu une légitimité pratique durant plusieurs décennies avant d’être justifiés par la théorie des valeurs extrêmes. Je ferais d’ailleurs bientôt un billet sur l’histoire des valeurs extrêmes en statistiques, en revenant en particulier sur les travaux de Gumbel ou de Fréchet.

Sinon le code utilisé dans le papier est en ligne ici, et la base
> sinpe = read.table("http://perso.univ-rennes1.fr/arthur.charpentier/sinpe.csv",header=TRUE,sep=";")
> head(sinpe)
      DSUR  MNTPE
1 19850206 240439
2 19851228 125674
3 19850504 488331
4 19851118 457347
5 19850220 990919
6 19851214 182939
> annee=as.numeric(substr(as.character(sinpe$DSUR),1,4))
> sinistres=sinpe$MNTPE[annee>1992]
> XS=sinistres/100000
On se limite ici aux sinistres survenus après 1992. Si l’on visualise ces montants de sinistres, on obtient
> datesur=as.Date(as.character(sinpe$DSUR),"%Y%m%d")
> jour=datesur[annee>1992]
> plot(jour,sinistres/100000,xlab="",ylab="Coût individuel",cex=.5,ylim=c(0,600))
> ded=50
> abline(h=ded)

On peut aussi faire le graphique Pareto-log-log,
> library(evir)
> n=length(X)
> plot(log(sort(X)),log((n:1)/(n+1)),
+ xlab="Coûts des sinistres (logarithme)",ylab="Fonction de survie (logarithme)",cex=.8)
> out <- gpd(XS, 15)
> XI=as.numeric(out$par.ests[1]); BETA=as.numeric(out$par.ests[2])
> x0=seq(2,8,.01)
> lines(x,-1/XI*(x-log(15)),col="red")

(l’ajustement de la loi de Pareto permettant de tracer la droite rouge) ou encore visualiser l’estimateur de Hill de l’indice de queue,
> hill(X)

Pour finir, le code suivant permet de calculer la prime pure (ou plus généralement une prime de Wang) soit de manière non-paramétrique (burning cost) ou bien en utilisant l’ajustement d’une loi de Pareto. Le papier étant un papier de vulgarisation, j’ai pris le seuil de manière arbitraire, sans aucune recherche de valeur "optimale",
> DEDUC = seq(10,50,by=5)
> lambda=0;  
seuil=5
> WG1=WG2=rep(NA,length(DEDUC))
> for(k in 1:length(DEDUC)){
+ deductible=DEDUC[k]
+ out <- gpd(XS, seuil)
+ XI=as.numeric(out$par.ests[1]); BETA=as.numeric(out$par.ests[2])
+ G0=function(x){1-pgpd(x+seuil, xi = XI, mu = seuil, beta = BETA)}
+ G=function(x){(G0(x+deductible-seuil))/(G0(deductible-seuil))}
+ F=function(x){pnorm(qnorm(G(x))+lambda)}
+ (wang1=integrate(F, 0, Inf))
+ X=XS[XS>deductible]
+ n=length(X)
+ FS= function(z){
+ m=rep(NA,length(z))
+ for(i in 1:length(m)){
+ m[i]=sum(X>z[i]+deductible)/n}
+ return(m)
+ }
+ G=function(x){pnorm(qnorm(FS(x))+lambda)}
+ (wang2=sum(G(seq(0,800,.01))*.01))
+ WG1[k]=as.numeric(wang1$value)
+ WG2[k]=wang2
+ }
> plot(DEDUC,WG2-DEDUC,type=’b’,xlab="Niveau de la priorité (’00 000 euros)",ylab="Prime pure par sinistres réassuré",ylim=c(0,50))
> lines(DEDUC,WG1-DEDUC,type=’b’,col="red",pch=4)
> legend(10,50,c("Aujustement d’une loi de Pareto","’Burning cost’"),
+ col=c("red","black"),lty=1,cex=.8,pch=c(4,1))

Econométrie de la finance de marché (1)

Le cours d’économétrie de la finance commence aujourd"hui. Les slides de cette semaine, et de la semaine prochaine, sont en ligne ici et . Le cours se terminera normalement en salle machine, pour une mise en oeuvre pratique sous R, mais je posterais un billet plus complet d’ici là.

lundi 1 février 2010

Calculer une provision mathématique (2)

Après un petit clin d’œil à Michel Rabagliati qui a gagné le prix du public hier à Angoulême.... Je vais poursuivre le billet commencé la semaine dernière (ici), histoire de montrer que ce n’est pas si compliqué que ça de calculer des PM avec R. Compte tenu des doubles indices dans toutes les fonctions, nous allons créer des matrices afin de stocker les différentes grandeurs.
> L=read.table("http://perso.univ-rennes1.fr/arthur.charpentier/TD8890.csv",
+ sep=";",header=TRUE)
> Lx=L$Lx
> m=length(Lx)
> p=matrix(0,m,m); d=p
> for(i in 1:m){
+ p[1:(m-i+1),i]=Lx[1+i:m]/Lx[i]
+ d[1:(m-i+1),i]=(Lx[i:(m)]-Lx[i:(m)+1])/Lx[i]}
> diag(d[m:1,])=0
> diag(p[m:1,])=0
> q=1-p
aussi, p[10,40] correspondra à http://perso.univ-rennes1.fr/arthur.charpentier/latex/PM80.png. Je stocke ainsi la probabilité de décéder au bout de k années, ou bien exactement dansk années. On peut alors calculer les différentes grandeurs,
> a=E=A=matrix(0,m,m)
> r=3/100
> for(i in 1:(m-1)){
+ a[,i]=cumsum(1/(1+r)^(0:(m-1))*c(1,p[1:(m-1),i]))
+ }
> for(i in 1:(m-1)){
+ #for(j in 1:(m-i)){
+ #A[j,i]=sum(1/(1+r)^(1:j)*p[1:j,i]*q[1,i-1+1:j])
+ A[,i]=cumsum(1/(1+r)^(1:m)*d[,i])
+ }
> for(i in 1:m){
+ E[,i]=(1/(1+r)^(1:m))*p[,i]
+ }

A partir de ces constantes (ou de ces matrices de constantes), on peut déjà calculer la prime annuelle du contrat décès que nous avions considéré,
> x=50; n=30
> prime=A[n,x]/a[n,x]
> sum(prime/(1+r)^(0:(n-1))*c(1,p[1:(n-1),x]))
[1] 0.3116454
> sum(1/(1+r)^(1:n)*d[1:n,x])
[1] 0.3116454

De plus, on peut calculer les provisions mathématiques à l’aide des différentes méthodes (et vérifier qu’elles donnent bien le même résultat).
> ### méthode rétrospective
> VR=(prime*a[1:n,x]-A[1:n,x])/E[1:n,x]
> plot(0:n,c(0,VR),col="red")
> ### méthode prospective
> VP=diag(A[n-(0:(n-1)),x+(0:(n-1))])-prime*diag(a[n-(0:(n-1)),x+(0:(n-1))])
> points(0:n,c(VP,0),pch=4,col="blue")
> # méthode itérative
> VI=0
> for(k in 1:n){
+ VI=c(VI,(VI[k]+prime-A[1,x+k-1])/E[1,x+k-1])
+ }
> points(0:n,VI,pch=5,col="green")
Pour un contrat signé par une personne de 50 ans, pour une couverture de 30 ans, au taux d’actualisation de 3%, on obtient les provisions mathématiques suivantes,

(les trois vecteurs coïncident, presque par miracle). Si la personne avait 45 ans, les provisions auraient l’allure suivante (en bleu),

Enfin, si on regarde la sensibilité au taux d’actualisation, un taux de 5% donnerait les provisions suivantes (en vert),

Ben voilà. La prochaine étape ça sera de distribuer les devoirs maison....

vendredi 29 janvier 2010

Calculer une provision mathématique (1)

Bon, je vais essayer de clarifier ici une partie du cours que, tous les ans, je rate: les méthodes de calculs de la provision mathématique. Après quelques éléments de droit, je pense que le plus simple est d’illustrer le calcul de la PM sur un cas simple, comme en cours, i.e. une temporaire décès avec prime annuelle. Cet exemple avait été repris ici, très (trop) brièvement. Sinon je peux renvoyer vers les slides de Pierre Devolder, ici,

  • Valorisation d’une assurance temporaire décès
Le "principe fondamental de valorisation", on doit avoir
http://perso.univ-rennes1.fr/arthur.charpentier/latex/PM01.png
en faisant une valorisation à la date 0, i.e. la date de souscription du contrat.
Pour l’assuré, il souhaite payer une prime annuelle constante http://perso.univ-rennes1.fr/arthur.charpentier/latex/PM02.png, noté plus simplement http://perso.univ-rennes1.fr/arthur.charpentier/latex/PM03.png, tant qu’il est en vie i.e.
http://perso.univ-rennes1.fr/arthur.charpentier/latex/PM04.png

http://perso.univ-rennes1.fr/arthur.charpentier/latex/PM05.png
(le http://perso.univ-rennes1.fr/arthur.charpentier/latex/PM06.png car le paiement se faisant ici en début de période).  De même,
http://perso.univ-rennes1.fr/arthur.charpentier/latex/PM07.png

http://perso.univ-rennes1.fr/arthur.charpentier/latex/PM08.png
http://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/.losser-france_s.jpg (l’indemnité étant versé par l’assureur à terme échu). J’utilie ici les notations officielles selons les normes françaises que l’on retrouve dans les ouvrages de Pierre Petauton et Christian Hess (ormis peut être pour le taux d’actualisation qui est généralement noté i). Il existe des nuances entre les notations françaises et les notations internationales. Je renvoie ici ou pour les standards internationaux. Parmi les notations internationales, on notera que ce que j’ai noté A ici est souvent noté \bar{A} dans la littérateure anglosaxonne (le A signifiant un paiement en milieu d’année).
Berf, pour revenir à notre prime annuelle, on en déduit que
http://perso.univ-rennes1.fr/arthur.charpentier/latex/PM09.png
  • Définition de la provision mathématique
Notons http://perso.univ-rennes1.fr/arthur.charpentier/latex/PM10.png la valeur actuelle probable, en http://perso.univ-rennes1.fr/arthur.charpentier/latex/PM11.png, des engagements de l’assurés pour la période http://perso.univ-rennes1.fr/arthur.charpentier/latex/PM12.png. Aussi, http://perso.univ-rennes1.fr/arthur.charpentier/latex/PM13.png sera la valeur actuelle probable, en 0, des k premières primes annuelles. Et on notera http://perso.univ-rennes1.fr/arthur.charpentier/latex/PM14.png la valeur actuelle probable, en 0, des engagements de l’assurés pour la période http://perso.univ-rennes1.fr/arthur.charpentier/latex/PM15.png, i.e. la valeur actuelle probable des n-k dernières primes annuelles.
De manière analogue, notons http://perso.univ-rennes1.fr/arthur.charpentier/latex/PM16.png la valeur actuelle probable, en http://perso.univ-rennes1.fr/arthur.charpentier/latex/PM11.png, des engagements de l’assureur pour la période http://perso.univ-rennes1.fr/arthur.charpentier/latex/PM12.png.
Comme tenu du principe fondamental de valorisation, pour un contrat arrivant à échéance au bout de n années, on doit avoir
http://perso.univ-rennes1.fr/arthur.charpentier/latex/PM20.png
pour un contrat soucrit à la date 0 et tel qu’il n’y a plus d’engagement de part et d’autre passé n années.
Aussi
http://perso.univ-rennes1.fr/arthur.charpentier/latex/PM21.png
avec, de manière générale
http://perso.univ-rennes1.fr/arthur.charpentier/latex/PM22.png
et
http://perso.univ-rennes1.fr/arthur.charpentier/latex/PM23.png
(d’où le principe d’inversion du cycle de production de l’assurance).
La provision mathématique (pures) de l’année k sera noté http://perso.univ-rennes1.fr/arthur.charpentier/latex/PM30.png si elles sont actualisée à la date t. La référence étant L'image “http://perso.univ-rennes1.fr/arthur.charpentier/latex/PM31.png” ne peut être affichée car elle contient des erreurs. (i.e. on actualise en k). C’est ce qui a été noté http://perso.univ-rennes1.fr/arthur.charpentier/latex/PM32.png dans le cours. On  définit http://perso.univ-rennes1.fr/arthur.charpentier/latex/PM33.png par
http://perso.univ-rennes1.fr/arthur.charpentier/latex/PM34.png
Cette définition sera dite rétrospective (car on se place sur la période antérieure à k).On peut aussi écrire, de manière équivalente (compte tenu du principe de valorisation)
http://perso.univ-rennes1.fr/arthur.charpentier/latex/MP000.png
Cette définition sera dite prospective (car on se place sur la période postérieure à k).
Enfin, il existe une dernière méthode, correspondant à une simple mise à jour, i.e.
http://perso.univ-rennes1.fr/arthur.charpentier/latex/PM36.png
Cette méthode sera dite itérative,voire en l’occurence itérative ascendante, car on initialise avec http://perso.univ-rennes1.fr/arthur.charpentier/latex/PM37.png. Mais il sera aussi possible de construire une méthode itérative descendante, récursive, commençant à la fin du contrat. Mais n’essayons pas de compliquer les choses (j’ai déjà suffisement de mal à m’y retrouver moi même).
  • Un peu de législation
Comme l’explique l’article R331-3 du Code des Assurance (ici), "Provision mathématique : différence entre les valeurs actuelles des engagements respectivement pris par l’assureur et par les assurés, à l’exception, pour les contrats mentionnés à l’article L. 142-1, des engagements relatifs à la provision de diversification" ou ici pour le lexique de l’ACAM. Le document ici reprend aussi une définition de ces provisions,
http://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/2008_11_30_Le_Monde_Gerner-1.jpg
  • La vision du Bowers et al.
Dans l’actuariat américain, le calcul des benefit reserves se fait en introduisant http://perso.univ-rennes1.fr/arthur.charpentier/latex/PM40.png, correspondant à la valeur actualisée à la date t des gains (ou des pertes) futurs de l’assureur (obtenus comme différence entre les engagements de part et d’autre). La provision mathématique est alors "l’espérance de la valeur actualisée à la date k des gains futurs de l’assureur conditionnelemnt au fait que (x) soit en vie". La forme concise est alors
http://perso.univ-rennes1.fr/arthur.charpentier/latex/PM41.png
C’est cette définition que nous avions utilisé en cours.
  • La méthode prospective
Nous avions commencé, en cours, par la méthode prospective, car c’est la plus naturelle (de mon point de vue) compte tenu de la législation. Nous avions posé
http://perso.univ-rennes1.fr/arthur.charpentier/latex/MP000.png
Notons que
http://perso.univ-rennes1.fr/arthur.charpentier/latex/PM51.png
où le terme de droite désigne la valeur actuelle probable d’un capital différé, relatif au versement d’un euro dans n année, conditionnée par la survie de (x), i.e.
http://perso.univ-rennes1.fr/arthur.charpentier/latex/PM52.png
Si l’on se place à la date k (car c’est le plus simple, mais l’assuré a alors l’âge x+k), notons que la différence entre les valeurs actuelles probables des engagements des deux parties donne, simplement
http://perso.univ-rennes1.fr/arthur.charpentier/latex/PM53.png
car d’un côté, on a une temporaire décès sur les n-k années restantes pour un assuré d’âge x+k, et de l’autre, l’assuré a pris l’engagement de verser sa prime (qui reste inchangée) pendant n-k années s’il vit. Aussi,
http://perso.univ-rennes1.fr/arthur.charpentier/latex/PM54.png
où l’on considère des assurances décès différées. On peut aussi écrire
http://perso.univ-rennes1.fr/arthur.charpentier/latex/PM55.png
  • La méthode rétrospective
La méthode rétrospective n’est pas beaucoup plus compliquée. On écrit simplement
http://perso.univ-rennes1.fr/arthur.charpentier/latex/PM60.png
i.e. http://perso.univ-rennes1.fr/arthur.charpentier/latex/PM61.png. Or http://perso.univ-rennes1.fr/arthur.charpentier/latex/PM62.png, et donc
http://perso.univ-rennes1.fr/arthur.charpentier/latex/PM63.png
  • La méthode itérative
L’idée est ici de décrire la variation de la provision mathématique entre deux dates en fonction des variation des engagements de part et d’autre. D’un côté il y a le paiement de la prime (en début de période, donc pas de problème d’actualisation et de non paiement), et de l’autre, une assurance décès sur un an. Aussi
http://perso.univ-rennes1.fr/arthur.charpentier/latex/PM70.png
Or http://perso.univ-rennes1.fr/arthur.charpentier/latex/PM72.png, ce qui donne, finalement
http://perso.univ-rennes1.fr/arthur.charpentier/latex/PM73.png
avec la convention que la première PM est nulle (de part notre principe fondamental de valorisation).
  • Mise en oeuvre pratique
Bon, c’est bien joli ces formules, mais il serait bon de voir si on peut faire des calculs. J’avais fait en cours les calculs dans une feuille excel (comme ici), mais avec R on devrait pouvoir y arriver facilement.... ça sera l’objet d’un billet qui arrivera rapidement (normalement).

Proposition de TER sur la vitesse des véhicules

J’avais fait, il y a quelque temps, un petit papier sur la régulation du marché de l’assurance (ici). Je rappelais ce que j’avais retenu de mes études, correspondant à la notion de "menace crédible": généralement la punition n’est pas nécessaire si la menace est suffisement forte. Dans le même ordre idée, je me suis longtemps demandé à quoi servent les radars ?
 http://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/.radar-automatique-fusille_m.jpg
Je vais essayer d’adopter ici un oeil d’économiste.... Si l’on regarde 2 minutes les radars classiques, la grande question que je me posais à leur sujet était  "faut-il faire des appels de phare lorsque l’on aperçoit des policiers avec un radar ? ". On va supposer que les agents ont une utilité qui soit fonction du "bien être collectif", autrement dit je suppose que les conducteurs souhaitent que le nombre d’accidents baisse, et donc ont intérêt à ce que les autres conducteurs baissent leur vitesse1. Comparons différents cas
  • si la personne arrivant en face dépasse la vitesse autorisée et que je ne fais pas d’appel de phare: la personne se fera arrêter, paiera une amende (et éventuellement roulera lentement pendant quelques dizaines de kilomètres, tout en pestant en son fort intérieur2)
  • si la personne arrivant en face dépasse la vitesse autorisée et que je fais un appel de phare: la personne en face a généralement un réflexe que l’on retrouve chez tous les conducteurs: elle ralentit. En fait, la personne avertie va souvent respecter la limitation de vitesse entre le moment où elle a reçu le premier appel de phare, et le moment où elle aperçoit la maréchaussée (voire davantage3). Bref, pendant quelques dizaines de kilomètres, elle respecte la limitation de vitesse. En fait l’idéal est même que la police soit cachée derrière un pont ou un bosquet,  et que la personne ne les voit pas.... elle peut  alors respecter  la limitation de vitesse pendant au moins 50 kilomètres (certes, en pestant sur la fin en se demandant "mais pourquoi donc ai-je reçu un appel de phare ?", mais peu importe).
Bref, il semble que l’on puisse arrive à la conclusion que pour aider la sécurité routière, faire des appels de phare n’est pas stupide. Bon, en fait la stratégie optimale est plutôt qu’il faut faire des appels de phare s’il n’y a pas de policiers...
Sinon il existe aussi les radars automatiques. Je passe sur le fait que certains prétendent que les radars automatiques ont été une très bonne chose pour la sécurité routière (je reviendrais dans un billet, bientôt, sur la différence entre de la corrélation et de la causalité). Je pensais que les radars automatiques pourraient être une mine d’or pour faire des stats... sauf que les radars n’enregistrent pas toutes les vitesses mesurées.... En tout les cas c’est ce qui m’a constament été dit. Malgré tout, il est possible de récupérer des données très intéresssantes sur les vitesses de véhicules et l’INRETS récupère ces données... Si un étudiant souhaite faire un TER (travail encadré de recherche) c’est possible ! J’ai les données !

1 aucun agent ne pense que si personne ne se fait contrôler, les impôts devront augmenter l’année suivante afin de combler le manque à gagner pour l’état. On exclura aussi les personnes sadiques qui aiment à penser que les autres conduisent plus mal qu’eux et que donc les autres doivent le payer à un moment donné.
2 ou au contraire, décidera d’accéler car elle a perdu pas mal de temps à discuter avec l’agent et se retrouve en retard, ou bien elle pense qu’elle est tombé sur le radar qui avait été installé aujourd’hui sur la route, et que ce n’est pas possible qu’il y ait un autre contrôle 12km plus loin !
3 car juste après avoir doublé les policiers, elle s’enflamme en disant qu’elle a eu de la chance, et qu’elle est heureuse de ne pas avoir payé d’amende, et que ça ne se joue à pas grand chose, et que finalement ce n’est pas si mal de rouler un peu moins vite, etc.

jeudi 28 janvier 2010

Optimal control, part 2

In the first part (here), we introduced Bellman’s idea of backward induction. But what if we consider now infinite time horizon ? Actually, the maths will be even more simple... and we will be able to use fixed pointed theorem to derive solutions.

  • The mathematical framework
Here, consider the following value function
http://perso.univ-rennes1.fr/arthur.charpentier/latex/opt-contr-01.png
and define 
http://perso.univ-rennes1.fr/arthur.charpentier/latex/opt-contr-02.png
A sequence http://perso.univ-rennes1.fr/arthur.charpentier/latex/opt-contr-04.png is said to be an admissible solution for starting point x,
http://perso.univ-rennes1.fr/arthur.charpentier/latex/opt-contr-05.png
if http://perso.univ-rennes1.fr/arthur.charpentier/latex/opt-contr-06.png and http://perso.univ-rennes1.fr/arthur.charpentier/latex/opt-contr-07.png. If we reformulate the dynamic programming idea, we obtain that if http://perso.univ-rennes1.fr/arthur.charpentier/latex/opt-contr-04.png is a solution to problem http://perso.univ-rennes1.fr/arthur.charpentier/latex/opt-contr-08.png, then for all http://perso.univ-rennes1.fr/arthur.charpentier/latex/opt-contr-09.png, sequence http://perso.univ-rennes1.fr/arthur.charpentier/latex/opt-contr-10.png is a solution to problem http://perso.univ-rennes1.fr/arthur.charpentier/latex/opt-contr-11.png. It comes that function v is a solution of Bellman’s equation
http://perso.univ-rennes1.fr/arthur.charpentier/latex/opt-contr-12.png
Note that is can be ssen as some fixed point resul, since
http://perso.univ-rennes1.fr/arthur.charpentier/latex/opt-contr-13.png
i.e.
http://perso.univ-rennes1.fr/arthur.charpentier/latex/opt-contr-14.png
So far, it shouldn’t be so hard....
  • Frank Ramsey’s model (discrete version)
In 1928, Frank Ramsey wanted to understand the amount of savings in a dynamic perspective (in how much of its income should a nation save). Consider the following infinite horizon problem, where some planifier wants to maximize
http://perso.univ-rennes1.fr/arthur.charpentier/latex/co-opt-01.png
subject to constraints http://perso.univ-rennes1.fr/arthur.charpentier/latex/co-opt-02.png http://perso.univ-rennes1.fr/arthur.charpentier/latex/co-opt-03.png and http://perso.univ-rennes1.fr/arthur.charpentier/latex/co-opt-04.png.
Before looking at dynamic programming answers, we might start with standard Lagrangian optimization techniques.
Assuming concavity of utility function and production function, we should look only for interior solutions. Define the Lagrangian as
http://perso.univ-rennes1.fr/arthur.charpentier/latex/co-opt-05.png
Thus, the first order conditions are then given by
http://perso.univ-rennes1.fr/arthur.charpentier/latex/co-opt-06.png
and
http://perso.univ-rennes1.fr/arthur.charpentier/latex/co-opt-088.png
Assume further some terminal condition, e.g.
http://perso.univ-rennes1.fr/arthur.charpentier/latex/co-opt-08.png
(also called transversality condition). If we combine those two conditions, and assume that the first constraint is saturated, we obtain the so-called Euler equation,
http://perso.univ-rennes1.fr/arthur.charpentier/latex/co-opt-09.png
It is also possible to use Bellman’s equation: given the dynamic of the capital
http://perso.univ-rennes1.fr/arthur.charpentier/latex/co-opt-10.png
http://perso.univ-rennes1.fr/arthur.charpentier/latex/co-opt-11.png
The first order condition states
http://perso.univ-rennes1.fr/arthur.charpentier/latex/co-opt-12.png
But since v is unknown, so is its derivative. But from the enveloppe theroem, we obtain something like
http://perso.univ-rennes1.fr/arthur.charpentier/latex/co-opt-13.png
where
http://perso.univ-rennes1.fr/arthur.charpentier/latex/co-opt-14.png
We can then write
http://perso.univ-rennes1.fr/arthur.charpentier/latex/co-opt-15.png
i.e.
http://perso.univ-rennes1.fr/arthur.charpentier/latex/co-opt-16.png
and finally
http://perso.univ-rennes1.fr/arthur.charpentier/latex/co-opt-17.png
which is, Euler’s equation.
  • A specified model, with calculations
As in the previous post (here), consider a log utility function, and a power production function,http://perso.univ-rennes1.fr/arthur.charpentier/latex/co-opt-18.png and http://perso.univ-rennes1.fr/arthur.charpentier/latex/co-opt-19.png. The dynamic is then
 http://perso.univ-rennes1.fr/arthur.charpentier/latex/co-opt-20.pngand http://perso.univ-rennes1.fr/arthur.charpentier/latex/co-opt-21.png
Note that fixed points are here
http://perso.univ-rennes1.fr/arthur.charpentier/latex/co-opt-22.png
and
http://perso.univ-rennes1.fr/arthur.charpentier/latex/co-opt-23.png
Recall that the value function is defined as
http://perso.univ-rennes1.fr/arthur.charpentier/latex/co-opt-24.png
A natural idea to derive the value function can be to iterate, i.e.
http://perso.univ-rennes1.fr/arthur.charpentier/latex/co-opt-25.png
starting with a simple function, e.g the null function, at step 0. At step n=1
http://perso.univ-rennes1.fr/arthur.charpentier/latex/co-opt-26.png
thus
http://perso.univ-rennes1.fr/arthur.charpentier/latex/co-opt-27.png
At step n=2,
http://perso.univ-rennes1.fr/arthur.charpentier/latex/co-opt-28.png
i.e.
http://perso.univ-rennes1.fr/arthur.charpentier/latex/co-opt-29.png
The first order condition is then
http://perso.univ-rennes1.fr/arthur.charpentier/latex/co-opt-30.png
and thus, we obtain
http://perso.univ-rennes1.fr/arthur.charpentier/latex/co-opt-31.png
that can be plugged in the previous equation, i.e.
http://perso.univ-rennes1.fr/arthur.charpentier/latex/co-opt-32.png
At step 3, we start from that new expression, derive the first order condition, and we get
http://perso.univ-rennes1.fr/arthur.charpentier/latex/co-opt-33.png
and
http://perso.univ-rennes1.fr/arthur.charpentier/latex/co-opt-34.png
and so on...
And finally, we can prove that
http://perso.univ-rennes1.fr/arthur.charpentier/latex/co-opt-35.png
i.e. http://perso.univ-rennes1.fr/arthur.charpentier/latex/co-opt-36.png. Assuming that
http://perso.univ-rennes1.fr/arthur.charpentier/latex/co-opt-37.png 
actually, we can prove that
http://perso.univ-rennes1.fr/arthur.charpentier/latex/co-opt-38.png
(and http://perso.univ-rennes1.fr/arthur.charpentier/latex/co-opt-39.png has a form that can be explicited).

vendredi 22 janvier 2010

Boire ou conduire ?... ou marcher ?

J’assistais en début de semaine à une présoutenance de thèse, où l’on apprenait (entre autres choses) que la consommation d’alcool augmentait fortement la probabilité d’avoir un accident de la route si on prendre le volant. Comme le disent la prévention routière et la FFSA (ici par exemple), le risque d’être impliqué dans un accident mortel est multiplié par 2 avec 0,5 g/l, par 10 à 0,8 g/l et 35 à 1,2 g/l (je n’ai pas eu accès à l’étude et aux données donc on prendra ces chiffres tels quels).
Ceci a fait remonter à la surface des débats interminables pendant les vacances de Noël  - où souvent ces débats n’ont lieu qu’entre personnes ayant déjà dépassé la dose limite pour prendre le volant (ce qui le rend d’autant plus intéressant (in vino veritas)): peut-on prendre le volant quand on a bu deux verres ? Sinon que doit-on faire ?
 Au même moment, je lisais Superfreakonomics de Steven Levitt et Stephen Dubner, qui proposait une question alternative (dans la première partie) "vaut-il mieux conduire ou rentrer à pied quand on trop bu ? ". Pour cela, ils utilisaient des statistiques que je résumerais de la manière suivante

  • les conducteurs (aux Etats-Unis) parcourent 5 milliards de miles, dont 1/140 sont conduits en état d’ivresse
  • les piétons (toujours aux Etats Unis) parcourent 43 milliards de miles, et Steven Levitt se dit que la proportion de miles parcourus en état d’ivresse n’a pas de raison d’être inférieure pour les piétons, soit 1/140. Soit un total de 307 millions de miles parcours en état d’ivresse.
Autrement dit, comme le conclue Steven Levitt, "on a per-miles basis, a drunk walker is eight times more likeliy to get killed than a drunk driver". Bon, le raisonnement est fallacieux, mais je trouvais la conclusion assez intéressante pour être creusée un peu (à partir de statistiques disponibles ici par exemple). Bon, comme toutes ces études, il est très dur de séparer les différents effets car les variables explicatives des accidents sont très fortement corrélées. Par exemple, si on se focalise sur les accidents de piétons (et que l’on regarde leur taux d’alcool) on note que la plupart des accidents où des piétons ont essentiellement lieux au milieu de la nuit

On peut d’ailleurs regarder le taux d’alcool de ces piétons victimes d’un accident,

Sur la figure ci-dessous, on voit pour les piétons à jeun (en vert) ou alcoolisés (en rouge) la répartition des moments des accidents,

ou de manière duale, la probabilité d’avoir un accident (piéton) saoul en fonction de l’heure de la journée (ou de la nuit),

Notons que l’on a des tableaux semblables pour 1999 (ici) par exemple.
Mais on peut aussi regarder si, par malheur, les gens ivres ne s’attirent pas entre eux ?

On peut faire un test du chi-deux pour voir si ces modalités sont corrélées, ou pas, sur les données observées sur 4 années, 2000 (ici), 1999 () et 1998 (),
> M01=t(matrix(c(2345,98,308,154,10,34,1094,87,260),3,3))
> M00=t(matrix(c(2323,80,276,149,10,32,1089,245,6),3,3))
> M99=t(matrix(c(2443,81,223,141,11,26,1211,72,238),3,3))
> M98=t(matrix(c(2583,92,284,189,11,33,1184,90,272),3,3))
> (M=M98+M99+M00+M01)
     [,1] [,2] [,3]
[1,] 9694  351 1091
[2,]  633   42  125
[3,] 4578  494  776
> chisq.test(M)$expected
          [,1]      [,2]       [,3]
[1,] 9333.2254 555.42240 1247.35223
[2,]  670.4903  39.90103   89.60864
[3,] 4901.2843 291.67656  655.03914
> chisq.test(M)
        Pearson’s Chi-squared test
data:  M
X-squared = 308.9695, df = 4, p-value < 2.2e-16

ce qui tend à montrer que les piétons saouls ont peut être la malchance de tomber également sur des chauffeurs saouls... Bref, si je reprends l’analyse de Lewitt, quand on a bu, il ne faut pas prendre le volant, et encore moins rentrer à pieds chez moi... Promis, la prochaine, je reste squatter chez le prochain qui m’offre un verre ! ou bien je vais arrêter de boire...

Présentation de l'économétrie en Licence

Je mets en ligne les slides de la courte intervention que j’avais faite devant les étudiants de Licence pour présenter le master de Statistique et d’Econométrie (ou plutôt de présenter l’intérêt de la modélisation statistique et économétrique), ici. Le desciptif de la forme du master risque de changer à moyen terme, mais les illustrations permettent de se faire une idée de ce qu’est la modélisation. L’essentiel des illustrations sont extraites de ce blog, donc je n’insisterais pas trop.

jeudi 21 janvier 2010

Modéliser des durées de vies humaines (1)

Dans le TD de ce matin, nous avons commencé par le modèle d’Abraham De Moivre, avant de présenter celui de Benjamin Gompertz. Je vais continuer à illustrer historiquement le cours d’assurance-vie, pour revenir sur les grands modèles paramétriques classiques...

  • Le modèle de De Moivre
Le modèle d’Abraham De Moivre peut paraître trop simple, mais il n’en est pas moins relativement riche, et intéressant. L’idée de De Moivre était que la durée de vie résiduelle pour un individu d’âge x était uniformément distribuée entre 0 et w-x, w étant la durée de vie limite d’un individu. Le taux de hasard est alors une fonction croissante de l’âge, ce qui est une hypothèse relativement générale (en particulier passé 10 ans). Rappelons qu’en 1724, c’était un blasphème de penser qu’il existait des lois (autres que divines) pour modéliser la vie humaine.
De Moivre a pu acquérir une légitimité sur la modélisation des durées de vies humaines après avoir quelque temps en Angleterre, à côtoyer en particulier Edmond Halley. L’autre légitimité a été acquise sur son lit de mort, en 1754, à Londres. La légende raconte que De Moivre avait remarqué qu’il dormait, chaque nuit, 15 minutes de plus. En faisant un petit calcul sur les suites arithmétique, il avait prédit qu’il mourrait le jour où il dormirait 24 heures. Et la légende prétend qu’il ne s’est pas trompé (ici).
  • Le modèle de Gompertz
Un siècle plus tard, Benjamin Gompertz proposait un autre modèle pour modéliser la durée de vie humaine. Si on retient souvent la forme simple sur le taux de hasard, notons qu’il est possible de caractériser cette loi en notant que
Il est alors possible d’ajuster cette loi de manière assez simple, en demandant à ce qu’elle passe par 3 points particuliers (comme il y a 3 paramètres). En cours, nous avions mis des contraintes sur les âges 50, 60 et 70. Le code est le suivant (pour une ajustement de la table TD88-90),
> L=read.table("http://perso.univ-rennes1.fr/arthur.charpentier/TD8890.csv",
+ sep=";",header=TRUE)
> h=10;  x1=50
> x2=x1+h;  x3=x2+h
> c10=log(L$Lx[L$Age==x2]/L$Lx[L$Age==x3])/log(L$Lx[L$Age==x1]/L$Lx[L$Age==x2])
> c=c10^(1/(x2-x1))
> g=exp(log(L$Lx[L$Age==x1]/L$Lx[L$Age==x2])/c^x1/(1-c^(x2-x1)))
> k=L$Lx[L$Age==x1]/(g^(c^x1))
> LG=k*g^(c^(L$Age))
> plot(L$Age,LG,type=’l’,col="blue",xlab="Age",ylab="Nombre de survivants Lx")
> lines(L$Age,L$Lx,col="red")
L’ajustement se visualise sur la figure ci-dessous,

avec l’erreur observée ci-dessous.

Mais on peut tenter un ajustement sur les âges 30, 50 et 70,

lundi 18 janvier 2010

Quand la formation continue se fourvoie

Le concept d’une formation professionnelle continue, i.e. "tout au long de la vie", remonte à Condorcet, qui partait de l’idée qu’"en continuant ainsi l’instruction pendant toute la durée de la vie, on empêchera les connaissances acquises dans les écoles de s’effacer trop promptement de la mémoire. On pourra montrer l’art de s’instruire par soi-même, comme à chercher des mots dans un dictionnaire, à se servir de la table d’un livre, à suivre sur une carte, sur un plan, sur un dessin, des narrations ou des descriptions, à faire des notes ou des extraits" (ici ou ). Mais pour être honnête, on retrouve des idées similaires chez Platon (ici par exemple). On retrouvera cette expression dans un certain nombre de références de la commission européenne en 2001, suivi d’une résolution du conseil européen en juin 2002. Bref, on y retrouve des idées que je ne peux que cautionner !
Mais forcément la vertue n’existe pas toujours, et c’est un peu ce que racontait Xavier Monnier il y a tout juste un mois dans Bakchich Hebdo, au sujet de la formation continue des avocats,

Il y a quelques jours, l’Institut des Actuaires a suivi le mouvement (même si l’idée était dans les têtes depuis quelques années), avec désormais l’obligation pour les actuaires de suivre l’équivalent de 3 jours de formation (ce qu’on peut rapprocher des 20 heures imposées aux avocats),   

Voilà quelques années que nous essayons avec Frédéric et Stéphane (ici) de monter des formations pour l’Institut.... Reste à espérer que l’effet pervers dénoncé dans l’article ne se retrouve pas chez les actuaires.

samedi 16 janvier 2010

Vraisemblance, quasi-vraisemblance, ou pseudo-vraisemblance ?

Un petit point technique sur les GLM (que j’ai passé rapidement en cours). Pour revenir un peu en arrière, il existe plusieurs méthodes usuelles d’estimation de paramètre. Supposons que l’on dispose d’un échantillon tiré suivant une loi http://perso.univ-rennes1.fr/arthur.charpentier/latex/vrais01.png
La méthode du maximum de vraisemblance consiste à maximiser la vraisemblance c’est à dire la loi jointe du n-uplet. C’est une vraie loi. Dans le cas discret
http://perso.univ-rennes1.fr/arthur.charpentier/latex/vrais02.png
soit, en supposant les observations indépendantes et de même loi
http://perso.univ-rennes1.fr/arthur.charpentier/latex/vrais03.png
Dans le cas continu,
http://perso.univ-rennes1.fr/arthur.charpentier/latex/vrais04.png
On suppose que l’on connaît ici parfaitement la famille de loi, et que connaître la valeur http://perso.univ-rennes1.fr/arthur.charpentier/latex/vrais0.png du paramètre suffit à caractériser complètement la loi. Le modèle est alors parfaitement spécifié. Parfois, on peut avoir du mal à spécifier complètement la loi (et en fait on n’en a pas vraiment besoin... On peut supposer que connaître le paramètre nous permet de caractériser les premiers moments (et pas complètement la loi). Ca sera l’idée de la pseudo-vraisemblance. Une pseudo vraisemblance adaptée à l’ordre 1 sera une fonction telle que l’espérance coïncide, i.e.
http://perso.univ-rennes1.fr/arthur.charpentier/latex/vrais05.png
Autrement dit la pseudo-vraisemblance n’est pas forcément la vraie loi des observations, mais elle nous garantie qu’ici, le premier moment est le bon. Sur la méthode du pseudo du pseudo-maximum de vraisemblance, je renvoie vers une très bon survey fait par Alain Trognon (ici) voilà une vingtaine d’année.
Une autre méthode est basée sur la quasi-vraisemblance. Ici, on ne fait pas toujours pas vraiment d’hypothèse de d’hypothèse de loi. Pour reprendre les notations des modèles GLM, on suppose juste que
http://perso.univ-rennes1.fr/arthur.charpentier/latex/vrais06.png
On suppose de plus qu’il existe un lien entre la variance et l’espérance, i.e.
http://perso.univ-rennes1.fr/arthur.charpentier/latex/vrais07.png
On introduit alors la loi de
http://perso.univ-rennes1.fr/arthur.charpentier/latex/vrais08.png
Alors
http://perso.univ-rennes1.fr/arthur.charpentier/latex/vrais09.png
optime une "quasi-vraisemblance", qui permet généralement d’obtenir un bon estimateur du paramètre (il est asymptotiquement sans biais et et optimal). Son avantage est de ne pas reposer sur la spécification d’une loi. Notons que l’on peut obtenir la variance de l’estimateur
http://perso.univ-rennes1.fr/arthur.charpentier/latex/vrais10.png
Pour aller un peu plus loin sur ces histoires de quasi-vraisemblance, on notera qu’il existe, sous R, des lois quasi (quasi poisson ou quasi binomiale) parmi les familles de lois exponentielles. Les lois de Poisson et binomiales sont des lois à un paramètre, c’est à dire que l’on force \phi à être égal à 1. Ce paramètre est simplement un modèle de surdispersion (dont on parlait tout à l’heure, ici). Mais si \phi n’est pas égal à 1, la loi exponentielle n’est plus nécessairement une loi de probabilité.
>  X=runif(100)
> Y=rpois(100,lambda=exp(1+X))
> summary(glm(Y~X,family=poisson))
Call:
glm(formula = Y ~ X, family = poisson)
Coefficients:
            Estimate Std. Error z value Pr(>|z|)  
(Intercept)   0.9263     0.1125   8.233  < 2e-16 ***
X             1.0987     0.1698   6.469 9.86e-11 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 138.25  on 99  degrees of freedom
Residual deviance:  94.62  on 98  degrees of freedom
AIC: 421.21
Number of Fisher Scoring iterations: 5

> summary(glm(Y~X,family=quasipoisson))
Call:
glm(formula = Y ~ X, family = quasipoisson)
Coefficients:
            Estimate Std. Error t value Pr(>|t|)  
(Intercept)   0.9263     0.1078   8.593 1.36e-13 ***
X             1.0987     0.1627   6.752 1.04e-09 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for quasipoisson family taken to be 0.917879)

    Null deviance: 138.25  on 99  degrees of freedom
Residual deviance:  94.62  on 98  degrees of freedom
AIC: NA
Number of Fisher Scoring iterations: 5

On note que dans ce dernier cas, le paramètre indique une sous-dispersion, avec une valeur de 0,9178.
Enfin, dans la terminologie, il existe aussi une vraisemblance empirique...  histoire d’embrouiller encore davantage. Mais je n’en parlerais pas....
En revanche je reviens très bientôt pour finir l’économétrie de l’assurance auto...

- page 1 de 37