9  Métodos de simulación

Podríamos definir la Simulación como una técnica que consiste en realizar experimentos de muestreo sobre el modelo de un sistema, con el objetivo de recopilar información bajo determinadas condiciones.

9.1 Pasos de una Simulación

  1. Construya un modelo para \(S\) que dependa de variables aleatorias \(X, Y, Z, \ldots\), donde se conozcan sus distribuciones y cualquier dependencia.

  2. Para \(j = 1, \ldots, n\) genera valores pseudoaleatorios \(x_{j}, y_{j}, z_{j}, \ldots\) y luego calcula \(s_{j}\) usando el modelo del paso 1

  3. El CDF de \(S\) puede ser aproximado por \(F_{n} (s)\), el CDF empírico basado en la muestra pseudoaleatoria \(s_{1}, \ldots, s_{n}\)

  4. Calcule las cantidades de interés, como la media, la varianza, los percentiles o las probabilidades, utilizando la función de distribución de datos empírica.

9.2 Generación de simulaciones

Método de inversión Un hecho que facilita la provisión de tales secuencias es que es suficiente poder generarlas para la distribución uniforme en el intervalo \((0,1)\). Esto se debe a que, si \(U\) tiene la distribución \((0,1)\) uniforme, entonces \(X = F_{X} ^ {- 1} (U)\) cuando exista la inversa tendrá \(F_{X} (x )\) como su cdf. Por lo tanto, simplemente obtenemos números pseudoaleatorios uniformes \(u_ {1} ^ {* *}, \ldots, u_ {n} ^ {* *}\) y luego dejamos \(x_{j} ^ {* *} = F_ {X} ^ {- 1} \left (u_{j} ^ {* *} \right).\)

Tamaño de la muestra Sabemos que cualquier estimulador consistente se acercará arbitrariamente al valor real con una alta probabilidad de que el tamaño de la muestra está incrementado.

El Teorema del Límite Central nos ayudara a conocer un valor apropiado de muestras para un determinado nivel de precisión.

Hagamos un caso sencillo con la pareto.

dist_pareto <- function(x) {
  1 - 1000 / ((1000 + x))^3
}

dist_pareto_inv <- function(u) {
  1000 * ((1 - u)^(-1 / 3) - 1)
}

u <- runif(1000)

plot(density(dist_pareto_inv(u)))

9.3 Casos especiales

9.3.1 Mezcla de distribuciones

La mezcla para una distribución discreta es \[\begin{equation*} F_Y(y)=a_1 F_{X_1}(y)+a_2 F_{X_2}(y) +\cdots+a_k F_{X_k}(y) \text {. } \end{equation*}\]

Puede ser difícil invertir esta función, pero puede ser fácil invertir 1. Simule un valor de la variable aleatoria discreta \(J\) donde \(\operatorname{Pr}(J=j)=a_j\). 2. Usar un método apropiado (generalmente inversión) para simular una observación de una variable aleatoria con función de distribución \(F_{X_J}(y)\).

Suponga que se tiene la siguiente mixtura \(F(x)=0.3\left(1-e^{-0.02 x}\right)+\) \(0.5\left(1-e^{-0.04 x}\right)+0.2\left(1-e^{-0.05 x}\right)\). Genere una secuencia aleatoria de esta.

La solución viene dato por lo siguiente

mixtura_exponenciales <- function(u) {
  if (u < 0.3) {
    qexp(p = runif(1), rate = 0.02)
  } else if (u < 0.8) {
    qexp(p = runif(1), rate = 0.04)
  } else {
    qexp(p = runif(1), rate = 0.05)
  }
}

u <- runif(1000)

plot(density(sapply(X = u, FUN = mixtura_exponenciales)))

9.3.2 Tiempo hasta fallo

En ocasiones hay distribuciones complejas de modelar, sobretodo cuando hay múltiples decrementos.

Se tiene un portafolio con 250 personas y las siguientes probabilidades en cierto momento

  1. Muerte: 0.01.
  2. Renuncia: 0.09.
  3. Invalidez: 0.03.
  4. Sano: 0.87.

Primero modele el número de muertes

n <- 1000
muertes <- qbinom(p = runif(n), size = 250, prob = 0.01)
renuncias <- qbinom(
  p = runif(n), size = 250 - muertes,
  prob = 0.09 / (1 - 0.01)
)
invalidez <- qbinom(
  p = runif(n),
  size = 250 - muertes - renuncias,
  prob = 0.03 / (1 - 0.01 - 0.09)
)

sanos <- 250 - muertes - renuncias - invalidez

plot(density(sanos))

9.4 Modulación de pérdidas agregadas

Los modelos recursivos que se vieron anteriormente para modelarse necesita 3 características fundamentales:

  1. El método es exacto según la aproximación de \(f_X(x)\). Esto se puede resolver agregando más puntos.
  2. Se asume la forma \(S = X_1 + \cdots + X_N\) con los \(X\)’s independientes e idénticamente distribuidos. Esto puede ser un problema para casos un poco más complejos.
  3. Se asume que la frecuencia es de clase \((a,b,0)\) o \((a,b,1)\). Si el patrón de frecuencias es diferente, entonces una distribución compuesta puede funcionar, pero si se necesita simulación para obtener la distribución.

Ejemplo 9.1 Suponga que hay un deducible, \(d\), sobre pérdidas individuales. Sin embargo, en el transcurso de un año, el titular de la póliza no pagará más de \(u\).

En este caso se tiene

  • \(X_j\) la cantidad de \(j\)-ésima pérdida. Los datos son i.i.d.
  • Defina \(W_j=X_j \wedge d\) la cantidad pagada por el asegurado debido al deducible y sea \(Y_j=X_j-W_j\) la cantidad pagada por la aseguradora.
  • En este caso \(R=W_1+\cdots+W_N\) es el monto total pagado por el titular de la póliza antes de imponer el desembolso máximo.
  • Entonces el monto realmente pagado por el asegurado \(R_u=R \wedge u\).
  • Sean \(S=X_1+\cdots+X_N\) las pérdidas totales, y el monto total pagado por el asegurador es \(T=S-R_u\).
  • Tenga en cuenta que las distribuciones de \(S\) y \(R_u\) se basan en distribuciones i.i.d. de severidad.
Advertencia

Los métodos analíticos descritos anteriormente se pueden utilizar para obtener sus distribuciones. Pero debido a que son dependientes, sus distribuciones individuales no pueden combinarse para producir la distribución de \(T\). Tampoco hay forma de escribir \(T\) como una suma aleatoria de i.i.d. variables Al comienzo del año, parece que \(T\) será la suma de i.i.d. \(Y_j\)’s, pero en algún momento los \(Y_j\) s pueden ser reemplazados por \(X_j\) s a medida que se alcanza el desembolso máximo.

Suponga que para el caso anterior se tiene que el deducible es \(d=250\) el máximo que pagaría el asegurado por año es \(u=1000\). Asuma que la frecuencia es binomial negativa con \(r=3\) y \(\beta=2\). Además las pérdidas individuales son Weibull \(\tau=2\) y \(\theta=600\).

  1. Simule la frecuencia con una binomial negativa. Use 1000 valores por ejemplo.
  2. Con los valores obtenidos simule las Weibull correspondientes. 3 En este paso considere cuánto paga realmente el asegurado de acuerdo al deducible.
  3. Considere el límite máximo de suma para determinar cuál sería el monto final a pagar por el asegurado.
  4. Con este último vector calcule el cuantil 95%. Ese el VaR.
  5. Del vector de sumas agregadas, filtre todas las pérdidas mayores al VaR. Tome el promedio de ese vector. Ese valor es el CVaR o TVaR.