Se eliminan los mensajes de precaución y de descarga de paquetes.

options(warn = -1)
suppressPackageStartupMessages({
  library(knitr)
  library(plotly)
  library(dplyr)
  library(ggplot2)
  library(plotly)
  library(moments)
  library(magrittr)
  library(stats)
  library(tidyverse)
})

Se analiza la serie diaria TRM desde 2000-01-31, hasta 2022-09-17

Datos_TRM = read.table("C:/Users/57314/Downloads/DB_TRM.txt", header = TRUE)

# Guardar el archivo como un archivo HTML
library(knitr)

# Convertir la tabla a HTML
html_TRM <- kable(Datos_TRM, format = "html")

# Guardar el código HTML en un archivo
write(html_TRM, "tabla_trm.html")

head(Datos_TRM)
##        Fecha     TRM
## 1 01/01/2000 1873.77
## 2 02/01/2000 1873.77
## 3 03/01/2000 1873.77
## 4 04/01/2000 1874.35
## 5 05/01/2000 1895.97
## 6 06/01/2000 1912.69

1. Primera Parte

1.1. Estadística Descriptiva de los Precios

install.packages("moments")
library(moments)
library(tidyverse)
library(magrittr)

TRM = Datos_TRM[,2]
resumen_trm <- Datos_TRM %>%
  mutate(media_trm = mean(TRM),
         mediana_trm = median(TRM),
         desviacion_estandar_trm = sd(TRM),
         rango_trm = max(TRM) - min(TRM),
         varianza_trm = var(TRM),
         coeficiente_variacion_trm = sd(TRM) / mean(TRM),
         rendimiento_anualizado = (mean(TRM) - 1) * 100,
         volatilidad_anual = sd(TRM) * sqrt(252),
         cuantiles = quantile(TRM, probs = 0.25:0.75),
         sesgo = skewness(TRM),
         curtosis = kurtosis(TRM),
         minimo = min(TRM),
         maximo = max(TRM)
  )

# Seleccionar solo las columnas que no son Fecha ni TRM
estadistica_trm <- resumen_trm %>%
  select(-c(Fecha, TRM))

fila_unica <- do.call(c, estadistica_trm)

  print(head(estadistica_trm, n = 1))
##   media_trm mediana_trm desviacion_estandar_trm rango_trm varianza_trm
## 1  2553.402     2345.47                 643.277   2975.05     413805.3
##   coeficiente_variacion_trm rendimiento_anualizado volatilidad_anual cuantiles
## 1                 0.2519294               255240.2          10211.71   1974.47
##       sesgo curtosis  minimo  maximo
## 1 0.7478894 2.685184 1652.41 4627.46
rendimientos_diarios <- diff(log(Datos_TRM$TRM))

Análisis Descriptivo

• Media: La TRM promedio durante el período 2000-2022 fue de $2660.966 COP por USD.

• Mediana: El valor de la TRM que representa el valor “del medio” en un conjunto de datos ordenados durante el período 2000-2022 fue de $2389.75 COP por USD.

• Desviación estándar: La desviación estándar de la TRM fue de $762.4161 COP por USD, lo que significa que en promedio el precio de la TRM puede subir o bajar $762.4161 en un día.

• Máximo: El valor máximo de la TRM fue de $5061.21 COP por USD.

• Mínimo: El valor mínimo de la TRM fue de $1652.41 COP por USD.

Análisis de Tendencias:

• La tendencia de la TRM ha sido creciente durante el período de tiempo analizado.

• El aumento de la TRM se ha visto afectada por varios factores tanto internos como externos:

Factores externos

• La economía estadounidense ha experimentado un crecimiento durante estos años, acompañado de que la Reserva Federal ha trabajado en que el dólar sea una moneda atractiva para invertir (aumento de las tasas de interés), por esto su precio frente a otras monedas ha aumentado.

• Debido a los diferentes conflictos en países productores de petróleo, su precio ha sido muy volátil, lo que puede conducir a una depreciación de la TRM.

• Los conflictos económicos globales como la guerra de Ucrania o la crisis financiera en USA del 2008, hacen que los inversores busquen activos de bajo riesgo para invertir como las divisas, por ejemplo, el dólar estadounidense.

• La pandemia de COVID-19 en 2020 generó incertidumbre económica global, lo que generó una depreciación del peso colombiano frente al dólar.

Factores internos

• En Colombia los bienes importados son mayores a los exportados, por esto existe un déficit en cuenta corriente, lo que genera un alza en la TRM.

• En los últimos años la inflación en Colombia ha sido superior a la de Estados Unidos esto hace que el peso colombiano pierda poder adquisitivo frente al dólar.

• Cuando el Banco de la República emite pesos para aumentar la oferta de dinero, puede conducir a una depreciación de la moneda.

Análisis de Autocorrelación:

La TRM tiene una autocorrelación positiva, lo que significa que los valores actuales de la TRM están correlacionados con los valores pasados.

• Sesgo: El sesgo de la TRM es de 0.8964715, lo que indica una ligera asimetría a la derecha.

• Curtosis:La curtosis de la TRM es de 3.040295, lo que indica una distribución ligeramente más aplanada que la normal.

#Graficar los datos

Se realiza la gráfica de la Serie Histórica TRM

library(ggplot2)
library(plotly)

if (!is.data.frame(Datos_TRM)) {
  stop("Datos_TRM no es un data frame. Por favor, proporcione datos de series de tiempo.")
} else {
  
  if (!("index" %in% names(Datos_TRM))) {
    Datos_TRM$index <- 1:nrow(Datos_TRM)
  }
}

ggplot_obj <- ggplot(Datos_TRM, aes(x = index, y = TRM)) +
  geom_line(color = "lightblue") +
  labs(title = "Serie Histórica TRM", x = "Fecha", y = "Valor") +
  theme_minimal() +
  theme(panel.border = element_blank())

p <- ggplotly(ggplot_obj)

p

Se realiza la gráfica de la Densidad Empírica de los Precios y de los Retornos

# Calcular la densidad empírica de los precios TRM
densidad_empirica <- density(TRM)

# Crear un data frame con los valores de la densidad
df_densidad <- data.frame(
  x = densidad_empirica$x,
  y = densidad_empirica$y
)

# Graficar la densidad empírica
ggplot_obj <- ggplot(data = df_densidad, aes(x = x, y = y)) +
  geom_line(color = "lightblue", size = 1.5) +
  geom_point(aes(tooltip = paste0("Valor TRM: ", x, "\nDensidad: ", round(y, 4))), color = "lightblue", size = 3) +
  geom_vline(xintercept = TRM[length(TRM)], color = "black", linetype = "dashed") +
  labs(title = "Densidad Empírica de Precios TRM",
       x = "Valor TRM",
       y = "Densidad") +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5))

plotly_densidad <- ggplotly(ggplot_obj)

plotly_densidad
# Calcular la densidad empírica de los retornos
densidad_retornos <- density(rendimientos_diarios)

# Crear un data frame con los valores de la densidad
df_densidad_retornos <- data.frame(
  x = densidad_retornos$x,
  y = densidad_retornos$y
)

ggplot_obj <- ggplot(data = df_densidad_retornos, aes(x = x, y = y)) +
  geom_line(color = "lightblue", size = 1.5) +
  geom_point(aes(tooltip = paste0("Retorno diario: ", x, "\nDensidad: ", round(y, 4))), color = "lightblue", size = 3) +
  geom_vline(xintercept = 0, color = "black", linetype = "dashed") +
  labs(title = "Densidad Empírica de Retornos Diarios",
       x = "Retorno diario",
       y = "Densidad") +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5))

plotly_densidad_retornos <- ggplotly(ggplot_obj)

plotly_densidad_retornos

Cuánto miden los precios y los retornos

length(TRM)
## [1] 8296
length(rendimientos_diarios)
## [1] 8295

Estadísticas Descriptivas de los Retornos

Media anualizada

Media_Rtos = mean(rendimientos_diarios)*252

Desviación Estandar anualizada

Desv_Rtos = sd(rendimientos_diarios)*sqrt(252)

Sesgo

Sesgo_Rtos = skewness(rendimientos_diarios)

Curtosis

Curtosis = kurtosis(rendimientos_diarios)

2. Simulación Movimiento Browniano Geométrico

Simulación de los precios a 5 meses

t = 5/12
delta_t = 1/252
N = t/delta_t
M = 10000
S = matrix(ncol = M, nrow = (N+1))
volatilidad = sd(diff(log(TRM)))
mu = mean(diff(log(TRM)))

Se toma el último valor de la serie de precios

S[1,] = 4435.84
S0 <- 4435.84

2.1. Se grafica la distribución empírica

# Simulación de la dinámica del precio
for (i in 1:M) {
  for (t in 2:(N + 1)) {
    S[t, i] <- S[(t - 1), i] * exp((Media_Rtos - Desv_Rtos^2 / 2) * delta_t + Desv_Rtos * sqrt(delta_t) * rnorm(1))
  }
}

matplot(S, type = "l",
        main = "Distribución Empírica",
        xlab = "Tiempo",
        ylab = "Precio")+
  theme_minimal()

## NULL

Movimiento Browniano Geométrico

 S <- matrix(NA, nrow = N + 1, ncol = M)
S[1, ] <- S0

for (i in 1:M) {
 for (j in 1:N) {
  epsilon <- rnorm(1)
  dW <- epsilon * sqrt(delta_t)
  S[j + 1, i] <- S[j, i] * exp((mu - 0.5 * volatilidad^2) * delta_t + volatilidad * dW)
 }
}

Proyeccion1 = matrix(nrow = (8296 + 105), ncol = 1)
Proyeccion1[1: 8296] = TRM
matplot(Proyeccion1, type="l", ylim = c(1500, 5000))
title("Proyección")
lines(Proyeccion1, col = "lightblue", type = "l")


Proyeccion2 = matrix(nrow = (8296 + 105), ncol = M)
Proyeccion2[8296:8401] = S
matlines(Proyeccion2, type="l")

2.2. Cálculo los cuantiles al 2.5% y 97.5% de los posibles precios

install.packages("stats")
library(stats)
alfa=0.05
q1 = quantile(S[(N+1),], (alfa/2))
q2 = quantile(S[(N+1),], (1-alfa/2))
q1
##     2.5% 
## 4403.987
q2
##    97.5% 
## 4469.147

2.3. Gráfico de los resultados

library(ggplot2)
library(plotly)

# Supongamos que Rtos_TRM contiene los rendimientos diarios
gg_qq <- ggplot(data.frame(rendimientos_diarios = rendimientos_diarios), aes(sample = rendimientos_diarios)) +
  stat_qq(color = "lightblue", size = 1.5, pch = 16) +
  labs(title = "Gráfico QQ Rendimientos Diarios TRM",
       x = "Cuantiles teóricos",
       y = "Cuantiles observados") +
  theme_minimal() +
  theme(plot.title = element_text(size = 14),
        axis.title = element_text(size = 12),
        axis.text = element_text(size = 10))

# Convierte la gráfica de ggplot a plotly
plotly_qq <- ggplotly(gg_qq)
plotly_qq

2.4. Análisis de los Precios Futuros

Se simuló la distribución empírica de los precios futuros del TRM colombiano utilizando el método MBG.Se define la matriz Proyeccion1 que contiene los valores históricos de la TRM y los 105 días de la proyección. También se define la matriz Proyeccion2 que contiene las simulaciones de la serie de precios a partir del día 8296 (actual).

• La distribución empírica de los precios futuros de la TRM presenta una forma asimétrica con una mayor probabilidad de valores superiores al precio actual.

• La volatilidad anualizada de la tasa de cambio es un factor importante que determina la dispersión de la distribución.

• La incertidumbre alrededor de la proyección es considerable, como se observa en la dispersión de las simulaciones.

3. Estimación de un intervalo de predicción del 95% sobre los posibles precios futuros de la TRM

mu = mean(S)
sigma = sd(S)

qnorm(c(0.05, 0.95), mu, sigma)
## [1] 4416.642 4455.481

2. segunda Parte

Contexto del problema

Una empresa exportadora que realizó crédito para ejercer una inversión empresarial el 17 de septiembre de 2022 con los siguientes criterios de pago por parte de su capitalizadora de inversión:

Monto de inversión: 1’790.000 USD Tasa de interés: 10% AMV Cuotas: 24 Frecuencia de pago: Mensual

fecha<- c("17/10/2022","17/11/2022","17/12/2022","17/01/2023","17/02/2023")
pagos<- c(91738,90992,90246,89500,88754)

tabla <- data.frame(fecha = fecha, pagos = pagos)

print(tabla)
##        fecha pagos
## 1 17/10/2022 91738
## 2 17/11/2022 90992
## 3 17/12/2022 90246
## 4 17/01/2023 89500
## 5 17/02/2023 88754

La tasa de cambio actual es de 4066 COP. Se sugiere una cobertura del 70% del total del monto de los flujos expuestos a riesgo cambiario. El intermediario financiero IF le ofrece contratos de futuros TRM__F (ver bolsa de valores de Colombia) sugiriendo el uso de los siguientes contratos que se encuentran en la bolsa de valores de Colombia.

2.1. Simulación

mu= mean(rendimientos_diarios)*252
sigma= sd(rendimientos_diarios)*sqrt(252)
T= 5/12 # 5 meses
delta_t=1/252
N=T/delta_t
M= 1000 # Numero de simulaciones
TRM= matrix(ncol=M, nrow = (N+1))
TRM[1,]=4435.84

for(i in 1:M){
  for (t in 2:(N+1)){
    TRM[t,i]=TRM[(t-1),i]*exp((mu-sigma^2/2)*delta_t+sigma*sqrt(delta_t)*rnorm(1))
  }
}

matplot(TRM, type = "l",
        main = "Simulación",
        xlab = "Tiempo",
        ylab = "Precio")+
  theme_minimal()

## NULL

2.2. Simulación de Montecarlo del Flujo de caja total esperado para cada cuota sin cobertura

Dada la proyección de la TRM, se hace una estimación del valor de las cuotas en COP para cada periodo.

Cuotas_USD = c(91738,90992,90246,89500,88754)
Cuotas_COP=TRM[2,]*Cuotas_USD[1]

for(j in 2:length(Cuotas_USD)){
  Cuotas_COP=rbind(Cuotas_COP,TRM[j+1,]*Cuotas_USD[j])
}

media_cuotas=vector()
volatilidad_cuotas=vector()
percentil_5_cuotas=vector()
percentil_95_cuotas=vector()

for(l in 1:5){
  media_cuotas[l]=mean(TRM[l+1,]*Cuotas_USD[l])
  volatilidad_cuotas[l]=sd(TRM[l+1,]*Cuotas_USD[l])
  percentil_5_cuotas[l]=quantile(TRM[l+1,]*Cuotas_USD[l],0.05)
  percentil_95_cuotas[l]=quantile(TRM[l+1,]*Cuotas_USD[l],0.95)
}

matplot(Cuotas_COP,type = "l", col="lightblue", main="Estimado Cuotas en COP", ylab = "Valor Cuota", xlab = "Cuota")
lines(percentil_5_cuotas, type = "l")
lines(percentil_95_cuotas, type = "l")
lines(media_cuotas, type = "l")

print("Volatilidad de cada cuota  en millones de COP")
## [1] "Volatilidad de cada cuota  en millones de COP"
print(volatilidad_cuotas/1000000)
## [1] 2.351191 3.270562 3.894329 4.475130 4.952336

2.3. Simulación de Montecarlo de los pagos de cuotas con cobertura

Dado que quien tiene una deuda en USD tiene una posición corta en la tasa de cambio, se debe de tomar una posición larga en los contratos Forward para hacer la cobertura.

Cobertura<-0.7
Precio_Entrega= 4500

Utilidad_POR_USD=TRM[-1,]-Precio_Entrega
matplot(Utilidad_POR_USD,type="l", col="lightblue", main="Utilidad por USD", ylab = "Utilidad (USD)", xlab = "Fecha")

Utilidad_total=Utilidad_POR_USD[1,]*Cuotas_USD[1]*0.7

for(j in 2:length(Cuotas_USD)){
  Utilidad_total=rbind(Utilidad_total,Utilidad_POR_USD[j,]*Cuotas_USD[j]*0.7)
}

matplot(Utilidad_total,type="l", col = "lightblue", main="Utilidad Total", ylab = "Utilidad (USD)", xlab = "Fecha")

Cuotas_COP_CON_C=Cuotas_COP-Utilidad_total

media_cuotas_COP_CON_C=vector()
volatilidad_cuotas_COP_CON_C=vector()
percentil_5_cuotas_COP_CON_C=vector()
percentil_95_cuotas_COP_CON_C=vector()

for(l in 1:5){
  media_cuotas_COP_CON_C[l]=mean(Cuotas_COP_CON_C[l,])
  volatilidad_cuotas_COP_CON_C[l]=sd(Cuotas_COP_CON_C[l,])
  percentil_5_cuotas_COP_CON_C[l]=quantile(Cuotas_COP_CON_C[l,],0.05)
  percentil_95_cuotas_COP_CON_C[l]=quantile(Cuotas_COP_CON_C[l,],0.95)
}

matplot(Cuotas_COP_CON_C,type = "l", col="lightblue", main="Estimado Cuotas en COP", ylab = "Valor Cuota", xlab = "cuota")
lines(percentil_5_cuotas_COP_CON_C, type = "l")
lines(percentil_95_cuotas_COP_CON_C, type = "l")
lines(media_cuotas_COP_CON_C, type = "l")

print("Volatilidad de cada cuota con cobertura  en millones de COP")
## [1] "Volatilidad de cada cuota con cobertura  en millones de COP"
print(volatilidad_cuotas_COP_CON_C/1000000)
## [1] 0.7053574 0.9811685 1.1682987 1.3425389 1.4857008

2.4. Beneficios reales de las series

2.4.1. Beneficio real sin cobertura

fecha <- c("17/10/2022", "17/11/2022", "17/12/2022", "17/01/2023", "17/02/2023")
valores_TRM <- c(4636.83, 4922.70, 4802.48, 4693.99, 4966.33)

# Crear un vector vacío para almacenar los valores de TRM_OBSERVADA
TRM_OBSERVADA <- numeric(length(fecha))

# Asignar los valores de TRM_OBSERVADA
for (j in 1:length(fecha)) {
  TRM_OBSERVADA[j] <- valores_TRM[j]
}

print(TRM_OBSERVADA)
## [1] 4636.83 4922.70 4802.48 4693.99 4966.33
Cuota_Real_COP=TRM_OBSERVADA*Cuotas_USD

matplot(Cuotas_COP,type = "l", col="lightblue", main="Estimado Cuotas en COP", ylab = "Valor Cuota", xlab = "Cuota")
lines(percentil_5_cuotas, type = "l")
lines(percentil_95_cuotas, type = "l")
lines(media_cuotas, type = "l")
lines(Cuota_Real_COP, type = "l",lwd=2)

"Cuota Real COP"
## [1] "Cuota Real COP"
print(Cuota_Real_COP)
## [1] 425373511 447926318 433404610 420112105 440781653
"Beneficio sin cobertura"
## [1] "Beneficio sin cobertura"
Utilidad_POR_USD_Observada=TRM_OBSERVADA-Precio_Entrega
Utilidad_total_observada= Cuotas_USD*Utilidad_POR_USD_Observada

2.4.2. Beneficio real con cobertura

Valor_Observado_Cuotas_con_cobertura=Cuota_Real_COP-Utilidad_total_observada

matplot(Cuotas_COP_CON_C,type = "l", col="lightblue", main="Estimado Cuotas en COP", ylab = "Valor Cuota", xlab = "Cuota")
lines(percentil_5_cuotas_COP_CON_C, type = "l")
lines(percentil_95_cuotas_COP_CON_C, type = "l")
lines(media_cuotas_COP_CON_C, type = "l")
lines(Valor_Observado_Cuotas_con_cobertura, type = "l", lwd=2)

2.5. Cuenta de margen con el flujo de caja operado

Se deben de definir las variables y las fórmulas a usar:

Numero_de_contratos <- 1
Nominal_del_contrato <- 1790000*0.7
Precio_inicial <- 4066 # TRM al momento 0
Garantia_inicial <- 0.063 #Futuros tasa de cambio BVC
Margen_mantenimiento <- 0.50

Exposicion_Total =Numero_de_contratos*Nominal_del_contrato*Precio_inicial
Valor_de_la_garantia_inicial = Exposicion_Total*Garantia_inicial
Valor_de_la_garantia_minima = Margen_mantenimiento*Valor_de_la_garantia_inicial
Apalancamiento = Exposicion_Total/Valor_de_la_garantia_inicial

Se liquidan los futuros y se nombran las columnas

TRM_Filas <- TRM[c(22, 43, 64, 85, 106), ]

Media_Filas <- rowMeans(TRM_Filas)

Media_TRM <- data.frame(Media_Futuro = Media_Filas)

Tabla_liquidacion <- Media_TRM
Media_TRM$fecha <- fecha
Media_TRM$TRM_OBSERVADA <- TRM_OBSERVADA
Media_TRM <- Media_TRM[, c("fecha", "TRM_OBSERVADA", "Media_Futuro")]
Tabla_liquidacion <- Media_TRM
Dif_Fut_TRM <- diff(Tabla_liquidacion[,3])*-1
Dif_Fut_TRM <- append(0,Dif_Fut_TRM)
Tabla_liquidacion <- cbind(Tabla_liquidacion,(Dif_Fut_TRM*Numero_de_contratos*Nominal_del_contrato))

colnames(Tabla_liquidacion) <- c("Fecha","TRM","Futuro","Liquidacion_Diaria")

Se calcula el llamado al margen

margen <- matrix()
tabaux <- matrix()
Llamado_al_margen <- matrix()

for (i in 1:(nrow(Tabla_liquidacion))) {
  # Tabla Auxiliar
  tabaux[1]=Valor_de_la_garantia_inicial
  tabaux[i+1] <- tabaux[i] + Tabla_liquidacion$Liquidacion_Diaria[i+1]
  # Llamado al Margen
  Llamado_al_margen[i] <- if(tabaux[i]<Valor_de_la_garantia_minima){Valor_de_la_garantia_inicial-tabaux[i]}
  else{0}
  # Margen 
  margen[i] <- tabaux[i]+Llamado_al_margen[i]
}

margen <- na.omit(margen)

Tabla_liquidacion <- cbind(Tabla_liquidacion,margen,Llamado_al_margen,Dif_Fut_TRM)

print(Tabla_liquidacion)
##        Fecha     TRM   Futuro Liquidacion_Diaria    margen Llamado_al_margen
## 1 17/10/2022 4636.83 4445.645                  0 320965974                 0
## 2 17/11/2022 4922.70 4457.082          -14330559 306635415                 0
## 3 17/12/2022 4802.48 4470.761          -17139927 289495488                 0
## 4 17/01/2023 4693.99 4479.008          -10334000 279161489                 0
## 5 17/02/2023 4966.33 4491.262          -15354463 263807026                 0
##   Dif_Fut_TRM
## 1    0.000000
## 2  -11.436999
## 3  -13.679111
## 4   -8.247406
## 5  -12.254160

2.6. ¿Cuál es la probabilidad de ser llamado al margen si solo se deposita en la cuenta el margen inicial?

Se calcula el número de periodos en los que se realiza una llamada al margen y luego se Calcula la probabilidad de ser llamado al margen

llamadas_al_margen <- sum(Llamado_al_margen > 0)

probabilidad_llamado_al_margen <- llamadas_al_margen / nrow(Tabla_liquidacion)

cat("La probabilidad de ser llamado al margen es:", probabilidad_llamado_al_margen)
## La probabilidad de ser llamado al margen es: 0

La probabilidad de ser llamado al margen si solo se deposita en la cuenta el margen inicial es del 0

2.7. ¿Cuánto debería de haber en la cuenta de margen al inicio para que la probabilidad de ser llamado al margen antes de cubrir el primer flujo sea menor al 1%?

Valor_de_la_garantia_minima
## [1] 160482987
Valor_garantia_minima <- 160482987

# Probabilidad máxima deseada
probabilidad_maxima <- 0.01

# Calcular el saldo inicial necesario usando la distribución binomial inversa
saldo_inicial <- qbinom(probabilidad_maxima, size = nrow(Tabla_liquidacion), prob = 0.5) * (-Valor_garantia_minima)

print(paste("El saldo inicial requerido para una probabilidad menor al 1% es:", saldo_inicial))
## [1] "El saldo inicial requerido para una probabilidad menor al 1% es: 0"

La garantía inicial necesaria para que la probabilidad de ser llamado al margen antes de cubrir el primer flujo sea menor al 1% es 0

3. Tercera Parte: Modelo Vasicek

3.1. Tasas simuladas usando una tasa de bono del mercado americano

Se definen los parámetros del modelo

r0_FED <- 0.0387 #19/09/2022 Tasa semestral tomada de: https://home.treasury.gov/resource-center/data-chart-center/interest-rates/TextView?type=daily_treasury_yield_curve&field_tdr_date_value=2022
theta <- 0.08
k <- 0.4
SD_FED <- 0.0094 #Desviación estándar de 3 años, calculada desde el 19/09/2019 hasta el 19/09/2022

Se simulan las trayectorias de Short-Rate

n <- 10    # MC simulation trials
T <- 10    # total time
m <- 252   # subintervals
dt <- T/m  # difference in time each subinterval

r <- matrix(0,m+1,n)  # matrix to hold short rate paths
r[1,] <- r0_FED

for(j in 1:n){
  for(i in 2:(m+1)){
    dr <- k*(theta-r[i-1,j])*dt + SD_FED*sqrt(dt)*rnorm(1,0,1)
    r[i,j] <- r[i-1,j] + dr
  }
}

Se grafican los resultados

t <- seq(0, T, dt)
rT.expected <- theta + (r0_FED-theta)*exp(-k*t)
rT.stdev <- sqrt( SD_FED^2/(2*k)*(1-exp(-2*k*t)))
matplot(t, r[,1:10], type="l", lty=1, main="Short Rate Paths", ylab="rt") 
abline(h=theta, col="lightblue", lty=2)
lines(t, rT.expected, lty=2) 
lines(t, rT.expected + 2*rT.stdev, lty=2) 
lines(t, rT.expected - 2*rT.stdev, lty=2) 
points(0,r0_FED)

3.2. Curvas Yield usando una tasa de bono del mercado americano

Función para encontrar el precio ZCB usando el modelo Vasicek

VasicekZCBprice <- 
  function(r0_FED, k, theta, SD_FED, T){
    b.vas <- (1/k)*(1-exp(-T*k)) 
    a.vas <- (theta-SD_FED^2/(2*k^2))*(T-b.vas)+(SD_FED^2)/(4*k)*b.vas^2
    return(exp(-a.vas-b.vas*r0_FED))
  }

Se definen los parámetros del modelo para trazar curvas de rendimiento

r0 <- seq(0.00, 0.20, 0.05)
n <- length(r0)
yield <- matrix(0, 10, n)
for(i in 1:n){
  for(T in 1:10){
    yield[T,i] <- -log(VasicekZCBprice(r0[i], k, theta, SD_FED, T))/T
  }
}

Se grafican los resultados

maturity <- seq(1, 10, 1)
matplot(maturity, yield, type="l", col="skyblue", lty=1, main="Yield Curve FED")
abline(h=theta, col="blue", lty=2)

3.3. Tasas simuladas usando la TES del mercado colombiano

Se definen los parámetros del modelo

r0_TES <- 0.0564 #19/09/2022 Tasa semestral tomada del Banco de la Rpública
theta <- 0.08
k <- 0.4
SD_TES <- 0.01274 #Desviación estándar de 3 años

Se simulan las trayectorias de Short-Rate

n <- 10    # MC simulation trials
T <- 10    # total time
m <- 252   # subintervals
dt <- T/m  # difference in time each subinterval

r <- matrix(0,m+1,n)  # matrix to hold short rate paths
r[1,] <- r0_TES

for(j in 1:n){
  for(i in 2:(m+1)){
    dr <- k*(theta-r[i-1,j])*dt + SD_TES*sqrt(dt)*rnorm(1,0,1)
    r[i,j] <- r[i-1,j] + dr
  }
}

Se grafican los resultados

t <- seq(0, T, dt)
rT.expected <- theta + (r0_TES-theta)*exp(-k*t)
rT.stdev <- sqrt( SD_TES^2/(2*k)*(1-exp(-2*k*t)))
matplot(t, r[,1:10], type="l", lty=1, main="Short Rate Paths", ylab="rt") 
abline(h=theta, col="lightblue", lty=2)
lines(t, rT.expected, lty=2) 
lines(t, rT.expected + 2*rT.stdev, lty=2) 
lines(t, rT.expected - 2*rT.stdev, lty=2) 
points(0,r0_TES)

3.4. Curvas Yield usando la TES del mercado colombiano

Función para encontrar el precio ZCB usando el modelo Vasicek

VasicekZCBprice <- 
  function(r0_TES, k, theta, SD_TES, T){
    b.vas <- (1/k)*(1-exp(-T*k)) 
    a.vas <- (theta-SD_TES^2/(2*k^2))*(T-b.vas)+(SD_TES^2)/(4*k)*b.vas^2
    return(exp(-a.vas-b.vas*r0_TES))
  }

Se definen los parámetros del modelo para trazar curvas de rendimiento

r0 <- seq(0.00, 0.20, 0.05)
n <- length(r0)
yield <- matrix(0, 10, n)
for(i in 1:n){
  for(T in 1:10){
    yield[T,i] <- -log(VasicekZCBprice(r0[i], k, theta, SD_TES, T))/T
  }
}

Se grafican los resultados

maturity <- seq(1, 10, 1)
matplot(maturity, yield, type="l", col="skyblue", lty=1, main="Yield Curve TES")
abline(h=theta, col="blue", lty=2)