Modélisation des montants de sinistres

Loi lognormale

Pour calibrer une loi lognormale, $X\sim \text{lognormal}(\mu, \sigma)$ de densité $$f_X(x) = \frac{1}{\sigma x\sqrt{2\pi}}\exp\left\{-\frac{(\ln(x) - \mu)^2}{2\sigma^2} \right\},\text{ }x> 0.$$ on infère les paramètres d’une loi normale sur les log données.

# Inference des paramètres de la loi lognormale
mu = mean(log(AutoBi$LOSS))
sigma = sd(log(AutoBi$LOSS))

# Quantile-quantile plot pour la loi lognormale
quantile_lnorm = data.frame(
   quantile_loss = sapply(seq(0.001,0.999,0.001), function(alpha) quantile(AutoBi$LOSS, alpha)),
   quantile_lnorm = sapply(seq(0.001,0.999,0.001), function(alpha) qlnorm(alpha, meanlog = mu, sdlog = sigma))
   )
ggplot(data = quantile_lnorm, mapping = aes(x  = quantile_lnorm , y = quantile_loss )) + geom_point()  + geom_abline(intercept = 0, slope = 1)

data.frame('mu'=mu, 'sigma'= sigma)

Loi de Weibull

Pour calibrer la loi de Weibull, de densité

$$f_X(x) =\left(\frac{\alpha}{\beta}\right)\left(\frac{x}{\beta}\right)^{\alpha-1}\exp\left\{-\left(\frac{x}{\beta}\right)^{\alpha}\right\},\text{ }x> 0.$$ on utilise la méthode du maximum de vraisemblance.

# Fonction qui donne l'oppposé de la log vraisemblance
log_like <- function(param){
  -sum(log(dweibull(AutoBi$LOSS, shape = param[1], scale = param[2])))
}
# La fonction nlm qui permet de minimiser la fonction en argument
mle_weib <- nlm(log_like, c(1/5, 5))$estimate

# QQ-plot pour la loi de Weibull
quantile_weibull = data.frame(
   quantile_loss = sapply(seq(0.001,0.999,0.001), function(alpha) quantile(AutoBi$LOSS, alpha)),
   quantile_weib = sapply(seq(0.001,0.999,0.001), function(alpha) qweibull(alpha, shape = mle_weib[1], scale = mle_weib[2]))
   )

ggplot(data = quantile_weibull, mapping = aes(x  = quantile_weib , y = quantile_loss )) + geom_point()  + geom_abline(intercept = 0, slope = 1)

data.frame('alpha' = mle_weib[1], 'beta' = mle_weib[2])

Loi Gamma

Pour calibrer la loi gamma $X\sim\text{Gamma}(\alpha,\beta)$ de densité

$$f_X(x) = \frac{e^{-x/\beta}x^{\alpha-1}}{\beta^\alpha\Gamma(\alpha)},\text{ }x>0,$$ on peut utiliser la méthode des moments. On résout le système d’équation

$$\begin{cases} \mathbb{E}(X) = \alpha\beta&\\ \mathbb{V}(X) = \alpha \beta^2&\\ \end{cases} \Leftrightarrow \widehat{\alpha} = \frac{\widehat{\sigma}^2}{\overline{X}}\text{ et }\widehat{\beta} =\frac{\overline{X}^2}{\widehat{\sigma}^2},$$

où

$$\overline{X} = \frac{1}{n}\sum_{i = 1}^{n}X_i\text{ et }\widehat{\sigma}^2 = \frac{1}{n-1}\sum_{i = 1}^{n}(X_i-\overline{X})^2.$$

# Inférence des paramètre de la loi gamma
alpha_g <- var(AutoBi$LOSS) / mean(AutoBi$LOSS)
beta_g <-  mean(AutoBi$LOSS)^2 / var(AutoBi$LOSS)

# QQ-plot pour la loi de gamma 
quantile_gamma = data.frame(
   quantile_loss = sapply(seq(0.001,0.999,0.001), function(alpha) quantile(AutoBi$LOSS, alpha)),
   quantile_gamma = sapply(seq(0.001,0.999,0.001), function(alpha) qgamma(alpha, shape = alpha_g, scale = beta_g))
   )

ggplot(data = quantile_gamma, mapping = aes(x  = quantile_gamma , y = quantile_loss )) + geom_point()  + geom_abline(intercept = 0, slope = 1)

data.frame('alpha' = alpha_g, 'beta_g'=beta_g)

Ce résultat plutôt décevant pourrait être mis sur le dos de la méthode d’inférence car utiliser la méthode des moments si la distribution des pertes est à queue lourde n’est pas forcément pertinent (si par exemple les pertes extrêmes suivent une loi de Pareto).

Comparaison des modèles

Visuellement, on peut superposer les densités des lois estimées à l’histogramme des données.

# Histogramme des données et des densités paramétriques
lnorm_pdf <- function(x) dlnorm(x, meanlog = mu, sdlog = sigma)
gamma_pdf <- function(x) dgamma(x, shape = alpha_g , scale = beta_g)
weibull_pdf <- function(x) dweibull(x, shape = mle_weib[1], scale = mle_weib[2])

hist_loss_AutoBi <- ggplot(data = AutoBi %>% filter(LOSS < 25)) + geom_histogram(mapping = aes(x = LOSS, y = ..density..), bg = 'lightgrey', binwidth = 0.5) + xlab('Montants de sinsitres') + ylab('Densité') + stat_function(fun = lnorm_pdf, linetype="solid") +
  stat_function(fun = gamma_pdf, linetype="dotted") +
  stat_function(fun = weibull_pdf, linetype="dashed") +
  theme_bw() +
  theme(axis.title.x = element_text(size = 12),
        axis.title.y = element_text(size = 12),
        axis.text = element_text(size = 12)
        )

hist_loss_AutoBi 

L’histogramme montre bien le manque d’adéquation des modèles par rapport aux données. L’estimations des paramètres est trop influencé par l’occurence de valeurs extrêmes. On peut tout de même calculer le Bayesian Information criterium

$$BIC = -2\ln\left[p(X|\theta)\right] + k\ln(n),$$ où $k$ désigne le nombre de paramètres

data.frame(
   Modele = c('lognormal', 'gamma', 'Weibull'),
   BIC = c(-2 * sum(log(lnorm_pdf(AutoBi$LOSS))) + 2 * log(length( AutoBi$LOSS)),
-2 * sum(log(gamma_pdf(AutoBi$LOSS))) + 2 * log(length( AutoBi$LOSS)),
-2 * sum(log(weibull_pdf(AutoBi$LOSS))) + 2 * log(length( AutoBi$LOSS))
)
)

Le modèle lognormale semble le plus adapté, on peut tester statistiquement l’adéquation en utilisant par exemple le test de Cramer-von Mises basé sur le critère suivant

$$\text{C-vM} = \int (F_n(x)-F_X(x))^2\text{d}F_X(x)$$ où $F_n(x) = \frac{1}{n}\sum_{i=1}^{n}\mathbb{I}_{X_i\leq x}$ désigne la fonction de répartition empirique et $F_X(x):=F_X(x;\theta)$ est la fonction de répartition du modèle théorique. La statistique de test est estimé via

$$T = n\times \text{C-vM}_n = \frac{1}{12n}+\sum_{i=1}^{n}\left(\frac{2i-1}{2n}-\widehat{F}_{X}(X_{i:n})\right)^2,$$ où $X_{i:n}$, $i=1,\ldots,n$ sont les statistiques d’ordre et $\widehat{F}_{X}(x):=F_{X}(x;\widehat{\theta})$. Le test rejette l’hypothèse

$$(H_0):\text{ }X\sim F_X$$ Lorsque $T$ dépasse la valeur critique de la statistique de test correspondant au quantile à $95\%$ de la distribution de $T$ lorsque $H_0$ est vérifiée. Pour calculer la valeur critique, on utilise un algorithme de bootstrap paramétrique.

Pour $k$ allant de $1$ à $B=10,000$,

1. Simuler un échantillon $X^{(k)}\sim F_X$
2. Estimer $\widehat{\theta}^{(k)}$ sur la base de $X^{(k)}$
3. Calcul de la statistique de test $T^{(k)}$

Si $T>\text{quantile}(T^{(k)}, k=1,\ldots, B\text{; }95\%)$.

X <- AutoBi$LOSS
# Fonction pour calculer la statistique de test
compute_CvM <- function(X){
  n <- length(X)
  mu = mean(log(X))
  sigma = sd(log(X))
  X_sort <- sort(X)
  CvM = 1/ 12 / n + sum((sapply(1:n, function(i) (2*i-1) / 2 / n - plnorm(X_sort[i], meanlog = mu, sdlog = sigma)))^2)
  return(CvM)
  }

B <- 1000

# Estimation des paramètres
n <- length(X)
mu = mean(log(X))
sigma = sd(log(X))
# Calcul de la statistiques 
T <- compute_CvM(X)
T_boot <- {}
for(k in 1:B){
  X_k <- rlnorm(n, meanlog = mu, sdlog = sigma)
  T_boot[k] <- compute_CvM(X_k)
}

data.frame(Test_stat = T, critical_value = quantile(T_boot, 0.95), Reject =  T>quantile(T_boot, 0.95))

Le test rejette clairement l’hypothèse nulle.