Tag Archives: Wüthrich

Incertitude(s) dans les triangles de paiements

Un petit billet aujourd’hui sur les méthodes de provisionnements, et sur les calculs d’incertitude à 1 an, ou à l’ultime, de montants de provisions pour sinistres restant à payer…. Ce sujet avait été évoqué ici (surtout dans les commentaires à l’époque).

  • Mack et incertitude à long terme 

Le mse dans le cadre du modèle de Mack est proposé dans la library(ChainLadder), via la fonction MackChainLadder(). Sur le triangle utilisé dans les exemples et formations en provisionnement, on obtient les mse estimés suivants

> source("https://perso.univ-rennes1.fr/arthur.charpentier/bases.R")
> library(ChainLadder)
> MackChainLadder(PAID)
MackChainLadder(Triangle = PAID)

  Latest Dev.To.Date Ultimate    IBNR Mack.S.E CV(IBNR)
1  4,456       1.000    4,456     0.0    0.000      NaN
2  4,730       0.995    4,752    22.4    0.639   0.0285
3  5,420       0.993    5,456    35.8    2.503   0.0699
4  6,020       0.989    6,086    66.1    5.046   0.0764
5  6,794       0.978    6,947   153.1   31.332   0.2047
6  5,217       0.708    7,367 2,149.7   68.449   0.0318

              Totals
Latest:    32,637.00
Ultimate:  35,063.99
IBNR:       2,426.99
Mack S.E.:     79.30
CV(IBNR):       0.03

Mais il s’agit fondamentalement d’une formule sur l’estimation ultime de la charge. Pour une année de survenance donnée, la formule dérivée par Thomas Mack ressemble à ça

(honteusement emprunté dans des slides de Mario Wüthrich, par paresse, ici), où pour rappel,

  • Muerz & Wüthrich et incertitude à court terme

Pour faire simple, il existe une formule fermée, un peu complexe, et une simplification, reposant sur l’idée que, si on note

alors, Michael Merz et Mario Wüthrich ont montré que

Le code sur le triangle mentionné ici donne le résultat suivant (ce code est parfois appelé dans divers documents MackMerzWuthrich.R)

> source("https://perso.univ-rennes1.fr/arthur.charpentier/merz-wuthrich-triangle.R")
> MSEP_Mack_MW(PAID,0)
   MSEP Mack MSEP observable approchée MSEP observable exacte
1  0.0000000                  0.000000               0.000000
2  0.6393379                  1.424131               1.315292
3  2.5025153                  2.543508               2.543508
4  5.0459004                  4.476698               4.476698
5 31.3319292                 30.915407              30.915407
6 68.4489667                 60.832875              60.832898
7 79.2954414                 72.574735              72.572700

La dernière ligne (notée ici 7 correspond au total – toutes années de survenance confondues). Ici, on regarde l’erreur à un an. Et Mario Wüthrich concluait – lors d’une jolie école d’été à Lyon il y a quelques années – en notant le point suivant

Sur le triangle étudié ici, on était à 91%. Mais supposer que l’on soit inférieur à 1 est un point qu’il me semble intéressant de discuter…

  • incertitude plus grande à court terme ?

Dans un billet datant un peu (ici), on laissait entendre qu’il était sûrement possible que l’incertitude à court terme excède l’incertitude à long terme.
Arnaud Lacoume a trouvé un vrai triangle dans lequel l’incertitude de court terme excède l’incertitude à long terme (détaillé ici dans une feuille xls). Pour information, il s’agit du risque GAV, garanties accidents de la vie, bref, rien de très exotique…

> Triangle = read.table("https://perso.univ-rennes1.fr/arthur.charpentier/GAV-triangle.csv",sep=";")

Si on utilise le code gentiment offert par Arnaud, on note qu’effectivement, l’incertitude à 1 an est grande, très grande…

> MSEP_Mack_MW(Triangle,0)
    MSEP Mack MSEP observable approchée MSEP observable exacte
1        0.00                       0.0                    0.0
2    12459.74                  129692.2               152605.9
3   209431.14                  214136.5               214419.6
4   258003.38                  198072.3               198773.0
5  3055297.40                 3048489.5              3065525.1
6 58429393.29                57056117.3             67375794.0
7 58669646.13                57301552.4             67586106.6

On note ici que le MSE exact à un an (car Merz et Wüthrich proposent aussi une formule approchée, basée sur une approximation qui contraint l’estimation à un an à être inférieure à l’estimation à horizon ultime) est – si on regarde toutes années de survenance confondues – 15% plus élevé que le MSE à ultime. Avec un résultat du même ordre pour la dernière année (ce n’est donc pas un problème d’agrégation). On peut jouer à simuler des triangles, via un processus de déroulé de paiements.

Regardons par exemple le code suivant pour simuler des triangles. L’idée est ici que l’on n’a pas d’incertude sur la valeur à l’ultime,

> M=matrix(NA,10,10)
> for(i in 1:10){
+ X=rexp(10)
+ Y=cumsum(X)/sum(X)
+ M[i,1:(11-i)]=Y[1:(11-i)]*100}
> Triangle = M

On tire des variables positives, car on veut un modèle à incréments positifs. Par exemple, une des simulations a donné le triangle suivant

> Triangle
      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
 [1,] 35.4 36.3 38.2 47.6 50.7 66.1 71.3 72.1 98.7   100
 [2,] 15.5 20.2 22.3 23.7 63.9 68.2 75.9 79.3 95.1    NA
 [3,]  4.1 30.3 32.3 33.3 49.4 53.1 62.8 77.1   NA    NA
 [4,]  1.1 19.4 38.9 41.9 48.0 66.5 68.3   NA   NA    NA
 [5,]  0.3  7.2 10.8 17.8 28.1 30.1   NA   NA   NA    NA
 [6,] 15.6 32.2 46.0 50.6 68.1   NA   NA   NA   NA    NA
 [7,]  4.0  9.9 10.6 15.8   NA   NA   NA   NA   NA    NA
 [8,]  6.1 24.0 29.4   NA   NA   NA   NA   NA   NA    NA
 [9,]  0.0  0.5   NA   NA   NA   NA   NA   NA   NA    NA
[10,]  2.3   NA   NA   NA   NA   NA   NA   NA   NA    NA

En faisant le calcul sur les mesures d’incertitude à la Merz & Würthrich, on obtient une très grande incertitude à court terme,

> MSEP_Mack_MW(Triangle,0)
   MSEP Mack MSEP observable approchée MSEP observable exacte
1   0.000000                  0.000000               0.000000
2   6.820602                 12.776729             118.050196
3  13.377536                 13.119557              13.126625
4  17.415678                 13.638880              13.651131
5  11.797752                  5.574916               5.580840
6  26.767233                 18.407916              18.429069
7  27.184076                 23.725348              23.773581
8  44.100071                 16.410374              16.444281
9   7.227724                  3.450343               3.460982
10 54.012270                 49.443633              49.593284
11 99.184167                 85.042393             145.008282

Ici je suis 50% plus élevé ! Bon, mais c’est beau de montrer pourquoi une idée partagée par beaucoup de monde n’est pas toujours vraie…mais il reste à comprendre plus en détails le pourquoi … à suivre donc !

En attendant, je mets ici le code utilisé pour quantifier cette incertitude à un

### MERZ & WUTHRICH

MWChainLadder <- function(Trian, bool = 0, warning = FALSE){    

##  variables
I <- nrow(Trian)
J <- ncol(Trian)
diag <- matrix(NA,I,1)
diag_inv <- matrix(NA,I,1)
S_I <- matrix(NA,1,J)
S_II <- matrix(NA,1,J)
Phi <- matrix(NA,I,1)
Psi <- matrix(NA,I,1) 
Delta <- matrix(NA,I,1)
Gamma <- matrix(NA,I,1)
Lambda <- matrix(NA,I,1)
Upsilon <- matrix(NA,I,1)
cov_obs <- matrix(0,I,J)
cov_reel <- matrix(0,I,J)
msep_obs <- matrix(NA,1,J+1)
msep_reel <- matrix(NA,1,J+1)
Delta[1] <- Phi[1] <- Psi[1] <- Upsilon[1] <- Lambda[1] <- 0

Mack<-MackChainLadder(Trian)
Mack$sigma[J-1] <- sqrt(min(Mack$sigma[J-2]^4/
Mack$sigma[J-3]^2,min(Mack$sigma[J-2]^2,Mack$sigma[J-3]^2)))
#plot(Mack)

for (i in 1:I){
    diag[i] = Trian[i,I-i+1]
    diag_inv[i] = Trian[I-i+1,i]
}

for (j in 1:J){
    S_I[j] <- sum(Trian[1:(I-j),j])
    S_II[j] <- sum(Trian[1:(I-j+1),j])
}
S_I[I] <- 0

Delta[2] <- Mack$sigma[I-1]^2/(S_I[I-1]*(Mack$f[I-1])^2)
Phi[2] <- 0
Psi[2] <- Mack$sigma[I-1]^2/(diag[2]*(Mack$f[I-1])^2)
Upsilon[2] <- Mack$sigma[I-1]^2/(S_II[I-1]*(Mack$f[I-1])^2)
Lambda[2] <- diag[2]*Mack$sigma[I-1]^2/((Mack$f[I-1]^2)*S_II[I-1]*S_I[I-1])
for (i in 3:I){
    Delta[i] <- Mack$sigma[I-i+1]^2/(S_I[I-i+1]*(Mack$f[I-i+1])^2) + sum(
            (diag_inv[(I-i+2):(J-1)]/S_II[(I-i+2):(J-1)])^2*
         Mack$sigma[(I-i+2):(J-1)]^2/(S_I[(I-i+2):(J-1)]*(Mack$f[(I-i+2):(J-1)])^2))        
    Phi[i] <- sum((diag_inv[(I-i+2):(J-1)]/S_II[(I-i+2):(J-1)])^2*
           Mack$sigma[(I-i+2):(J-1)]^2/(diag_inv[(I-i+2):(J-1)]*(Mack$f[(I-i+2):(J-1)])^2))    
    Psi[i] <- Mack$sigma[I-i+1]^2/(diag[i]*(Mack$f[I-i+1])^2)
    Upsilon[i] <- Phi[i] + Mack$sigma[I-i+1]^2/(S_II[I-i+1]*(Mack$f[I-i+1])^2)

Lambda[i] <- diag[i]*Mack$sigma[I-i+1]^2/((Mack$f[I-i+1]^2)*S_II[I-i+1]*S_I[I-i+1]) + sum(
            (diag_inv[(I-i+2):(J-1)]/S_II[(I-i+2):(J-1)])^2*
           Mack$sigma[(I-i+2):(J-1)]^2/(S_I[(I-i+2):(J-1)]*(Mack$f[(I-i+2):(J-1)])^2))
}
Gamma = Phi + Psi

#MSEP par acc. year
for (i in 1:I){
    msep_obs[i] = Mack$FullTriangle[i,J]^2 * (Gamma[i] + Delta[i])
    msep_reel[i] = Mack$FullTriangle[i,J]^2 * (Phi[i] + Delta[i])
}

#Covariance
for (i in 2:(I-1)){ 
        for (k in (i+1):I){
        cov_obs[i,k] <- Mack$FullTriangle[i,J]*Mack$FullTriangle[k,J]*(Upsilon[i] + Lambda[i])
        cov_reel[i,k] <- Mack$FullTriangle[i,J]*Mack$FullTriangle[k,J]*(Phi[i] + Lambda[i])    
    }
}

#MSEP aggregated
msep_obs[I+1] = sum(msep_obs[2:I]) + 2*sum(cov_obs)
msep_reel[I+1] = sum(msep_reel[2:I]) + 2*sum(cov_reel)

##################################################################################################################
## EXACTES
##################################################################################################################
#DÈclaration des variables
facteur <- matrix(1,I,1)
Phi_exact <- matrix(NA,I,1)
Psi_exact <- matrix(NA,I,1)
Gamma_exact <- matrix(NA,I,1)
Upsilon_exact <- matrix(NA,I,1)
cov_obs_exact <- matrix(0,I,J)
cov_reel_exact <- matrix(0,I,J)
msep_obs_exact <- matrix(NA,1,J+1)
msep_reel_exact <- matrix(NA,1,J+1)
Phi_exact[1] <- Upsilon_exact[1]<- Gamma_exact[1] <- 0

# i=2
Phi_exact[2] <- Psi_exact[2] <- 0
Gamma_exact[2] <- Mack$sigma[I-1]^2/(diag[I-1]*(Mack$f[I-1])^2)
Upsilon_exact[2] <- Mack$sigma[I-1]^2/(S_II[I-1]*(Mack$f[I-1])^2)

# loop for prior
for (i in 3:I){
    facteur[i] <-prod(1 + diag_inv[(I-i+2):(J-1)]*(Mack$sigma[(I-i+2):(J-1)]^2)/((S_II[(I-i+2):(J-1)]*
            Mack$f[(I-i+2):(J-1)])^2))
    Phi_exact[i] <- (1 + Mack$sigma[I-i+1]^2/(diag_inv[I-i+1]*(Mack$f[I-i+1])^2))* (facteur[i]-1)    
    Psi_exact[i] <- (1 + Mack$sigma[I-i+1]^2/(S_II[i]*(Mack$f[I-i+1])^2)) * Phi[i] /
                (1 + Mack$sigma[I-i+1]^2/(diag_inv[I-i+1]*(Mack$f[I-i+1])^2))
    Upsilon_exact[i] <- (1 + Mack$sigma[I-i+1]^2/(S_II[I-i+1]*(Mack$f[I-i+1])^2))* facteur[i]-1    
    Gamma_exact[i] <- (1 + Mack$sigma[I-i+1]^2/(diag_inv[I-i+1]*(Mack$f[I-i+1])^2))* facteur[i]-1    
}

#MSEP per accident year
for (i in 1:I){
    msep_obs_exact[i] = Mack$FullTriangle[i,J]^2 * (Gamma_exact[i] + Delta[i])
    msep_reel_exact[i] = Mack$FullTriangle[i,J]^2 * (Phi_exact[i] + Delta[i])
}

#Covariance
for (i in 2:(I-1)){ 
    for (k in (i+1):I){
        cov_obs_exact[i,k] <- Mack$FullTriangle[i,J]*Mack$FullTriangle[k,J]*(Upsilon_exact[i] + Lambda[i])
        cov_reel_exact[i,k] <- Mack$FullTriangle[i,J]*Mack$FullTriangle[k,J]*(Psi_exact[i] + Lambda[i])    
    }
}

#MSEP ag.
msep_obs_exact[I+1] <- sum(msep_obs_exact[2:I]) + 2*sum(cov_obs_exact)
msep_reel_exact[I+1] <- sum(msep_reel_exact[2:I]) + 2*sum(cov_reel_exact)

##################################################################################################################
##Outputs Fonction 
##################################################################################################################
msep_Mack <- array(0,c(1,I+1))
msep_Mack[1:I] <- Mack$Mack.S.E[,I]
msep_Mack[I+1] <- Mack$Total.Mack.S.E
Vari <- array(0,c(1,I+1))
for (i in 1:I){
    Vari[i] <- Mack$FullTriangle[i,J]^2 * Psi[i]
}
Vari[I+1] <- sum(Vari[1:I])

# complete vision
if(bool==1){    
    result <- cbind(t(msep_Mack), t(sqrt(Vari)), t(sqrt(msep_obs)), t(sqrt(msep_obs_exact)), t(sqrt(msep_reel)), t(sqrt(msep_reel_exact)))
    result <- as.data.frame(result)
    l <- list("Variance CDR", "MSEP Mack","MSEP obs. approx.","MSEP obs. exact", "MSEP real approx.","MSEP real exact")
    names(result) <- l
    return(result)
    }

#si l'on souhaite uniquement l'aspect solvabilité
if(bool==0){
    result <- cbind(t(msep_Mack), t(sqrt(msep_obs)), t(sqrt(msep_obs_exact)))
    result <- as.data.frame(result)
    l <- list("MSEP Mack","MSEP obs. approx","MSEP obs. exact")
    names(result) <- l
    return(result)
    }
}