Freakonometrics

Nous avions obtenu, ici et là, des modèles pour les lois des fréquence et pour les coûts individuels. Nous avions supposé

• suit une loi binomiale négative (ici)
• suit une loi mélange de 3 lois (exponentielle tronquée, Dirac et lognormale translatée) (là)

• Utilisation des simulations paramétriques

Simuler une loi binomiale négative peut se faire en se souvenant que la loi binomiale négative est une loi Poisson-Gamma (mélange Poisson avec un paramètre Gamma). On peut alors simplement simuler une loi Gamma, puis simuler une loi de Poisson.  Mais faisons plus simple, via la fonction
> Ns=rnbinom(nsimul,.3,mu=.12)
Avant de simuler la charge totale, il nous faut quelques lignes de codes pour générer une loi exponentielle tronquée,
> rexpt=function(n,t,m){
+ ze=rexp(2*n,m)
+ z=ze[ze<t]
+ while(length(z)<n){
+ ze=rexp(2*n,m)
+ z=c(z,ze[ze<t])
+ }
+ return(z[1:n])
+ }
On peut alors enfin simuler la charge totale simplement via
> Ns=rnbinom(nsimul,.3,mu=.12)
> Z =sample(c("bleu","vert","rouge"),sum(Ns)+1,replace=TRUE,prob=c(.25,.4,.35))
> X = (Z=="vert")*1170 +
+     (Z=="bleu")*(1170+rlnorm(sum(Ns)+1,7.07,1.55)) +
+     (Z=="rouge")*rexpt(sum(Ns)+1,1170,1/400)
On commence par regarder, par police, combien il y a de sinistres, puis on simule autant de sinistres qu'il en faut. Comme on a un mélange, on simule le paramètre du mélange, i.e. le type de sinistre.
Enfin, le vecteur de la charge totale est ici obtenu via le code suivant,
> S=rep(0,nsimul)
> s=0
> for(i in (1:nsimul)[Ns>0]){
+ S[i]=sum(X[(s+1):(s+Ns[i])])
+ s=s+Ns[i]
+ }
On commence par vérifier que les ordres de grandeurs coïncident avec ce que nous dit la théorie
> Nemp=read.table("D:\ -data\\base-nb.txt",header=TRUE)$x
> Xemp=read.table("D:\ -data\\base-cout.txt",header=TRUE)$x
> mean(Ns)
[1] 0.120989
> mean(Nemp)
[1] 0.1216735
> mean(X)
[1] 1831.396
> mean(Xemp)
[1] 2881.173
> mean(S)
[1] 221.5804
> mean(Ns)* mean(X)
[1] 221.5788
On notera une différence significative sur les coûts moyens: on rate ici les très gros couts de sinistres.
Si l'on suppose l'ajustement valide, on en déduit alors une estimation d'un quantile élevé simplement
> quantile(S,.99)
     99%
3984.903
ou encore la probabilité qu'une police coute plus de 15000 euros sur un an,
> sum(S>15000)/length(S)
[1]  0.00179

• Utilisation des simulations paramétriques et nonparamétriques

Comme tenu de la différence observée sur le coût moyen des sinistres individuels, on peut aussi tenter de bootrapper dans les charges observées,
> Sb=rep(0,nsimul)
> for(i in (1:nsimul)[Ns>0]){
+ Sb[i]=sum(sample(Xemp,size=Ns[i],replace=TRUE))
+ }
>  mean(Sb)
[1] 336.3616
>  mean(Nemp)*mean(Xemp)
[1] 350.5623

Cette fois, l'ordre de grandeur paraît plus acceptable (on retrouve la prime pure mentionnée ici). Et si on regarde les quantiles et les probabilités d'évnènement rare,
> quantile(Sb,.99)
     99%
3602.766
> sum(Sb>15000)/length(Sb)
[1] 0.001403
On peut visualiser ces deux distributions ci-dessous, avec en rouge la fonction de survie de variable simulée, suivant des lois paramétriques, et en bleu la fonction de survie lorsqu'on boostrappe les coûts des sinistres,

• Approximation Normal Power

A l'aide de l'algorithme précédant, on peut récuper les premiers moments de la charge totale,
> library(fBasics)
Le chargement a nécessité le package : timeDate
Le chargement a nécessité le package : timeSeries
> mean(S)
[1] 221.5804
> var(S)
[1] 5781070
> skewness(S)
[1] 125.4405
attr(,"method")
[1] "moment"
> kurtosis(S)
[1] 32134.53
attr(,"method")
[1] "excess"
La distribution de S est alors
>  m=mean(S)
>  s=sd(S)
>  (g=as.numeric(skewness(S)))
[1] 125.4405
>  mean((S-mean(S))^3)/sd(S)^3
[1] 125.4405
> x <- seq(0,30000,length=500)
> z=(x-m)/s
> y=1-pnorm(sqrt(9/g^2 + 6*z/g + 1) - 3/g) 
qui n'est pas vraiment réaliste...
On voit que pour une police sur 20, la charge sinistre dépasse les 100000 euros. Cela n'est pas réaliste ! La moyenne devrait alors forcément excéder 5000.

• Algorithme de Panjer
  

Histoire de conclure, on peut aussi regarder ce qu'aurait donné l'algorithme de Panjer. Commençons par discrétiser la distribution de la charge individuelle.
> F=function(x){
+ y=NA
+ if(x<1170){y=pexp(x,1/400)/pexp(1170,1/400)*.35}
+ if(x==1170){y=.75}
+ if(x>1170){y=.75+plnorm(x-1170,7.07,1.55)*.25}
+ return(y)
+ }
> x  = seq(0,100000,length=1000)
> Vf = rep(NA,length(x)-1)
> for(i in 1:(length(x)-1)){
+   Vf[i] = F(x[i+1])-F(x[i])
+ }
> sum(Vf)
[1] 0.9994685
Il faut faire attention lors de la censure, car la discrétisation sera la même pour la charge totale. Une fois discrétisée la loi de la charge individuelle, on travaille sur la loi de la fréquence. Avec les notations de la classe de Sundt, on a
> size=.3
> mu=.12
> a=1-size/(size+mu)
> b=(size-1)*a
> a
[1] 0.2857143
> b
[1] -0.2
L'algorithme de Panjer s'écrit alors
>  G=rep(NA,length(Vf))
>  G[1]=dnbinom(0,size=.3,mu=.12)
>  for(i in 1:(length(G)-1)){
+  G[i+1]=sum((a+b*(1:i)/i)*Vf[1:i]*G[i-(1:i)+1])
+  }
>  sum(G)
[1] 0.9999346
et l'on peut visualiser la fonction de survie avec le code suivant
> lines(x[2:1000],1-cumsum(G),lwd=4)
Fort logiquement, ces méthodes (simulations et Panjer) donnent des grandeurs très proches en terme de quantiles.