Distribución lognormal

Usando los datos de seguros de automotor para privados de la data set CASdatasets y tomando el número de cédula del compañero Ronald Palencia, calculamos una muestra de la base de datos y de la función.

Ajuste distribucional

Lo primero que se debe hacer es identificar el tipo de distribución que tiene la variable, para esto usaremos paquetes existentes en R como :

• univariateML: estimación de máxima verosimilitud para densidades univariadas.

• fitdistrplus: ayuda para ajustar una distribución paramétrica a datos censurados o no censurados.

• gamlss: modelos aditivos generalizados para ubicación, escala y forma.

Estos paquetes nos ayudaran a ajustar las distrubuciones.

Ahora ajustando la variable IDclaim como una variable númerica obtenemos una tabla resumen de los datos.

# Ajuste variable numerica
data$IDclaim <- data$IDclaim %>%  as.numeric() 
summary(data)  %>%
  kbl() %>%
  kable_paper("hover", 
              full_width = F)

Ahora filtramos los valores de la variable payment mayores a cero para hacer un mejor manejo de la data.

data <- data %>% dplyr::filter(Payment > 0)

summary(data) %>%
  kbl() %>%
  kable_paper("hover", 
              full_width = F)

p1 <- ggplot(data = data, aes(x = data$Payment)) +
     geom_histogram(fill = "steelblue") +
     labs(title = "Distribución del precio de Payment",
          x = "Precio") +
     theme_bw() +
     theme(plot.title = element_text(face = "bold"))

p2 <- ggplot(data = data, aes(x = data$Payment)) +
      stat_ecdf(geom = "step", color = "steelblue", size = 1) +
      labs(title = "Función de distribución empírica",
           x = "Precio",
           y = "CDF") +
      theme_bw() +
      theme(plot.title = element_text(face = "bold"))

ggpubr::ggarrange(plotlist = list(p1, p2), ncol = 2)

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

# Se comparan únicamente las distribuciones con un dominio [0, +inf)
comparacion_aic <- AIC(
                    mlbetapr(data$Payment),
                    mlexp(data$Payment),
                    mlinvgamma(data$Payment),
                    mlgamma(data$Payment),
                    mllnorm(data$Payment),
                    mlrayleigh(data$Payment),
                    mlinvgauss(data$Payment),
                    mlweibull(data$Payment),
                    mlinvweibull(data$Payment),
                    mllgamma(data$Payment)
                   )

comparacion_aic %>% rownames_to_column(var = "distribucion") %>% arrange(AIC) 

##                  distribucion df      AIC
## 1       mllnorm(data$Payment)  2 2286.346
## 2      mllgamma(data$Payment)  2 2289.562
## 3    mlinvgauss(data$Payment)  2 2298.129
## 4  mlinvweibull(data$Payment)  2 2306.168
## 5      mlbetapr(data$Payment)  2 2313.264
## 6    mlinvgamma(data$Payment)  2 2313.994
## 7     mlweibull(data$Payment)  2 2315.688
## 8       mlgamma(data$Payment)  2 2329.234
## 9         mlexp(data$Payment)  1 2340.433
## 10   mlrayleigh(data$Payment)  1 2865.163

# Se comparan únicamente las distribuciones con un dominio [0, +inf)
comparacion_bic <- BIC(
                    mlbetapr(data$Payment),
                    mlexp(data$Payment),
                    mlinvgamma(data$Payment),
                    mlgamma(data$Payment),
                    mllnorm(data$Payment),
                    mlrayleigh(data$Payment),
                    mlinvgauss(data$Payment),
                    mlweibull(data$Payment),
                    mlinvweibull(data$Payment),
                    mllgamma(data$Payment)
                    )
comparacion_bic %>% rownames_to_column(var = "distribucion") %>% arrange(BIC) 

##                  distribucion df      BIC
## 1       mllnorm(data$Payment)  2 2292.327
## 2      mllgamma(data$Payment)  2 2295.543
## 3    mlinvgauss(data$Payment)  2 2304.109
## 4  mlinvweibull(data$Payment)  2 2312.149
## 5      mlbetapr(data$Payment)  2 2319.245
## 6    mlinvgamma(data$Payment)  2 2319.975
## 7     mlweibull(data$Payment)  2 2321.669
## 8       mlgamma(data$Payment)  2 2335.215
## 9         mlexp(data$Payment)  1 2343.423
## 10   mlrayleigh(data$Payment)  1 2868.153

Ajuste de la distribución log-normal

# Ajuste de una distribución log-normal
distribucion1 <- mllnorm(x = data$Payment)

distribucion <- fitdist(data$Payment, distr = "lnorm")
summary(distribucion)

## Fitting of the distribution ' lnorm ' by maximum likelihood 
## Parameters : 
##         estimate Std. Error
## meanlog 6.102583 0.10501440
## sdlog   1.273232 0.07425619
## Loglikelihood:  -1141.173   AIC:  2286.346   BIC:  2292.327 
## Correlation matrix:
##         meanlog sdlog
## meanlog       1     0
## sdlog         0     1

Por lo tanto, tenemos que la media estimada para la lognormal usando maxima verosimilitud es de \(6.10\) con una desviacón estadar de \(1.2\), también podemos ver que lo errores de estimación son relativamente bajos.

plot(distribucion)

hist(data$Payment,
     main = "Distribución precio Payment",
     freq = FALSE,
     ylim = c(0, 0.002))
lines(mllnorm(data$Payment), lwd = 2, lty = 1, col = "blue")
lines(mlinvgauss(data$Payment), lwd = 2, lty = 2, col = "red")
legend(x = 15000, y = 0.001, legend = c("lnorm", "invgauss"),
       col = c("blue", "red"), lty = 1:2)
rug(data$Payment)

ggplot(data = data) +
  geom_histogram(aes(x = Payment, y =  after_stat(density)),
                 bins = 40,
                 alpha = 0.3, color = "black") +
  geom_rug(aes(x = Payment)) +
  stat_function(fun = function(.x){dml(x = .x, obj = mllnorm(data$Payment))},
                aes(color = "log-normal"),
                size = 1) +
  stat_function(fun = function(.x){dml(x = .x, obj = mlinvgauss(data$Payment))},
                aes(color = "inverse-normal"),
                size = 1) +
  scale_color_manual(breaks = c("log-normal", "inverse-normal"),
                     values = c("log-normal" = "red", "inverse-normal" = "blue")) +
  labs(title = "Distribución precio Payment",
       color = "Distribución") +
  theme_bw() +
  theme(legend.position = "bottom")

Valor en riesgo VaR

Calculamos VaR(\(\alpha\)) con \(\alpha = 0.05\)

VaR <- qlnorm(0.95,meanlog=distribucion$estimate['meanlog'],sdlog =distribucion$estimate['sdlog'])

Por tanto tenemos que \[ VaR = 3629.5487 \]

Deficit esperado ES

Calculamos ES(\(\alpha\)) con \(\alpha = 0.05\)

ES <- mean(data$Payment[data$Payment>=VaR])

Por tanto tenemos que \[ES = 8777\]

Reservas necesarias al inicio de un periodo UV

Con un portafolio de \(20\) pólizas, una tasa de reclamos por poliza de \(0.015\) y un cargo de seguridad \(\lambda = 0.05\) en la prima, entonces para que la aseguradora pueda honrar los contratos con una confianza del \(95\%\) el número de reservas debe ser de:

n <- 20*0.015
m <- mean(data$Payment)
lambda <- 0.05
(UV <- VaR-(1+lambda)*n*m)

## [1] 3299.688

Es decir las reservas tiene un valor de \[U_V = 3299.69 \]

Reservas más confiables de una forma sensible UE

Ahora calculemos las reservas más confiables de forma más sensible

(UE <- ES-(1+lambda)*n*m)

## [1] 8447.139

Por lo tanto, el valor de las reservas necesarias más confiables de forma más sensible para que la aseguradora pueda honrar los contratos con una confianza del \(95\%\) es de:

\[U_E = 8447.139\]

Ajuste una distribución de Pareto

Para la solución de este punto primero verificamos que tennga solución el sistema de ecuaciones, que se definen de la siguiente manera:

\[\gamma = \frac{-e^{\mu + \frac{\sigma^2}{2}}}{\beta - -e^{\mu + \frac{\sigma^2}{2}}}\\ \beta = exp(\mu + z_{1-\alpha}\sigma)\alpha^{\frac{1}{\gamma}}\]

La ecuaciones anteriores se muestran en el siguiente gráfico

modelo_LN = distribucion
mu_<-modelo_LN$estimate[1]/sqrt(1+(modelo_LN$estimate[2]/modelo_LN$estimate[1])^2)
sigma_<-sqrt(log(1+modelo_LN$estimate[2]^2/modelo_LN$estimate[1]^2 ))
significancia<-0.05
ecu_1<-exp(mu_+sigma_^2/2)
ecu_2<-exp(mu_+qnorm(1-significancia)*sigma_)

gamma_beta<-function(beta_){
  -ecu_1/(beta_-ecu_1)
}

beta_gamma<-function(gamma_){
  ecu_2*significancia^(1/gamma_)
}
beta_1<-seq(0.01,ecu_1*0.9,length.out=1000)
gamma_1<-gamma_beta(beta_1)
gamma_2<-seq(0.01,max(gamma_1), length.out=1000)
beta_2<-beta_gamma(gamma_2)

#plot(beta_1,gamma_1,type="l",col="red")
#points(beta_2,gamma_2,type="l",col="blue")

datos<-data.frame(gamma=c(gamma_1,gamma_2), beta=c(beta_1,beta_2),ecuacion=c(rep("1",1000),rep("2",1000))  ) 

fig <- plot_ly(data=datos, x = ~beta, y = ~gamma, type = 'scatter', mode = 'lines',color = ~ecuacion)
fig

Del gráfico anterio podemos notar que el sistema si tiene solución, ademas logramos ver que no solo tiene una, sino que presenta dos puntos de corte, dichos puntos son \((36.18, 1.099)\) y \((334.17, 5.97)\)

beta1.1 = 36.18
beta1.2 = 334.17
gamma2.1 = 1.099
gamma2.2=5.977

par(mfrow = c(1,2))
u <-rpareto(1000, gamma2.1, beta1.1)

hist(u, freq = FALSE, breaks = 100,
     main = "Histograma de aleatorios de paretou(5.97, 334.17)", ylab = "Densidad", col = "blue")

u1 <-rpareto(1000, gamma2.2, beta1.2)

hist(u1, freq = FALSE, breaks = 100,
     main = "Histograma de aleatorios de paretou(1.099,36.18)", ylab = "Densidad",col="blue")

Haciendo la comparación de los dos gráficos anteriores y el siguiente notamos que la mejor distribución pareto se da con los siguientes parametros:

\[\beta = 334.17\]

\[\gamma = 5.97\]

Ya que es la que toma mas información de las colas de la log normal y ademas tiene un mejor ajuste.

ggplot(data = data) +
  geom_histogram(aes(x = Payment, y =  after_stat(density)),
                 bins = 40,
                 alpha = 0.3, color = "black") +
  geom_rug(aes(x = Payment)) +
  stat_function(fun = function(.x){dml(x = .x, obj = mllnorm(data$Payment))},
                aes(color = "log-normal"),
                size = 1) +
  stat_function(fun = function(.x){dml(x = .x, obj = mlinvgauss(data$Payment))},
                aes(color = "inverse-normal"),
                size = 1) +
  scale_color_manual(breaks = c("log-normal", "inverse-normal"),
                     values = c("log-normal" = "red", "inverse-normal" = "blue")) +
  labs(title = "Distribución precio Payment",
       color = "Distribución") +
  theme_bw() +
  theme(legend.position = "bottom")

Con base en la gráfica cual cree Ud que sería la utilidad de la distribución de Pareto para una aseguradora o Reaseguradora?.

La importancia de la distribución pareto para la las aseguradoras es que se dan cuenta de los precios que ellos estaria asumiedo, ya que muestra las colas con más valores elevados, por lo que le seria de mucha ayda para identificar cuales serian esas polizas de las que no se haría carga la aseguradora y si la reaseguradora.