Graficando el Movimiento Browniano (caminata aleatoria aproximante)

El siguiente código genera trayectorias como las de las imágenes

MB <- function(m,n,t){

        #m <- 5; n <- 20; t <- 1

        pasos <- 2*rbinom(n*m,1,0.5)-1

        r <- matrix(pasos, nrow=n, ncol=m)

        delta <- t/n

        x <- (0:n)*delta ##Valores del eje tiempo.

        B <- matrix(rep(0,(n+1)*m), nrow=n+1, ncol=m)

        for(k in 1:m){

                B[1,k] <- 0

                for(j in 1:n){

                        aux <- r[j,k]*sqrt(delta)

                        B[j+1,k] <- B[j,k]+aux

                }

        }

        par(bg='lightyellow',mar=c(2,4,2,1)) 

        matplot(x,B, type="l",lwd=2,lty = 1,

                main='Movimiento Browniano Est?ndar',

                ylab='Espacio de estados',xlab='')

        grid(5)

        abline(v=0,h=0)

}

##########################

#Ejemplo##################

##########################

par(mfrow=c(3,1),oma=c(1,0,1,0))

MB(20,10,5)

MB(20,100,5)

MB(20,1000,5)

par(mfrow=c(1,1))

MB(50,1000,5)

MB(100,1000,5)

set.seed(321)

MB(100,1000,6)

#################################

x <- seq(-6,6,0.01)

y <- dnorm(x)

lines(1-y,x,lwd=10,col='orange')

lim_sup <- qnorm(0.975)

lim_sup

abline(h=c(lim_sup,-lim_sup),lwd=2)

#################################

y <- dnorm(x,mean=0,sd=sqrt(1/2))

lines(1/2-y,x,lwd=10,col=2)

lim_sup <- qnorm(0.975,0,sqrt(1/2))

lim_sup

abline(h=c(lim_sup,-lim_sup),col=2,lwd=2)

#################################

y <- dnorm(x,mean=0,sd=sqrt(2))

lines(2-y,x,lwd=10,col=3)

lim_sup <- qnorm(0.975,0,sqrt(2))

lim_sup

abline(h=c(lim_sup,-lim_sup),col=3,lwd=2)

#################################

y <- dnorm(x,mean=0,sd=sqrt(3))

lines(3-y,x,lwd=10,col=4)

lim_sup <- qnorm(0.975,0,sqrt(3))

lim_sup

abline(h=c(lim_sup,-lim_sup),col=4,lwd=2)

#################################

y <- dnorm(x,mean=0,sd=sqrt(4))

lines(4-y,x,lwd=10,col=5)

lim_sup <- qnorm(0.975,0,sqrt(4))

lim_sup

abline(h=c(lim_sup,-lim_sup),col=5,lwd=2)

#################################

y <- dnorm(x,mean=0,sd=sqrt(5))

lines(5-y,x,lwd=10,col=6)

lim_sup <- qnorm(0.975,0,sqrt(5))

lim_sup

abline(h=c(lim_sup,-lim_sup),col=6,lwd=2)

#################################

y <- dnorm(x,mean=0,sd=sqrt(6))

lines(6-y,x,lwd=10,col=1)

lim_sup <- qnorm(0.975,0,sqrt(6))

lim_sup

abline(h=c(lim_sup,-lim_sup),col=1,lwd=2)

#########################################

#########################################

set.seed(321)

m <- 100

n <- 1000

t <- 6

pasos <- 2*rbinom(n*m,1,0.5)-1

r <- matrix(pasos, nrow=n, ncol=m)

delta <- t/n

x <- (0:n)*delta ##Valores del eje tiempo.

B <- matrix(rep(0,(n+1)*m), nrow=n+1, ncol=m)

for(k in 1:m){

        B[1,k] <- 0

        for(j in 1:n){

                aux <- r[j,k]*sqrt(delta)

                B[j+1,k] <- B[j,k]+aux

        }

}

par(bg='black',mar=c(2,4,2,1)) 

matplot(x,B, type="l",lwd=2,lty = 1,

        col=heat.colors(m),

        col.main=7,

        main='Movimiento Browniano Est?ndar',

        ylab='Espacio de estados',xlab='')

grid(6,col='green')

abline(v=0,h=0,col='white')

La fórmula de Panjer 2017

Fórmula de Panjer

David J. Santana

2 de noviembre de 2017

#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