Tarea 3
José Roberto del Valle
2022-08-16
if("moments" %in% rownames(installed.packages()) == FALSE) {install.packages("moments",repos="http://cran.r-project.org")}
try(suppressPackageStartupMessages(library(moments, quietly = TRUE,warn.conflicts = FALSE)),silent = TRUE)Pregunta 1
#a) Utiliza el método Monte Carlo para simular una distribución normal con media μ = 0.1 y = 0.20. Reporta la media, la desviación estándar, el sesgo y la curtosis de la muestra simulada. Utiliza 1,000, 10,000 y 100,000 simulaciones (N).
moments <- function(mu,sigma,n){
simulaciones <- rnorm(n,mu,sigma)
mean <- mean(simulaciones)
sd <- sd(simulaciones)
bias <- skewness(simulaciones)
kurtosis <- kurtosis(simulaciones)
mom_df <- data.frame(c(mean,sd,bias,kurtosis))
row.names(mom_df) <- c("Media","Desviación Estándar", "Sesgo","Curtosis")
return(mom_df)
}sim1 <- moments(0.1,0.2,1000)
sim2 <- moments(0.1,0.2,10000)
sim3 <- moments(0.1,0.2,100000)
moments_df <- cbind(sim1,sim2,sim3)
colnames(moments_df) <- c("1000 sim", "10,000 sim", "100,000 sim")
moments_df## 1000 sim 10,000 sim 100,000 sim
## Media 0.09239689 0.0986058373 0.100301998
## Desviación Estándar 0.20591761 0.1992585057 0.200849753
## Sesgo -0.02027549 0.0002710352 -0.002516025
## Curtosis 2.89837732 2.9515019458 2.995976432b) Interpreta los resultados: Conforme aumentamos simulaciones, los valores se acercan más a los valores poblaciones. Sesgo a cero y curtosis 3
c) Para cada N, repite la simulación del inciso a) 100 veces y reporta la raíz del error cuadrático medio
de cada uno de los estimadores que obtuviste (media, desviación estándar, sesgo y curtósis).
Parámetros
media <- 0.1
sd <- 0.2
sesgo <- 0
kurt <- 3Simulaciones = 1000
sim_rep1 <- rep(1000,100)
sim_rep1_list <- lapply(sim_rep1, function (x) moments (0.1,0.2,x))
sim_rep1_mom <- matrix(unlist(sim_rep1_list),nrow =4,byrow=F)
row.names(sim_rep1_mom) <-c("Media","Desviación Estándar", "Sesgo","Curtosis")n <- 100
ECM_media_1 <- sum((sim_rep1_mom[1,]-media)^2)/n
ECM_sd_1 <- sum((sim_rep1_mom[2,]-sd)^2)/n
ECM_sesgo_1 <- sum((sim_rep1_mom[3,]-sesgo)^2)/n
ECM_kurt_1 <- sum((sim_rep1_mom[2,]-kurt)^2)/n
ECM_1 <- rbind(ECM_media_1,ECM_sd_1,ECM_sesgo_1,ECM_kurt_1)
row.names(ECM_1) <- c("Media","Desviación Estándar", "Sesgo","Curtosis")
ECM_1## [,1]
## Media 4.098004e-05
## Desviación Estándar 1.824402e-05
## Sesgo 7.586456e-03
## Curtosis 7.837722e+00Simulaciones 10,000
sim_rep2 <- rep(10000,100)
sim_rep2_list <- lapply(sim_rep2, function (x) moments (0.1,0.2,x))
sim_rep2_mom <- matrix(unlist(sim_rep2_list),nrow =4,byrow=F)
row.names(sim_rep2_mom) <-c("Media","Desviación Estándar", "Sesgo","Curtosis")
n <- 100
ECM_media_2 <- sum((sim_rep2_mom[1,]-media)^2)/n
ECM_sd_2 <- sum((sim_rep2_mom[2,]-sd)^2)/n
ECM_sesgo_2 <- sum((sim_rep2_mom[3,]-sesgo)^2)/n
ECM_kurt_2 <- sum((sim_rep2_mom[2,]-kurt)^2)/n
ECM_2 <- rbind(ECM_media_2,ECM_sd_2,ECM_sesgo_2,ECM_kurt_2)
row.names(ECM_2) <- c("Media","Desviación Estándar", "Sesgo","Curtosis")
ECM_2## [,1]
## Media 3.594292e-06
## Desviación Estándar 2.167906e-06
## Sesgo 6.107227e-04
## Curtosis 7.839504e+00Simulaciones 100,000
sim_rep3 <- rep(100000,100)
sim_rep3_list <- lapply(sim_rep3, function (x) moments (0.1,0.2,x))
sim_rep3_mom <- matrix(unlist(sim_rep3_list),nrow =4,byrow=F)
row.names(sim_rep3_mom) <-c("Media","Desviación Estándar", "Sesgo","Curtosis")
n <- 100
ECM_media_3 <- sum((sim_rep3_mom[1,]-media)^2)/n
ECM_sd_3 <- sum((sim_rep3_mom[2,]-sd)^2)/n
ECM_sesgo_3 <- sum((sim_rep3_mom[3,]-sesgo)^2)/n
ECM_kurt_3 <- sum((sim_rep3_mom[2,]-kurt)^2)/n
ECM_3 <- rbind(ECM_media_3,ECM_sd_3,ECM_sesgo_3,ECM_kurt_3)
row.names(ECM_3) <- c("Media","Desviación Estándar", "Sesgo","Curtosis")
ECM_3## [,1]
## Media 4.631797e-07
## Desviación Estándar 2.010545e-07
## Sesgo 7.331753e-05
## Curtosis 7.840263e+00# Pregunta 2
Simula el precio de la acción a diferentes horizontes
- Calibra los paramétros u, d y q para un “time-step” diario
rf <- 0.03
S_0 <- 100
media <- 0.1
sd <- 0.2
T <- 1
n <- 252
dt <- T/n
u <- exp(sd*sqrt(dt))
d <- 1/u
q <- (exp(rf*dt)-d)/(u-d)
params <- cbind(u,d,q)
colnames(params) <- c("u","d","q")
params## u d q
## [1,] 1.012679 0.9874802 0.501575- simula el precio de la acción dentro de una semana
S_n <- function(T,n,S_0){
dt <- T/n
u <- exp(sd*sqrt(dt))
d <- 1/u
q <- (exp(rf*dt)-d)/(u-d)
S_sim <- runif(n)
sign <- sign(q-S_sim)
ret_d <- u*(sign>=0)+d*(sign<0)
S_n <- S_0*prod(ret_d)
return(S_n)
}simulaciones semanales
S0_sem_1000 <- rep(S_0,1000)
Sn_sem_1000 <- lapply(S0_sem_1000,function (x) S_n(1/52,5,x))
Sn_sem_vector<- unlist(Sn_sem_1000)para aplicar a 10,000 y 100,000
S0_sem_10000 <- rep(S_0,10000)
Sn_sem_10000 <- lapply(S0_sem_10000,function (x) S_n(1/52,5,x))
Sn_sem_vector2<- unlist(Sn_sem_10000)
S0_sem_100000 <- rep(S_0,100000)
Sn_sem_100000 <- lapply(S0_sem_100000,function (x) S_n(1/52,5,x))
Sn_sem_vector3<- unlist(Sn_sem_100000)Momentos Semanales
#N = 1,000
Sn_sem_mom <- cbind(mean(Sn_sem_vector),sd(Sn_sem_vector),skewness(Sn_sem_vector),kurtosis(Sn_sem_vector))
colnames(Sn_sem_mom) <-c("Media","Desviación Estándar", "Sesgo","Curtosis")
Sn_sem_mom## Media Desviación Estándar Sesgo Curtosis
## [1,] 100.1714 2.824482 -0.01169812 2.518002#N=10,000
Sn_sem_mom2 <- cbind(mean(Sn_sem_vector2),sd(Sn_sem_vector2),skewness(Sn_sem_vector2),kurtosis(Sn_sem_vector2))
colnames(Sn_sem_mom2) <-c("Media","Desviación Estándar", "Sesgo","Curtosis")
Sn_sem_mom2## Media Desviación Estándar Sesgo Curtosis
## [1,] 100.0486 2.757265 0.06234758 2.628591#N = 100,000
Sn_sem_mom3 <- cbind(mean(Sn_sem_vector3),sd(Sn_sem_vector3),skewness(Sn_sem_vector3),kurtosis(Sn_sem_vector3))
colnames(Sn_sem_mom3) <-c("Media","Desviación Estándar", "Sesgo","Curtosis")
Sn_sem_mom3## Media Desviación Estándar Sesgo Curtosis
## [1,] 100.0553 2.769969 0.06210704 2.593623Histogramas
hist(Sn_sem_vector,main="Histograma precios terminales a 1 semana (N=1,000)")
hist(Sn_sem_vector2,main="Histograma precios terminales a 1 semana (N=10,000)")
hist(Sn_sem_vector3,main="Histograma precios terminales a 1 semana (N=100,000)")
#ylab para cambiar los nombres de los ejes- Log precios
Log_sn_sem <- log(Sn_sem_vector)
Log_sn_sem_mom <- cbind(mean(Log_sn_sem),sd(Log_sn_sem),skewness(Log_sn_sem),kurtosis(Log_sn_sem))
colnames(Log_sn_sem_mom) <-c("Media","Desviación Estándar", "Sesgo","Curtosis")
Log_sn_sem_mom## Media Desviación Estándar Sesgo Curtosis
## [1,] 4.606485 0.02822422 -0.07583054 2.520618Log_sn_sem2 <- log(Sn_sem_vector2)
Log_sn_sem_mom2 <- cbind(mean(Log_sn_sem2),sd(Log_sn_sem2),skewness(Log_sn_sem2),kurtosis(Log_sn_sem2))
colnames(Log_sn_sem_mom2) <-c("Media","Desviación Estándar", "Sesgo","Curtosis")
Log_sn_sem_mom2## Media Desviación Estándar Sesgo Curtosis
## [1,] 4.605277 0.02755811 -0.004700617 2.621164Log_sn_sem3 <- log(Sn_sem_vector3)
Log_sn_sem_mom3 <- cbind(mean(Log_sn_sem3),sd(Log_sn_sem3),skewness(Log_sn_sem3),kurtosis(Log_sn_sem3))
colnames(Log_sn_sem_mom3) <-c("Media","Desviación Estándar", "Sesgo","Curtosis")
Log_sn_sem_mom3## Media Desviación Estándar Sesgo Curtosis
## [1,] 4.60534 0.0276831 -0.00385996 2.588783Diario, Mensual y Anual
rf <- 0.03
S_0 <- 100
media <- 0.1
sd <- 0.2
T_m <- 1
n_m <- 12
dt_m <- T/n_m
u_m <- exp(sd*sqrt(dt_m))
d_m <- 1/u
q_m <- (exp(rf*dt_m)-d_m)/(u_m-d_m)
params_m <- cbind(u_m,d_m,q_m)
colnames(params_m) <- c("u mes","d mes","q mes")
T_a <- 1
n_a <- 1
dt_a <- T/n_a
u_a <- exp(sd*sqrt(dt_a))
d_a <- 1/u
q_a <- (exp(rf*dt_a)-d_a)/(u_a-d_a)
params_a <- cbind(u_a,d_a,q_a)
colnames(params_a) <- c("u año","d año","q año")
params## u d q
## [1,] 1.012679 0.9874802 0.501575params_m## u mes d mes q mes
## [1,] 1.059434 0.9874802 0.2087849params_a## u año d año q año
## [1,] 1.221403 0.9874802 0.1837117Valuación de Derivados: 40 puntos
Utilizando las distribuciones de precios del ejercicio
- El valor de un call europeo con un strike
K <- 100
call<-mean((Sn_sem_vector>K)*(Sn_sem_vector-K))*exp(-rf)
#BSM
d1 <- (log(S_0/K)+(rf+sd^2/2)*dt)/(sd*sqrt(dt))
d2 <- d1-sd * sqrt(dt)
call_BSM <- S_0*pnorm(d1)-K*exp(-rf*dt)*pnorm(d2)
print(list("Montecarlo"=call,"Black-Scholes-Merton"= call_BSM))## $Montecarlo
## [1] 1.238694
##
## $`Black-Scholes-Merton`
## [1] 0.5085613- obten el precio de un call asiático con un strike
S_path <- function(T,n,S_0){
dt <- T/n
u <- exp(sd*sqrt(dt))
d <- 1/u
q <- (exp(rf*dt)-d)/(u-d)
S_sim <- runif(n)
sign <- sign(q-S_sim)
ret_d <- u*(sign>=0)+d*(sign<0)
S_path <- S_0*cumprod(ret_d)
return(S_path)
}simulaciones paths semanales
T <- 1/252
n <- 5
S0_sem_1000 <- rep(S_0,1000)
Sn_path_1000 <- lapply(S0_sem_1000,function (x) S_path(1/52,5,x))
Sn_path_matrix<- matrix(unlist(Sn_path_1000),nrow =n,byrow=F)
S_promedio <- apply(Sn_path_matrix,2,mean)
call_asiatica<- mean(exp(-rf)*(S_promedio>K)*(S_promedio-K))
print(paste ("Call asiática =" ,call_asiatica))## [1] "Call asiática = 0.808968391676552"- para un inversionista que se va largo este contrato
p_pos <- sum(Sn_sem_vector>=140)/length(Sn_sem_vector>=140)
p_neg <- sum(Sn_sem_vector<=80)/length(Sn_sem_vector>=140)
print(list("Probabilidad que le paguen" = p_pos," Probabilidad de pagar" = p_neg))## $`Probabilidad que le paguen`
## [1] 0
##
## $` Probabilidad de pagar`
## [1] 0- calcula el precio del contrato
pago_pos <- 1000000
pago_neg <- -25000000
vp <- exp(-rf)*(sum(Sn_sem_vector>=140)/length(Sn_sem_vector)*pago_pos+
sum(Sn_sem_vector<=80)/length(Sn_sem_vector)*pago_neg)
print( paste ("Precio del contrato", vp))## [1] "Precio del contrato 0"- calcula el rendimiento esperado promedio de tu contrato
rdto <- (sum(Sn_sem_vector>=140)/length(Sn_sem_vector)*pago_pos+
sum(Sn_sem_vector<=80)/length(Sn_sem_vector)*pago_neg)/vp
print( paste ("Rendimiento esperado del contrato", rdto))## [1] "Rendimiento esperado del contrato NaN"