7 Aplicación práctica de modelos de frecuencia y severidad

set.seed(123)
S <- rcompound(
  n = 500,
  model.freq = rpois(lambda = 1.5),
  model.sev = rpareto(shape = 7, scale = 2000),
  SIMPLIFY = FALSE
)

hist(S$severity)

hist(S$frequency)

hist(S$aggregate)

severity <- S$severity
severity[severity <= 100] <- NA
severity[severity >= 1000] <- 1000
severity <- na.omit(severity)
hist(severity)

frequency <- S$frequency

delta <- severity != 1000
surv <- Surv(
  time = rep(100, length(severity)),
  time2 = severity,
  event = delta
)
head(surv)

[1] (100,320.2094] (100,308.4212] (100,390.3076] (100,869.1446] (100,309.3110]
[6] (100,559.2225]

descdist(as.numeric(severity))

summary statistics
------
min:  100.4654   max:  1000 
median:  311.0787 
mean:  416.6421 
estimated sd:  281.5671 
estimated skewness:  0.9175211 
estimated kurtosis:  2.572411

custom_pareto <- list(
  name = "pareto",
  pars = c("shape", "scale"),
  location = c("scale"),
  transforms = c(log, log),
  inv.transforms = c(exp, exp)
)
fit_severity <- flexsurvreg(
  formula = surv ~ 1,
  dist = custom_pareto,
  inits = c(1, 1)
)

Forming integrated rmst function...

Forming integrated mean function...

plot(fit_severity)

fit_severity <- flexsurvreg(
  formula = surv ~ 1,
  dist = "llogis"
)
plot(fit_severity)

descdist(frequency, discrete = TRUE)

summary statistics
------
min:  0   max:  7 
median:  1 
mean:  1.5 
estimated sd:  1.243017 
estimated skewness:  0.9334414 
estimated kurtosis:  3.929327

fit_frequency <- fitdist(data = frequency, dist = "pois")

x <- seq.int(0, 900, 10)
severity_discrete <- discretize(
  cdf = pllogis(x,
    shape = fit_severity$coefficients["shape"],
    scale = fit_severity$coefficients["scale"]
  ),
  from = 0,
  to = 900,
  by = 10,
  method = "rounding"
)
severity_discrete

 [1] 0.4838405802 0.1551435104 0.0650334216 0.0388972701 0.0267776369
 [6] 0.0199341817 0.0156063898 0.0126575627 0.0105386768 0.0089541702
[11] 0.0077318154 0.0067649801 0.0059843953 0.0053432702 0.0048089562
[16] 0.0043580400 0.0039733329 0.0036419570 0.0033540917 0.0031021290
[21] 0.0028800926 0.0026832292 0.0025077149 0.0023504427 0.0022088642
[26] 0.0020808715 0.0019647078 0.0018588979 0.0017621952 0.0016735398
[31] 0.0015920253 0.0015168725 0.0014474083 0.0013830483 0.0013232834
[36] 0.0012676677 0.0012158095 0.0011673635 0.0011220239 0.0010795195
[41] 0.0010396087 0.0010020756 0.0009667271 0.0009333896 0.0009019068
[46] 0.0008721380 0.0008439554 0.0008172437 0.0007918978 0.0007678221
[51] 0.0007449297 0.0007231409 0.0007023828 0.0006825886 0.0006636971
[56] 0.0006456519 0.0006284009 0.0006118962 0.0005960936 0.0005809520
[61] 0.0005664337 0.0005525032 0.0005391281 0.0005262778 0.0005139241
[66] 0.0005020406 0.0004906028 0.0004795877 0.0004689738 0.0004587410
[71] 0.0004488706 0.0004393449 0.0004301472 0.0004212622 0.0004126752
[76] 0.0004043723 0.0003963408 0.0003885683 0.0003810435 0.0003737553
[81] 0.0003666937 0.0003598488 0.0003532116 0.0003467733 0.0003405257
[86] 0.0003344611 0.0003285721 0.0003228516 0.0003172931 0.0003118902

frequency_discrete <- dpois(x = 1:10, lambda = fit_frequency$estimate)
frequency_discrete

 [1] 3.346952e-01 2.510214e-01 1.255107e-01 4.706652e-02 1.411996e-02
 [6] 3.529989e-03 7.564262e-04 1.418299e-04 2.363832e-05 3.545748e-06

fX <- severity_discrete
pn <- frequency_discrete

ncolumnas <- length(pn) - 1
nfilas <- 50
tabla_convolucion <- matrix(0, nrow = nfilas + 1, ncol = ncolumnas + 1)

tabla_convolucion[1, 1] <- 1

for (nc in 1:ncolumnas) {
  for (nf in 1:nfilas) {
    s <- 0
    for (k in seq_len(min(nf, length(fX)))) {
      s <- s + tabla_convolucion[nf - k + 1, nc] * fX[k]
    }
    tabla_convolucion[nf + 1, nc + 1] <- s
  }
}

pn_tabla <- kronecker(matrix(1, nrow = nfilas + 1), t(pn))

fS <- rowSums(pn_tabla * tabla_convolucion)

plot(0:nfilas, fS, type = "l", xlab = "payment", ylab = "aggregate density")

lambda <- fit_frequency$estimate["lambda"]
fS_recursive <- exp(-lambda)

for (x in 1:nfilas) {
  s <- 0
  for (y in 1:20) {
    if (x - y >= 0) {
      s <- s + y * fX[y] * fS_recursive[x - y + 1]
    }
  }
  fS_recursive[x + 1] <-  lambda / x * s
}

plot(fS_recursive, type = "l", xlab = "payment", ylab = "aggregate density")

Usando la funciones de recursión construidas en capítulos anteriores, trate de encontrar la función de distribución de la S.

Además revise la función aggregateDist del paquete actuar y trate de reconstruir la S.

FS_convolution <- aggregateDist(
  method = "convolution",
  model.freq = frequency_discrete,
  model.sev = severity_discrete
)

plot(FS_convolution, xlim = c(0, 50))

plot(diff(FS_convolution), xlim = c(0, 50), type = "l")

FS_recursive <- aggregateDist(
  method = "recursive",
  model.freq = "poisson",
  lambda = fit_frequency$estimate["lambda"],
  model.sev = severity_discrete
)

Warning in panjer(fx = model.sev, dist = dist, p0 = p0, x.scale = x.scale, :
maximum number of recursions reached before the probability distribution was
complete

plot(FS_recursive, xlim = c(0, 50))

plot(diff(FS_recursive), xlim = c(0, 50), type = "l")

Use la siguiente tabla para encontrar la distribución de la S.

seguros <- read.table("http://www.karlin.mff.cuni.cz/~pesta/NMFM402/moped.txt",
header=T)

Aquí nos interesa las variables

severity: Monto de reclamos (en moneda sueca SEK)
number: Número de reclamos.