Petit article sur la modélisation du système bonus-malus français dans le prochain numéro de la revue Risques (le numéro 83), intitulé "le bonus-malus français a-t-il encore un avenir ?" (en ligne ici). L'article est inspirée d'un mini-cours qu'avait donné Jean Lemaire il y a une dizaines d'années, alors qu'il était invité par AON à Hong Kong sur "markov chains in insurance". En particulier, il avait insisté sur la vitesse de convergence des systèmes bonus-malus, que j'ai repris en partie ici, sur le cas français.
Au niveau du code, pour étudier cette vitesse de convergence (au moins graphiquement), on peut reprendre le code suivant. On se donne une distribution initiale (ainsi que les niveaux de prime), et une matrice de transition

>U=c(1,0,0); P=c(121.8,.9*121.8,.81*121.8)
>p=.1
> M=matrix(c(p,p,0,1-p,0,p,0,1-p,1-p),3,3)

Ensuite, on garde la prime moyenne collectée, année après année, et la répartition dans les 3 classes,

>  MP=t(U)%*%P
>  MU=U
>  for(i in 1:10){
+  V=U%*%M
+  U=V
+  MU=rbind(MU,U)
+  MP=c(MP,U%*%P)}

Le code suivant (il n'est pas forcément optimal, mais il tourne) permet alors de faire un graphique montrant l'évolution par fréquences, ainsi que la prime moyenne collectée,

>  ylim <- c(0,100)
>  plot(0:1,0:1, col="white", type="p", ylim=ylim,
+ axes=FALSE, ylab="",xlab="Années", xlim=c(0,10))
> COULEUR=heat.colors(length(P))
> MU2=cbind(rep(0,nrow(MU)),MU[,length(P):1])
> MU3=t(apply(MU2,1,cumsum))
> for(i in 0:9){
+  for(j in 1:length(P)){
+  polygon(c(i-.4,i-.4,i+.4,i+.4),
+  100*c(MU3[i+1,j],MU3[i+1,j+1],MU3[i+1,j+1],
+ MU3[i+1,j]),col=COULEUR[j])
+ }}
>  axis(1);
>  ylim <- c(90,130)
>  axis(2)
>  par(new=TRUE)
>  plot(0:1,0:1, col="white", type="p",
+ ylim=ylim,axes=FALSE, ylab="",
+  xlab="",xlim=c(0,9))
>  axis(4)
> points(0:8,MP[1:9],type="b",
+ lwd=4,col="blue")
> abline(h=100,col="blue",lty=2)
>  mtext("Proportion dans les classes", 2, line=2)
>  mtext("     Prime moyenne", 4, line=-1.5 )

On peut alors adapter le code à n'importe quel système de bonus malus (à classes), comme ceux décrit dans le livre de Jean Lemaire (dont j'avais parlé ici). Par exemple pour
http://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/.BM-hong-kong-B-9_m.jpg
il suffit d'entrer les paramètres (je décris un peu mieux la matrice de transition, pour que les choses soient claires, en particulier il est souvent plus simple de rentrer la matrice par ligne, en calculant les probabilités d'arriver dans telle ou telle classe, en fonction de là où on est)
> U=c(1,0,0,0,0,0); P=c(100,80,70,60,50,40)
> p=.1
> M=t(matrix(
+ c(p,1-p,0,0,0,0,
+   p,0,1-p,0,0,0,
+   p,0,0,1-p,0,0,
+   p,0,0,0,1-p,0,
+   0,0,p,0,0,1-p,
+   0,0,0,p,0,1-p),
+ length(P),length(P)))

Sinon côté références, le livre de Jean Lemaire est en partie dans google books, ici, et pour l'article d'Edouard Franckx, il est là. Pour aller plus loin, il y a le livre de Michel Denuit (et al.), là, et sinon l'article de Magali Kelle (paru dans le BFA) qui parle du système bonus-malus français est en ligne ici.