6 Estimación paramétrica de modelos de frecuencia y severidad

6.1 Modelos continuos

En los modelos paramétricos continuos, se asume que los datos siguen una distribución de probabilidad con un conjunto fijo de parámetros.\(\\\)

Algunos modelos asumen que las variables de tiempo de falla o pérdida siguen un cierto familia de distribuciones.\(\\\)

Los principales métodos de cálculo:

Método de los momentos
Máxima verosimilitud
Estimación bayesiana

6.1.1 Método de los momentos

En la estimación del método de momentos \(\hat{\theta}\) es la solución de \(\theta\) en las ecuaciones.

\[ \mu_{r}^{\prime}(\theta)=\hat{\mu}_{r}^{\prime}, \quad \text { para } r=1, \cdots, k \]

Por lo tanto, tenemos un conjunto de \(k\) ecuaciones que involucran \(k\) incógnitas \(\theta_1, \cdots, \theta_k\).

Asumimos que existe una solución a las ecuaciones sin embargo, puede haber varias , en cuyo caso la estimación del método de momentos no es único.

6.1.2 Método de máxima verosimilitud

Suponga que tenemos una muestra aleatoria de \(n\) observaciones de \(X\), denotado por \(x = (x_1, · · ·, x_n)\). Dado la pdf o pf de \(X\), \(f (\cdot, \theta)\), definimos la función de verosimilitud de la muestra como el producto de \(f (x_i; \theta)\), denotado por \(L (\theta, x)\). Por lo tanto, tenemos

\[ L(\theta ; \boldsymbol{x})=\prod_{i=1}^{n} f\left(x_{i} ; \theta\right) \]

que se toma en función de \(\theta\) dado \(x\). Como las observaciones son independientes, \(L(\theta, x)\) es la probabilidad conjunta o la densidad conjunta de las observaciones.

Definimos la función logarítmica de verosimilitud como el logaritmo de \(L(\theta, x)\), es decir: \[ \log L(\theta ; \boldsymbol{x})=\sum_{i=1}^{n} \log f\left(x_{i} ; \theta\right) \]

El valor de \(\theta\), denotado por \(\hat{\theta}\), que maximiza la función de verosimilitud se llama estimador de máxima verosimilitud (MLE) de \(\theta\).\(\\\)

Aplicamos la función de verosimilitud a los datos:

\[ L(\theta)=\prod_{j=1}^{n} f_{X_{j}}\left(x_{j} \mid \theta\right), \quad l(\theta)=\sum_{j=1}^{n} \ln f_{X_{j}}\left(x_{j} \mid \theta\right) \]

6.1.3 Máxima verosimilitud para datos agrupados

Considere el conjunto de números \(c_0 < c_1 < \cdots < c_k\) y sea \(n_j\) número de observaciones en el intervalo .\((c_{j-1},c_j]\) para tales datos, la función de verosimilitud es:

\[ L(\theta)=\prod_{j=1}^{k}\left[F\left(c_{j} \mid \theta\right)-F\left(c_{j-1} \mid \theta\right)\right]^{n_{j}} \] y su logaritmo es: \[ l(\theta)=\sum_{j=1}^{k} n_{j} \ln \left[F\left(c_{j} \mid \theta\right)-F\left(c_{j-1} \mid \theta\right)\right] \]

Determine el estimador de máxima verosimilitud de la exponencial para los siguientes datos:

6.1.4 Máxima verosimilitud para datos censurados

Ahora consideramos el caso en el que tenemos observaciones individuales que son derecho censurado. Hay un límite de póliza de \(u\), solo reclamos de montos en el intervalo \((0, u]\) son posibles. Las pérdidas que excedan el valor de \(u\) tienen probabilidad \(1 - F(u,\theta )\).

Tenemos las observaciones \(x = (x_1, \ldots , x_{n_1})\), donde \(0 < x_1,\cdots, x_{n_1} < u,\) y \(n_2\) son los reclamos de monto \(u\), con \(n = n_1 +n_2\), entonces la función de verosimilitud es:

\[ L\left(\theta ; \boldsymbol{x}, n_{2}\right)=\left[\prod_{i=1}^{n_{1}} f\left(x_{i} ; \theta\right)\right][1-F(u ; \theta)]^{n_{2}} \] y el logaritmo es: \[ \log L\left(\theta ; \boldsymbol{x}, n_{2}\right)=n_{2} \log [1-F(u ; \theta)]+\sum_{i=1}^{n_{1}} \log f\left(x_{i} ; \theta\right) \]

6.1.5 Máxima verosimilitud con datos truncados

Si la póliza de seguro tiene un deducible de \(d\), los datos de los pagos de la reclamación se muestrean de una población con truncamiento, es decir, solo pérdidas con cantidades exceden \(d\) son muestreados.

La función de verosimilitud es

\[ L(\theta ; \boldsymbol{x})=\prod_{i=1}^{n} \frac{f\left(x_{i} ; \theta\right)}{1-F(d ; \theta)}=\frac{1}{[1-F(d ; \theta)]^{n}} \prod_{i=1}^{n} f\left(x_{i} ; \theta\right) \]

y el logaritmo es:

\[ \log L(\theta ; \boldsymbol{x})=-n \log [1-F(d ; \theta)]+\sum_{i=1}^{n} \log f\left(x_{i} ; \theta\right) \]

6.1.6 Ejemplos

library(flexsurv)

Loading required package: survival


Attaching package: 'flexsurv'

The following objects are masked from 'package:actuar':

    dllogis, pllogis, qllogis, rllogis

library(actuar)

tabla_d <- read.csv(file = "./data/tabla_d.csv")
tabla_d <- tabla_d %>%
  mutate(delta = case_when(
    event == "d" ~ 1,
    event %in% c("s", "e") ~ 0
  ))

surv1 <-
  Surv(
    time = tabla_d$d,
    time2 = tabla_d$obs,
    event = tabla_d$delta
  )




custom_pareto <-
  list(
    name = "pareto",
    pars = c("shape", "scale"),
    location = c("scale"),
    transforms = c(log, log),
    inv.transforms = c(exp, exp)
  )

parametric_pareto <- flexsurvreg(
  surv1 ~ 1,
  dist = custom_pareto,
  inits = c(1, 1)
)

Forming integrated rmst function...

Forming integrated mean function...

summary(parametric_pareto)

 
   time       est       lcl ucl
1   0.1 0.9939674 0.9892973   1
2   0.5 0.9701999 0.9481873   1
3   0.8 0.9527492 0.9195127   1
4   1.8 0.8968222 0.8275895   1
5   2.1 0.8806953 0.8023068   1
6   2.5 0.8596443 0.7721217   1
7   2.8 0.8441880 0.7490045   1
8   2.9 0.8390981 0.7377343   1
9   3.1 0.8290104 0.7224203   1
10  3.9 0.7898617 0.6616148   1
11  4.0 0.7851004 0.6522493   1
12  4.1 0.7803678 0.6418248   1
13  4.8 0.7480315 0.6041672   1
14  5.0 0.7390419 0.5916442   1

parametric_pareto

Call:
flexsurvreg(formula = surv1 ~ 1, dist = custom_pareto, inits = c(1, 
    1))

Estimates: 
       est       L95%      U95%      se      
shape  3.15e+02  4.02e-20  2.47e+24  8.11e+03
scale  5.21e+03  6.53e-19  4.15e+25  1.34e+05

N = 40,  Events: 8,  Censored: 32
Total time at risk: 132.1
Log-likelihood = -30.43446, df = 2
AIC = 64.86891

plot(parametric_pareto)

parametric_gamma <- flexsurvreg(
  surv1 ~ 1,
  dist = "gamma",
  inits = c(1, 1)
)

summary(parametric_gamma)

 
   time       est       lcl       ucl
1   0.1 0.9999728 0.9946796 1.0000000
2   0.5 0.9983177 0.9691321 0.9999939
3   0.8 0.9946061 0.9488375 0.9998864
4   1.8 0.9636197 0.8661713 0.9939274
5   2.1 0.9488755 0.8398372 0.9893672
6   2.5 0.9257877 0.7989214 0.9785396
7   2.8 0.9061905 0.7683855 0.9686592
8   2.9 0.8992701 0.7598907 0.9651045
9   3.1 0.8849010 0.7417467 0.9560545
10  3.9 0.8216590 0.6344052 0.9130853
11  4.0 0.8132444 0.6199456 0.9082173
12  4.1 0.8047432 0.6044254 0.9019412
13  4.8 0.7434278 0.4932738 0.8591167
14  5.0 0.7255315 0.4658636 0.8436135

parametric_gamma

Call:
flexsurvreg(formula = surv1 ~ 1, dist = "gamma", inits = c(1, 
    1))

Estimates: 
       est     L95%    U95%    se    
shape  2.6167  1.0755  6.3662  1.1870
rate   0.3020  0.0795  1.1475  0.2057

N = 40,  Events: 8,  Censored: 32
Total time at risk: 132.1
Log-likelihood = -28.52685, df = 2
AIC = 61.0537

plot(parametric_gamma)

La base dataCarse basa en pólizas de seguro de vehículos de un año contratadas en 2004 o 2005. Hay 67856 pólizas. La variable claimcst0 tiene los reclamos incluidos con el 0.

library(insuranceData)
library(fitdistrplus)

Loading required package: MASS


Attaching package: 'MASS'

The following object is masked from 'package:dplyr':

    select

library(survminer)

Loading required package: ggpubr


Attaching package: 'survminer'

The following object is masked from 'package:survival':

    myeloma

data("dataCar")

summary(dataCar)

   veh_value         exposure             clm            numclaims      
 Min.   : 0.000   Min.   :0.002738   Min.   :0.00000   Min.   :0.00000  
 1st Qu.: 1.010   1st Qu.:0.219028   1st Qu.:0.00000   1st Qu.:0.00000  
 Median : 1.500   Median :0.446270   Median :0.00000   Median :0.00000  
 Mean   : 1.777   Mean   :0.468651   Mean   :0.06814   Mean   :0.07276  
 3rd Qu.: 2.150   3rd Qu.:0.709103   3rd Qu.:0.00000   3rd Qu.:0.00000  
 Max.   :34.560   Max.   :0.999316   Max.   :1.00000   Max.   :4.00000  
                                                                        
   claimcst0          veh_body        veh_age      gender    area     
 Min.   :    0.0   SEDAN  :22233   Min.   :1.000   F:38603   A:16312  
 1st Qu.:    0.0   HBACK  :18915   1st Qu.:2.000   M:29253   B:13341  
 Median :    0.0   STNWG  :16261   Median :3.000             C:20540  
 Mean   :  137.3   UTE    : 4586   Mean   :2.674             D: 8173  
 3rd Qu.:    0.0   TRUCK  : 1750   3rd Qu.:4.000             E: 5912  
 Max.   :55922.1   HDTOP  : 1579   Max.   :4.000             F: 3578  
                   (Other): 2532                                      
     agecat                     X_OBSTAT_    
 Min.   :1.000   01101    0    0    0:67856  
 1st Qu.:2.000                               
 Median :3.000                               
 Mean   :3.485                               
 3rd Qu.:5.000                               
 Max.   :6.000

descdist(dataCar$claimcst0, discrete = FALSE)

summary statistics
------
min:  0   max:  55922.13 
median:  0 
mean:  137.2702 
estimated sd:  1056.298 
estimated skewness:  17.5025 
estimated kurtosis:  482.9257

descdist(dataCar$claimcst0[dataCar$clm == 1], discrete = FALSE)

summary statistics
------
min:  200   max:  55922.13 
median:  761.565 
mean:  2014.404 
estimated sd:  3548.907 
estimated skewness:  5.04178 
estimated kurtosis:  43.25942

surv2 <- Surv(time = dataCar$claimcst0)


km_policy <- survfit(surv2 ~ 1)


ggsurvplot(km_policy, data = dataCar)

surv3 <- Surv(time = dataCar$claimcst0[dataCar$clm == 1])
km_policy <- survfit(surv3 ~ 1)
ggsurvplot(km_policy, data = dataCar)

data_car_gamma <- flexsurvreg(surv3 ~ 1, dist = "gamma", inits = c(1, 1))

plot(data_car_gamma)

Con la base de datos dataCar suponga que se impone una límite en la poliza de 10,000. ¿Cómo se verían los estimadores anteriores con este cambio?

6.2 Modelos discretos

Para todos los modelos discretos de tipo \((a,b,0)\) se estiman los parámetros de la misma forma que ustedes saben desde estadística I.

La diferencia reside en los modelos \(a,b,1\), donde se debe considerar \(p_0^{M}\) o \(p_0^{T}\) como parámetros adicionales a ajustar.

Si quisiéramos modelar una clase \((a,b,1)\) con modificación en 0, entonces tendríamos lo siguiente

\[\begin{equation*} L=\left(p_0^M\right)^{n_0} \prod_{k=1}^{\infty}\left(p_k^M\right)^{n_k}=\left(p_0^M\right)^{n_0} \prod_{k=1}^{\infty}\left[\left(1-p_0^M\right) p_k^T\right]^{n_k} \end{equation*}\] The loglikelihood is \[\begin{equation*} \begin{aligned} l &=n_0 \ln p_0^M+\sum_{k=1}^{\infty} n_k\left[\ln \left(1-p_0^M\right)+\ln p_k^T\right] \\ &=n_0 \ln p_0^M+\sum_{k=1}^{\infty} n_k \ln \left(1-p_0^M\right)+\sum_{k=1}^{\infty} n_k\left[\ln p_k-\ln \left(1-p_0\right)\right] \end{aligned} \end{equation*}\]

La última expresión es dado que \(p_k^{T}= \frac{p_k}{1-p_k}\).

Para mayor facilidad, \[\begin{equation*} l=l_0+l_1 \end{equation*}\] con \[\begin{equation*} \begin{aligned} &l_0=n_0 \ln p_0^M+\sum_{k=1}^{\infty} n_k \ln \left(1-p_0^M\right) \\ &l_1=\sum_{k=1}^{\infty} n_k\left[\ln p_k-\ln \left(1-p_0\right)\right] \end{aligned} \end{equation*}\]

\[\begin{equation*} \frac{\partial l}{\partial p_0^M}=\frac{n_0}{p_0^M}-\sum_{k=1}^{\infty} \frac{n_k}{1-p_0^M}=\frac{n_0}{p_0^M}-\frac{n-n_0}{1-p_0^M}, \end{equation*}\]

Resultando en \[\begin{equation*} \hat{p}_0^{M} = \frac{n_0}{n}. \end{equation*}\]

Suponga que se tiene una binomial, y defina como \(m\) como el punto con observaciones de la tabla.

\[\begin{align*} l &=\sum_{k=0}^m n_k \ln p_k \\ &=\sum_{k=0}^m n_k\left[\ln \left(\begin{array}{l} m \\ k \end{array}\right)+k \ln q+(m-k) \ln (1-q)\right] . \end{align*}\]

Si \(m\) es fijo entonces

\[\begin{equation*} \hat{q}=\frac{1}{m} \frac{\sum_{k-0}^m k n_k}{\sum_{k=0}^m n_k} \end{equation*}\]

Si \(m\) es desconocido entonces

\[\begin{equation*} \hat{q}=\frac{1}{\hat{m}} \frac{\sum_{k-0}^{\infty} k n_k}{\sum_{k=0}^x n_k} \end{equation*}\] donde \(\hat{m}\) es el estimador de máxima verosimilitud de \(m\).

Para el caso de una modificación en cero, se tiene que

\[\begin{equation*} \begin{aligned} l_1=& \sum_{k=1}^m n_k\left\{\ln \left[\left(\begin{array}{c} m \\ k \end{array}\right) q^k(1-q)^{m-k}\right]-\ln \left[1-(1-q)^m\right]\right\} \\ =&\left(\sum_{k=1}^m k n_k\right) \ln q+\sum_{k=1}^m(m-k) n_k \ln (1-q) \\ &-\sum_{k=1}^m n_k \ln \left[1-(1-q)^m\right]+c \\ =& n x \ln q+m\left(n-n_0\right) \ln (1-q)-n \bar{x} \ln (1-q) \\ &-\left(n-n_0\right) \ln \left[1-(1-q)^m\right]+c \end{aligned} \end{equation*}\]

\[\begin{equation*} \frac{\partial l_1}{\partial q}=\frac{n x}{q}-\frac{m\left(n-n_0\right)}{1-q}+\frac{n x}{1-q}-\frac{\left(n-n_0\right) m(1-q)^{m-1}}{1-(1-q)^m} . \end{equation*}\] Igualando a cero se tiene que \[\begin{equation*} \bar{x}=\frac{1-\hat{p}_0^M}{1-p_0} m q, \end{equation*}\]

Hagamos un ejemplo para entender el proceso. Se tienen la siguiente tabla (14.5 del libro)

library(fitdistrplus)
library(actuar)

tabla_14_6 <- read.csv(file = "./data/tabla_14_6.csv")
tabla_14_6_flat <- rep(tabla_14_6$claims, tabla_14_6$policies)

hist(tabla_14_6_flat)

hist(rbinom(1000, size = 10, prob = 0.09))

binomial_loglik <- data.frame(
  size = numeric(0),
  prob = numeric(0),
  loglik = numeric(0)
)

for (s in 7:12) {
  fit_binom <- fitdist(
    data = tabla_14_6_flat,
    dist = "binom",
    fix.arg = list(size = s),
    start = list(prob = 0.3)
  )
  binomial_loglik <- rbind(
    binomial_loglik,
    data.frame(
      size = s,
      prob = fit_binom$estimate,
      fit = dbinom(
        x = s,
        size = s,
        prob = fit_binom$estimate["prob"]
      ),
      loglik = fit_binom$loglik
    )
  )
}


knitr::kable(binomial_loglik, row.names = FALSE)

size	prob	fit	loglik
7	0.1407857	1.1e-06	-19273.56
8	0.1231788	1.0e-07	-19265.37
9	0.1094940	0.0e+00	-19262.02
10	0.0985422	0.0e+00	-19260.98
11	0.0895806	0.0e+00	-19261.11
12	0.0821181	0.0e+00	-19261.84

zmbinomial_loglik <- data.frame(
  size = numeric(0),
  prob = numeric(0),
  p0 = numeric(0),
  loglik = numeric(0)
)

for (s in 7:12) {
  fit_zmbinom <- fitdist(
    data = tabla_14_6_flat,
    dist = "zmbinom",
    fix.arg = list(size = s),
    start = list(prob = 0.3, p0 = 0.1)
  )
  zmbinomial_loglik <- rbind(
    zmbinomial_loglik,
    data.frame(
      size = s,
      prob = fit_zmbinom$estimate["prob"],
      p0 = fit_zmbinom$estimate["p0"],
      loglik = fit_zmbinom$loglik
    )
  )
}


knitr::kable(zmbinomial_loglik, row.names = FALSE)

size	prob	p0	loglik
7	0.1453357	0.3539867	-19267.72
8	0.1253778	0.3540334	-19263.55
9	0.1102720	0.3538646	-19261.73
10	0.0984373	0.3540363	-19260.97
11	0.0888397	0.3541315	-19260.74
12	0.0809597	0.3539984	-19260.77

Use la siguiente tabla para ajustar distribuciones binomial, geometrica, poisson y sus contrapartes modificadas. Indique cuál sería la mejor ajustada entre esas 6.

tabla_14_7 <- read.csv("./data/tabla_14_7.csv")