Fórmula recursiva de Panjer

El objetivo de la fórmula recursiva de Panjer es calcular la densidad exacta de las pérdidas agregadas \(S= \sum_{i=1}^N X_i\) del modelo de riesgo colectivo que tiene las siguientes características:

\(N\) es una v.a. discreta con distribución perteneciente a la clase \((a,b,0)\).

Recordatorio: La distribución de un v.a. discreta \(M\) pertenece a la clase \((a,b,0)\) si \[P(M=m)=\left(a+\frac{b}{m}\right)P(M=m-1)\] para \(m=1,2,3,...\) donde \(a,b\) son constantes. Se asume que \(P(M=0)>0\). Las distribuciones más conocidas que pertenecen a esta clase son la Poisson, Binomial y Binomial Negativa.
Las \(X_i\)’s son v.a.i.i.d. discretas no negativas.

Cuando esto se cumple, entonces \[f_S(s) = \frac{1}{1-af_X(0)}\sum_{x=1}^s \left(a+\frac{bx}{s}\right)f_X(x)f_S(s-x)\] para \(s=1,2,...\).

Para \(s=0\), \[f_S(0) = P(N=0)\] si \(X_i\) son estrictamente mayores a 0 y \[f_S(0) = \sum_{j=n}^\infty P(N=n)P(X=0)^n\] cuando \(X_i\) toma el valor 0 con probabilidad positiva.

Algoritmo para implementar la fórmula recursiva de Panjer

Para el caso en el que \(X\) tiene soporte finito:

Input: s (punto hasta el cual evaluar), a,b,g0\(=P(S=0)\),fx (vector con los valores de la densidad de X)
Inicio:
1. Crear el vector g para almacenar los valores de \(f_S(i), i=0,1,...,s\).
2. Asignar g[1]=g0.
3. Hacer i=1.
4. Mientras i<=1:
- g[i+1]=0.
- Para j desde 1 hasta el valor mínimo entre el valor máximo que puede tomar \(X\) y s:
  - g[i+1]=g[i+1] + (a+b*j/i)* fx[j+1]*g[i-j+1].
- Aumentar i en 1.
Fin. Regresar el vector g.

Implementación

#Fórmula recursiva de Panjer para X con soporte finito
panjer <- function(s,a,b,g0,fx){
  g <- rep(NA,s+1)
  g[1] <- g0
  i=1
  while(i<=s) {
    g[i+1] <- 0
    for (j in 1:min(i,length(fx)-1)) {
      g[i+1] <- g[i+1] + (1/(1-a*fx[1]))*(a+b*j/i)* fx[j+1]*g[i+1-j]
    }
  i=i+1
  }
  return(g)
}

Ejemplo 1:

La frecuencia de reclamaciones \(N\) sigue una distribución Poisson con media 100, mientras que el monto de cada reclamación solo toma los valores enteros entre 0 y 10 con probabilidad dada por la siguiente tabla:

x	fx	x	fx
0	0.2212	6	0.0252
1	0.3064	7	0.0153
2	0.1859	8	0.0093
3	0.1127	9	0.0056
4	0.0684	10	0.0087
5	0.0415	.	.

a <- 0
b <- 100
g0 <- exp(-77.88)
fx <- c(0.2212, 0.3064,0.1859,0.1127,0.0684,0.0415, 0.0252,0.0153,0.0093,0.0056,0.0087)
s <- 250

gs <- panjer(s,a,b,g0,fx)
gs[1:10]

 [1] 1.503647e-34 4.607173e-33 7.337717e-32 8.082183e-31 6.912964e-30
 [6] 4.889456e-29 2.974308e-28 1.598373e-27 7.736604e-27 3.422527e-26

Para el caso en el que \(X\) tiene soporte infinito:

En este caso, solo cambiamos el vector de probabilidades de \(X\) (fx) por su función de densidad con soporte infinito.

panjer_f <- function(s,a,b,g0,fX){
  g <- rep(NA,s+1)
  g[1] <- g0
  i=1
  while(i<=s) {
    g[i+1] <- 0
    for (j in 1:i) {
      g[i+1] <- g[i+1] + (1/(1-a*fX(0)))*(a+b*j/i)* fX(j)*g[i+1-j]
    }
    i=i+1
  }
  return(g)
}

Ejemplo 2:

La frecuencia de reclamaciones tiene distribución Poisson con media 2, mientras que el monto de cada reclamación puede tomar cualquier valor entero positivo con probabilidad \(f_X(x) = (0.6)(0.4^{x-1}), x=1,2,...\).

b2 <- 2
(g02 <- dpois(0,2))

[1] 0.1353353

f <- function(x){
  if(x>0){
    return(0.6*(0.4^(x-1))) 
  }else{
    return(0)
  }
}
s <- 20

gs2 <- panjer_f(s,a,b2,g02,f)
gs2[1:10]

 [1] 0.13533528 0.16240234 0.16240234 0.14291406 0.11563047 0.08803506
 [7] 0.06396314 0.04476440 0.03037416 0.02007861

Fórmula recursiva de De Pril

Esta fórmula sirve para calcular la función densidad de las pérdidas agregadas \(S=\sum_{j=1}^n Y_j\) en un modelo de riesgo individual que cumpla con las siguientes características:

El portafolio contiene \(J\) grupos de riesgo determinados por la probabilidad de reclamación de la póliza. El grupo j tiene probabilidad de reclamación \(\theta_j\) para \(j=1,2,...,J\).
Cada póliza puede tener una suma asegurada de \(i\in \{1,2,...,I\}\) unidades monetarias.
Tenemos \(n_{ij}\) pólizas independientes con suma asegurada (SA) \(i\) y probabilidad de reclamo \(\theta_j\), para \(i=1,2,...I\) y \(j=1,2,...,J\).
Denotamos como \(n=\sum_{i=1}^I\sum_{j=1}^J i\cdot n_{ij}\) al monto máximo de las pérdidas agregadas.

En consecuencia la función de densidad de \(S\) para este modelo, denotada como \(f_S(s)\), tiene soporte en \(\{1,2,...,n\}\).

Para este modelo de riesgo individual, se puede calcular \(f_S(s)\) como:

\[f_S(s)=\frac{1}{s}\sum_{i=1}^{min(s,J)}\sum_{k=1}^{[s/i]}f_S(s-ik)h(i,k)\] para \(s=1,2,...,n\) con \[h(i,k)=i(-1)^{k-1}\sum_{j=1}^J n_{ij}\left(\frac{\theta_j}{1-\theta_j}\right)^k\] para \(k=1,2,...,I\) y \(h(i,j)=0\) en otro caso.

El valor inicial para esta fórmula recursiva es \[f_S(0)=\prod_{i=1}^I\prod_{j=1}^J (1-\theta_j)^{n_{ij}}.\]

Implementación

La función recibe s (punto hasta el cual evaluar), n (matriz de dimensión \(I\times J\) con los valores \(n_ij\)) y q (vector de probabilidades \(\theta_j\) con longitud \(J\)). Regresa el vector fs que contiene los valores de \(f_S(r), r=1,...,s\).

#Función auxiliar
h <- function(i,k,n){
  I <- nrow(n)
  if(i>=1 && i<=I){
    return(i*(-1)^(k-1)*sum(n[i,]*(q/(1-q))^k))
  }else{
    return(0)
  }
}
#Fórmula de De Pril
de_pril <- function(s,n,q){
  I <- nrow(n)
  J <- length(q)
  fs <- rep(NA,s+1)
  fs[1] <- 1
  for (i in 1:I) {
    fs[1] <- fs[1]*prod((1-q)^n[i,])
  }
  r=1
  while(r<=s){
    fs[r+1] <- 0
    for (i in 1:min(r,I)) {
      for (k in 1:floor(r/i)) {
        fs[r+1] <- fs[r+1] + fs[r+1-i*k]*h(i,k,n)
      }
    }
    fs[r+1] <- fs[r+1]/r
    r=r+1
  }
  return(fs)
}

Ejemplo 3:

Un portafolio contiene pólizas separadas en grupos según su suma asegurada que puede ser de 1, 2 o 3 unidades monetarias y su probabilidad de reclamación que puede ser 0.0015, 0.0024, 0.0031, o 0.0088. El número de pólizas de cada grupo y su correspondiente SA y probabilidad de reclamación se muestra en la siguiente tabla:

	j=1	j=2	j=3	j=4
SA=1	30	40	50	60
SA=2	40	50	70	80
SA=3	70	60	80	90

[1] 0.0015 0.0024 0.0031 0.0088

s<-30
fs1 <- de_pril(s,n,q)
#Función de densidad en s=1,...,10
fs1[1:10]

 [1] 0.03977521 0.03299227 0.05765641 0.09185897 0.08287425 0.09915229
 [7] 0.10352584 0.08935326 0.08578208 0.07421027

Ejemplo 4:

La siguiente tabla muestra el número de pólizas de seguro de vida agrupadas según la suma asegurada y la tasa de mortalidad del correspondiente grupo de edad.

#La matriz n contiene ceros en los renglones correspondientes 
# a las sumas aseguradas no definidas.
n

      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
 [1,]  400    0    0    0    0    0    0    0    0     0
 [2,]    0  400    0    0    0    0    0    0    0     0
 [3,]    0    0    0    0    0    0    0    0    0     0
 [4,]    0    0  400    0    0    0    0    0    0     0
 [5,]    0    0    0    0    0    0    0    0    0     0
 [6,]    0    0    0  400    0    0    0    0    0     0
 [7,]    0    0    0    0    0    0    0    0    0     0
 [8,]    0    0    0    0  400    0    0    0    0     0
 [9,]    0    0    0    0    0    0    0    0    0     0
[10,]    0    0    0    0    0  400    0    0    0     0
[11,]    0    0    0    0    0    0  400    0    0     0
[12,]    0    0    0    0    0    0    0  400    0     0
[13,]    0    0    0    0    0    0    0    0    0     0
[14,]    0    0    0    0    0    0    0    0  600     0
[15,]    0    0    0    0    0    0    0    0    0   600

s<-250

f <- de_pril(s,n,q)
#Función de densidad en s=1,...,10
f[1:10]

 [1] 1.800951e-07 4.395768e-07 9.269649e-07 1.389615e-06 2.200847e-06
 [6] 2.958986e-06 4.298655e-06 5.515409e-06 7.611877e-06 9.455523e-06

#Cálculo de la función de distribución.
x <- seq(26,251,by=25)
Fx <- sapply(x, function(a){sum(f[1:a])})

Método de Kornya

Este método provee una aproximación de la función de densidad de las pérdidas agregadas \(S=\sum_{j=1}^n Y_j\) para un modelo de riesgo individual con las mismas características que se requieren para usar la fórmula de De Pril. Usando la notación de la sección anterior, la aproximación de Kornya es: \[g_S^{(K)}(s)=\frac{1}{s}\sum_{j=1}^{min(x,IK)}jb_j^{(K)}g_S^{(K)}(s-j)\] donde \[b_x^{(K)} = \sum_{k=\{x/I\},k|x}^{min(K,x)}\frac{(-1)^{(k+1)}}{k}\sum_{j=1}^J n_{x/k,j}\left(\frac{\theta_j}{1-\theta_j}\right)^k\], \[g_S^{(K)}(0) = \exp{(b_0^{(K)})}\] y \[b_0^{(K)} = \sum_{k=1}^K\frac{(-1)^{k}}{k}\sum_{i=1}^I \sum_{j=1}^J n_{ij}\left(\frac{\theta_j}{1-\theta_j}\right)^k\] Mientras mayor sea \(K\) mejor será la aproximación de \(g_S^{(K)}(s)\) a \(f_S(s)\), sin embargo, con \(K=4\) ya suele tenerse una aproximación suficientemente buena.

Implementación

La función kornya() recibe el punto s hasta el cual se calcula la aproximación, el grado K de aproximación, la matrix n con los grupos de pólizas y el vector de probabilidades de reclamación q; estos dos últimos definidos de igual manera que para la función de de_pril.

#Función auxiliar
b <- function(x,K,I,J){
  if(x==0){
    bx <- 0
    for (k in 1:K) {
      for (i in 1:I) {
        bx <- bx + ((-1)^k/k)*sum(n[i,]*(q/(1-q))^k)
      }
    }
  }else{
    k <- seq(floor(x/I) + 1,min(x,K))
    k <- subset(k,x%%k==0)
    bx <- 0
    for(i in k){
      bx <- bx + (((-1)^(i+1))/i)*sum(n[x/i,]*(q/(1-q))^i)
    } 
  }
 return(bx) 
}
# Fórmula de Kornya
kornya <- function(s,K,n,q){
  I <- nrow(n)
  J <- ncol(n)
  d <- rep(NA,s+1)
  d[1] <- exp(b(0,K,I,J))
  m=1
  while (m<=s) {
    d[m+1] <- 0
    for (j in 1:min(m,I*K)) {
      d[m+1] <- d[m+1] + j*b(j,K,I,J)*d[m+1-j]
    }
    d[m+1] <- d[m+1]/m
    m=m+1
  }
  return(d)
}

Ejemplo 5

Usando el mismo plantemiento del Ejemplo 4.

s<-250
#Con K=2
kor <- kornya(s,2,n,q)
x <- seq(26,251,by=25)
Fx_K2 <- sapply(x, function(a){sum(kor[1:a])})
#Aproximación de la función de densidad en s=1,...,10
kor[1:10]

 [1] 1.801129e-07 4.396201e-07 9.270562e-07 1.389747e-06 2.201050e-06
 [6] 2.959250e-06 4.299033e-06 5.515876e-06 7.612518e-06 9.456294e-06

Fórmulas recursivas para la densidad de las pérdidas agregadas

Fórmula recursiva de Panjer

Algoritmo para implementar la fórmula recursiva de Panjer

Implementación

Ejemplo 1:

Ejemplo 2:

Fórmula recursiva de De Pril

Implementación

Ejemplo 3:

Ejemplo 4:

Método de Kornya

Implementación

Ejemplo 5