Un modelo de credibilidad tiene el objetivo de establecer una prima para un grupo heterogéneo de pólizas, es decir, en el que se tienen distintos niveles de riesgo. Existen varios enfoques desde los que se aborda el tema de credibilidad. En esta nota mostramos el modelo de credibilidad de Bühlmann-Straub y el modelo de credibilidad bayesiana con distintos ejemplos usando los datos de la paquetería insuranceData.

Sean \(X_1,X_2,...,X_n\) observaciones pasadas de que representan el monto de reclamación y supongamos que el nivel de riesgo de una póliza está caracterizado por el parámetro \(\theta\), el cual varia para cada póliza, y es una realización de la v.a. \(\Theta\). Al variar entre cada póliza, asumimos que \(\theta\) tiene función de densidad de probabilidad \(\pi(\theta)\) y distribución \(\Pi(\theta)\).

Credibilidad de Bühlmann-Straub

El modelo de credibilidad de Bühlmann-Straub supone que tenemos un elemento que representa la exposición (\(m_i, i = 1,...,n\)), mientras que las observaciones \(X_1,X_2,...,X_n\) representan el promedio de pérdida por unidad de exposición.

A \[\mu(\theta)=E(X_j|\Theta=\theta)\] se le conoce como la media hipotética mientras que \[v(\theta) = Var(X_j|\Theta=\theta)\] es el proceso de varianza.

A la esperanza del proceso de varianza (\(E(v(\theta))\)) también se le conoce como within risk variance y a la varianza de la media hipotética (\(Var(\mu(\theta))\)) también se le llama between risk variance.

La prima de credibilidad de Bühlmann-Straub es \[\hat{X} = Z\sum_{i}^n \frac{m_i}{m}X_i + (1-Z)\hat{\mu}\] donde \(m=\sum_{i=1}^n m_i\) y \(Z\) es el factor de credibilidad \(Z_=\frac{m}{k+m}\) con \(k=\frac{E(v(\theta))}{Var(\mu(\theta))}\).

Ejemplo

En la siguiente tabla se muestra el costo promedio observado durante tres años para cuatro territorios.

Territory	Year	Risk_count	Avg_cost
A	2016	700	900
A	2017	750	800
A	2018	875	1000
B	2016	350	200
B	2017	400	550
B	2018	425	625
C	2016	100	1300
C	2017	125	1800
C	2018	175	2000
D	2016	675	750
D	2017	700	800
D	2018	725	925

Para usar la función cm() del paquete actuar debemos transformar el data.frame a una forma wider en la que tengamos solo una observación por territorio.

# Usamos la función pivot_wider del paquete tidyr.
datos_wider <- datos %>% pivot_wider(
  names_from = Year, 
  values_from = c(Avg_cost, Risk_count),
  names_sep = ".",
)
kable(datos_wider)

Territory	Avg_cost.2016	Avg_cost.2017	Avg_cost.2018	Risk_count.2016	Risk_count.2017	Risk_count.2018
A	900	800	1000	700	750	875
B	200	550	625	350	400	425
C	1300	1800	2000	100	125	175
D	750	800	925	675	700	725

Ajustamos el modelo y calculamos las primas de Bühlmann-Straub para cada territorio.

modelo_bs <- cm(~Territory, datos_wider, ratios = Avg_cost.2016:Avg_cost.2018, 
   weights = Risk_count.2016:Risk_count.2018)

summary(modelo_bs)

Call:
cm(formula = ~Territory, data = datos_wider, ratios = Avg_cost.2016:Avg_cost.2018, 
    weights = Risk_count.2016:Risk_count.2018)

Structure Parameters Estimators

  Collective premium: 962.4466 

  Between Territory variance: 114891.9 
  Within Territory variance: 12171436 

Detailed premiums

    Territory Indiv. mean Weight Cred. factor Cred. premium
    A          905.3763   2325   0.9564209     907.8634    
    B          472.8723   1175   0.9172964     513.3619    
    C         1762.5000    400   0.7906104    1594.9772    
    D          827.0833   2100   0.9519759     833.5840

En el resumen del modelo podemos ver que la esperanza de proceso de varianza es \[E(v(\theta)) = 12,171,436 \] y la varianza de la media hipotética es \[Var(\mu(\theta)) = 11,4891.9\] En el resumen también podemos ver la media del costo como Indiv.mean, el total de unidades de riesgo (Risk_count) como Weight y el factor de credibilidad (Cred.factor) para cada territorio así como la prima colectiva que es 962.4466. Con estos elementos y la fórmula de la prima B-S se calculan las primas para cada territorio, las cuales también están en dicha tabla.

Credibilidad Bayesiana

La prima Bayesiana es \[B_{n+1}=E(\mu(\Theta)|X_1,...,X_n)\] y es la mejor aproximación (usando mínimos cuadrados) a la prima de riesgo desconocida \(\mu(\theta)=E(X_t|\Theta=\theta)\) basandonos en las observaciones pasadas \(X_1,...,X_n\)

Puede demostrarse que \[B_{n+1}=z\overline{X} + (1-z)\mu_X\] donde \(\mu_X=E(\mu(\Theta))\) y \(z=\frac{n}{n+k}\) para alguna constante \(k\).

Para calcular la prima Bayesiana, la paquetería actuar provee la función cm() la cual recibe el parámetro "bayes" y la familia paramétrica de la verosimilitud (likelihood) que puede ser una de las siguientes:

"poisson", "exponencial", "gamma", "normal", "bernoulli", "binomial", "geometric", "negative binomial" o "pareto"

Dependiendo del dicha familia paramétrica la distribución a priori es la distribución conjugada correspondiente por lo que se deben pasar los parámetros para esta distribución.

Las familias paramétricas para la verosimilitud y su correspondiente familia paramétrica de distribución inicial conjugada son las siguientes:

\(X|\Theta = \theta \sim Poisson(\theta), \Theta \sim Gamma(\alpha,\lambda)\);
\(X|\Theta = \theta \sim Exponencial(\theta), \Theta \sim Gamma(\alpha,\lambda)\);
\(X|\Theta = \theta \sim Normal(\theta,\sigma_2^2), \Theta \sim Normal(\mu,\sigma_1^2)\);
\(X|\Theta = \theta \sim Bernoulli(\theta), \Theta \sim Beta(a,b)\);
\(X|\Theta = \theta \sim Geométrica(\theta), \Theta \sim Beta(a,b)\);
\(X|\Theta = \theta \sim Gamma(\tau,\theta), \Theta \sim Gamma(\alpha,\lambda)\);
\(X|\Theta = \theta \sim Binomial(v,\theta), \Theta \sim Beta(a,b)\);
\(X|\Theta = \theta \sim Binomial Negativa(r,\theta), \Theta \sim Beta(a,b)\)

Ejemplos

Severidad

Usando los datos de WorkersComp nos interesa la variable que representa el monto de compensaciones por incapacidad parcial permanente (LOSS). Para cada clase de trabajadores (CL) se tienen observaciones de 7 años, por lo que usaremos las observaciones de los primeros 6 años para calcular \(B_{7}\) y compararla con la observación del año 7.

A continuación se presenta cada modelo usado y la predicción obtenida con cada uno.

# SEVERIDAD -------------
# Cargamos los datos.
library(insuranceData)
data("WorkersComp")

# Seleccionamos las observaciones de la clase 1.
X <- WorkersComp %>% filter(CL==1) %>% select(LOSS)
# Monto de pérdida anual observado para 7 años
head(X,n=7)

# Tenemos 7 observaciones solo usaremos 6 para ajustar los modelos y compararemos las predicciones de cada modelo con la observación 7.
X_n <- X[7,] # observaciones que usaremos para el cálculo.
X <- X[1:6,] # monto que queremos predecir.

## MODELOS

# Exponencial - Gamma (POCO INFORMATIVO)
exp_pinf <- cm("bayes",X,likelihood = "exponential",shape=0.01,rate=0.01)
summary(exp_pinf)

Call:
cm(formula = "bayes", data = X, likelihood = "exponential", shape = 0.01, 
    rate = 0.01)

Structure Parameters Estimators

  Collective premium: -0.01010101 

  Between variance: -0.00005127156 
  Within variance: 0.00005075884 

Detailed premiums

    Indiv. mean Weight Cred. factor Bayes premium
    783331.7    6      1.197605     938121.8

predict(exp_pinf)

[1] 938121.8

# Exponencial - Gamma (INFORMATIVO)
exp_inf <- cm("bayes",X,likelihood = "exponential",shape=10,rate=800)
summary(exp_inf)

Call:
cm(formula = "bayes", data = X, likelihood = "exponential", shape = 10, 
    rate = 800)

Structure Parameters Estimators

  Collective premium: 88.88889 

  Between variance: 987.6543 
  Within variance: 8888.889 

Detailed premiums

    Indiv. mean Weight Cred. factor Bayes premium
    783331.7    6      0.4          313386

predict(exp_inf)

[1] 313386

# Normal - Normal
norm_pinf <- cm("bayes",X,likelihood = "normal",sd.lik=10000,mean=0,sd=10000)
summary(norm_pinf)

Call:
cm(formula = "bayes", data = X, likelihood = "normal", sd.lik = 10000, 
    mean = 0, sd = 10000)

Structure Parameters Estimators

  Collective premium: 0 

  Between variance: 100000000 
  Within variance: 100000000 

Detailed premiums

    Indiv. mean Weight Cred. factor Bayes premium
    783331.7    6      0.8571429    671427.1

predict(norm_pinf)

[1] 671427.1

# Gamma - Gamma
gm_pinf <- cm("bayes",X,likelihood = "gamma",shape.lik=1000,shape=0.01,rate=0.01)
summary(gm_pinf)

Call:
cm(formula = "bayes", data = X, likelihood = "gamma", shape.lik = 1000, 
    shape = 0.01, rate = 0.01)

Structure Parameters Estimators

  Collective premium: -10.10101 

  Between variance: -51.27156 
  Within variance: 0.05075884 

Detailed premiums

    Indiv. mean Weight Cred. factor Bayes premium
    783331.7    6      1.000165     783460.9

predict(gm_pinf)

[1] 783460.9

En la siguiente tabla se muestran las estimaciones obtenidas con cada modelo y el valor observado para el año 7.

Observación	609833.0
Exponencial (poco informativa)	938121.8
Exponencial (informativa)	313386.0
Normal	671427.1
Gamma	783460.9

Frecuencia

Usaremos los datos norfire de la paquetería CASdatasets. En este caso la variable de interés es el número de reclamaciones anuales por incendio. Se tienen observaciones en el periodo de 1972 a 1992 provenientes de una aseguradora noruega.

Al igual que en el ejemplo anterior, no tomaremos en cuenta la observación del último año con la finalidad de comparar con las estimaciones obtenidas de cada modelo.

# FRECUENCIA ---------
# Cargamos los datos.
library(CASdatasets)
data("norfire")
# Cada observación es una reclamación y cada una tiene asociado el año de ocurrencia.
head(norfire)

  Year Loss Loss2012
1 1972  520 3345.106
2 1972  529 3403.002
3 1972  530 3409.435
4 1972  530 3409.435
5 1972  544 3499.495
6 1972  545 3505.928

# Obtenemos el número de reclamaciones por año.
norfire_ncl <- norfire %>% group_by(Year) %>% summarize(Num_claims=n())
# Primeros 6 años de observación y su correspondiente número de reclamaciones.
head(norfire_ncl)

# A tibble: 6 x 2
   Year Num_claims
  <dbl>      <int>
1  1972         97
2  1973        109
3  1974        110
4  1975        142
5  1976        207
6  1977        235

# Tenemos 21 observaciones solo usaremos 20 para ajustar los modelos y compararemos las predicciones de cada modelo con la observación 21.
X_ncl <- norfire_ncl$Num_claims[-21] # observaciones que usaremos para el cálculo.
X_ncl_n <- norfire_ncl$Num_claims[21] # Núm. de reclamaciones que queremos predecir.

# MODELOS
# Poisson - Gamma
poi_pinf <- cm("bayes",X_ncl,likelihood = "poisson",shape=0.01,rate=0.01)
summary(poi_pinf)

Call:
cm(formula = "bayes", data = X_ncl, likelihood = "poisson", shape = 0.01, 
    rate = 0.01)

Structure Parameters Estimators

  Collective premium: 1 

  Between variance: 100 
  Within variance: 1 

Detailed premiums

    Indiv. mean Weight Cred. factor Bayes premium
    428.3       20     0.9995002    428.0865

predict(poi_pinf)

[1] 428.0865

# Binomial Negativa - Beta
bn_pinf <- cm("bayes",X_ncl,likelihood = "negative binomial",size=10,shape1=0.1,shape2=1)
summary(bn_pinf)

Call:
cm(formula = "bayes", data = X_ncl, likelihood = "negative binomial", 
    size = 10, shape1 = 0.1, shape2 = 1)

Structure Parameters Estimators

  Collective premium: -11.11111 

  Between variance: -6.497726 
  Within variance: 0.5847953 

Detailed premiums

    Indiv. mean Weight Cred. factor Bayes premium
    428.3       20     1.00452      430.2863

predict(bn_pinf)

[1] 430.2863

En la siguiente tabla se muestran las estimaciones obtenidas con cada modelo y el valor observado.

Observación	615.0000
Poisson	428.0865
Binomial Negativa	430.2863

Los modelos usados en estos ejemplos son para ilustrar el uso de la función cm(), sin embargo, para obtener una buena estimación se debe tomar en cuenta la forma de elegir los parámetros de la distribución a priori.