###############################################################
# Script en R (no RMarkdown): Convolución, Modelos Individual/
# Colectivo, Estimación de severidad y frecuencia, VaR y Primas
# — Totalmente comentado para que puedas copiar/pegar y correr —
###############################################################
#--------------------------------------------------------------
# 0) Preparación: semillas y paquetes
#--------------------------------------------------------------
set.seed(123)
need <- c("ggplot2","dplyr","tibble","fitdistrplus","MASS")
to_install <- setdiff(need, rownames(installed.packages()))
if (length(to_install)) install.packages(to_install, quiet = TRUE)
invisible(lapply(need, library, character.only = TRUE))
##
## Adjuntando el paquete: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
## Warning: package 'fitdistrplus' was built under R version 4.4.3
## Cargando paquete requerido: MASS
## Warning: package 'MASS' was built under R version 4.4.3
##
## Adjuntando el paquete: 'MASS'
## The following object is masked from 'package:dplyr':
##
## select
## Cargando paquete requerido: survival
#--------------------------------------------------------------
# 1) Teorema de Convolución (demostración por simulación)
# Idea: Si X,Y ~ Exp(λ) iid, entonces S=X+Y ~ Gamma(k=2, rate=λ)
# f_S = f_X * f_Y (convolución); lo verificamos con hist + curva teórica
#--------------------------------------------------------------
lambda <- 0.8
n_sim_conv <- 1e5
x1 <- rexp(n_sim_conv, rate=lambda)
x2 <- rexp(n_sim_conv, rate=lambda)
s <- x1 + x2
# Comparación empírica vs teórica Gamma(k=2, rate=lambda)
hist(s, prob=TRUE, breaks=100, col="gray90", border="white",
main="Convolución: Exp(λ)+Exp(λ) ≈ Gamma(2, λ)", xlab="Suma")
curve(dgamma(x, shape=2, rate=lambda), add=TRUE, lwd=2)
# (Opcional) chequeo por integral numérica de la convolución:
f_conv_num <- function(s, lambda){
if (s <= 0) return(0)
integrand <- function(x) lambda*exp(-lambda*x) * lambda*exp(-lambda*(s-x))
integrate(integrand, lower=0, upper=s)$value
}
grid <- seq(0.05, 6, by=0.05)
num <- sapply(grid, f_conv_num, lambda=lambda)
theo <- dgamma(grid, shape=2, rate=lambda)
max_abs_err <- max(abs(num - theo))
cat("\n[Convolución] Máx. error abs (integral vs Gamma teórica):", round(max_abs_err, 6), "\n")
##
## [Convolución] Máx. error abs (integral vs Gamma teórica): 0
#--------------------------------------------------------------
# 2) Datos "juguete" para riesgo:
# - Severidad: cars$dist (positivos)
# - Frecuencia: InsectSprays$count (conteos)
# Nota: Son datasets de R usados solo con fines pedagógicos.
#--------------------------------------------------------------
sev_raw <- datasets::cars$dist
freq_raw <- datasets::InsectSprays$count
# Limpieza simple
sev <- sev_raw[is.finite(sev_raw) & sev_raw > 0]
freq <- freq_raw[is.finite(freq_raw) & freq_raw >= 0]
cat("\nResumen severidad (cars$dist):\n"); print(summary(sev))
##
## Resumen severidad (cars$dist):
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 2.00 26.00 36.00 42.98 56.00 120.00
cat("\nResumen frecuencia (InsectSprays$count):\n"); print(summary(freq))
##
## Resumen frecuencia (InsectSprays$count):
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 3.00 7.00 9.50 14.25 26.00
#--------------------------------------------------------------
# 3) Estimación de severidad (Gamma, Lognormal, Weibull) por MLE
# Selección por AIC/BIC. Guardamos la mejor para usar luego.
#--------------------------------------------------------------
fits_sev <- list(
gamma = fitdistrplus::fitdist(sev, "gamma"),
lnorm = fitdistrplus::fitdist(sev, "lnorm"),
weibull = fitdistrplus::fitdist(sev, "weibull")
)
aic_tbl <- tibble::tibble(
modelo = names(fits_sev),
AIC = sapply(fits_sev, AIC),
BIC = sapply(fits_sev, BIC)
)
aic_tbl <- dplyr::arrange(aic_tbl, AIC)
cat("\nComparación de severidad por AIC/BIC:\n"); print(aic_tbl)
##
## Comparación de severidad por AIC/BIC:
## # A tibble: 3 × 3
## modelo AIC BIC
## <chr> <dbl> <dbl>
## 1 weibull 461. 465.
## 2 gamma 463. 467.
## 3 lnorm 473. 477.
best_sev_name <- aic_tbl$modelo[1]
best_sev <- fits_sev[[best_sev_name]]
cat("\nMejor modelo de severidad por AIC:", best_sev_name, "\n")
##
## Mejor modelo de severidad por AIC: weibull
print(best_sev)
## Fitting of the distribution ' weibull ' by maximum likelihood
## Parameters:
## estimate Std. Error
## shape 1.723925 0.1908716
## scale 48.148193 4.1547159
# Momentos de la severidad estimada: E[X], Var[X]
EX <- switch(best_sev_name,
"gamma" = with(as.list(best_sev$estimate), shape/rate),
"lnorm" = with(as.list(best_sev$estimate), exp(meanlog + 0.5*sdlog^2)),
"weibull" = with(as.list(best_sev$estimate), scale*gamma(1 + 1/shape))
)
VarX <- switch(best_sev_name,
"gamma" = with(as.list(best_sev$estimate), shape/(rate^2)),
"lnorm" = with(as.list(best_sev$estimate), (exp(sdlog^2)-1)*exp(2*meanlog + sdlog^2)),
"weibull" = with(as.list(best_sev$estimate), scale^2*(gamma(1+2/shape) - (gamma(1+1/shape))^2))
)
cat("\nMomentos de severidad estimada: E[X] y Var[X]\n")
##
## Momentos de severidad estimada: E[X] y Var[X]
print(c(EX=EX, VarX=VarX, SDX=sqrt(VarX)))
## EX VarX SDX
## 42.92058 658.40777 25.65946
#--------------------------------------------------------------
# 4) Estimación de frecuencia (Poisson vs NegBin) por MLE
# Selección por AIC. Calculamos E[N] y Var[N] del modelo elegido.
#--------------------------------------------------------------
# Poisson
fit_pois <- MASS::fitdistr(freq, "Poisson")
loglik_pois <- fit_pois$loglik
k_pois <- length(fit_pois$estimate)
AIC_pois <- 2*k_pois - 2*loglik_pois
# NegBin (param: size y mu)
fit_nb <- MASS::fitdistr(freq, "Negative Binomial")
## Warning in densfun(x, parm[1], parm[2], ...): Se han producido NaNs
loglik_nb <- fit_nb$loglik
k_nb <- length(fit_nb$estimate)
AIC_nb <- 2*k_nb - 2*loglik_nb
freq_comp <- tibble::tibble(
modelo=c("poisson","negbin"),
AIC=c(AIC_pois, AIC_nb),
params = c(paste(names(fit_pois$estimate), round(fit_pois$estimate,3), collapse=", "),
paste(names(fit_nb$estimate), round(fit_nb$estimate,3), collapse=", "))
) |> dplyr::arrange(AIC)
cat("\nComparación de frecuencia por AIC:\n"); print(freq_comp)
##
## Comparación de frecuencia por AIC:
## # A tibble: 2 × 3
## modelo AIC params
## <chr> <dbl> <chr>
## 1 negbin 472. size 1.736, mu 9.5
## 2 poisson 677. lambda 9.5
if (AIC_pois <= AIC_nb){
freq_model <- "poisson"
lambda_hat <- as.numeric(fit_pois$estimate) # E[N], Var[N]=λ
EN <- lambda_hat
VarN <- lambda_hat
} else {
freq_model <- "negbin"
size_hat <- fit_nb$estimate[["size"]]
mu_hat <- fit_nb$estimate[["mu"]]
# Para NegBin parametrizada con (size, mu): E[N]=mu; Var[N]=mu + mu^2/size
EN <- mu_hat
VarN <- mu_hat + mu_hat^2/size_hat
}
cat("\nModelo de frecuencia elegido:", freq_model, "\n")
##
## Modelo de frecuencia elegido: negbin
print(c(EN=EN, VarN=VarN))
## EN VarN
## 9.50000 61.48668
#--------------------------------------------------------------
# 5) Modelo Individual: S = X1 + ... + Xm (iid)
# E[S]=m E[X], Var[S]=m Var[X]
# Estimamos también VaR por simulación.
#--------------------------------------------------------------
m <- 5L
ES_ind <- m * EX
VarS_ind <- m * VarX
SDS_ind <- sqrt(VarS_ind)
# Función generadora de severidad según el modelo elegido
rsev <- switch(best_sev_name,
"gamma" = function(n) rgamma(n, shape = best_sev$estimate["shape"], rate = best_sev$estimate["rate"]),
"lnorm" = function(n) rlnorm(n, meanlog = best_sev$estimate["meanlog"], sdlog = best_sev$estimate["sdlog"]),
"weibull" = function(n) rweibull(n, shape = best_sev$estimate["shape"], scale = best_sev$estimate["scale"])
)
# VaR por simulación
nsim_ind <- 2e5
S_ind_sim <- replicate(nsim_ind, sum(rsev(m)))
VaR_ind_95 <- unname(quantile(S_ind_sim, 0.95))
VaR_ind_99 <- unname(quantile(S_ind_sim, 0.99))
cat("\nModelo INDIVIDUAL (m=5): métricas\n")
##
## Modelo INDIVIDUAL (m=5): métricas
print(c(E_S=ES_ind, Var_S=VarS_ind, SD_S=SDS_ind, VaR_95=VaR_ind_95, VaR_99=VaR_ind_99))
## E_S Var_S SD_S VaR_95 VaR_99
## 214.60288 3292.03883 57.37629 314.95238 363.19191
#--------------------------------------------------------------
# 6) Modelo Colectivo: S = sum_{i=1}^N X_i
# Fórmulas: E[S]=E[N]E[X]
# Var[S]=E[N]Var[X] + Var[N]*(E[X])^2
# Validamos con simulación y estimamos VaR.
#--------------------------------------------------------------
ES_col <- EN * EX
VarS_col <- EN * VarX + VarN * EX^2
SDS_col <- sqrt(VarS_col)
# Simulación de N según frecuencia elegida
rfreq <- if (freq_model=="poisson") {
function(n) rpois(n, lambda=lambda_hat)
} else {
function(n) rnbinom(n, size=size_hat, mu=mu_hat)
}
simulate_aggregate <- function(nsim, rfreq, rsev){
# Devuelve nsim simulaciones de S = sum_{i=1}^N X_i
S <- numeric(nsim)
N <- rfreq(nsim)
for (i in seq_len(nsim)){
ni <- N[i]
S[i] <- if (ni==0) 0 else sum(rsev(ni))
}
S
}
nsim_col <- 2e5
S_col_sim <- simulate_aggregate(nsim_col, rfreq, rsev)
VaR_col_95 <- unname(quantile(S_col_sim, 0.95))
VaR_col_99 <- unname(quantile(S_col_sim, 0.99))
cat("\nModelo COLECTIVO: métricas teóricas y VaR simulado\n")
##
## Modelo COLECTIVO: métricas teóricas y VaR simulado
print(c(E_S=ES_col, SD_S=SDS_col, VaR_95=VaR_col_95, VaR_99=VaR_col_99))
## E_S SD_S VaR_95 VaR_99
## 407.7455 345.7226 1082.2848 1561.9209
#--------------------------------------------------------------
# 7) Primas técnicas bajo distintos principios
# - Valor esperado: Π_EV = (1+θ) E[S]
# - Varianza: Π_VAR = E[S] + β Var[S]
# - Desviación estándar: Π_SD = E[S] + k SD[S]
# - Exponencial (Esscher aprox por MC): (1/α) log E[e^{αS}]
# - VaR (cuantil): Π_VaR_q = VaR_q(S) (como referencia de cola)
#--------------------------------------------------------------
theta <- 0.30 # recargo del 30% sobre el valor esperado
beta <- 0.001 # peso para la varianza
k <- 0.8 # peso para la desviación típica
alpha <- 0.002 # aversión al riesgo (prima exponencial)
Pi_EV <- (1+theta) * ES_col
Pi_VAR <- ES_col + beta * VarS_col
Pi_SD <- ES_col + k * SDS_col
# Prima "exponencial" estable numéricamente (centrando en media)
Pi_EXP <- (1/alpha) * log(mean(exp(alpha*(S_col_sim - mean(S_col_sim))))) + mean(S_col_sim)
# Prima basada en VaR al 99% (solo en referencia)
Pi_VaR99 <- VaR_col_99
cat("\nPRIMAS (modelo COLECTIVO):\n")
##
## PRIMAS (modelo COLECTIVO):
print(
tibble::tibble(
Principio = c("Valor esperado","Varianza","Desv. estándar","Exponencial (MC)","VaR 99%"),
Parametros = c(paste0("theta=",theta),
paste0("beta=",beta),
paste0("k=",k),
paste0("alpha=",alpha),
"q=0.99"),
Prima = c(Pi_EV, Pi_VAR, Pi_SD, Pi_EXP, Pi_VaR99)
) |> dplyr::arrange(Prima)
)
## # A tibble: 5 × 3
## Principio Parametros Prima
## <chr> <chr> <dbl>
## 1 Varianza beta=0.001 527.
## 2 Valor esperado theta=0.3 530.
## 3 Exponencial (MC) alpha=0.002 595.
## 4 Desv. estándar k=0.8 684.
## 5 VaR 99% q=0.99 1562.
#--------------------------------------------------------------
# 8) Comparación rápida Individual vs Colectivo (métricas clave)
#--------------------------------------------------------------
cat("\nComparación INDIVIDUAL vs COLECTIVO:\n")
##
## Comparación INDIVIDUAL vs COLECTIVO:
print(
tibble::tibble(
Modelo = c("Individual (m=5)","Colectivo"),
E_S = c(ES_ind, ES_col),
SD_S = c(SDS_ind, SDS_col),
VaR_95 = c(VaR_ind_95, VaR_col_95),
VaR_99 = c(VaR_ind_99, VaR_col_99)
)
)
## # A tibble: 2 × 5
## Modelo E_S SD_S VaR_95 VaR_99
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 Individual (m=5) 215. 57.4 315. 363.
## 2 Colectivo 408. 346. 1082. 1562.
#--------------------------------------------------------------
# 9) Paso a paso (comentarios):
# - Ajustamos severidad (Gamma/LN/Weibull) y elegimos por AIC.
# - Ajustamos frecuencia (Poisson/NegBin) y elegimos por AIC.
# - Calculamos E[X], Var[X], E[N], Var[N].
# - Modelo COLECTIVO: E[S], Var[S]; validamos VaR por simulación.
# - Calculamos primas con distintos principios (parámetros editables).
# - Ajusta theta, beta, k, alpha según política y riesgo deseado.
#--------------------------------------------------------------