Proyecto de Simulación
Procesos de Renovación
1. Código de una función que simule trayectorias de un proceso de renovación.
# renewal: función que modela un proceso de renovación
# tf = tiempo final; a,b = parámetros de Ti's gamma(a,b)
renewal <- function(tf, a, b){
Ti <- c()
Wn <- c()
W <- 0
aux<-0
repeat{
Ti[aux+1] <- rgamma(1, a, b)
W <- sum(Ti)
if(W >= tf){
break
}
aux <- aux + 1
}
for(i in 1:length(Ti)){
Wn[i] <- sum(Ti[1:i])
}
T_inter <- Ti[-length(Ti)]
# La función imprimirá lo gráfica de la trayectoria apartir del sig. data frame
df_renov <- data.frame(x = c(0,cumsum(T_inter)), y = 0:length(T_inter), xend = c(cumsum(T_inter),tf),
yend = 0:length(T_inter))
graf <- ggplot(df_renov, aes(x = x, y = y, xend = xend, yend = yend)) + geom_segment(color = "red")+
ggtitle("Trayectoria de un proceso de Renovación")+
xlab("Tiempo")+
ylab("Número de arribos")
return(list(Nt = length(Wn[-length(Wn)]), Ti_s = Ti, Wn_s = Wn[-length(Wn)], graf))
}La función regresa una lista con el número de llegadas (Nt), los tiempos de interarribo (Ti), los tiempos de arribo (Wn) y la gráfica de la trayectoria.
Prueba
Tomando los siguientes parametros: tf = 1000, a = 25, b = 5. Verificaremos que se cumpla el teorema elemental de renovación y se graficará la trayectoria
Sea:
\(c_1 = \frac{1}{\mu} = \frac{1}{(\frac{a}{b})}\)
\(d_1 = \frac{Nt}{tf}\)
## [1] 0.2## [1] 0.201Podemos ver que sí se cumple el teorema elemental de renovación. Finalmente se muestra la gráfica de la trayectoria.
2. Código de una función que simule la trayectoria de un proceso Poisson compuesto.
# poic: función que modela un proceso Poisson compuesto
# tf = tiempo final; l = intensidad del poisson; a,b = parámetros de Yn's gamma(a,b)
poic <- function(tf, l, a, b){
# Como en un p. Poisson las Ti's, son exp(l) podemos usar la función definida
# en 1) estableciendo "a = 1" pues gamma(1,l) = exp(l)
Yn <- c()
p <- renewal(tf, 1, l)
for(i in 1:p$Nt){
Yn[i] = rgamma(1, a, b)
}
df_renov <- data.frame(x = c(1:p$Nt), y = c(cumsum(Yn)))
graf <- ggplot(df_renov, aes(x = x , y = y, xend = c(1:p$Nt), yend = c(cumsum(Yn)))) +
geom_step(color = "red")+
ggtitle("Trayectoria de un proceso Poisson Compuesto")+
xlab("Número de saltos")+
ylab("Valor")
return(list(numero_saltos = p$Nt, tiempo_saltos = p$Wn_s, tamaño_salto = Yn, Xt = cumsum(Yn), graf ))
}La función regresa una lista con el número de saltos (numero_saltos), los tiempos de salto (tiempo_saltos), el tamaño de cada salto (tamaño_salto), el valor del proceso después de cada salto (Xt), y la gráfica de la trayectoria.
Prueba
Tomando los siguientes parametros: tf = 1000, l = 2, a = 10, b = 2. Verificaremos que se cumpla el teorema elemental de renovación y se graficará la trayectoria
Sea:
\(c_2 = \frac{1}{\frac{1}{\lambda}} = \frac{1}{\frac{1}{2}}\)
\(d_2 = \frac{Nt}{tf}\)
## [1] 2## [1] 1.936Podemos ver que sí se cumple el teorema elemental de renovación. Finalmente se muestra la gráfica de la trayectoria.
3. Código de una función que simula la trayectoria del proceso de Cramer Lundberg.
Comenzamos creando una función que modele una trayectoria de CL.
# CL: función que modela un proceso de Cramer Lundberg
# tf = tiempo final; u = capital inicial; c = prima; l = intensidad del poisson;
# a,b = parámetros de Yn's gamma(a,b)
CL <- function(tf, u, c, l, a, b){
p <- poic(tf, l, a, b)
Rt <- c()
for (i in 1:p$numero_saltos){
Rt[i] <- u + c*p$tiempo_saltos[i] - sum(p$tamaño_salto[1:i])
}
df_renov <- data.frame(x = c(0,p$tiempo_saltos), y = c(u,Rt))
graf <- ggplot(df_renov, aes(x = x , y = y, xend = c(0, p$tiempo_saltos), yend = c(u, Rt))) +
geom_line(color = "red")+
ggtitle("Trayectoria de un Proceso de Cramer Lundberg")+
xlab("Tiempo")+
ylab("Valor")
return(list(numero_saltos = p$numero_saltos, tiempo_saltos = p$tiempo_saltos,
tamaño_salto = p$tamaño_salto, trayectoria = Rt, graf ))
}La función regresa una lista con el número de saltos (numero_saltos), los tiempos de salto (tiempo_saltos), el tamaño de cada salto (tamaño_salto), el valor del proceso después de cada salto (Rt), y la gráfica de la trayectoria.
Una vez que ya tenemos esta función creamos otra función que simula N trayectorias de un proceso de CL con el fin de calcular la probabilidad de ruina.
# SCL: función que simula N procesos de Cramer Lundberg para mostrar la probablidad de ruina
# N = número de simulaciones; tf = tiempo final; u = capital inicial;
# c = prima; l = intensidad del poisson; a,b = parámetros de Yn's gamma(a,b)
SCL <- function(N, tf, u, c, l, a, b){
vt <- c()
for (i in 1:N){
vt[i] <- sum(CL(tf, u, c, l, a, b)$trayectoria < 0)
}
return(sum(vt > 0)/N)
}Con ambas funciones creadas pasamos a la parte de la comprobación.
Prueba
Dado que nos piden que la las Yn sean exp(\(\alpha\)) escribimos “a = 1” y “b = \(\alpha\)” y tomando los siguientes parametros: tf = 100, u = 0, c = 5, l = 2, a = 1, b = 3/5. Verificaremos que se cumpla el teorema elemental de renovación y se graficará la trayectoria
Sea:
\(c_3 = \frac{\lambda}{\alpha c}e^{-(\alpha-\frac{\lambda}{c})u} = \frac{2}{\frac{3}{5}(5)}\)
\(d_3 = \mathbb{P}[R_t < 0]\)
La probabilidad de ruina se calculara con la función SCL.
## [1] 0.6666667## [1] 0.681Vemos que la probabilidad de ruina sí se aproxima. Finalmente se muestra la gráfica de la trayectoria, podemos ver que llega a ruina casi inmediatamente.
4. Cramer Lundberg general.
Al igual que en el inciso anterior comenzamos por crear una función que modele un proceso de CL pero está vez general.
# CLG: función que modela un proceso de Cramer Lundberg con un proceso "general"
# tf = tiempo final; u = capital inicial; c = prima; ta,tb = parámetros del proceso
# a,b = parámetros de Yn's gamma(a,b)
CLG <- function(tf, u, c, ta, tb, a, b){
# Como cambiaremos el proceso poisson por un general, usaremos otra función
poicg <- function(tf, ta, tb, a, b){
Yn <- c()
p <- renewal(tf, ta, tb)
for(i in 1:p$Nt){
Yn[i] = rgamma(1, a, b)
}
return(list(num_llegadas = p$Nt, tiempo_salto = p$Wn_s, tamaño_salto = Yn))
}
p <- poicg(tf, ta, tb, a, b)
Rt <- c()
for (i in 1:p$num_llegadas){
Rt[i] <- u + c*p$tiempo_salto[i] - sum(p$tamaño_salto[1:i])
}
df_renov <- data.frame(x = c(0,p$tiempo_salto), y = c(u,Rt))
graf <- ggplot(df_renov, aes(x = x , y = y, xend = c(0, p$tiempo_salto), yend = c(u, Rt))) +
geom_line(color = "red")+
ggtitle("Trayectoria de un Proceso de Cramer Lundberg")+
xlab("Tiempo")+
ylab("Valor")
return(list(numero_saltos = p$num_llegadas, tiempo_salto = p$tiempo_salto,
tamaño_salto = p$tamaño_salto, trayectoria = Rt, graf ))
}Con está función creamos una nueva que simule N procesos CL para calcular la probabilidad de ruina.
# SCLG: función que simula N procesos de Cramer Lundberg General para mostrar la probablidad de ruina
# N = número de simulaciones; tf = tiempo final; u = capital inicial;
# c = prima; ta,tb = parámetros del proceso; a,b = parámetros de Yn's gamma(a,b)
SCLG <- function(N, tf, u, c, ta, tb, a, b){
vt <- c()
for (i in 1:N){
vt[i] <- sum(CLG(tf, u, c, ta, tb, a, b)$trayectoria < 0)
}
return(sum(vt > 0)/N)
}Prueba
Tomando los siguientes parametros: tf = 50000, u = 1000, c = 5, ta = 10, tb = 2, a = 1, b = 1/20.
Sea:
\(d_4 = \mathbb{P}[R_t < 0]\)
\(d_4 \hspace{1mm}\) la probabilidad de ruina, se calculará con la funución que acabamos de crear
## [1] 0Gráfica del proceso.
Podemos ver que la probablidad de ruina, \(\hspace{1mm} d_4 = 0 \hspace{1mm}\) tienen sentido pues el proceso nunca toma valores menores que 0.
Cadenas de Markov a tiempo continuo
1. Genere una trayectoria de la cadena de Markov a tiempo continuo.
Apartir de \(\hspace{1mm} T, \tau \hspace{1mm}\) obtenemos la matriz de transición así como de parámetros.
mt <- matrix(c(-2, 1, 1, 0, 0, 0, -4, 1, 0, 2, 1, 1, -5, 1, 0, 1, 0, 0, -3, 0,
0, 0, 1, 1, -3), nrow = 5, byrow = TRUE)
pi <- c(1/5, 1/5, 1/5, 1/5, 1/5, 0)
e <- matrix(rep(1, 5), nrow = 5, byrow = TRUE)
t <- -mt %*% e
Q <- cbind(mt, t)
Q <- rbind(Q, rep(0,6))
l <- -diag(Q)
p = matrix(0, nrow = nrow(Q), ncol = ncol(Q))
for(i in 1:nrow(Q)){
for(j in 1:ncol(Q)){
if(i != j & Q[i,i] != 0){
p[i,j] <- -Q[i,j]/Q[i,i]
}
}
}Después se creará la función que genera la trayectoria de la cadena.
PH <- function(p, pi, l){
tr <- c()
ti <- c()
tr[1] <- sample(c(1:nrow(p)), 1, prob = pi)
ti[1] <- rexp(1, l[tr[1]])
i <- 2
cond <- FALSE
while(cond == FALSE){
tr[i] <- sample(c(1:nrow(p)), 1, prob = p[tr[i-1], ])
tr[i]
if(tr[i] != 6){
ti[i] <- rexp(1,l[tr[i]])
i <- i + 1
}
else{
cond <- TRUE
}
}
df_renov <- data.frame(Estado= tr, Tiempo = c(0,cumsum(ti)))
graf <- ggplot(df_renov, aes(x = Tiempo , y = Estado, xend = c(cumsum(ti), sum(ti) + 1), yend = tr)) +
geom_segment(color = "red")+
ggtitle("Trayectoria Cadena de Markov")+
xlab("Tiempo")+
ylab("Estados")
return(list('Tayectoria' = df_renov, "tau" = sum(ti), graf))
}La función devuelve un data frame con los tiempos en el que se está en alguno de los estados y la gráfica de la trayectoria.
2. Esperanza y varianza
A partir de la función anterior crearemos una nueva que simule N trayectorias y nos devuelva la esperanza y varianza de las \(\hspace{1mm} \tau \hspace{1mm}\) simuladas.
SPH <- function(N, p, pi, l){
tau <- c()
for(i in 1:N){
tau[i] <- PH(p, pi, l)$tau
}
return(cat('Esperanza:', mean(tau), ', Varianza:', var(tau)))
}Obtendremos primero la esperanza y varianza teóricas.
\[\mathbb{E}(t) = \pi(-T)^{-1}e = 0.932\] \[Var(\tau) = 2 \pi (-T)^{-2}e-\mathbb{E}(\tau)^2 = 0.950\]
Ahora las estimaremos mediante 100, 000 simulaciones.
## Esperanza: 0.9331064 , Varianza: 0.939557Como podemos ver ambas son buenas estimaciones.
Browniano, cálculo estocástico, derivados
1. Código de una función que devuelve los precios teóricos a partir de la fórmula de Black-Scholes.
# BS: función que calcula los precios teóricos de put y call
# tf = tiempo final, s0 = precio inicial, K = precio strike, R = tasa, sigma = volatilidad
BS <- function(tf, s0, K, R, sigma){
d1 <- (log(s0/K) + (R + ((sigma^2)/2))*tf)/(sigma*sqrt(tf))
d2 <- d1 - sigma*sqrt(tf)
c0 <- s0*pnorm(d1) - K*exp(-(R*tf))*pnorm(d2)
p0 <- K*exp(-(R*tf))*pnorm(-(d2)) - s0*pnorm(-(d1))
return(cat('Call:', c0, ', Put:', p0))
}2. Código de una función que estima los precios usando S_CRR
# S_CRR: función que simula los precios de calls y puts
# s = número de simulaciones, tf = tiempo final, s0 = precio inicial
# K = precio strike, R = tasa, N = periodos en el árbol, sigma = volatilidad
SCRR <- function(s, tf, s0, K, R, sigma, N){
# Funión de Alex modificada
S_CRR<-function(s0,tf,R,N){
r <- R*tf/N
d <- (1+r)*exp(-sigma*sqrt(tf/N))
u <- (1+r)*exp(sigma*sqrt(tf/N))
q <- ((1+r)-d)/(u-d)
updowns<-rbinom(N,1,q)
incre<-ifelse(updowns==0,d,u)
return(c(s0,s0*cumprod(incre)))
}
c <- c()
p <- c()
for(i in 1:s){
c[i] <- S_CRR(s0,tf,R,N)[N+1] - K
p[i] <- K - S_CRR(s0,tf,R,N)[N+1]
}
cp <- ifelse(c > 0, c, 0)
pp <- ifelse(p > 0, p, 0)
c0 <- exp(-(R*tf))*mean(cp)
p0 <- exp(-(R*tf))*mean(pp)
return(cat('Call:', c0, ', Put:', p0))
}3. Código de una función que estima los precios usando S_BS1
# SBS1: función que simula los precios de calls y puts
# s = número de simulaciones, tf = tiempo final, s0 = precio inicial
# K = precio strike, R = tasa, N = discretización del tiempo, sigma = volatilidad
SBS1 <- function(s, tf, s0, K, R, sigma, N){
# Función de Alex
S_BS1<-function(s0,R,sigma,tf,N){
dt<-tf/N
S_t <- rep(s0,(N+1))
for(i in 2:(N+1)){
S_t[i]<- S_t[i-1] + S_t[i-1]*(R*dt + sigma * sqrt(dt)*rnorm(1))
}
return(S_t)
}
c <- c()
p <- c()
for(i in 1:s){
c[i] <- S_BS1(s0,R,sigma,tf,N)[N+1] - K
p[i] <- K - S_BS1(s0,R,sigma,tf,N)[N+1]
}
cp <- ifelse(c > 0, c, 0)
pp <- ifelse(p > 0, p, 0)
c0 <- exp(-(R*tf))*mean(cp)
p0 <- exp(-(R*tf))*mean(pp)
return(cat('Call:', c0, ', Put:', p0))
}Comparación.
A continuación se muestra la comparación entre los resultados obtenidos por las simulaciones contra el valor teórico.