jueves, 2 de noviembre de 2017

La fórmula de Panjer 2017

Fórmula de Panjer
#Se comienza estableciendo los valores de a, b y g0=P(S=0).
#Estos valores dependen de la distribución de la frecuencia. 
#Se usará para tal fin la función abg0, donde
#dist = 'bin', 'binN' o 'poi' según la distribución de N.
#param = vector de parámetros de la distribución de N.
###########################################################
abg0 <- function(dist,param){
        if(dist=='bin'){
                p <- param[2]
                n <- param[1]        
                a <- -p/(1-p) 
                b <- -(n+1)*a
                g0 <- (1-p+p*f(0))^n
        }
        if(dist=='binN'){
                p <- param[2]
                n <- param[1]        
                a <- (1-p) 
                b <- (n-1)*a
                g0 <- (p/(1-(1-p)*f(0)))^n
        }
        if(dist=='poi'){
                a <- 0 
                b <- param
                g0 <- exp(-param*(1-f(0)))
        }
        c(a,b,g0)
}
#La función panjer calcula la probabilidad siguiente:
#    g_r = P(S=r), de forma recursiva.
# Debe estar definida previamente la función f(x)
# que corresponde a la función de probabilidad de las
# severidades.
###########################################################
panjer <- function(r,dist,param){
        info <- abg0(dist,param)
        a <- info[1]; b <- info[2]; g0 <- info[3]
        g <- g0
        C <- 1/(1-a*f(0))
        for(i in 1:r){
                aux <- 0
                for(j in 1:i){
                        aux <- aux + C*(a+(i-j+1)*b/i)*f(i-j+1)*g[j]
                }
                g[i+1] <- aux
        }
        g[r+1]
}

#Ejemplo 1
###########################################################
#Severidades uniformes discretas con parámetro n = 10
#Frecuencia geométrica con parámetro p=1/4
f <- function(x){(1/10)*(x>=1 && x<=10)} 
r <- 20
panjer(r,'binN',c(1,1/4))
## [1] 0.01557006
#Comprobación usando simulación
m <- 10000
i <- 0
uno_o_cero <- numeric(m) 
set.seed(123)
N <- rgeom(m,1/4)
while(i<m){
        i <- i + 1
        if(N[i] == 0){
                S <- 0
        }else{
                Y <- sample(1:10,N[i],replace = T)
                S <- sum(Y)
        }
        uno_o_cero[i] <- 1*(S == r)
}
mean(uno_o_cero)
## [1] 0.0159
#Ejemplo 2
###########################################################
#Severidades Poisson con parámetro lambda = 7
#Frecuencia binomial con parámetros n= 10, p=1/3
f <- function(x){dpois(x,7)} 
r <- 20
panjer(r,'bin',c(10,1/3))
## [1] 0.03408564
#Comprobación usando simulación
m <- 10000
i <- 0
uno_o_cero <- numeric(m) 
set.seed(123)
N <- rbinom(m,10,1/3)
while(i<m){
        i <- i + 1
        if(N[i] == 0){
                S <- 0
        }else{
                Y <- rpois(N[i],7)
                S <- sum(Y)
        }
        uno_o_cero[i] <- 1*(S == r)
}
mean(uno_o_cero)
## [1] 0.034
#Ejemplo 3
###########################################################
#Severidades binomiales con parámetros n= 10, p=1/3
#Frecuencia Poisson con parámetro lambda = 7
f <- function(x){dbinom(x,10,1/3)} 
r <- 20
panjer(r,'poi',7)
## [1] 0.04153745
#Comprobación usando simulación
m <- 10000
i <- 0
uno_o_cero <- numeric(m) 
set.seed(123)
N <- rpois(m,7)
while(i<m){
        i <- i + 1
        if(N[i] == 0){
                S <- 0
        }else{
                Y <- rbinom(N[i],10,1/3)
                S <- sum(Y)
        }
        uno_o_cero[i] <- 1*(S == r)
}
mean(uno_o_cero)
## [1] 0.0404

No hay comentarios:

Publicar un comentario