Nous utiliserons des données disponibles dans le package ìnsuranceData.

library(insuranceData)
library(VGAM)
data(AutoBi)
AutoBi

1340 montants de sinistres auto incluant des dommages corporels, exprimés en milliers de Dollars. L’analyse de la distribution des montants démarre avec des statistiques descriptives, un histogramme et un boxplot

library(tidyverse)
library(ggplot2)
t(AutoBi %>% summarise(
  moyenne = mean(LOSS),
  ecart_type = sd(LOSS),
  minimum = min(LOSS),
  maximum = max(LOSS),
  Q1 = quantile(LOSS, 0.25),
  Q2 = quantile(LOSS, 0.5),
  Q3 = quantile(LOSS, 0.75)
))
##                   [,1]
## moyenne       5.953461
## ecart_type   33.136205
## minimum       0.005000
## maximum    1067.697000
## Q1            0.640000
## Q2            2.331000
## Q3            3.994750
ggplot(data = AutoBi) + geom_histogram(mapping = aes(x= LOSS), binwidth = 10)

ggplot(data = AutoBi) + geom_boxplot(aes(y = LOSS))

Il s’agit d’un cas classique avec une forte occurence de sinsitres de faible intensité et une occurence non négligeable de sinistres de plus forte intensité. Les graphiqueS qunatile-quantile sont très utile pour exhiber la présence de valeurs extrêmes, on compare les quantiles de

\[ \frac{x_i}{\frac{1}{n}x_i \sum_{i=1}^n},\text{ }i=1,\ldots,n \] et \[ \frac{\ln(x_i / \underset{i=1,\ldots, n}{min} x_i)}{\frac{1}{n}\sum_{i=1}^{n}\ln(x_i/ \underset{i=1,\ldots, n}{min} x_i))} \]

data_quantile = data.frame(
   quantile_loss = sapply(seq(0.001,0.999,0.001), function(alpha) quantile(AutoBi$LOSS, alpha)),
   quantile_exp = sapply(seq(0.001,0.999,0.001), function(alpha) qexp(alpha, rate = 1 / mean(AutoBi$LOSS))),
   quantile_log_exp = sapply(seq(0.001,0.999,0.001), function(alpha) qexp(alpha, rate = 1 / mean(log(AutoBi$LOSS / min(AutoBi$LOSS))))),
   quantile_log_loss = sapply(seq(0.001,0.999,0.001), function(alpha) quantile(log(AutoBi$LOSS / min(AutoBi$LOSS)), alpha))
)

ggplot(data = data_quantile, mapping = aes(x  = quantile_loss , y = quantile_exp )) + geom_point()  

ggplot(data = data_quantile, mapping = aes(x  = quantile_log_loss , y = quantile_exp )) + geom_point()

Calibration des lois classiques

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.

Analyse de valeurs extrêmes

Une technique classique consiste à grapher l’estimateur de Hill en faisant varier le seuil pour la loi de Pareto avec

\[ \widehat{\alpha}_j = \frac{n-j}{\sum_{i=j}^{n}\ln\left(x_{n-i:n}/x_{n-j:n}\right)}. \]

# On trie par ordre croissant les observations 
X_sort <- sort(X)
n <- length(X)


alpha_hat = sapply(1:(n-1), function(j) (n-j) / sum(sapply((j+1):n, function(i) log(X_sort[i] / X_sort[j]))))
plot(rev(alpha_hat))

Le seuil et le paramètre semble être autour de

j <- n - 150
data.frame(seuil = X_sort[j], alpha = (n-j) / sum(sapply((j+1):n, function(i) log(X_sort[i] / X_sort[j]))))

L’analyse menée ici est très approximative, il existe bien d’autres outils d’aide au diagnostique en ce qui concerne les valeurs extrêmes.

Calibration du modèle composite

La méthode précédente ne suppose pas une distribution pour le ventre et la queue de distribution. Le choix de telle ou telle loi paramétrique pour le ventre de la distribution va influencer la valeur estimer du seuil et du paramètre de la loi de Pareto sous jacente. Nous allons calibrer un modèle composite Weibull-Pareto de densité

\[ f(x) = \begin{cases} r\frac{f_1(x)}{F_1(\theta)}&\text{ pour }0\leq x\leq \theta,\\ (1-r)\frac{f_2(x)}{1-F_2(\theta)}&\text{ pour } x> \theta. \end{cases} \] avec

\[ f_{1}(x)=\frac{k}{\beta}\left(\frac{x}{\beta}\right)^{k-1}e^{-(x/\beta)^k},\text{ }x>0. \] et

\[ f_2(x)= \begin{cases} \frac{\alpha \theta^\alpha}{x^{\alpha+1}},&\text{ }x\geq \theta,\\ 0,&\text{ }x< \theta, \end{cases} \] La continuité et dérivabilité en \(\theta\) impose que

\[ r = \frac{\frac{\alpha}{\theta}\left[1-\exp\left(-\frac{k+\alpha}{k}\right)\right]}{\frac{\alpha}{\theta}+\frac{k}{\theta}\exp\left(-\frac{k+\alpha}{k}\right)}\text{ et }\beta = \left(\frac{k}{k+\alpha}\right)^{1/k}\theta \] Pour inférer les paramètres \(k,\alpha\) et \(\theta\), nous utilisons le maximum de vraisemblance

X <- AutoBi$LOSS
n <- length(X)
# Densité du modèle composite
pdf_weibull_pareto <- function(x, alpha, theta, k){
  # The scale parameter of the Weibull and the mixing parameter are fixed to ensure smoothness of the resulting pdf
  lam <- (k / (k + alpha))^(1 / k) * theta
  r <- VGAM::dpareto(theta, scale = theta, shape = alpha) * pweibull(theta, shape = k, scale = lam) /
    (VGAM::dpareto(theta, scale = theta, shape = alpha) * pweibull(theta, shape = k, scale = lam) +
       dweibull(theta, shape = k, scale = lam))

  ifelse(x < theta, r * dweibull(x, shape = k, scale = lam) / pweibull(theta, shape = k, scale = lam),
         (1 - r) * VGAM::dpareto(x, shape = alpha, scale = theta))
}


# Fonction qui donne l'oppposé de la log vraisemblance
log_like <- function(param){
  -sum(log(pdf_weibull_pareto(X, param[1], param[2], param[3])))
}

# La fonction nlm qui permet de minimiser la fonction en argument
mle_weib_pareto <- optim(c(1, 3, 1.2), log_like)$par
#estimate

data.frame(k = mle_weib_pareto[1],
           beta = (mle_weib_pareto[1] / (mle_weib_pareto[3] + mle_weib_pareto[1]))^(1 / mle_weib_pareto[1]) * mle_weib_pareto[2],
           r = mle_weib_pareto[3] / mle_weib_pareto[2] * (1-exp(-(mle_weib_pareto[3] + mle_weib_pareto[1]) / mle_weib_pareto[1])) / (mle_weib_pareto[3] / mle_weib_pareto[2] + mle_weib_pareto[1] / mle_weib_pareto[2] * exp(-(mle_weib_pareto[3] + mle_weib_pareto[1]) / mle_weib_pareto[1])),
           theta = mle_weib_pareto[2],
           alpha = mle_weib_pareto[3],
           BIC = -2 * sum(log(pdf_weibull_pareto(X, mle_weib_pareto[1], mle_weib_pareto[2], mle_weib_pareto[3]))) + 3 * log(n)
           )

Nous gagnons quelques points de BIC par rapport aux modèles classiques.

Calibration d’une loi mélange de Erlang

Nous allons tenter d’améliorer l’adéquation du modèle aux données en utilisant la loi Erlang mélange de densité

\[ f_{X}(x) = \sum_{i=1}^{m}\pi_{i}\frac{e^{-x/\beta}x^{z_i-1}}{\beta^{z_i}(z_i-1)!}. \]

Il s’agit d’un modèle de mélange avec un variable aléatoire latente \(Z\sim\pi\). On simplifie ici en considérat des mélange de Erlang avec de plus en plus de composante tel que \(z_i = i,\text{ }i =1,\ldots, m\). On rappelle qu’à partir de valeurs initiale (ou précédente), l’algorithme EM met à jour la valeur des paramètres de la façon suivante

\[ \pi^{(1)}_{i}=\frac{\sum_{k=1}^{n}\mathbb{P}\left(Z_i=z_i|x_k;\theta^{(0)}\right)}{n}\text{ et } \beta^{(1)}=\frac{\frac{1}{n}\sum_{k=1}^n x_k}{\sum_{i=1}^{m}\pi_i^{(1)}z_i}, \] avec

\[ \mathbb{P}\left(Z=z_i|X;\theta^{(0)}\right)=\frac{\pi^{ (0)}_{i}\frac{e^{-x/\beta^{(0)}}x^{z_i-1}}{\left(\beta^{(0)}\right)^{z_i}(z_i-1)!}}{\sum_{i=1}^{m}\pi^{ (0)}_{i}\frac{e^{-x/\beta^{(0)}}x^{z_i-1}}{\left(\beta^{(0)}\right)^{z_i}(z_i-1)!}}. \]

# Densité d'un mélange de Erlang
dMErl <- function(x, pi, beta){
  nb_erlang <- length(pi)
  sum(sapply(1:nb_erlang, function(i) pi[i] * dgamma(x, shape = i, scale = beta)))
}

# Fonction qui calcule la loi de proba de Z sachant un point x et la valeur courante des paramètres
p_z <- function(x, pi, beta, nb_erlang){
  sapply(1:nb_erlang, function(i) pi[i] * dgamma(x, shape = i, scale = beta)) / dMErl(x, pi, beta)
}

# Fonction pour inférer les paramètre d'un mélange de Erlang avec un nombre donnée de composantes
infer_MErl<- function(X, pi_0, beta_0, nb_erlang, eps_pi, eps_beta){
  n <- length(X)
  pi <- list()
  beta <- {}
  pi[[1]] <- pi_0
  beta[1] <- beta_0
  pi[[2]] <- apply(sapply(X, function(x) p_z(x, pi[[1]], beta[1], nb_erlang)), 1, sum) / n
  beta[2] <- mean(X) / (pi[[2]] %*% (1:nb_erlang))
  k <- 2
  while( sum(abs(pi[[k]]-pi[[k-1]])) > eps_pi & abs(beta[k] - beta[k-1]) > eps_beta){
    k <- k + 1
    pi[[k]] <- apply(sapply(X, function(x) p_z(x, pi[[k-1]], beta[k-1], nb_erlang)), 1, sum) / n
    beta[k] <- mean(X) / (pi[[k]] %*% (1:nb_erlang))
  }
  return(list('pi' = pi[[k]], 'beta' = beta[k]))
}

# Initialisation des paramètres
beta_0 <- 10 
eps_pi <- 0.05 
eps_beta <- 0.05
# Nombre de 
nb_max_erlang <- 20
params_MErl <- list()
BIC <- {}
for(nb_erlang in 2:nb_max_erlang){
  params_MErl[[nb_erlang - 1]] <- infer_MErl(X, rep(1/nb_erlang, nb_erlang), beta_0 = 10, nb_erlang, 
                                             eps_pi, eps_beta)
  BIC[nb_erlang-1] <- - 2 * sum(log(sapply(X, function(x) dMErl(x, pi = params_MErl[[nb_erlang - 1]]$pi, params_MErl[[nb_erlang - 1]]$beta)))) + length(params_MErl[[nb_erlang - 1]]$pi) * log(n)
}

plot(BIC)

params_MErl[[nb_max_erlang-1]]
## $pi
##  [1] 0.8700524425 0.0823127111 0.0117333549 0.0035322990 0.0022758300
##  [6] 0.0021349978 0.0021850514 0.0021657738 0.0020075964 0.0017432004
## [11] 0.0014494841 0.0011960970 0.0010181275 0.0009171549 0.0008784476
## [16] 0.0008913068 0.0009724391 0.0012253203 0.0021832110 0.0091251547
## 
## $beta
## [1] 3.909805

Le résultat est assez décevant, on constate bien une diminution du BIC au fur et à mesure qu’on ajoute des composantes au mélange.