Cargamos las librerías necesarias y fijamos un tiempo preliminar pequeño para que se vean bonitas las gráficas, también inicializamos algunos parámetros.
library(ggplot2)
library(tidyverse)
library(plotly)
source('C:/Users/luisa/OneDrive/Documentos/own_themes.R')
theme_set(theme_luis())
T_f <- 50Simularemos trayectorias de un proceso de renovación cuyos tiempos interarribo distribuyen \(Gamma(\alpha, \beta)\), es decir, tienen la siguiente función de distribución: \[ f(x) = \frac{\beta^\alpha}{\Gamma(\alpha)}~x^{\alpha-1}~e^{-\beta x} \quad \quad \text{donde $shape = \alpha~$ y $~rate = \beta$ } \]
# Function that simulates the trajectories
# Receives:
# t: final time
# alpha: shape parameter
# beta: rate parameter
# Returns: Arrival times (excluding the initial 0)
ren_traj <- function(t, alpha, beta){
i <- 0
Tn <- i
Ti <- rgamma(1, shape = alpha, rate = beta)
Wn <- cumsum(Tn)
while (Wn[length(Wn)] <= t){
Tn <- c(Tn, Ti)
Wn <- cumsum(Tn)
Ti <- rgamma(1, shape = alpha, rate = beta)
}
Wn <- Wn[-c(1,length(Wn))]
return(Wn)
}Veamos la gráfica de una trayectoria
alpha <- 2
beta <- 1
example <- ren_traj(T_f, alpha, beta)
steps <- tibble::tibble(w = c(0, example), Arrivals = seq_along(c(0, example)) - 1)
df <- steps %>%
mutate(type = "Arrivals") %>%
bind_rows(steps %>%
mutate(type = "Prior",
Arrivals = lag(Arrivals))) %>%
drop_na() %>%
arrange(w, desc(type))
aux <- as.numeric(df[length(df$w),1])
auy <- as.numeric(df[length(df$w),2])
plot <- ggplot(df) + geom_point(aes(w, Arrivals, fill = type), shape = 21, size = 1) +
scale_fill_manual(values = c("black", "white")) +
geom_segment(aes(x = lag(w), y = lag(Arrivals), xend = w, yend = Arrivals, lty = type)) +
scale_linetype_manual(values = c("dashed", "solid")) + xlab("Time") +
theme(legend.position = "none") + geom_segment(aes(x = aux, y = auy, xend = T_f, yend = auy))
ggplotly(plot)Simularemos trayectorias de un proceso Poisson compuesto \((X_t)\) donde la distribución de los saltos es \(Gamma(\alpha, \beta)\), es decir: \[ X_t = \sum\limits_{n = 1}^{N_t} Y_n \quad \quad \text{donde $N_t$ es un proceso Poisson($\lambda$)} \]
# Function that simulates the trajectories of a Compound Poisson Process
# Receives:
# t: final time
# lambda: intensity of the Poisson Process
# alpha: shape parameter
# beta: rate parameter
# Returns: List of 2
# - Arrival times (excluding the initial 0)
# - Size of the jumps
compound_poisson <- function(t, lambda, alpha, beta){
i = 0
Tn <- i
Ti <- rexp(1, lambda)
Wn <- cumsum(Tn)
while (Wn[length(Wn)] <= t){
Tn <- c(Tn, Ti)
Wn <- cumsum(Tn)
Ti <- rexp(1, lambda)
}
Wn <- Wn[-c(1,length(Wn))]
jumps <- length(Wn)
Yn <- rgamma(jumps, shape = alpha, rate = beta)
final <- list(time_jumps = Wn, size_jumps = Yn)
}Veamos la gráfica de una trayectoria
T_f <- 25
lambda <- 2
example <- compound_poisson(T_f, lambda, alpha, beta)
steps <- tibble::tibble(t = c(0, example$time_jumps),
Arrivals = c(0, cumsum(example$size_jumps)))
df <- steps %>%
mutate(type = "Arrivals") %>%
bind_rows(steps %>%
mutate(type = "Prior",
Arrivals = lag(Arrivals))) %>%
drop_na() %>%
arrange(t, desc(type))
aux <- as.numeric(df[length(df$t),1])
auy <- as.numeric(df[length(df$t),2])
plot <- ggplot(df) + geom_point(aes(t, Arrivals, fill = type), shape = 21, size = 1) +
scale_fill_manual(values = c("black", "white")) +
geom_segment(aes(x = lag(t), y = lag(Arrivals), xend = t, yend = Arrivals, lty = type)) +
scale_linetype_manual(values = c("dashed", "solid")) + xlab("Time") +
theme(legend.position = "none") + geom_segment(aes(x = aux, y = auy, xend = T_f, yend = auy))
ggplotly(plot)El siguiente código simula las trayectorias de un proceso de Cramer-Lundberg, donde usamos nuestra función del proceso Poisson compuesto. El proceso de Cramer-Lundberg se define como: \[ R_t = u + ct - \sum\limits_{n=1}^{N_t} Y_n \]
# Function that simulates the trajectories of a Cramer-Lundberg Process
# Receives:
# t: final time
# u: initial surplus
# c: premium
# lambda: intensity of the Poisson Process
# alpha: shape parameter
# beta: rate parameter
# Returns: List of 2
# - Arrival times (including the initial 0)
# - Size of the jumps (Rt, including initial surplus)
cramer_lundberg <- function(t, u, c, lambda, alpha, beta){
aux <- compound_poisson(t, lambda, alpha, beta)
times <- aux$time_jumps
jumps <- aux$size_jumps
Rt <- c(u, u + c*times - cumsum(jumps))
final <- list(time_jumps = c(0,times), size_jumps = Rt)
return(final)
}Veamos una trayectoria del modelo de Cramer-Lundberg
u <- 0
c <- 7
alpha <- 1
beta <- 1/3
example <- cramer_lundberg(T_f, u, c, lambda, alpha, beta)
steps <- tibble::tibble(t = example$time_jumps,
Rt = example$size_jumps)
df <- steps %>%
mutate(type = "Rt") %>%
bind_rows(steps %>%
mutate(type = "t",
Rt = lag(Rt))) %>%
drop_na() %>%
arrange(t, desc(type))
#Sacamos las Ti's (incluyendo el 0) para multiplicarlo por la prima
dif <- diff(example$time_jumps)
sum_ct <- c*dif
# Como el data frame (df) tiene en las posiciones pares el valor final de Rt a
# tiempo t-1, le sumaremos la prima*t en ese intervalo para luego caer al valor
# de Rt a tiempo t
even <- seq(2, dim(df[,1])[1] - 1, 2)
df$Rt[even] <- df$Rt[even] + sum_ct
aux <- as.numeric(df[length(df$t),1])
auy <- as.numeric(df[length(df$t),2])
plot <- ggplot(df) + geom_point(aes(t, Rt, fill = type), shape = 21, size = 1) +
scale_fill_manual(values = c("black", "white")) +
geom_segment(aes(x = lag(t), y = lag(Rt), xend = t, yend = Rt, lty = type)) +
scale_linetype_manual(values = c("dashed", "solid")) +
theme(legend.position = "none") + geom_segment(aes(x = aux, y = auy, xend = T_f, yend = auy + c*(T_f-example$time_jumps[length((example$time_jumps))])))
ggplotly(plot)Por el teorema elemental de renovación, sabemos que para un proceso \(X_t\):
\[ \mathbb{E}[X_t] = \frac{1}{\mathbb{E}[T_i]} = \frac{1}{\alpha/\beta} \]
El código es:
T_f <- 10000
alpha <- 1
beta <- 5
renewal <- ren_traj(T_f, alpha, beta)
N <- length(renewal)
(Muestral <- N/T_f)## [1] 5.0149## [1] 5
Este proceso es un proceso de renovación con premio, entonces:
\[ \mathbb{E}[X_t] = \frac{\mathbb{E}[Y_i]}{\mathbb{E}[T_i]} = \frac{\alpha/\beta}{1/\lambda} \]
El código es:
T_f <- 10000
lambda <- 2
alpha <- 1
beta <- 5
c_poiss <- compound_poisson(T_f, lambda, alpha, beta)
mu <- 1/lambda
ri <- alpha/beta
(TERP <- ri/mu)## [1] 0.4## [1] 0.400391
Suponemos que \(Y_n \sim Exp(\nu)\), sabemos que la probabilidad de ruina se acerca a:
\[ \psi(u) = \frac{\lambda}{\nu c}~e^{-(\nu - \frac{\lambda}{c})u} \]
Vamos a simular 1000 trayectorias de un Proceso de Cramer-Lundberg hasta un \(T_f = 200\) y contar en cuentas de ellas hubo ruina.
#Recordando que si X ~ Gamma(alpha = 1, beta = nu) entonces X ~ Exp(nu)
T_f <- 200
u <- 0
c <- 5
lambda <- 2
alpha <- 1
nu <- 0.8
nSims <- 10000
#Prob teórica
(TEO <- (lambda/(nu*c))*exp(-(nu - (lambda/c))*u))## [1] 0.5# Simulación
ruin <- 0
for(i in 1:nSims){
cr_lbg <- cramer_lundberg(T_f, u, c, lambda, alpha, nu)
if(sum(cr_lbg$size_jumps < 0) > 0){
ruin <- ruin +1
}
}
ruin/nSims## [1] 0.4903Aquí adecuamos la función del proceso de Cramer-Lundberg, ahora con \(T_n \sim \Gamma(10,2)\) y las \(Y_n \sim Exp(\frac{1}{20})\). Estimaremos la probabilidad de ruina con \(u = 1000\) y \(c = 5\).
La nueva función es:
# Function that simulates the trajectories of a Cramer-Lundberg Process
# Receives:
# t: final time
# u: initial surplus
# c: premium
# alpha_r: shape parameter of interarrival distribution
# beta_r: rate parameter of interarrival distribution
# lambda: parameter of the reclaim distribution (in this case, exponential)
# Returns: List of 2
# - Arrival times (including the initial 0)
# - Size of the jumps (Rt, including initial surplus)
cram_lund_renew <- function(t, u, c, alpha_r, beta_r, lambda){
times <- ren_traj(t, alpha_r, beta_r)
jumps <- rexp(length(times), lambda)
Rt <- c(u, u + c*times - cumsum(jumps))
final <- list(time_jumps = c(0, times), size_jumps = Rt)
return(final)
}Y la probabilidad de ruina es:
T_f <- 500
nSims <- 1000
ruin <- 0
ruin <- 0
for(i in 1:nSims){
clr <- cram_lund_renew(T_f, u = 1000, c = 5, 10, 2, 1/20)
if(sum(clr$size_jumps < 0) > 0){
ruin <- ruin +1
}
}
ruin/nSims## [1] 0Que la probabilidad de ruina de \(0\) es lógico ya que en promedio las reclamaciones llegan cada 5 unidades de tiempo, lo cual multiplicado por la prima nos da que recabamos 25 unidades monetarias durante ese tiempo, lo que nos daría hasta este momento una ganancia de 25 y como las reclamaciones son en promedio de 20 unidades monetarias, la ganancia neta esperada por reclamación sería de 5 unidades monetarias, y al tener un colchón de 1,000 unidades monetarias como capital inicial, suena casi imposible ir a la ruina.
Por curiosidad podemos ver una trayectoria con estos valores.
T_f <- 50
example <- cram_lund_renew(T_f, u = 1000, c = 5, 10, 2, 1/20)
steps <- tibble::tibble(t = example$time_jumps,
Rt = example$size_jumps)
df <- steps %>%
mutate(type = "Rt") %>%
bind_rows(steps %>%
mutate(type = "t",
Rt = lag(Rt))) %>%
drop_na() %>%
arrange(t, desc(type))
#Sacamos las Ti's (incluyendo el 0) para multiplicarlo por la prima
dif <- diff(example$time_jumps)
sum_ct <- c*dif
# Como el data frame (df) tiene en las posiciones pares el valor final de Rt a
# tiempo t-1, le sumaremos la prima*t en ese intervalo para luego caer al valor
# de Rt a tiempo t
even <- seq(2, dim(df[,1])[1] - 1, 2)
df$Rt[even] <- df$Rt[even] + sum_ct
aux <- as.numeric(df[length(df$t),1])
auy <- as.numeric(df[length(df$t),2])
plot <- ggplot(df) + geom_point(aes(t, Rt, fill = type), shape = 21, size = 1) +
scale_fill_manual(values = c("black", "white")) +
geom_segment(aes(x = lag(t), y = lag(Rt), xend = t, yend = Rt, lty = type)) +
scale_linetype_manual(values = c("dashed", "solid")) +
theme(legend.position = "none") + geom_segment(aes(x = aux, y = auy, xend = T_f, yend = auy + c*(T_f-example$time_jumps[length((example$time_jumps))])))
ggplotly(plot)Consideremos \(\tau \sim PH(\pi, T)\), con \(\pi = (\frac{1}{5},\frac{1}{5},\frac{1}{5},\frac{1}{5},\frac{1}{5})\) y \[ T = \begin{pmatrix} -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 \end{pmatrix} \]
Declarando nuestras variables:
pi <- matrix(rep(1/5,5), nrow = 1)
T_matrix <- 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 = T)
#Auxiliares para simular una sola trayectoria
t <- matrix(-rowSums(T_matrix), ncol = 1)
zeros <- matrix(rep(0,6), nrow = 1)
Lambda <- as.matrix(rbind(cbind(T_matrix, t),zeros))
diagonal <- c(diag(T_matrix)^(-1), 0)
###El código que genera una trayectoria de la cadena de Markov asociada a \(\tau\). Consiste en seleccionar aleatoriamente un estado inicial \(i\) de acuerdo al vector de probabilidades \(\pi\). El siguiente estado es seleccionado por un nuevo vector de probabilidad dado por la i-ésima fila de la cadena asociada a \(\Lambda\), esto es:
\[ S = I - diag(\Lambda)^{-1}~\Lambda \]
Entonces usando \(S_i\) como nuevo vector de probabilidades se selecciona un nuevo estado \(j\), se simula el tiempo que se tardó en llegar a este nuevo estado con una \(Exp(-\lambda_{ii})\) y se guarda en \(x\). Si el nuevo estado \(j\) es el estado absorbente, devolvemos \(x\), sino, se vuelve a elegir otro estado \(k\), con el vector de probabilidades \(S_j\) y se suma a \(x\) el tiempo en que llegó de \(j\) a \(k\) y así sucesivamente.
Primero veamos una trayectoria de esta cadena de Markov:
x <- c()
state <- c()
y <- 1
name <- as.character(1:6)
state[y] <- sample(name, size = 1, prob = c(pi,0))
while(state[y] != name[length(name)]){
S <- diag(6) - (`diag<-`(matrix(0, 6, 6), diagonal)) %*% Lambda
aux <- as.integer(state[y])
a_ii <- Lambda[aux, aux]
x[y] <- rexp(1, -a_ii)
y <- y + 1
S_i <- S[aux,]
state[y] <- sample(name, size = 1, prob = S_i)
}
steps <- tibble::tibble(Time = c(0, cumsum(x)),
States = as.integer(state))
df <- steps %>%
mutate(type = "States") %>%
bind_rows(steps %>%
mutate(type = "Time",
States = lag(States))) %>%
drop_na() %>%
arrange(Time, desc(type))
plot <- ggplot(df) + geom_point(aes(Time, States, fill = type), shape = 21, size = 1) +
scale_fill_manual(values = c("black", "white")) +
geom_segment(aes(x = lag(Time), y = lag(States), xend = Time, yend = States, lty = type)) +
scale_linetype_manual(values = c("dashed", "solid")) +
theme(legend.position = "none") + scale_y_continuous(breaks = seq(1,6,1))
plotly::ggplotly(plot)Utilizando la lógica de nuestro código anterior, crearemos una función que nos de el tiempo de absorción en cada realización o simulación.
ph_type <- function(n, pi, T_matrix){
times <- rep(0, n)
size <- length(pi)
name <- as.character(1:(size + 1))
t <- matrix(-rowSums(T_matrix), ncol = 1)
zeros <- matrix(rep(0, size + 1), nrow = 1)
Lambda <- as.matrix(rbind(cbind(T_matrix, t), zeros))
diagonal <- c(diag(T_matrix)^(-1), 0)
for(i in 1:n){
x <- 0
state <- sample(name, size = 1, prob = c(pi,0))
while(state != name[length(name)]){
S <- diag(size + 1) - (`diag<-`(matrix(0, size + 1, size + 1), diagonal)) %*% Lambda
aux <- as.integer(state)
S_i <- S[aux,]
state <- sample(name, size = 1, prob = S_i)
a_ii <- Lambda[aux, aux]
x <- x + rexp(1, -a_ii)
}
times[i] <- x
}
return(times)
}Estimaremos la esperanza al realizar 100,000 realizaciones de \(\tau\) y la comparamos con la media teórica (\(\mathbb{E}[\tau] = \pi(-T)^{-1}e\))
## [1] 0.9290841## [1] 0.932Podemos ver que son casi iguales.
Ahora estimamos la varianza con el mismo número de realizaciones (la varianza teórica está dada por \(Var(\tau) = 2\pi(-T)^{-2}e~ - (\mathbb{E}[\tau])^2\)).
## [1] 0.9465479## [1] 0.950576También son casi iguales.
Aquí creamos una función que nos de los precios de opciones europeas usando la fórmula de Black-Scholes.
# Función que calcula el valor una opción europea utilizando Black Scholes
# Parámetros:
# s: Precio del subyacente a tiempo t
# k: Precio strike
# r: Tasa libre de riesgo, tiene que ponerse en decimales
# sigma: Volatilidad anual, tiene que ponerse en decimales
# t: Tiempo de valuación, t está entre 0 y el vencimiento (años)
# V: Vencimiento de la opción en años (T)
# tipo: "c" para call, "p" para put
# q: Tasa de dividendos anual compuesta continuamente en decimales
BlackScholes_t <- function(s, k, r, sigma, t = 0, V, tipo="c", q = 0){
if(tipo == "c") j <- 1
if(tipo == "p") j <- (-1)
d1 <- (log(s/k)+(r-q+(sigma^2)/2)*(V-t))/(sigma*sqrt(V-t))
d2 <- (log(s/k)+(r-q-(sigma^2)/2)*(V-t))/(sigma*sqrt(V-t))
f_t <- j*s*exp(-q*(V-t))*pnorm(j*d1) - j*k*exp(-r*(V-t))*pnorm(j*d2)
return(f_t)
} Veamos el precio (a tiempo 0) que devuelve esta función de un activo con las siguientes características:
sigma <- 0.18
R<- 0.06
K <- 25
S0<- 20
t_F <- 2
(callBS <- BlackScholes_t(S0, K, R, sigma, 0, t_F, "c"))## [1] 1.220978## [1] 3.393989Ahora simularemos el precio a partir de dos funciones, una usando el modelo de Cox-Ross-Rubinstein y otra usando el método de Euler para generar la trayectoria del browniano geométrico. Éstas son las funciones:
# Cox-Ross-Rubinstein
S_CRR <- function(S_0, r, u, d, N){
#calculamos q
q <- ((1+r)-d)/(u-d)
#generamos las bernoullis de cada paso
updowns <- rbinom(N,1,q)
#la convertimos al un vector de u's y d's
incre <- ifelse(updowns==0,d,u)
return(c(S_0, S_0*cumprod(incre)))
}
#Forma 1, Euler
S_BS1 <- function(S_0, R, sigma, t_F, N){
dt <- t_F/N
#generamos la trayectoria del browniano geométrico por el método de Euler
S_t <- rep(S_0,(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)
}Y las primas para el call y el put las calcularemos usando las siguientes funciones (hechas para simular las trayectorias, promediar las primas y traerlas a valor presente):
#Simulación de Cox-Ross-Rubinstein
sim_CRR <- function(n_sims, T_f, k, S0, R, sigma, periodos, tipo){
r <- R*t_F/periodos
d <- (1+r)*exp(-sigma*sqrt(t_F/periodos))
u <- (1+r)*exp(sigma*sqrt(t_F/periodos))
payoff <- 0
payoff_call <- function(s,k){max(s-k,0)}
payoff_put <- function(s,k){max(k-s,0)}
for(i in 1:n_sims){
precio <- S_CRR(S0, r, u, d, periodos)[periodos+1]
if(tipo == "c"){
payoff <- payoff + payoff_call(precio, k)
}
if(tipo == "p"){
payoff <- payoff + payoff_put(precio, k)
}
}
prom <- (payoff/n_sims)*exp(-R*T_f)
return(prom)
}
#Simulación bajo Euler
sim_BS1 <- function(n_sims, T_f, k, S0, R, sigma, periodos, tipo){
payoff <- 0
payoff_call <- function(s,k){max(s-k,0)}
payoff_put <- function(s,k){max(k-s,0)}
for(i in 1:n_sims){
precio <- S_BS1(S0, R, sigma, T_f, periodos)[periodos+1]
if(tipo == "c"){
payoff <- payoff + payoff_call(precio, k)
}
if(tipo == "p"){
payoff <- payoff + payoff_put(precio, k)
}
}
prom <- (payoff/n_sims)*exp(-R*T_f)
return(prom)
}Ahora las primas con los datos que tenemos son:
## [1] 1.220978## [1] 1.229409## [1] 1.234873## [1] 3.393989## [1] 3.375643## [1] 3.367148
En una tabla quedan así:
| Method | Call | Put |
|---|---|---|
| Black-Scholes | 1.220978 | 3.393989 |
| Cox-Ross-Rubinstein | 1.229409 | 3.375643 |
| Euler | 1.234873 | 3.367148 |