Cálculo de Black-Litterman en R - ¿dónde me equivoco?

Question

Estoy intentando calcular un pequeño modelo de Black Litterman en R. Estoy siguiendo un vídeo de Youtube y traduciendo la implementación de Excel en R.

Tengo una matriz var cov S

              INTC          AEP          AMZN          MRK           XOM        ^GSPC
INTC  0.0119535151 0.0005721887  0.0072352418 0.0016447926  0.0005925077 0.0024795274
AEP   0.0005721887 0.0042225253  0.0008231236 0.0011854049  0.0010758889 0.0011941026
AMZN  0.0072352418 0.0008231236  0.0191091776 0.0009086193 -0.0002442391 0.0017836173
MRK   0.0016447926 0.0011854049  0.0009086193 0.0063486415  0.0009187387 0.0009943984
XOM   0.0005925077 0.0010758889 -0.0002442391 0.0009187387  0.0027747986 0.0009486789
^GSPC 0.0024795274 0.0011941026  0.0017836173 0.0009943984  0.0009486789 0.0012362303

Datos:

S_cov_var <- structure(c(0.0119535151035911, 0.000572188710022071, 0.00723524182537011, 
0.00164479256302833, 0.000592507747499871, 0.00247952741729956, 
0.000572188710022071, 0.00422252526205478, 0.000823123610432928, 
0.00118540486616208, 0.00107588894445389, 0.00119410264013768, 
0.00723524182537011, 0.000823123610432928, 0.0191091775682989, 
0.000908619322530227, -0.000244239135715373, 0.00178361731695959, 
0.00164479256302833, 0.00118540486616208, 0.000908619322530227, 
0.00634864154256473, 0.000918738733973792, 0.00099439837734023, 
0.000592507747499871, 0.00107588894445389, -0.000244239135715373, 
0.000918738733973792, 0.00277479857981738, 0.000948678870995285, 
0.00247952741729956, 0.00119410264013768, 0.00178361731695959, 
0.00099439837734023, 0.000948678870995285, 0.00123623026419288
), .Dim = c(6L, 6L), .Dimnames = list(c("INTC", "AEP", "AMZN", 
"MRK", "XOM", "^GSPC"), c("INTC", "AEP", "AMZN", "MRK", "XOM", 
"^GSPC")))

Tengo una matriz de enlaces P:

       INTC AEP AMZN MRK XOM
View 1    0   1    0   0  -1
View 2    1   0   -1   0   0

Datos:

P <- structure(c(0, 1, 1, 0, 0, -1, 0, 0, -1, 0), .Dim = c(2L, 5L), .Dimnames = list(
    c("View 1", "View 2"), c("INTC", "AEP", "AMZN", "MRK", "XOM"
    )))

Calculo Omega como:

tau = 1
Omega = tau * P %*% S_cov_var[1:5 ,1:5] %*% t(P)

Calculo la primera parte de la fórmula como:

$$((\tau S)^{-1} + P^{T}\Omega^{-1}P)^{-1}$$

first <- ((tau * S_cov_var[1:5 ,1:5])^(-1) + (t(P) %*% Omega^(-1) %*% P))^(-1)

Entonces la segunda parte de la fórmula como:

$$(\tau S)^{-1}\pi + P^{T}\Omega^{-1}Q$$

Datos:

Q <- c(0.01, 0.0175)   # uncertainty about my views

implied_equilib_excess_rets <- structure(c(0.00933950373355221, 0.0031834850342374, 0.00648459638783838, 
0.00560398430525973, 0.00578609504932214), .Dim = c(5L, 1L), .Dimnames = list(
    c("INTC", "AEP", "AMZN", "MRK", "XOM"), NULL))

Cálculo:

second <- (tau * S_cov_var[1:5 ,1:5])^(-1) %*% implied_equilib_excess_rets[,1] + (t(P) %*% (Omega^(-1)) %*% Q)

Lo que me da el resultado (para la segunda parte):

          [,1]
INTC 12.274655
AEP  21.034321
AMZN -3.885805
MRK  22.681126
XOM  14.381804

Lo cual es completamente erróneo.

He comprobado todas mis cifras hasta este momento y casi coinciden con las del vídeo que estoy siguiendo (él usa precios ajustados de Yahoo yo uso precios de Closing ya que el vídeo es de hace unos años). Espero que los resultados no coincidan pero no coinciden ni de lejos. Por ejemplo el resultado esperado debería ser (para la segunda parte)

INTC  1.175
AEP   2.304
AMZN -1.074
MRK   0.448
XOM  -0.431

Minuto 11:27 aquí muestra cómo debería ser la segunda parte de la fórmula.

Adicional:

Aquí hay un volcado del código R que tengo del video de excel (obtengo resultados bastante cercanos basados en la salida del video hasta el second parte del código):

library(tsibble)
library(tidyverse)
library(tidyquant)

start_date <- "2002-01-01"
end_date <- "2007-08-01"
symbols <- c("INTC", "AEP", "AMZN", "MRK",  "XOM", "^GSPC")

portfolio_prices <- tq_get(
  symbols,
  from = start_date,
  to = end_date,
) %>% 
  select(symbol, date, close)

portfolio_monthly_prices <- portfolio_prices %>% 
  group_by(symbol) %>% 
  tq_transmute(
    select = close,
    mutate_fun = to.period,
    period = "months"
  ) %>% 
  pivot_wider(names_from = symbol, values_from = close) %>% 
  tk_xts(., date_var = date)

portfolio_monthly_returns <- portfolio_prices %>% 
  group_by(symbol) %>% 
  tq_transmute(
    select = close,
    mutate_fun = periodReturn,
    period = "monthly",                           
    type = "log",
  ) %>% 
  pivot_wider(names_from = symbol, values_from = monthly.returns) %>% 
  tk_xts(., date_var = date)

portfolio_monthly_returns[,1:5]

Asset_Ave_Rets <- colMeans(portfolio_monthly_returns[, 1:5])
Market_Ave_Rets <- colMeans(portfolio_monthly_returns[, 6])
Market_variance <- var(portfolio_monthly_returns[, 6])

obs <- nrow(portfolio_monthly_returns) - 1

S_cov_var <- as.matrix(cov(portfolio_monthly_returns))

Variance <- diag(S_cov_var)
StandardDev <- sqrt(Variance)

lambda = c(1.5, 1.5, 1.5, 1.5, 1.5)

Market_caps <- data.frame(
  stock = c("INTC", "AEP", "AMZN", "MRK", "XOM"),
  mkt_cap = c(153.42, 19.2, 36.62, 125.5, 505.49)
) %>% 
mutate(
    market_weights = mkt_cap / sum(mkt_cap)
  )

weights <- as.vector(Market_caps$market_weights)

implied_equilib_excess_rets <- 2*c(lambda) * (S_cov_var[1:5, 1:5] %*% weights[1:5]) # AKA pi
implied_equilib_excess_rets

#VIEW 1: AP outperforms exxon mobile by 1% per month
#VIEW 2: Intel outperforms Amazon by 1.75 % per month
Q <- c(0.01, 0.0175)

VIEWS = matrix(data = 0, nrow = 2, ncol = ncol(S_cov_var[,1:5]))
rownames(VIEWS) = c(paste("View", seq(1:2)))  
colnames(VIEWS) = colnames(S_cov_var[, 1:5])
# Fill out the link matrix

VIEWS[1, 2] <- 1 
VIEWS[1, 5] <- -1

VIEWS[2, 1] <- 1
VIEWS[2, 3] <- -1

P = as.matrix(VIEWS) # link matrix

tau = 1
Omega = tau * P %*% S_cov_var[1:5 ,1:5] %*% t(P)  # uncertainty associated with our views
Omega

# black litterman formula

# part 1:
# expected returns calculation

first <- ((tau * S_cov_var[1:5 ,1:5])^(-1) + (t(P) %*% Omega^(-1) %*% P))^(-1)
first
# part 2:
second <- (tau * S_cov_var[1:5 ,1:5])^(-1) %*% implied_equilib_excess_rets[,1] + (t(P) %*% (Omega^(-1)) %*% Q)
second

Answer 1

He echado un vistazo a tu código y parece que estás aplicando el inverse en sus cálculos. Por ejemplo, en el second ecuación hice los siguientes cambios:

sub1 = tau * S_cov_var[1:5 ,1:5]
Isub1 = solve(sub1)

IOmega = solve(Omega)

second <- (Isub1 %*% implied_equilib_excess_rets[,1] + (t(P) %*% IOmega %*% Q))

Que devuelve lo siguiente:

> second
           [,1]
INTC  1.7555493
AEP   2.4034124
AMZN -1.0770208
MRK   0.4480916
XOM  -0.5300325

Espero que esto sea lo que está buscando. De lo contrario, hágamelo saber.