BLACK LITTERMAN
Genjis Ossa
2026-05-03
library(tidyverse)## Warning: package 'tidyverse' was built under R version 4.4.3## Warning: package 'ggplot2' was built under R version 4.4.3## Warning: package 'tibble' was built under R version 4.4.3## Warning: package 'readr' was built under R version 4.4.3## Warning: package 'dplyr' was built under R version 4.4.3## Warning: package 'lubridate' was built under R version 4.4.2## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.5.2 ✔ tibble 3.3.0
## ✔ lubridate 1.9.4 ✔ tidyr 1.3.1
## ✔ purrr 1.0.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errorslibrary(quantmod)## Warning: package 'quantmod' was built under R version 4.4.3## Cargando paquete requerido: xts
## Cargando paquete requerido: zoo
##
## Adjuntando el paquete: 'zoo'
##
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
##
##
## ######################### Warning from 'xts' package ##########################
## # #
## # The dplyr lag() function breaks how base R's lag() function is supposed to #
## # work, which breaks lag(my_xts). Calls to lag(my_xts) that you type or #
## # source() into this session won't work correctly. #
## # #
## # Use stats::lag() to make sure you're not using dplyr::lag(), or you can add #
## # conflictRules('dplyr', exclude = 'lag') to your .Rprofile to stop #
## # dplyr from breaking base R's lag() function. #
## # #
## # Code in packages is not affected. It's protected by R's namespace mechanism #
## # Set `options(xts.warn_dplyr_breaks_lag = FALSE)` to suppress this warning. #
## # #
## ###############################################################################
##
## Adjuntando el paquete: 'xts'
##
## The following objects are masked from 'package:dplyr':
##
## first, last
##
## Cargando paquete requerido: TTR
## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoolibrary(DataExplorer)## Warning: package 'DataExplorer' was built under R version 4.4.3library(corrplot)## corrplot 0.92 loadedlibrary(scales)## Warning: package 'scales' was built under R version 4.4.3##
## Adjuntando el paquete: 'scales'
##
## The following object is masked from 'package:purrr':
##
## discard
##
## The following object is masked from 'package:readr':
##
## col_factorlibrary(ggplot2)
library(quadprog)Definición de activos del portafolio.
activos = c("VTI", "VEA", "VWO", "BND", "GLD")
fecha_inicio = "2021-01-01"
fecha_final = "2023-08-31"getSymbols(activos,from = fecha_inicio, to = fecha_final, auto.assign = TRUE)## [1] "VTI" "VEA" "VWO" "BND" "GLD"# extraer precios ajustados y convertirlos a frecuencia mensual
precios_mensuales = lapply(activos, function(activo) {
to.monthly(Ad(get(activo)),indexAt = "lastof",OHLC = FALSE)
})
precios_activos = do.call(cbind, precios_mensuales) ## una sola matriz de precios
colnames(precios_activos) = activos ## cambio de nombres de columnas
precios_activos = na.omit(precios_activos) # eliminar datos faltantesCalcular rendimientos mensuales
rendimientos_aritmeticos = diff(precios_activos) / lag(precios_activos)
rendimientos_aritmeticos = na.omit(rendimientos_aritmeticos) ## eliminar primera fila con NA
head(rendimientos_aritmeticos, n = 3)## VTI VEA VWO BND GLD
## 2021-02-28 0.03139328 0.02432264 0.015673616 -0.015476956 -0.06256881
## 2021-03-31 0.03647936 0.02768455 -0.007104038 -0.012713823 -0.01143311
## 2021-04-30 0.05036514 0.03054375 0.017867430 0.008659562 0.03563389Analisis exploratorio
plot_intro(as.data.frame(rendimientos_aritmeticos))
rendimientos_acumulados = cumprod(1 + rendimientos_aritmeticos)
autoplot(rendimientos_acumulados, facets = NULL) +
ggtitle("Rendimiento acumulado del universo de inversión") +
ylab("Rendimiento acumulado") +
xlab("Año")
summary(rendimientos_aritmeticos)## Index VTI VEA
## Min. :2021-02-28 Min. :-0.092305 Min. :-0.098600
## 1st Qu.:2021-10-15 1st Qu.:-0.024442 1st Qu.:-0.035534
## Median :2022-05-31 Median : 0.017368 Median : 0.013133
## Mean :2022-05-31 Mean : 0.007138 Mean : 0.002975
## 3rd Qu.:2023-01-15 3rd Qu.: 0.037283 3rd Qu.: 0.031841
## Max. :2023-08-31 Max. : 0.093468 Max. : 0.125486
## VWO BND GLD
## Min. :-0.100823 Min. :-0.041830 Min. :-0.071477
## 1st Qu.:-0.033086 1st Qu.:-0.014095 1st Qu.:-0.024048
## Median :-0.004598 Median :-0.003103 Median :-0.011132
## Mean :-0.004063 Mean :-0.004210 Mean : 0.002203
## 3rd Qu.: 0.016331 3rd Qu.: 0.007088 3rd Qu.: 0.027302
## Max. : 0.143019 Max. : 0.036659 Max. : 0.084919corrplot(
cor(rendimientos_aritmeticos),type = "lower",method = "shade", tl.col = "black",
cl.pos = "r",
tl.cex = 1)
## Definición de parametros
numero_activos = ncol(rendimientos_aritmeticos)
numero_observaciones = nrow(rendimientos_aritmeticos)
vector_unos = rep(1, numero_activos)
# Como los datos son mensuales, se anualizan con 12
periodicidad = 12
tasa_libre_riesgo = 0.04 ##risk freeMatriz de covarianzas
# Matriz de covarianzas anualizada
matriz_covarianzas = cov(rendimientos_aritmeticos) *
(numero_observaciones - 1) /
(numero_observaciones - numero_activos - 2) *
periodicidad
parametro_regularizacion = 0.0001 ## parametro de regularización
matriz_covarianzas = matriz_covarianzas +
diag(parametro_regularizacion, ncol(matriz_covarianzas))
# Retornos esperados anualizados en exceso de la tasa libre de riesgo
retornos_esperados = colMeans(rendimientos_aritmeticos) *
periodicidad -
tasa_libre_riesgo
retornos_esperados = as.matrix(retornos_esperados)Portafolio de retornos historicos
constante_historica = t(vector_unos) %*%
solve(matriz_covarianzas) %*%
retornos_esperados
pesos_historicos = as.numeric(1 / constante_historica) *
solve(matriz_covarianzas) %*%
retornos_esperados
retorno_portafolio_historico = t(retornos_esperados) %*%
pesos_historicos
print("Pesos del portafolio histórico:")## [1] "Pesos del portafolio histórico:"print(pesos_historicos)## [,1]
## VTI -0.34648917
## VEA -0.08973644
## VWO 0.03657776
## BND 1.64074384
## GLD -0.24109600print("Retorno esperado del portafolio histórico:")## [1] "Retorno esperado del portafolio histórico:"print(retorno_portafolio_historico)## [,1]
## [1,] -0.1639312Opiniones del inversionista
# Estas son opiniones subjetivas sobre el rendimiento esperado anual de cada activo.
## Entre mas se aleje de cero (0) mas positiva es.
opiniones_inversionista = numeric(numero_activos)
opiniones_inversionista[1] = 0.06
opiniones_inversionista[2] = 0.04
opiniones_inversionista[3] = 0.03
opiniones_inversionista[4] = 0.02
opiniones_inversionista[5] = 0.05
opiniones_inversionista = as.matrix(opiniones_inversionista)
print(opiniones_inversionista)## [,1]
## [1,] 0.06
## [2,] 0.04
## [3,] 0.03
## [4,] 0.02
## [5,] 0.05tau = 0.95 ## parámetro tau del modelo Black-Litterman
## matriz Omega que se asocia a incertidumbre de las opiniones
matriz_incertidumbre_opiniones = diag(
diag(matriz_covarianzas),
numero_activos,
numero_activos)
print(matriz_incertidumbre_opiniones)## [,1] [,2] [,3] [,4] [,5]
## [1,] 0.03879556 0.00000000 0.00000000 0.000000000 0.00000000
## [2,] 0.00000000 0.03888318 0.00000000 0.000000000 0.00000000
## [3,] 0.00000000 0.00000000 0.03368497 0.000000000 0.00000000
## [4,] 0.00000000 0.00000000 0.00000000 0.005557468 0.00000000
## [5,] 0.00000000 0.00000000 0.00000000 0.000000000 0.02497601Estimación Black Litterman
retornos_ajustados_BL = solve(
solve(tau * matriz_covarianzas) +
solve(matriz_incertidumbre_opiniones)
) %*% (
solve(tau * matriz_covarianzas) %*% retornos_esperados +
solve(matriz_incertidumbre_opiniones) %*% opiniones_inversionista
)
print("Retornos esperados ajustados por Black-Litterman:")## [1] "Retornos esperados ajustados por Black-Litterman:"print(retornos_ajustados_BL)## [,1]
## VTI 0.09972105
## VEA 0.07285422
## VWO -0.00132202
## BND -0.04112773
## GLD 0.04616245Portafolio con solo opiniones
constante_opiniones = t(vector_unos) %*%
solve(matriz_covarianzas) %*%
opiniones_inversionista
pesos_opiniones = as.numeric(1 / constante_opiniones) *
solve(matriz_covarianzas) %*%
opiniones_inversionista
print("Pesos del portafolio usando solo opiniones:")## [1] "Pesos del portafolio usando solo opiniones:"print(pesos_opiniones)## [,1]
## VTI 1.0919437
## VEA -1.1337302
## VWO 0.1431666
## BND 0.2023850
## GLD 0.6962349Optimización con pesos no negativos
matriz_D = 2 * matriz_covarianzas
vector_d = as.vector(retornos_ajustados_BL)
# Restricciones:
# 1. La suma de los pesos debe ser igual a 1
# 2. Cada peso debe ser mayor o igual a 0
matriz_restricciones = cbind(
rep(1, numero_activos),
diag(numero_activos)
)
vector_restricciones = c(1,rep(0, numero_activos))
resultado_optimizacion = solve.QP(
Dmat = matriz_D,
dvec = vector_d,
Amat = matriz_restricciones,
bvec = vector_restricciones,
meq = 1)
# Pesos finales del portafolio Black-Litterman
pesos_finales_BL = resultado_optimizacion$solution
names(pesos_finales_BL) = activos
print("Pesos finales del portafolio Black-Litterman:")## [1] "Pesos finales del portafolio Black-Litterman:"print(pesos_finales_BL)## VTI VEA VWO BND GLD
## 9.028157e-01 -9.398333e-17 5.725852e-18 -1.918577e-16 9.718425e-02Grafico
barplot(
pesos_finales_BL,
col = "blue",
names.arg = activos,
ylim = c(0, max(pesos_finales_BL) + 0.05),
main = "Pesos del portafolio Black-Litterman",
ylab = "Peso",
xlab = "Activos"
)
## Retorno y riesgo del portafolio final
retorno_portafolio_BL = t(pesos_finales_BL) %*%
retornos_ajustados_BL
riesgo_portafolio_BL = sqrt(
t(pesos_finales_BL) %*%
matriz_covarianzas %*%
pesos_finales_BL
)
print("Retorno esperado del portafolio Black-Litterman:")## [1] "Retorno esperado del portafolio Black-Litterman:"print(retorno_portafolio_BL)## [,1]
## [1,] 0.094516print("Riesgo esperado del portafolio Black-Litterman:")## [1] "Riesgo esperado del portafolio Black-Litterman:"print(riesgo_portafolio_BL)## [,1]
## [1,] 0.1820017