Creación de Arboles binomiales para cobertura de un portafolio de inversión

Introducción

Con este proceso vamos a llevar a cabo un análisis de cobertura europea para las acciones NU, ONON y NVDA, utilizando datos recopilados desde el 01/06/2022 al 06/04/2025.

Al final, vamos a considerar una opción de estás 3 acciones, que liquide de forma trimestral hasta la fecha de inversión máxima de la opción (o máximo dos años en el futuro), donde el precio strike se determinará conforme los criterios de liquidez (bid-ask), volatilidad implícita y open interest. Se realizará un “rolling” al finalizar cada contrato, ajustando el precio de prima cada trimestre. La tasa libre de riesgo se establecerá a un trimestre, basada en los bonos del Tesoro (Tbond).

Instalación y carga de paquetes

Antes de iniciar con el proceso de analisis, primero vamos a instalar y cargar los paquetes necesarios

install.packages(c(“tidyquant”,“plotly”,“timetk”,“tibble”, “tidyr”, “quantmod”, “rvest”, “dplyr”, “glue”))

library(tidyquant) # Para descargar los datos
library(plotly) # Para crear gr?ficos interactivos
library(timetk) # Para manipular las series de datos
library(tibble) #Para crear y manejar data frames modernos
library(tidyr) #Para organizar y transformar datos
library(forcats) #Para manipular variables categ?ricas (factores) de forma eficiente.
library(dplyr) #Para manipula datos: filtrar, resumir, transformar, ordenar, agrupar.
library(tidyverse) #Es un conjunto de paquetes integrados para an?lisis de datos limpio y estructurado.
library(quantmod) #Para descargar, analizar y visualizar datos financieros y series de tiempo del mercado.
library(rvest) #Sirve para hacer web scraping (extraer datos de páginas web HTML de forma estructurada)
library(glue) #Sirve para crear strings dinámicos fácilmente, insertando variables dentro de texto
library(stringr) #Sirve para manipular cadenas de texto (buscar, reemplazar, dividir, extraer patrones).
library(ggplot2) #Sirve para crear gráficos elegantes y personalizables con una sintaxis coherente basada en capas.
options(scipen=999)

set.seed(1234)

Descarga y conversión de la información de las acciones

Luego procedemos a descargar los precios historicos de las acciones NU, ONON y NVDA, utilizando datos recopilados desde el 01/06/2022 al 06/04/2025

tick <- c('NU','ONON','NVDA')
price_data <- tq_get(tick,from='2022-06-01',to='2025-04-06',get='stock.prices')
head(price_data)

## # A tibble: 6 × 8
##   symbol date        open  high   low close   volume adjusted
##   <chr>  <date>     <dbl> <dbl> <dbl> <dbl>    <dbl>    <dbl>
## 1 NU     2022-06-01  3.84  3.95  3.66  3.82 31535600     3.82
## 2 NU     2022-06-02  3.82  4.5   3.80  4.44 79498000     4.44
## 3 NU     2022-06-03  4.31  4.54  4.22  4.5  54708600     4.5 
## 4 NU     2022-06-06  4.60  4.75  4.26  4.29 24604400     4.29
## 5 NU     2022-06-07  4.07  4.32  4     4.32 31965900     4.32
## 6 NU     2022-06-08  4.23  4.52  4.23  4.33 17304400     4.33

Para calcular el rendimiento diario de estas acciones

log_ret_tidy <- price_data%>%group_by(symbol)%>%tq_transmute(select=adjusted,mutate_fun=periodReturn,period='daily',col_rename='ret',type='log')
head(log_ret_tidy)

## # A tibble: 6 × 3
## # Groups:   symbol [1]
##   symbol date            ret
##   <chr>  <date>        <dbl>
## 1 NU     2022-06-01  0      
## 2 NU     2022-06-02  0.150  
## 3 NU     2022-06-03  0.0134 
## 4 NU     2022-06-06 -0.0478 
## 5 NU     2022-06-07  0.00697
## 6 NU     2022-06-08  0.00231

lo convertimos a un formato ancho y en un objeto de serie temporal

log_ret_xts <- log_ret_tidy%>%tidyr::pivot_wider(names_from=symbol, values_from=ret)%>%tk_xts()

## Warning: Non-numeric columns being dropped: date

## Using column `date` for date_var.

head(log_ret_xts)

##                      NU         ONON         NVDA
## 2022-06-01  0.000000000  0.000000000  0.000000000
## 2022-06-02  0.150403984  0.071262994  0.067127945
## 2022-06-03  0.013423007 -0.034540346 -0.045528668
## 2022-06-06 -0.047790673 -0.045281916  0.003519312
## 2022-06-07  0.006968718  0.010451369  0.007424809
## 2022-06-08  0.002312082  0.002831501 -0.014586369

A continuación, calculamos los rendimientos diarios medios de cada activo.

mean_ret <- colMeans(log_ret_xts)
print(round(mean_ret, 5))

##      NU    ONON    NVDA 
## 0.00129 0.00088 0.00230

Para calcular la matriz de covarianza para las 3 acciones; La anualizaremos multiplicandola por 252.

cov_mat <- cov(log_ret_xts) * 252
print(round(cov_mat,4))

##          NU   ONON   NVDA
## NU   0.2813 0.0908 0.1089
## ONON 0.0908 0.2727 0.1072
## NVDA 0.1089 0.1072 0.2906

Indice de Sharpe

Antes de aplicar nuestros métodos a miles de carteras aleatorias, demostramos los pasos en una sola cartera. Para calcular la rentabilidad y el riesgo (desviación estándar) de la cartera, necesitaremos:

-Rendimiento medio de los activos -Ponderaciones de la cartera -Matriz de covarianza de todos los activos -Ponderaciones aleatorias

Primero, vamos a crear ponderaciones aleatorias

wts <- runif(n=length(tick))
wts <- wts/sum(wts)

A continuación calculamos la rentabilidad anualizada de la cartera.

port_returns <- (sum(wts*mean_ret)+1)^252-1

además, calculamos el riesgo de la cartera (desviación estándar). Esta será la desviación estándar anualizada de la cartera.

port_risk <- sqrt(t(wts)%*%(cov_mat%*%wts))

Finalemente calculamos el ratio de Sharpe

sharpe_ratio <- port_returns/port_risk
print(wts)

## [1] 0.08452041 0.46258068 0.45289891

print(port_returns)

## [1] 0.4799842

print(port_risk)

##           [,1]
## [1,] 0.4246365

print(sharpe_ratio)

##          [,1]
## [1,] 1.130341

Simulación de 5000 carteras aleatorias

Tenemos todo lo necesario para realizar la optimización. Ahora solo necesitamos ejecutar este código en 5000 carteras aleatorias. Para ello, usaremos un bucle llamado for. Antes de hacerlo, necesitamos crear vectores y una matriz de vacíos para almacenar nuestros valores.

num_port <- 5000

Creamos una matriz para almacenar las ponderaciones

all_wts <- matrix(nrow =num_port,ncol=length(tick))

Creamos un vector vacío para almacenar Rendimientos de la cartera

port_returns <- vector('numeric',length=num_port)

Creamos un vector vacío para almacenar Desviación estándar de la cartera

port_risk <- vector('numeric',length=num_port)

Creamos un vector vacío para almacenar índice de Sharpe de portafolio

sharpe_ratio <- vector('numeric',length=num_port)

Procedemos a ejecutar el bucle [for] 5000 veces.

for (i in seq_along(port_returns)) {
  wts <- runif(length(tick))
  wts <- wts/sum(wts)
  # Almacenar ponderaci?n en la matriz
  all_wts[i,] <- wts
  # Rendimientos de la cartera
  port_ret <- sum(wts*mean_ret)
  port_ret <- ((port_ret+1)^252)-1
  # Almacenamiento de valores de retorno de cartera
  port_returns[i] <- port_ret
  # Creaci?n y almacenamiento de riesgos de cartera
  port_sd <- sqrt(t(wts)%*%(cov_mat%*%wts))
  port_risk[i] <- port_sd
  # Creaci?n y almacenamiento de ratios de Sharpe de cartera
  # Suponiendo una tasa libre de riesgo del 0 %
  sr <- port_ret/port_sd
  sharpe_ratio[i] <- sr
}

Ahora creamos una tabla de datos para almacenar todos los valores juntos.

portfolio_values <- tibble(Return=port_returns,Risk=port_risk,SharpeRatio=sharpe_ratio)

Convertimos una matriz en un tibble y cambiamos los nombres de las columnas

all_wts <- tk_tbl(all_wts)

## Warning in tk_tbl.data.frame(as.data.frame(data), preserve_index, rename_index,
## : Warning: No index to preserve. Object otherwise converted to tibble
## successfully.

colnames(all_wts) <- colnames(log_ret_xts)

Combinamos todos los valores juntos

portfolio_values <- tk_tbl(cbind(all_wts, portfolio_values))

## Warning in tk_tbl.data.frame(cbind(all_wts, portfolio_values)): Warning: No
## index to preserve. Object otherwise converted to tibble successfully.

Tenemos las ponderaciones de cada activo, el riesgo y la rentabilidad, junto con el ratio de Sharpe de cada cartera.

La cartera de mínima varianza

Aunque algunos activos tengan alta volatilidad individual, al combinarlos con otros que estén poco correlacionados, el riesgo total puede reducirse.

min_var <- portfolio_values[which.min(portfolio_values$Risk),]

Graficamos las ponderaciones de la cartera de mínima varianza

p <- min_var%>%
  pivot_longer(cols=NU:NVDA, names_to="Asset", values_to="Weights")%>%
  mutate(Asset=as.factor(Asset))%>%
  ggplot(aes(x=fct_reorder(Asset, Weights), y=Weights, fill=Asset))+
  geom_bar(stat='identity')+
  theme_minimal()+
  labs(x='Assets', y='Weights', title="Minimum Variance Portfolio Weights")+
  scale_y_continuous(labels=scales::percent)
ggplotly(p)

ANALISIS:

Bajo el portafolio de minima varianza, asumiendo una posicion conservadora, el peso de los 3 activos seleccionados son 29.49% para NVDA, 34.24% para NU y de 36.26% para ONON.

Asignamos USD 1.000.000 según las ponderaciones Extraemos las ponderaciones de min_var

weights_row <- min_var %>% select(NU, ONON, NVDA)

Calculamos la asignación en dólares (1 millón)

allocation_row <- weights_row * 1e6

Creamos y mostramos la tabla usando [rbind]

min_var_allocation <- rbind(Weights = weights_row,AllocationUSD = allocation_row)
print(min_var_allocation)

## # A tibble: 2 × 3
##           NU       ONON       NVDA
## *      <dbl>      <dbl>      <dbl>
## 1      0.340      0.366      0.293
## 2 340497.    366031.    293473.

Traemos los últimos precios historicos

ultimos_precios <- price_data %>%
  group_by(symbol) %>%
  filter(date == max(date)) %>%
  select(symbol, adjusted) %>%
  rename(Precio_Final = adjusted)

Tabla resumen con cantidades iniciales

resumen_inicial <- tibble(
  Activo = c("NU", "ONON", "NVDA"),USD_Asignado = as.numeric(allocation_row))%>%
  left_join(ultimos_precios, by = c("Activo" = "symbol")) %>%
  mutate(Cantidad_Acciones = floor(USD_Asignado / Precio_Final),Valor_Real = Cantidad_Acciones * Precio_Final)

Valor inicial del portafolio

valor_inicial <- sum(resumen_inicial$Valor_Real)

SIMULACIÓN NU

nu_data <- price_data %>%
  filter(symbol == "NU") %>%
  select(date, adjusted)

price_nu=nu_data[,2]
log_ret_nu <- coredata(log_ret_xts[, "NU"])[-1]

plot.ts(price_nu)

plot.ts(log_ret_nu)

factores iniciales

year1=252
s01=as.numeric(nu_data[NROW(nu_data),2])
T1=504
r1=exp(mean(log_ret_nu)*year1)-1
sigma1= as.numeric(sqrt(var(log_ret_nu)*year1))
N1=5000 #N?mero de simulaciones
dt1=1/year1
stx1=as.numeric(s01*exp(r1*(T1/year1)))

vec1=rep(s01,N1)
mb1=matrix(ncol=N1, nrow=T1)
mb1[1,]=vec1

for (i in 1:N1) {
  for (t in 2:T1) {
    mb1[t,i]=mb1[(t-1),i]*exp((r1-(0.5*(sigma1^2)))*(dt1)+sigma1*(dt1^(1/2))*qnorm(runif(1, min=0, max=1)))                                                                      
  }
}

validacion antes de grafico

print(paste("Último precio de NU:", s01))

## [1] "Último precio de NU: 9.60000038146973"

print(paste("Volatilidad anualizada (sigma1):", sigma1))

## [1] "Volatilidad anualizada (sigma1): 0.530790771707883"

print(paste("Rendimiento anual esperado (r1):", r1))

## [1] "Rendimiento anual esperado (r1): 0.384994190335734"

final_prices <- mb1[nrow(mb1), ]
quantile(final_prices, probs = c(0.01, 0.05, 0.5, 0.95, 0.99))

##        1%        5%       50%       95%       99% 
##  2.789032  4.585439 15.578939 51.583031 84.813421

max_value <- max(mb1)
print(max_value)

## [1] 354.8852

which_max <- which(mb1 == max_value, arr.ind = TRUE)
print(which_max)

##      row  col
## [1,] 498 1278

ANALISIS:

Dado que bajo el Movimiento Browniano Geométrico (MBG) los precios se simulan de forma exponencial, se observa que aproximadamente el 1% de las trayectorias presentan comportamientos explosivos. Esto se debe a la combinacion de un retorno anual esperado (drift) y una volatilidad anualizada relativamente altos, lo cual amplifica el crecimiento acumulado en el tiempo.

Estos valores extremos distorsionan el analisis general de la distribucion de precios simulados. Por esta razon, se decide excluir las trayectorias por encima del percentil 99% para representar de forma mas realista el comportamiento esperado del precio.

Extrae los precios finales

final_prices <- mb1[nrow(mb1), ]

Umbral del percentil 99

p99 <- quantile(final_prices, 0.99)

Filtramos trayectorias que terminan debajo del p99

filtered_mb1 <- mb1[, final_prices <= p99]

Calculamos percentiles a lo largo del tiempo

p_mat_nu <- apply(filtered_mb1, 1, function(x) quantile(x, probs = c(0.05, 0.5, 0.95)))

Transponemos para graficar

p_mat_nu <- t(p_mat_nu)

Grafica de trayectorias filtradas

matplot(filtered_mb1[, 1:100], type = "l", col = "gray", lty = 1, 
        main = "Simulaciones NU (trayectorias filtradas)", ylab = "Precio", xlab = "Días")
# Agrega lineas de percentiles
lines(p_mat_nu[, 1], col = "blue", lwd = 2, lty = 2)  # P5
lines(p_mat_nu[, 2], col = "black", lwd = 2)          # P50 (mediana)
lines(p_mat_nu[, 3], col = "red", lwd = 2, lty = 2)   # P95
legend("topleft", legend = c("P5", "P50", "P95"),
       col = c("blue", "black", "red"), lty = c(2,1,2), lwd = 2, bty = "n")

P5 (azul punteado): Muestra el límite inferior de los precios simulados, interpretado como una zona de riesgo de pérdida significativa. Es decir, un 5% de los caminos se encuentran por debajo de esa línea en cada día. Este trazo representa el VaR visual diario. P50 (negro sólido): Representa la mediana de las simulaciones. Es la mejor estimación no sesgada del precio futuro. Da una visión realista del comportamiento central. P95 (rojo punteado): Es el límite superior, marcando el rango superior del 5% más optimista. Ideal para ver zonas potenciales de ganancia extraordinaria.

SIMULACION ONON

onon_data <- price_data %>%
  filter(symbol == "ONON") %>%
  select(date, adjusted)

price_onon=onon_data[,2]
log_ret_onon <- coredata(log_ret_xts[, "ONON"])[-1]

plot.ts(price_onon)

plot.ts(log_ret_onon)

factores iniciales

year2=252
s02=as.numeric(onon_data[NROW(onon_data),2])
T2=504
r2=exp(mean(log_ret_onon)*year2)-1
sigma2= as.numeric(sqrt(var(log_ret_onon)*year2))
N2=5000 #numero de simulaciones
dt2=1/year2
stx2=as.numeric(s02*exp(r2*(T2/year2)))

vec2=rep(s02,N2)
mb2=matrix(ncol=N2, nrow=T2)
mb2[1,]=vec2

for (j in 1:N2) {
  for (u in 2:T2) {
    mb2[u,j]=mb2[(u-1),j]*exp((r2-(0.5*(sigma2^2)))*(dt2)+sigma2*(dt2^(1/2))*qnorm(runif(1, min=0, max=1)))                                                                      
  }
}

validación antes de grafico

print(paste("ultimo precio de ONON:", s02))

## [1] "ultimo precio de ONON: 39.6100006103516"

print(paste("Volatilidad anualizada (sigma2):", sigma2))

## [1] "Volatilidad anualizada (sigma2): 0.522581747655572"

print(paste("Rendimiento anual esperado (r2):", r2))

## [1] "Rendimiento anual esperado (r2): 0.248898893896643"

final_prices2 <- mb2[nrow(mb2), ]
quantile(final_prices2, probs = c(0.01, 0.05, 0.5, 0.95, 0.99))

##         1%         5%        50%        95%        99% 
##   9.065471  14.491122  49.448533 168.537093 277.576746

max_value2 <- max(mb2)
print(max_value2)

## [1] 902.331

which_max2 <- which(mb2 == max_value2, arr.ind = TRUE)
print(which_max2)

##      row  col
## [1,] 504 2667

###ANALISIS: Dado que bajo el Movimiento Browniano Geométrico (MBG) los precios se simulan de forma exponencial, se observa que aproximadamente el 5% de las trayectorias presentan comportamientos explosivos. Esto se debe a la combinacion de un retorno anual esperado (drift) y una volatilidad anualizada relativamente altos, lo cual amplifica el crecimiento acumulado en el tiempo.

Estos valores extremos distorsionan el analisis general de la distribucion de precios simulados. Por esta razon, se decide excluir las trayectorias por encima del percentil 95% para representar de forma mas realista el comportamiento esperado del precio.

Extraemos los precios finales

final_prices2 <- mb2[nrow(mb2), ]

Umbral del percentil 95

p95 <- quantile(final_prices2, 0.95)

Filtramos trayectorias que terminan debajo del p95

filtered_mb2 <- mb2[, final_prices2 <= p95]

Calculamos percentiles a lo largo del tiempo

p_mat_onon <- apply(filtered_mb2, 1, function(x) quantile(x, probs = c(0.05, 0.5, 0.95)))

Transponemos para graficar

p_mat_onon <- t(p_mat_onon)

Graficamos las trayectorias filtradas

matplot(filtered_mb2[, 1:100], type = "l", col = "gray", lty = 1, 
        main = "Simulaciones ONON (trayectorias filtradas)", ylab = "Precio", xlab = "Dias")
# Agrega lineas de percentiles
lines(p_mat_onon[, 1], col = "blue", lwd = 2, lty = 2)  # P5
lines(p_mat_onon[, 2], col = "black", lwd = 2)          # P50 (mediana)
lines(p_mat_onon[, 3], col = "red", lwd = 2, lty = 2)   # P95
legend("topleft", legend = c("P5", "P50", "P95"),
       col = c("blue", "black", "red"), lty = c(2,1,2), lwd = 2, bty = "n")

SIMULACION NVDA

nvda_data <- price_data %>%
  filter(symbol == "NVDA") %>%
  select(date, adjusted)

price_nvda=nvda_data[,2]
log_ret_nvda <- coredata(log_ret_xts[, "NVDA"])[-1]

plot.ts(price_nvda)

plot.ts(log_ret_nvda)

factores iniciales

year3=252
s03=as.numeric(nvda_data[NROW(nvda_data),2])
T3=504
r3=exp(mean(log_ret_nvda)*year3)-1
sigma3=as.numeric(sqrt(var(log_ret_nvda)*year3))
r3 <- 0.4
N3=5000 #numero de simulaciones
dt3=1/year3
stx3=as.numeric(s03*exp(r3*(T3/year3)))

vec3=rep(s03,N3)
mb3=matrix(ncol=N3, nrow=T3)
mb3[1,]=vec3

for (K in 1:N3) {
  for (V in 2:T3) {
    mb3[V,K]=mb3[(V-1),K]*exp((r3-(0.5*(sigma3^2)))*(dt3)+sigma3*(dt3^(1/2))*qnorm(runif(1, min=0, max=1)))                                                                      
  }
}

validación antes de grafico

print(paste("Ultimo precio de NVDA:", s03))

## [1] "Ultimo precio de NVDA: 94.3099975585938"

print(paste("Volatilidad anualizada (sigma3):", sigma3))

## [1] "Volatilidad anualizada (sigma3): 0.539449779245821"

print(paste("Rendimiento anual esperado (r3):", r3))

## [1] "Rendimiento anual esperado (r3): 0.4"

final_prices3 <- mb3[nrow(mb3), ]
quantile(final_prices3, probs = c(0.01, 0.05, 0.5, 0.95, 0.99))

##       1%       5%      50%      95%      99% 
##  26.2168  45.3470 153.7164 534.8470 931.4145

max_value3 <- max(mb3)
print(max_value3)

## [1] 2585.748

which_max3 <- which(mb3 == max_value3, arr.ind = TRUE)
print(which_max3)

##      row col
## [1,] 469 546

Extraemos los precios finales

final_prices3 <- mb3[nrow(mb3), ]

Umbral del percentil 95

p95_2 <- quantile(final_prices3, 0.95)

Filtramos trayectorias que terminan debajo del p95

filtered_mb3 <- mb3[, final_prices3 <= p95_2]

Calculamos percentiles a lo largo del tiempo

p_mat_nvda <- apply(filtered_mb3, 1, function(x) quantile(x, probs = c(0.05, 0.5, 0.95)))

Transponemos para graficar

p_mat_nvda <- t(p_mat_nvda)

Grafica de las trayectorias filtradas

matplot(filtered_mb3[, 1:100], type = "l", col = "gray", lty = 1, 
        main = "Simulaciones NVDA (trayectorias filtradas)", ylab = "Precio", xlab = "Dias")
# Agrega lineas de percentiles
lines(p_mat_nvda[, 1], col = "blue", lwd = 2, lty = 2)  # P5
lines(p_mat_nvda[, 2], col = "black", lwd = 2)          # P50 (mediana)
lines(p_mat_nvda[, 3], col = "red", lwd = 2, lty = 2)   # P95
legend("topleft", legend = c("P5", "P50", "P95"),
       col = c("blue", "black", "red"), lty = c(2,1,2), lwd = 2, bty = "n")

COMPARACION PORTAFOLIO REAL VS PORTAFOLIOS SIMULADOS

Aplicamos una función para seleccionar una trayectoria decreciente

elegir_trayectoria_decreciente <- function(sim_matrix) {
  indices <- which(sim_matrix[nrow(sim_matrix), ] < sim_matrix[1, ])
  if (length(indices) == 0) stop("No hay trayectorias en caida")
  sim_matrix[, sample(indices, 1)]
}

Elegimos trayectorias en caida

trayectoria_NU <- elegir_trayectoria_decreciente(mb1)
trayectoria_ONON <- elegir_trayectoria_decreciente(mb2)
trayectoria_NVDA <- elegir_trayectoria_decreciente(mb3)

Seguimiento trismestral del portafolio Dias a observar: cada 63 dias desde el inicio hasta 504 (8 trimestres + dia inicial)

dias_seguimiento <- c(1, seq(63, 504, by = 63))

Extraemos cantidades para cada acción

cantidad_NU <- resumen_inicial %>% filter(Activo == "NU") %>% pull(Cantidad_Acciones)
cantidad_ONON <- resumen_inicial %>% filter(Activo == "ONON") %>% pull(Cantidad_Acciones)
cantidad_NVDA <- resumen_inicial %>% filter(Activo == "NVDA") %>% pull(Cantidad_Acciones)

Creamos la tabla para registrar evolución del portafolio

seguimiento_portafolio <- tibble(
  Trimestre = 0:(length(dias_seguimiento) - 1),
  Dia_Simulacion = dias_seguimiento,
  Cantidad_NU = cantidad_NU,
  Cantidad_ONON = cantidad_ONON,
  Cantidad_NVDA = cantidad_NVDA,
  Precio_NU = trayectoria_NU[dias_seguimiento],
  Precio_ONON = trayectoria_ONON[dias_seguimiento],
  Precio_NVDA = trayectoria_NVDA[dias_seguimiento],
  Valor_NU = Precio_NU * cantidad_NU,
  Valor_ONON = Precio_ONON * cantidad_ONON,
  Valor_NVDA = Precio_NVDA * cantidad_NVDA,
  Valor_Portafolio = Valor_NU + Valor_ONON + Valor_NVDA,
  Diferencia_Valor = Valor_Portafolio - valor_inicial,
)

Descargamos la tasa libre de riesgo de FRED (3-Month Treasury Bill: ^IRX o DGS3MO)

getSymbols("DGS3MO", src = "FRED")  # Tasa a 3 meses

## [1] "DGS3MO"

rf_data <- na.omit(DGS3MO)
rf_trimestral <- tail(as.numeric(rf_data), 1) / 100  # Ultimo valor y lo llevamos a decimal

Calculamos retornos del portafolio (log-return trimestral)

seguimiento_portafolio <- seguimiento_portafolio %>%
  mutate(
    Retorno = c(NA, diff(log(Valor_Portafolio))),
    Volatilidad = zoo::rollapply(Retorno, width = 2, FUN = sd, fill = NA, align = "right"),
    Sharpe = ifelse(!is.na(Retorno), (Retorno - (rf_trimestral / 4)) / Volatilidad, NA)
  )
print(seguimiento_portafolio)

## # A tibble: 9 × 16
##   Trimestre Dia_Simulacion Cantidad_NU Cantidad_ONON Cantidad_NVDA Precio_NU
##       <int>          <dbl>       <dbl>         <dbl>         <dbl>     <dbl>
## 1         0              1       35468          9240          3111      9.60
## 2         1             63       35468          9240          3111      9.83
## 3         2            126       35468          9240          3111      6.72
## 4         3            189       35468          9240          3111      5.09
## 5         4            252       35468          9240          3111      4.16
## 6         5            315       35468          9240          3111      2.40
## 7         6            378       35468          9240          3111      2.82
## 8         7            441       35468          9240          3111      2.62
## 9         8            504       35468          9240          3111      2.90
## # ℹ 10 more variables: Precio_ONON <dbl>, Precio_NVDA <dbl>, Valor_NU <dbl>,
## #   Valor_ONON <dbl>, Valor_NVDA <dbl>, Valor_Portafolio <dbl>,
## #   Diferencia_Valor <dbl>, Retorno <dbl>, Volatilidad <dbl>, Sharpe <dbl>

VaR 1% y 5%

Calculamos VaR historico (al 1% y 5%) basado en retorno Se usa rolling de 2 para cada trimestre

seguimiento_portafolio <- seguimiento_portafolio %>%
  mutate(
    VaR_1 = zoo::rollapply(Retorno, width = 2, FUN = function(x) quantile(x, 0.01, na.rm = TRUE), fill = NA, align = "right"),
    VaR_5 = zoo::rollapply(Retorno, width = 2, FUN = function(x) quantile(x, 0.05, na.rm = TRUE), fill = NA, align = "right"),
    VaR_1 = pmin(VaR_1, 0),
    VaR_5 = pmin(VaR_5, 0),
    VVaR_1_USD = VaR_1 * lag(Valor_Portafolio),
    VaR_5_USD = VaR_5 * lag(Valor_Portafolio)
  )

Resultado final con todas las metricas por trimestre

print(seguimiento_portafolio)

## # A tibble: 9 × 20
##   Trimestre Dia_Simulacion Cantidad_NU Cantidad_ONON Cantidad_NVDA Precio_NU
##       <int>          <dbl>       <dbl>         <dbl>         <dbl>     <dbl>
## 1         0              1       35468          9240          3111      9.60
## 2         1             63       35468          9240          3111      9.83
## 3         2            126       35468          9240          3111      6.72
## 4         3            189       35468          9240          3111      5.09
## 5         4            252       35468          9240          3111      4.16
## 6         5            315       35468          9240          3111      2.40
## 7         6            378       35468          9240          3111      2.82
## 8         7            441       35468          9240          3111      2.62
## 9         8            504       35468          9240          3111      2.90
## # ℹ 14 more variables: Precio_ONON <dbl>, Precio_NVDA <dbl>, Valor_NU <dbl>,
## #   Valor_ONON <dbl>, Valor_NVDA <dbl>, Valor_Portafolio <dbl>,
## #   Diferencia_Valor <dbl>, Retorno <dbl>, Volatilidad <dbl>, Sharpe <dbl>,
## #   VaR_1 <dbl>, VaR_5 <dbl>, VVaR_1_USD <dbl>, VaR_5_USD <dbl>

COBERTURA

Ya que las opciones son una herramienta poderosa para cubrirse de los riesgos asociados, especialmente si el mercado es volátil o incierto, como en nuestro caso donde las acciones tienen una tendencia bajista, si bien no eliminan completamente el riesgo, proporcionan una forma efectiva de gestionarlo, lo que puede ser crucial para evitar pérdidas significativas. Sin embargo, es importante recordar que las opciones también implican costos (como la prima de la opción), por lo tanto se debe analizar cuales opciones realmente nos ayudan a compensar más las pérdidas.

Para nuestro análisis de cobertura con opciones elegimos realizarlo unicamente con NVDA ya que evaluando los valores de las opciones, y sus primas en el mercado para las 3 acciones elegidas (NVDA-NU-ONON), pudimos encontrar que NVDA presente una alta liquidez, multiples estrategias disponibles y un bajo costo; por el contrario para NU encontramos algunas opciones viables pero con spreads grandes; asi como para ONON no encontramos coberturas recomendables ya que hay pocos strikes activos.

Primero procedemos a definir un precio Strike

Precio actual de NVDA (ya calculado previamente)

nvda_price <- s03

URL de opciones para NVDA

url <- "https://finance.yahoo.com/quote/NVDA/options?p=NVDA"

Leer la página con cabecera de navegador (para evitar bloqueos)

page <- read_html(httr::GET(url, httr::add_headers(`User-Agent` = "Mozilla/5.0")))

Extraer la tabla de PUTs (segunda tabla)

puts_table <- page %>%
  html_elements("table") %>%
  .[[2]] %>%
  html_table(fill = TRUE)

Limpiar nombres de columnas

names(puts_table) <- gsub("\\s+", "_", names(puts_table))

Convertir columnas clave y calcular el spread

puts_clean <- puts_table %>%
  mutate(
    Strike = as.numeric(Strike),
    Bid = as.numeric(Bid),
    Ask = as.numeric(Ask),
    Last_Price = as.numeric(Last_Price),
    Implied_Volatility = as.numeric(gsub("%", "", Implied_Volatility)) / 100,
    spread = Ask - Bid
  ) %>%
  filter(
    !is.na(Strike), !is.na(Bid), !is.na(Ask),
    Strike > nvda_price  # <-- solo ITM
  ) %>%
  arrange(spread)

Convertimos también Open Interest (por si es character con comas)

puts_clean <- puts_clean %>%
  mutate(
    Open_Interest = as.numeric(gsub(",", "", Open_Interest)),
    # Score que premia alto interés abierto y bajo spread
    score = Open_Interest / (spread + 1e-6)  # evitamos división por 0
  ) %>%
  arrange(desc(score))  # Mejor score = mayor liquidez y menor spread

Seleccionamos la mejor opción

best_put <- puts_clean %>% slice(1)
selected_strike <- best_put$Strike
selected_prima <- best_put$Last_Price

Mostramos resultados

glue::glue("Strike seleccionado (PUT ITM con alto interés y bajo spread): {selected_strike}")

## Strike seleccionado (PUT ITM con alto interés y bajo spread): 101

print(best_put)

## # A tibble: 1 × 13
##   Contract_Name      Last_Trade_Date_(EDT…¹ Strike Last_Price   Bid   Ask Change
##   <chr>              <chr>                   <dbl>      <dbl> <dbl> <dbl>  <dbl>
## 1 NVDA250411P001010… 4/4/2025  3:59 PM         101       8.25  6.55  6.55      5
## # ℹ abbreviated name: ¹`Last_Trade_Date_(EDT)`
## # ℹ 6 more variables: `%_Change` <chr>, Volume <chr>, Open_Interest <dbl>,
## #   Implied_Volatility <dbl>, spread <dbl>, score <dbl>

ARBOL BINOMIAL

PARÁMETROS INICIALES

precio_inicial <- s03       # Precio actual de NVDA
strike <- selected_strike   # Precio de ejercicio (strike)
horizonte <- 2              # Horizonte de tiempo (años)
pasos <- 8                  # Número de pasos (8 trimestres)
volatilidad <- sigma3       # Volatilidad anual
tasa_rf <- rf_trimestral    # Tasa libre de riesgo anual

CÁLCULOS BASE

delta_t <- horizonte / pasos
sube <- exp(volatilidad * sqrt(delta_t))
baja <- 1 / sube
probabilidad <- (exp(tasa_rf * delta_t) - baja) / (sube - baja)

Árbol de precios del subyacente

precios_arbol <- matrix(0, nrow = pasos+1, ncol = pasos+1)
for (columna in 0:pasos) {
  for (fila in 0:columna) {
    precios_arbol[fila+1, columna+1] <- precio_inicial * sube^(columna-fila) * baja^fila
  }
}

Árbol de valoración de la put

valor_opcion <- matrix(0, nrow = pasos+1, ncol = pasos+1)
valor_opcion[, pasos+1] <- pmax(strike - precios_arbol[, pasos+1], 0)

Valoración hacia atrás

for (columna in (pasos-1):0) {
  for (fila in 0:columna) {
    valor_opcion[fila+1, columna+1] <- exp(-tasa_rf * delta_t) * (
      probabilidad * valor_opcion[fila+1, columna+2] + 
        (1 - probabilidad) * valor_opcion[fila+2, columna+2]
    )
  }
}

Preparamos DataFrame para graficar

datos_grafico <- data.frame()
for (columna in 0:pasos) {
  for (fila in 0:columna) {
    datos_grafico <- rbind(datos_grafico, data.frame(
      Paso = columna,
      Nodo = fila,
      Precio = precios_arbol[fila+1, columna+1],
      Valor = round(valor_opcion[fila+1, columna+1], 2)
    ))
  }
}

Ajustamos posiciones para graficar tipo árbol

datos_grafico <- datos_grafico %>%
  mutate(x = Paso,
         y = -Nodo + Paso / 2)

Graficamos árbol binomial

ggplot(datos_grafico, aes(x = x, y = y)) +
  geom_point(aes(color = Valor > 0), size = 5) +
  geom_text(aes(label = paste0("S=", round(Precio,1), "\nP=", Valor)), size = 3) +
  scale_color_manual(values = c("FALSE" = "gray70", "TRUE" = "red")) +
  labs(title = "Árbol Binomial de una Put Europea sobre NVDA",
       subtitle = "Precios del activo y valor de la opción en cada nodo",
       x = "Paso (Trimestre)", y = "Nodo") +
  theme_minimal() +
  theme(legend.position = "none")

ANALISIS COBERTURA

Parámetros

porcentaje_cubierto <- 0.85
strike_price <- selected_strike  # Asumimos que ya está definida
num_opciones <- ceiling(porcentaje_cubierto * seguimiento_portafolio$Cantidad_NVDA[1] / 100)
costo_prima <- num_opciones * selected_prima * 100

Función para encontrar valor de opción desde el árbol en un trimestre y precio dado

obtener_valor_opcion_desde_arbol <- function(trimestre, precio_simulado, precios_arbol, valor_opcion) {
  if (trimestre == 0) return(0)  # No hay cobertura en t=0
  # Extraer precios y valores para el trimestre correspondiente
  precios_trimestre <- precios_arbol[1:(trimestre+1), trimestre+1]
  valores_trimestre <- valor_opcion[1:(trimestre+1), trimestre+1]
  # Buscar el nodo más cercano al precio simulado
  idx_cercano <- which.min(abs(precios_trimestre - precio_simulado))
  return(valores_trimestre[idx_cercano])
}

Aplicamos la función por cada fila

seguimiento_portafolio$Valor_Put <- mapply(
  obtener_valor_opcion_desde_arbol,
  seguimiento_portafolio$Trimestre,
  seguimiento_portafolio$Precio_NVDA,
  MoreArgs = list(precios_arbol = precios_arbol, valor_opcion = valor_opcion)
)

Calculamos el ingreso por cobertura (valor PUT * 100 * opciones) por trimestre

seguimiento_portafolio$Ingreso_Cobertura_NVDA <- seguimiento_portafolio$Valor_Put * num_opciones * 100 - costo_prima
seguimiento_portafolio$Ingreso_Cobertura_NVDA[1]=0
seguimiento_portafolio$Costo_Prima[1]=0

## Warning: Unknown or uninitialised column: `Costo_Prima`.

Asi como también encontramos el valor del portafolio con cobertura

seguimiento_portafolio$Valor_Portafolio_Cubierto <- seguimiento_portafolio$Valor_Portafolio +
  seguimiento_portafolio$Ingreso_Cobertura_NVDA - cumsum(seguimiento_portafolio$Costo_Prima)

Las Diferencias respecto al valor inicial

valor_inicial <- seguimiento_portafolio$Valor_Portafolio[1]
seguimiento_portafolio$Diferencia_Sin_Cobertura <- seguimiento_portafolio$Valor_Portafolio - valor_inicial
seguimiento_portafolio$Diferencia_Con_Cobertura <- seguimiento_portafolio$Valor_Portafolio_Cubierto - valor_inicial

Gráfico

Convertimos el dataframe a formato largo para ggplot

df_plot <- seguimiento_portafolio[, c("Trimestre", "Valor_Portafolio", "Valor_Portafolio_Cubierto")]
df_plot <- reshape2::melt(df_plot, id.vars = "Trimestre")

Cambiar nombres para presentación

df_plot$variable <- factor(df_plot$variable,
                           levels = c("Valor_Portafolio", "Valor_Portafolio_Cubierto"),
                           labels = c("Sin Cobertura", "Con Cobertura"))

Graficar

ggplot(df_plot, aes(x = factor(Trimestre), y = value, fill = variable)) +
  geom_bar(stat = "identity", position = "dodge", color = "black") +
  labs(
    title = "Valor del Portafolio por Trimestre: Con y Sin Cobertura",
    x = "Trimestre",
    y = "Valor del Portafolio (USD)",
    fill = "Estrategia"
  ) +
  scale_fill_manual(values = c("steelblue", "darkorange")) +
  theme_minimal(base_size = 14) +
  theme(
    legend.position = "top",
    plot.title = element_text(face = "bold", hjust = 0.5)
  ) +
  geom_text(aes(label = round(value, 0)), 
            position = position_dodge(width = 0.9), 
            vjust = -0.5, size = 4)

A partir del trimestre 1 en adelante, el portafolio con cobertura mantiene consistentemente un valor mayor que el portafolio sin cobertura, lo que indica que la estrategia de cobertura si es efectiva para mitigar las pérdidas, especialmente durante los trimestres en que el valor del portafolio cae más fuertemente como en los trimestres 5 y 6. Asi mismo,la cobertura contribuye a reducir la volatilidad del portafolio, proporcionando mayor estabilidad en escenarios adversos.

Finalmente, esta gráfica justifica el uso de estrategias de cobertura dinámica, especialmente en portafolios expuestos a activos volátiles como acciones tecnológicas, para nuestra caso NVIDIA. Además, se evidencia que cuando aplicamos la cobertura mejoramos el perfil riesgo-retorno, ya que se respaldar cuantitativamente con indicadores como el índice de Sharpe y el VaR.

Creación de Arboles binomiales para cobertura de un portafolio de inversión

Realizado por Cesar Zapata - Oscar Ramirez - John Cano

6/4/2025

Introducción

Instalación y carga de paquetes

Descarga y conversión de la información de las acciones

Indice de Sharpe

Simulación de 5000 carteras aleatorias

La cartera de mínima varianza

ANALISIS:

SIMULACIÓN NU

ANALISIS:

SIMULACION ONON

SIMULACION NVDA

COMPARACION PORTAFOLIO REAL VS PORTAFOLIOS SIMULADOS

VaR 1% y 5%

COBERTURA

ARBOL BINOMIAL

ANALISIS COBERTURA

Gráfico