1 Il modello di Black-Litterman

Il primo passaggio nell’implementare il modello di BL è quello di estrarre i rendimenti impliciti scontati “a priori” dalle condizioni correnti di mercato. Chiaramente, questo passaggio comporta delle ipotesi, prima tra queste che gli investitori agendo per massimizzare la propria utilità, in base ad una funzione come quella riportata di seguito

\[U = w^\top \mu - \frac{\lambda}{2}w^\top \Sigma w\]

utility <- function(weights, returns, var_covar, lambda){
  w <- weights
  mu <- returns
  S <- var_covar
  l <- lambda
  utility <- t(w) %*% mu - (l/2 * t(w)) %*% S %*% w
  return(utility)
}

dove:

• \(w\) è il vettore (colonna) dei pesi di mercato delle \(n\) attività finanziarie scelte dato dalle capitalizzazioni di mercato note agli investitori;
• \(\mu\) si intende il vettore (colonna) dei rendimenti attesi delle \(n\) attività finanziarie;
• \(\lambda\) è il coefficiente di avversione al rischio e indica il modo in cui gli investitori scambiano rischio con rendimento attesi; un coefficiente \(\lambda\) elevato indica che gli investitori sono molto avversi al rischio —domandano molto rendimento in più per accettare maggior rischio— e viceversa;
• \(\Sigma\) è la matrice varianza-covarianza che rappresenta il modo in cui gli investitori misurano il rischio.

Si possono ricavare i rendimenti impliciti \(\mu\) denominandoli \(\Pi\) e, imponendo che la derivata prima della funzione di utilità rispetto ai pesi \(w\) (le allocazioni di portafoglio date dal mercato) sia pari a zero (nel suo punto di massimo), otteniamo quindi \[\mu = \Pi = \lambda \Sigma w\]

Pi <- function(lambda, Sigma, weight){
  lambda * Sigma %*% w
}

Il modello di BL si compone di fatto di due formule: una per stimare i rendimenti attesi tenendo conto “verosimilmente” delle views di portafoglio dell’investitore e dei rendimenti attesi impliciti “a priori” nelle condizioni correnti di mercato e una per stimare le varianze-covarianze date le views dell’investitore. Il modello di BL esprime quindi la distribuzione dei rendimenti delle attivitià finanziarie in portafoglio in base ad una \(N\left(\bar{\mu}, \bar{\Sigma} \right)\) i cui parametri² media (\(\bar{\mu}\)) e varianza (\(\bar{\Sigma}\)) sono dati da:

\[\bar{\mu} = \bigg[ \big( \tau \Sigma \big)^{-1} + P^\top \Omega^{-1} P\bigg]^{-1} \bigg[ \big( \tau \Sigma \big)^{-1} \Pi + P^\top \Omega^{-1} Q\bigg]\]

mu_bl <- function(tau, S, P, O, Pi, Q) {
  solve(solve(tau * S) + t(P) %*% solve(O) %*% P) %*% 
    (solve(tau * S) %*% Pi + t(P) %*% solve(O) %*% Q)
}

Mentre le varianze-covarianze \(\bar{\Sigma}\), applicando il modello BL sono date da:

\[\bar{\Sigma} = \Sigma + \bigg[ \big( \tau \Sigma \big)^{-1} + P^\top \Omega^{-1} P\bigg]^{-1}\]

s_bl <- function(tau, S, P, O) {
  S + solve(solve(tau * S) + t(P) %*% solve(O) %*% P)
}

dove con:

• \(\bar{\mu}\) si intendono i rendimenti attesi calcolati applicando la formula di BL;
• \(\tau\) è uno scalare, il cui valore è impostato in base alla calibratura del modello (di seguito);
• \(\Sigma\) è la matrice varianza-covarianza \(n \times n\) (\(n\) sono le asset class selezionate) dei rendimenti periodali;
• \(P\) è una matrice di “servizio” o di collegamento \(k \times n\) dove \(k\) rappresenta il numero di views dell’investitore; serve ad inserire le views dell’investitore nel modello per intrecciarle con i rendimenti impliciti attesi a priori scontati nelle condizioni correnti di mercato;
• \(\Omega\) è la matrice \(n \times n\) che contiene le varianze (e anche le covarianze) relative alle views dell’investitore;
• \(\Pi\) è il vettore colonna \(n \times 1\) dei rendimenti attei impliciti già scontati nelle condizioni correnti di mercato;
• \(Q\) è il vettore colonna \(k \times 1\) contenente le views dell’investitore.

Fondamentalmente, la formula per il calcolo dei rendimenti attesi di BL rappresenta una media ponderata complessa \(\hat{w}_{\left(\tau, V, \Omega\right)}\) dei rendimenti attesi impliciti già scontati dalle condizioni correnti di mercato, \(\Pi\), e dei rendimenti attesi contenuti nelle views dell’investitore, \(Q\).

\[\bar{\mu} = \hat{w}_{\left(\tau, \Sigma, \Omega\right)} \Pi + \left(1 - \hat{w}_{\left(\tau, \Sigma, \Omega\right)}\right) Q\]

Si può pensare ai pesi \(\hat{w}_{\left(\tau, V, \Omega\right)}\) come una funzione che dipende dalla verosimiglianza implicita nelle condizioni correnti di mercato data dalle varianze e covarianze contenute nella matrice \(V\) e dalla fiducia a priori, con cui l’investitore esprime le proprie attese nelle views, contenuta nella matrice \(\Omega\).

2 Un’applicazione con alcuni ETF

In questa parte dimostriamo un’applicazione concreta e interamente riproducibile del modello BL realizzata in ambiente R. Di seguito sono riportate, in percentuale, parte delle serie storiche dei rendimenti mensili utilizzate (scaricabili da qui, il pacchetto forecast di Hyndman et al. (2019) è molto utile per trattare le serie storiche in R–si tratta del periodo dal 1/1/2007 al 30/6/2019) per stimare i dati con cui alimentare lo sviluppo del modello e la sua applicazione.

library(tidyverse)
ret <-  read_csv("data/ETFs.csv")
ret <- ret %>% select(1:4)
library(forecast)
subset(ts(round(ret, 4), 
    start = c(2007, 1),
    end = c(2019, 6),
    frequency = 12,
    names = colnames(ret)
    ), start = 138
) * 100

#>            SPY   IEF    XLF    IWM
#> Jun 2018  0.13  0.00  -2.17   0.61
#> Jul 2018  3.70 -0.71   5.11   1.28
#> Aug 2018  3.19  0.82   1.36   4.31
#> Sep 2018  0.14 -1.40  -2.65  -2.58
#> Oct 2018 -6.91 -0.50  -4.71 -10.99
#> Nov 2018  1.85  1.11   2.63   1.73
#> Dec 2018 -9.33  2.38 -11.68 -12.27
#> Jan 2019  8.01  0.65   8.90  11.32
#> Feb 2019  3.24 -0.73   2.24   5.18
#> Mar 2019  1.36  2.46  -3.05  -2.35
#> Apr 2019  4.09 -0.72   8.98   3.40
#> May 2019 -6.38  2.84  -7.17  -7.85
#> Jun 2019  4.98 -0.20   5.15   4.04

Per applicare il modello di BL sono necessari —assieme ai parametri \(\lambda\) e \(\tau\) — come dati di input:

• la matrice-varianza covarianza \(\Sigma\);
• i pesi (dati dalle capitalizzazioni) del portafoglio di mercato;
• le views dell’investitore

(S <- cov(ret))

#>           SPY       IEF       XLF      IWM
#> SPY  0.001817 -0.000236  0.002504  0.00213
#> IEF -0.000236  0.000328 -0.000417 -0.00037
#> XLF  0.002504 -0.000417  0.004565  0.00310
#> IWM  0.002131 -0.000370  0.003098  0.00306

Il coefficiente di avversione al rischio³ \(\lambda\) gioca un ruolo importante: facendo delle ipotesi supportate dalla ricerca applicato ipotizziamo che sia pari a 3.

l <- 3.0

Per poter ottenere i rendimenti “a priori” impliciti attesi e scontati dalle condizioni correnti di mercato è necessario assegnare i pesi \(w\) derivanti dalle capitalizzazioni.

w <- c(0.22, 0.33, 0.25, 0.20)
names(w) <- colnames(ret)
w

#>  SPY  IEF  XLF  IWM 
#> 0.22 0.33 0.25 0.20

Ricaviamo il vettore dei rendimenti impliciti attesi “a priori” nelle condizioni correnti di mercato — dati i parametri impostati e le covarianze storiche — applicando l’equazione che massimizza la funzione di utilità dell’investitore —sono i pesi \(w\) che rendono nulla la derivata prima.

(Pi <- Pi(lambda = l, Sigma = S, weight = w))

#>          [,1]
#> SPY  0.004122
#> IEF -0.000365
#> XLF  0.006523
#> IWM  0.005200

In altre parole si può dire che il mercato sconti un rendimento implicito a priori per il mese successivo —è l’orizzonte temporale sottinteso nell’uso che facciamo delle serie storiche mensili— sull’SPY pari a 0,4122 (%), per esempio, con una volatilità (rappresentata dalla varianza e contenuta nella matrice-varianza covarianza \(\Sigma\)) pari a 0,001817.

A questo punto è necessario introdurre le views e ipotizziamo che l’investitore abbia delle prospettive circa le asset class (gli ETF in questo caso) sull’orizzonte temporale ipotizzato (in questo caso mensile), espresse come rendimenti attesi che, a loro volta, possono essere trattati come views assolute o views relative; vediamo meglio ipotizzando per i 4 ETF che l’investitore si attenda che:

(Q <- c(0.01, 0.005, 0.005))

#> [1] 0.010 0.005 0.005

Ossia:

1. l’SPY sovraperformi IEF di 1,00% (view relativa);
2. IEF abbia un rendimento di 0,50% (view assoluta);
3. XLF sovraperformi IWM di 0,50% (view relativa).

Costruiamo ora la matrice di collegamento (o di servizio) \(P\) necessaria per inserire le views nel modello.

P <- matrix(c(1, -1,  0,  0,   
              0,  1,  0,  0,   
              0,  0,  1, -1),
            nrow = 3, 
            ncol = 4, 
            byrow = TRUE)
rownames(P) <- seq_along(Q)
colnames(P) <- colnames(ret)
# lo schema di pesatura ignora le capitalizzazioni relative, ahimé!!!
P

#>   SPY IEF XLF IWM
#> 1   1  -1   0   0
#> 2   0   1   0   0
#> 3   0   0   1  -1

Impostiamo lo scalare \(\tau = 0{,}025\) —di fatto, il modo in cui viene implementato il modello rende i rendimenti attesi in base alle views indipendenti dal valore dello scalare \(\tau\) (Idzorek 2007).

tau <- 0.025

Calcoliamo la matrice \(\Omega\) contenente l’incertezza riguardante le views — dipende dalle varianze-covarianze storiche. Normalmente, la matrice \(\Omega\) viene calcolata come matrice diagonale — anche se non è strettmanente necessario. Il calcolo di \(\Omega\) è probabilmente il passaggio più dibattutto nell’implementazione del modello di BL —assieme alla definizione dello scalare \(\tau\) nel caso in cui si voglia tenere conto delle covarianze tra le views dell’investitore.

La matrice \(\Omega\) può essere costruita tenendo conto delle covarianze nel modo seguente.

O <- P %*% S %*% t(P) * tau

Optiamo per un’implementazione del modello BL in cui la matrice \(\Omega\) è una matrice diagonale.

(O <- diag(diag(P %*% S %*% t(P) * tau)))

#>          [,1]     [,2]     [,3]
#> [1,] 6.54e-05 0.00e+00 0.00e+00
#> [2,] 0.00e+00 8.21e-06 0.00e+00
#> [3,] 0.00e+00 0.00e+00 3.57e-05

Prima di applicare le equazioni per calcolare i rendimenti attesi e le varianze può essere utile spezzare in alcuni passaggi il calcolo dei primi. Come si è visto i rendimenti attesi nel modello di BL sono il risultato del prodotto tra due fattori. Il primo fattore è una matrice \(n \times n\) data dall’inversa della somma di due addendi:

\[\bigg[ \big( \tau \Sigma \big)^{-1} + P^\top \Omega^{-1} P\bigg]^{-1} \] Il primo addendo è dato da

(i <- solve(tau * S))

#>        SPY    IEF    XLF    IWM
#> SPY 163382 -18251 -40075 -75454
#> IEF -18251 143720   8655  21325
#> XLF -40075   8655  37982  -9500
#> IWM -75454  21325  -9500  77849

Mentre il secondo è dato da:

(ii <- t(P) %*% solve(O) %*% P)

#>        SPY    IEF    XLF    IWM
#> SPY  15280 -15280      0      0
#> IEF -15280 137130      0      0
#> XLF      0      0  28019 -28019
#> IWM      0      0 -28019  28019

La somma dei due addendi è data da:

i + ii

#>        SPY    IEF    XLF    IWM
#> SPY 178662 -33530 -40075 -75454
#> IEF -33530 280850   8655  21325
#> XLF -40075   8655  66001 -37519
#> IWM -75454  21325 -37519 105868

Il primo fattore è la matrice inversa della somma dei due addendi, ossia:

(I <- solve(i + ii))

#>           SPY       IEF       XLF       IWM
#> SPY  2.49e-05 -2.04e-07  3.16e-05  2.90e-05
#> IEF -2.04e-07  3.68e-06 -1.39e-06 -1.38e-06
#> XLF  3.16e-05 -1.39e-06  5.95e-05  4.39e-05
#> IWM  2.90e-05 -1.38e-06  4.39e-05  4.60e-05

Il secondo fattore nella formula per il calcolo delle medie nei rendimenti attesi è un vettore colonna \(n \times 1\), anch’esso dato dalla somma di due addendi:

\[\bigg[ \big( \tau \Sigma \big)^{-1} \Pi + P^\top \Omega^{-1} Q\bigg]\]

Il primo addendo è dato da:

(j <- solve(tau * S) %*% Pi)

#>     [,1]
#> SPY 26.4
#> IEF 39.6
#> XLF 30.0
#> IWM 24.0

Il secondo addendo è dato da:

(jj <- t(P) %*% solve(O) %*% Q)

#>     [,1]
#> SPY  153
#> IEF  456
#> XLF  140
#> IWM -140

Il secondo fattore è dunque dato dalla somma di:

(J <- j + jj)

#>     [,1]
#> SPY  179
#> IEF  496
#> XLF  170
#> IWM -116

Con il prodotto dei due fattori otteniamo:

I %*% J

#>        [,1]
#> SPY 0.00637
#> IEF 0.00171
#> XLF 0.01000
#> IWM 0.00665

Che è, ovviamente, lo stesso risultato di applicare la formula BL per il calcolo dei rendimenti attesi

(r_bl <- mu_bl(tau, S, P, O, Pi, Q))

#>        [,1]
#> SPY 0.00637
#> IEF 0.00171
#> XLF 0.01000
#> IWM 0.00665

Facciamo la stessa cosa per le varianze BL

(V_bl <- s_bl(tau, S, P, O))

#>           SPY       IEF       XLF       IWM
#> SPY  0.001842 -0.000237  0.002536  0.002160
#> IEF -0.000237  0.000332 -0.000418 -0.000371
#> XLF  0.002536 -0.000418  0.004625  0.003142
#> IWM  0.002160 -0.000371  0.003142  0.003105

(w_opt <- solve(l * V_bl) %*% r_bl)

#>       [,1]
#> SPY  0.888
#> IEF  2.933
#> XLF  0.627
#> IWM -0.188

arggs <- list(decimal.mark = ',',
              big.mark = ".",
              digits = 2,
              nsmall = 2)
knitr::kable(tibble(w_bl = w_opt, 
       w_pi = w , 
       diff = w_bl - w_pi), 
       format.args = arggs
)

f_utility <- function(x){
  1/utility(weights = x, 
          returns = r_bl, 
          var_covar = V_bl, 
          lambda = l)
}

max_u_etf <- optim(
  par = w,
  fn = f_utility,
  gr = NULL, 
  method = "L-BFGS-B",
  lower = rep(0.0001, length(w)),
  upper = rep(0.9999, length(w))
)

(w_bl <- max_u_etf$par/sum(max_u_etf$par)) * 100

#>  SPY  IEF  XLF  IWM 
#> 23.9 31.2 23.8 21.1

bl_1_out <- tibble(
  asset = names(S[1,]), 
  r_bl = r_bl[ , 1] * 100, 
  r_pi = Pi[ , 1] * 100, 
  s_bl = diag(V_bl)^0.5 * 100,
  s_pi = diag(S)^0.5 * 100, 
  w_bl = w_bl * 100, 
  w_pi = w * 100
)
knitr::kable(bl_1_out,
             format.args = arggs)

ggplot(data = data.frame(x = c(-10, 10)), aes(x = x)) +
  stat_function(fun = dnorm, 
                n = 101, 
                args = list(mean = bl_1_out$r_pi[2], sd = bl_1_out$s_pi[2]), 
                aes(linetype = "a priori")) +
  stat_function(fun = dnorm, 
                n = 101, 
                args = list(mean = Q[2]*100, sd = bl_1_out$s_pi[2]), 
                aes(linetype = "verosimile")) + 
  stat_function(fun = dnorm, 
                n = 101, 
                args = list(mean = bl_1_out$r_bl[2], sd = bl_1_out$s_bl[2]), 
                aes(linetype = "a posteriori")) + 
  labs(x = "IEF") + ylab("") +
  scale_y_continuous(breaks = NULL) +
  scale_linetype_manual("Distribuzione", values = c(1, 2, 3)) + 
  theme(legend.position = "bottom", 
        legend.box = "vertical")

bl <- bl_1_out %>% select(asset, r_bl, s_bl, w_bl) %>% 
  mutate(model = "modello Black-Litterman") %>% 
  rename(s = s_bl) %>% 
  rename(r = r_bl) %>% 
  rename(w = w_bl)

mkt_impl <- bl_1_out %>% select(asset, r_pi, s_pi, w_pi) %>% 
  mutate(model = "modello implicito nel mercato") %>% 
  rename(s = s_pi) %>% 
  rename(r = r_pi) %>%
  rename(w = w_pi)

3 Un’applicazione con i titoli dell’S&P 500

Un contesto in cui l’applicazione del modello di BL è sicuramente più complesso è quello in cui le asset class sono rappresentate da singoli titoli le cui serie storiche di rendimento sono rilevate con frequenza giornaliera. La complessità deriva dal fatto che si tende a lavorare con portafogli composti da alcune decine di titoli e questo può comportare da un punto di vista computazionale grossi problemi.

Accenniamo —vale la pena di farlo— al caso dell’applicazione del modello BL ad un indice come ad esempio l’S&P 500 —i tickers sono scaricabili da qui mentre le serie storiche dei rendimenti degli stessi titoli ret_sp500 sono scaricabili da qui.

Nella tabella che segue sono rappresentati i primi 10 titoli per capitalizzazione con dati aggiornati al 15 luglio 2019.

# library(tidyquant)
# tickers <- tq_index("SP500")
# write_csv(tickers, "data/tickers.csv")
tickers <- read_csv("data/tickers.csv")
knitr::kable(tickers[1:10, -5] %>% 
  mutate(weight = weight * 100), 
  format.args = arggs
)

Un ulteriore esempio di difficoltà che si incontrano nell’applicazione del modello di BL riguarda la completezza dei dati. Nel caso dell’indice S&P 500 è facile che una delle serie storiche dei rendimenti di un titolo azionario sia incompleta e, naturalmente, per questa motivo è necessario tenerne conto.

ret_sp500 <- read_csv("data/sp500.csv")
ret_sp500 <- ret_sp500  %>% 
  gather(security, return, 2:(ncol(ret_sp500) - 1), na.rm = TRUE) %>%
  spread(key = security, value = return)

L’output della matrice var-covar dei primi 4 titoli è riportato di seguito. Uno dei problemi che si può manifestare nell’implementazione del modello di BL è quello della non invertibilità della matrice varianza-covarianza (\(\Sigma\)).

V_sp500 <- ret_sp500 %>% 
  select(-date) %>% 
  cov(use ="complete.obs")
V_sp500[1:4, 1:4]

#>          ZTS        A      AAL      AAP
#> ZTS 1.48e-04 1.73e-05 1.36e-04 4.88e-05
#> A   1.73e-05 7.68e-05 7.86e-05 3.28e-05
#> AAL 1.36e-04 7.86e-05 6.86e-04 8.18e-05
#> AAP 4.88e-05 3.28e-05 8.18e-05 1.76e-04

Con un numero elevato di attività finanziarie questa eventualità diviene molto probabile che si manifesti: dipende dalle correlazioni esistenti tra due o più attività finanziarie. Nel caso di poche attività finanziarie è poco probabile che questo si accada, ma la probabilità aumenta all’aumentare della loro numerosità.

Nel caso specifico la matrice varianza-covarianza dei rendimenti giornalieri dell’S&P500 non è invertibile —il determinante è nullo perché la matrice var-covar è singolare— e per questo motivo non è possibile ricavare il vettore dei rendimenti attesi BL \(\mu\). La ragione è legata alla numerosità di titoli: si tratta di 500 titoli (503 per l’esattezza), e di conseguenza la possibilità che due o più titoli siano significativamente correlati è piuttosto elevata. Questa situazione causa la non invertibilità delle matrice varianza-covarianza.

det(V_sp500)

#> [1] 0

Può essere utile sviluppare una funzione che calcoli il determinante della matrice \(\Sigma\) per poi utilizzarla in una procedura che reiterando il calcolo del determinante trovi il numero massimo di attività finanziarie in cui il determinante è diverso da zero.

n_smpl <- function(ret_sp500, tickers, n) {
  n <- n
  ccc <- colnames(ret_sp500)
  sm <- intersect(ccc, tickers$symbol[1:n])
  ret_sp_sm <- ret_sp500 %>% select(sm)
  V_sp_sm <- ret_sp_sm %>% cov(use ="complete.obs")
  det(V_sp_sm)
}

n_smpl(ret_sp500 = ret_sp500, tickers = tickers, n = 400)

#> [1] 0

Con 400 attività finanziarie il determinante è ancora pari a zero. Proviamo quindi a reiterare il calcolo del determinante per definire quel numero massimo di attività finanziarie n per cui il determinante non è zero.

det <- 0
n <- 500
while(det == 0){
  n <- n - 1
  det <- n_smpl(ret_sp500 = ret_sp500, tickers = tickers, n = n)
}
n

#> [1] 78

Nel caso specifico, il numero massimo di titoli con cui è implementabile il modello BL è 78, superando questo limite la matrice \(\Sigma\) diventa non invertibile.

Calcoliamo quindi la matrice varianza-covarianza \(\Sigma\) con i primi 78 titoli.

ccc <- colnames(ret_sp500)
sm <- intersect(ccc, tickers$symbol[1:n])
ret_sp_sm <- ret_sp500 %>% select(sm)
V_sp_sm <- ret_sp_sm %>% cov(use ="complete.obs")
w_sp_sm <- pull(
  tickers %>% filter(symbol %in% sm) %>% select(weight)
)

Si procede estraendo i rendimenti attesi “a priori” e impliciti nelle condizioni di mercato.

Pi_sp_sm <- l * V_sp_sm %*% (w_sp_sm/sum(w_sp_sm)) * 252 * 100
(Pi_sm <- tibble(symbol = row.names(Pi_sp_sm), 
                 Pi = as.vector(Pi_sp_sm)
                 )
  )

#> # A tibble: 77 x 2
#>    symbol    Pi
#>    <chr>  <dbl>
#>  1 AAPL    9.25
#>  2 ABBV    7.12
#>  3 ABT     5.92
#>  4 ACN     4.82
#>  5 ADBE    7.96
#>  6 ADP     5.65
#>  7 AMGN    4.83
#>  8 AMT     1.03
#>  9 AMZN    8.00
#> 10 ANTM    5.41
#> # … with 67 more rows

I risultati —relativi ai primi 10 titoli per capitalizzazione di mercato— sono riportati di seguito nella tabella che segue, indicando con s_pi, r_pi e con w_pi rispettivamente la deviazione standard dei rendimenti storici, i rendimenti impliciti attesi, entrambi espressi su base annua⁴, e i pesi di mercato.

bl_sm <-  arrange(left_join(Pi_sm, 
                  tickers %>% select(-shares_held)
                  ), 
        desc(weight*100))

bl_sm_dt <- left_join(bl_sm, 
          tibble(symbol = row.names(V_sp_sm), 
                 sd = sqrt(diag(V_sp_sm)) * sqrt(252) * 100
                 )
          ) %>% 
    mutate(weight = weight * 100) %>% 
    select(symbol, company, Pi, sd, weight, sector) 

knitr::kable(bl_sm_dt %>% 
               select(-sector, -company) %>% 
               filter(weight >= 1.25) %>% 
               rename( 
                      r_pi = Pi, 
                      s_pi = sd, 
                      w_pi  = weight), 
             format.args = arggs)

bl_sm_dt %>% ggplot(aes(sd, Pi, size = weight, color = sector)) + 
  geom_point() + 
  theme(legend.position = "bottom", 
        legend.box = "vertical", 
        legend.title=element_text(size=9), 
        legend.text=element_text(size=8)) + 
  scale_colour_brewer(palette = "Paired") +
  ggrepel::geom_label_repel(data = filter(bl_sm_dt, weight >= 1.25),
                            aes(label = symbol), 
                            show.legend = FALSE, 
                            size = 4.0) + 
  labs(x = "deviazione standard (%, giorn. su base annua)", 
       y = "rendimento atteso (%, giorn. su base annua)"
       ) + 
  guides(size = guide_legend(title = "pesi (%)"), 
         color = guide_legend(title = "settore"))

V_sp_sm <- V_sp_sm[bl_sm_dt$symbol, bl_sm_dt$symbol]

Q_sp_500 <- c(25.0, 15.0, 10.0)

P_sp_500 <-  matrix(c(1,  0, rep(0, nrow(V_sp_sm) - 2), 
               0,  1, -1, 0, rep(0, nrow(V_sp_sm) - 4), 
               0,  0, 0,  -1, 1, 0, rep(0, nrow(V_sp_sm) - 6)), 
              nrow = 3, 
             ncol = nrow(V_sp_sm), byrow = TRUE)
colnames(P_sp_500) <- bl_sm_dt$symbol
P_sp_500[ , 1:5]

#>      MSFT AAPL AMZN FB JNJ
#> [1,]    1    0    0  0   0
#> [2,]    0    1   -1  0   0
#> [3,]    0    0    0 -1   1

(O_sp_500 <- diag(diag(P_sp_500 %*% V_sp_sm %*% t(P_sp_500) * tau)))

#>          [,1]    [,2]     [,3]
#> [1,] 4.78e-06 0.0e+00 0.00e+00
#> [2,] 0.00e+00 6.2e-06 0.00e+00
#> [3,] 0.00e+00 0.0e+00 1.01e-05

r_bl_sp_sm <- mu_bl(tau = tau, 
                S = V_sp_sm, 
                P = P_sp_500, 
                O = O_sp_500, 
                Pi = bl_sm$Pi, 
                Q = Q_sp_500)
bl_sm_dt$r_bl <- r_bl_sp_sm[ , 1]

V_bl_sp_sm <- s_bl(tau = tau, 
                                S = V_sp_sm, 
                                P = P_sp_500, 
                                O = O_sp_500)

bl_sm_dt %>% select(symbol, Pi, r_bl)

#> # A tibble: 77 x 3
#>    symbol    Pi  r_bl
#>    <chr>  <dbl> <dbl>
#>  1 MSFT    7.27 15.2 
#>  2 AAPL    9.25 19.7 
#>  3 AMZN    8.00 11.1 
#>  4 FB      6.85  4.91
#>  5 JNJ     3.04  6.42
#>  6 JPM     4.31  7.74
#>  7 GOOG    6.78  9.85
#>  8 GOOGL   6.74  9.66
#>  9 XOM     4.10  6.67
#> 10 V       5.37 10.7 
#> # … with 67 more rows

f_utility_sp <- function(x){
  1/utility(weights = x, 
          returns = bl_sm_dt$r_bl, 
          var_covar = V_bl_sp_sm, 
          lambda = l)
}

f_utility_sp(bl_sm_dt$weight)

#>         [,1]
#> [1,] 0.00176

max_u_sp <- optim(par = bl_sm_dt$weight, 
               fn = f_utility_sp)

(bl_temp <- bl_sm_dt %>% mutate(weight_bl = max_u_sp$par)  %>% 
  select(symbol, sd, Pi, weight, r_bl, weight_bl, everything()) %>% 
  rename(implied_ret = Pi, bl_ret = r_bl, mkt_weight = weight, bl_weight = weight_bl) %>%
  mutate(diff_ret = bl_ret - implied_ret) %>% 
  mutate(diff_weight = bl_weight - mkt_weight) %>% 
  select(company, everything())
)

#> # A tibble: 77 x 10
#>    company symbol    sd implied_ret mkt_weight bl_ret bl_weight sector
#>    <chr>   <chr>  <dbl>       <dbl>      <dbl>  <dbl>     <dbl> <chr> 
#>  1 Micros… MSFT    22.0        7.27       4.26  15.2      4.38  Infor…
#>  2 Apple … AAPL    29.5        9.25       3.54  19.7      3.67  Infor…
#>  3 Amazon… AMZN    26.6        8.00       3.33  11.1      3.39  Consu…
#>  4 Facebo… FB      32.2        6.85       1.94   4.91     1.87  Commu…
#>  5 Johnso… JNJ     14.3        3.04       1.50   6.42     1.44  Healt…
#>  6 JPMorg… JPM     18.6        4.31       1.49   7.74     1.52  Finan…
#>  7 Alphab… GOOG    25.5        6.78       1.41   9.85     1.46  Commu…
#>  8 Alphab… GOOGL   25.4        6.74       1.38   9.66     1.42  Commu…
#>  9 Exxon … XOM     17.2        4.10       1.32   6.67     0.485 Energy
#> 10 Visa I… V       17.2        5.37       1.26  10.7      1.32  Infor…
#> # … with 67 more rows, and 2 more variables: diff_ret <dbl>,
#> #   diff_weight <dbl>

bl_2_out <- tibble(
  stock = bl_temp$company, 
  symbol = bl_temp$symbol, 
  r_bl = bl_temp$bl_ret, 
  r_pi = bl_temp$implied_ret, 
  s_bl = diag(V_bl_sp_sm)^0.5 * sqrt(252) * 100,
  s_pi = bl_temp$sd,
  w_bl = bl_temp$bl_weight, 
  w_pi = bl_temp$mkt_weight, 
  sector = bl_temp$sector
)

knitr::kable(bl_2_out[1:10, ] %>% 
               select(symbol, r_bl, r_pi, s_bl, s_pi, w_bl, w_pi), 
             format.args = arggs)

ggplot(data = data.frame(x = c(-75, 100)), aes(x = x)) +
  stat_function(fun = dnorm, 
                n = 101, 
                args = list(mean = bl_2_out$r_pi[1], 
                            sd = bl_2_out$s_pi[1]), 
                aes(linetype = "a priori")) +
  stat_function(fun = dnorm, 
                n = 101, 
                args = list(mean = Q_sp_500[1], 
                            sd = bl_2_out$s_pi[1]), 
                aes(linetype = "verosimile")) + 
  stat_function(fun = dnorm, 
                n = 101, 
                args = list(mean = bl_2_out$r_bl[1], 
                            sd = bl_2_out$s_bl[1]), 
                aes(linetype = "a posteriori")) + 
  labs(x = "MSFT") + ylab("") +
  scale_y_continuous(breaks = NULL) +
  scale_linetype_manual("Distribuzione", values = c(1, 2, 3)) + 
  theme(legend.position = "bottom", 
        legend.box = "vertical")

implied <- bl_2_out %>% 
  select(stock, symbol, s_pi, r_pi, w_pi, sector) %>% 
  mutate(model = "modello implicito nel mercato") %>% 
  rename(return = r_pi, devst = s_pi, weight = w_pi)

bl <- bl_2_out %>% 
  select(stock, symbol, s_bl, r_bl, w_bl, sector) %>% 
  mutate(model = "Black-Litterman") %>% 
  rename(return = r_bl, devst = s_bl, weight = w_bl)

sp500 <- bind_rows(implied, bl)

sp500 %>% ggplot(aes(devst, return, size = weight, color = sector)) + 
  geom_point(alpha = 1/3) +  theme(legend.position = "bottom", 
        legend.box = "vertical", 
        legend.title=element_text(size=9), 
        legend.text=element_text(size=8)) + 
  scale_colour_brewer(palette = "Paired") +
  facet_wrap( ~ model) + 
  labs(x = "deviazione standard (%, giorn. su base annua)", 
       y = "rendimento atteso\n(%, giorn. su base annua)") + 
  geom_point(data = sp500 %>% filter(weight >=1.25)) +
  ggrepel::geom_label_repel(data = filter(sp500, weight >= 1.25), 
                            aes(label = symbol), 
                            show.legend = FALSE, 
                            size = 4.0) + 
  guides(color = guide_legend(title = "settore"), 
         size = guide_legend(title = "pesi (%)"))