This R session will illustrate the incorporation of prior information in the estimation of \(\boldsymbol{\mu}\) and \(\boldsymbol{\Sigma}\).

(Useful R links: Cookbook R, Quick-R, R documentation, CRAN, METACRAN.)

Shrinkage

Suppose we have an \(N\)-dimensional i.i.d. Gaussian time series \[ \mathbf{x}_t \sim \mathcal{N}(\boldsymbol{\mu}, \boldsymbol{\Sigma}). \] The sample estimators for the mean and covariance matrix are, respectively, the sample mean \[ \hat{\boldsymbol{\mu}} = \frac{1}{T}\sum_{t=1}^T \mathbf{x}_t \] and the sample covariance matrix \[ \hat{\boldsymbol{\Sigma}} = \frac{1}{T-1}\sum_{t=1}^T (\mathbf{x}_t - \hat{\boldsymbol{\mu}})^T(\mathbf{x}_t - \hat{\boldsymbol{\mu}}). \] We will next explore the idea of shrinkage in order to improve those sample estimators when \(T\) is not large enough.

Shrinkage in \(\boldsymbol{\mu}\)

The James-Stein estimator is \[ \hat{\boldsymbol{\mu}}^\textsf{JS} = (1-\rho) \hat{\boldsymbol{\mu}} + \rho \mathbf{t} \] where \(\mathbf{t}\) is the shrinkage target and \(0\le\rho\le1\) is the amount of shrinkage. Common choices for the target \(\mathbf{t}\) are:

  • \(\mathbf{t} = \mathbf{0}\)
  • \(\mathbf{t} = 0.1 \times \mathbf{1}\)
  • \(\mathbf{t} = \frac{\mathbf{1}^T \hat{\boldsymbol{\mu}}}{N} \times \mathbf{1}\) (grand mean target)
  • \(\mathbf{t} = \frac{\mathbf{1}^T\hat{\boldsymbol{\Sigma}}^{-1}\hat{\boldsymbol{\mu}}}{\mathbf{1}^T\hat{\boldsymbol{\Sigma}}^{-1}\mathbf{1}} \times \mathbf{1}\) (volatility-weighted mean)

A choice of \(\rho\) that improves the estimatior error is \[ \rho = \frac{1}{T}\frac{N\bar{\boldsymbol{\lambda}}-2\boldsymbol{\lambda}_\textsf{max}}{\|\hat{\boldsymbol{\mu}} - \mathbf{t}\|^2} \] where \(\bar{\boldsymbol{\lambda}}\) and \(\boldsymbol{\lambda}_\textsf{max}\) are the average and maximum values, respectively, of the eigenvalues of \(\boldsymbol{\Sigma}\) (or \(\hat{\boldsymbol{\Sigma}}\)).

We will first explore this estimator using synthetic data. Let’s start generating the data:

Now let’s try the James-Stein estimator with the four targets:

Let’s now compare the estimation error measured in terms of squared error and PRIAL (percentage relative improvement in average loss): \[ \begin{align} \textsf{SE} &= \|\hat{\boldsymbol{\mu}}^\textsf{JS} - \boldsymbol{\mu}\|^2\\ \textsf{PRIAL} &= \frac{\|\hat{\boldsymbol{\mu}} - \boldsymbol{\mu}\|^2 - \|\hat{\boldsymbol{\mu}}^\textsf{JS} - \boldsymbol{\mu}\|^2}{\|\hat{\boldsymbol{\mu}} - \boldsymbol{\mu}\|^2} \end{align} \]

The improvement in performance is quite clear.

We can now plot the error as a function of the number of samples to get a more complete picture:

# first generate all the data
set.seed(357)
N <- 100
T_max <- 300
mu <- runif(N)
U <- t(rmvnorm(n = round(0.7*N), sigma = 0.1*diag(N)))
Sigma <- U %*% t(U) + diag(N)
X <- rmvnorm(n = T_max, mean = mu, sigma = Sigma)

# now loop over subsets of the samples
SE_mu_vs_T <- NULL
T_sweep <- ceiling(seq(1.01*N, T_max, length.out = 20))
for (T in T_sweep) {
  X_ <- X[1:T, ]
  # sample estimates
  mu_sm <- colMeans(X_)
  Sigma_scm <- cov(X_)
  # James-Stein estimators
  t_3 <- rep(mean(mu_sm), N)
  t_4 <- rep(sum(solve(Sigma_scm, mu_sm))/sum(solve(Sigma_scm, rep(1, N))), N)
  lambdas <- eigen(Sigma_scm)$values
  lmd_mean <- mean(lambdas)
  lmd_max <- max(lambdas)
  rho_1 <- (1/T)*(N*lmd_mean - 2*lmd_max)/norm(mu_sm - t_1, "2")^2
  rho_2 <- (1/T)*(N*lmd_mean - 2*lmd_max)/norm(mu_sm - t_2, "2")^2
  rho_3 <- (1/T)*(N*lmd_mean - 2*lmd_max)/norm(mu_sm - t_3, "2")^2
  rho_4 <- (1/T)*(N*lmd_mean - 2*lmd_max)/norm(mu_sm - t_4, "2")^2
  mu_JS_1 <- (1-rho_1)*mu_sm + rho_1*t_1
  mu_JS_2 <- (1-rho_2)*mu_sm + rho_2*t_2
  mu_JS_3 <- (1-rho_3)*mu_sm + rho_3*t_3
  mu_JS_4 <- (1-rho_4)*mu_sm + rho_4*t_4
  # compute errors
  SE_mu_vs_T <- rbind(SE_mu_vs_T, c(norm(mu_sm - mu, "2")^2,
                                    norm(mu_JS_1 - mu, "2")^2,
                                    norm(mu_JS_2 - mu, "2")^2,
                                    norm(mu_JS_3 - mu, "2")^2,
                                    norm(mu_JS_4 - mu, "2")^2))
}
colnames(SE_mu_vs_T) <- names(SE_mu)
rownames(SE_mu_vs_T) <- paste("T =", T_sweep)
# compute PRIAL
PRIAL_mu_vs_T <- (SE_mu_vs_T[, 1] - SE_mu_vs_T[, -1])/SE_mu_vs_T[, 1]

# plots
matplot(T_sweep, SE_mu_vs_T, 
        main = "Squared error in estimation of mu", xlab = "T", ylab = "squared error",
        type = "b", pch = 20, col = rainbow(5))
legend("topright", inset=0.01, legend = colnames(SE_mu_vs_T), pch = 20, col = rainbow(5))

Indeed, we can clearly conclude the huge improvement from using shrinkage in the mean estimation. In particular, the James-Stein estimator using as target the grand mean and the volatility weighted mean are the best by a big difference.

Shrinkage in \(\boldsymbol{\Sigma}\)

The shrinkage estimator for the covariance matrix is \[ \hat{\boldsymbol{\Sigma}}^\textsf{sh} = (1-\rho) \hat{\boldsymbol{\Sigma}} + \rho \mathbf{T} \] where \(\mathbf{T}\) is the shrinkage target and \(0\le\rho\le1\) is the amount of shrinkage.

Common choices for the target \(\mathbf{T}\) are:

  • \(\mathbf{T} = \frac{1}{N}\textsf{Tr}\left(\hat{\boldsymbol{\Sigma}}\right) \times \mathbf{I}\) (scaled identity)
  • \(\mathbf{T} = \textsf{Diag}(\hat{\boldsymbol{\Sigma}})\) (diagonal)

To determine \(\rho\) we will use three different criteria via random matrix theory (RMT):

  • MSE of covariance matrix
  • quadratic loss of precision matrix
  • Sharpe ratio.

Let’s start by generating some synthetic data:

The two targets are easily computed:

We next compute the Ledoit-Wolf estimator: \[ \rho =\min\left(1,\frac{\frac{1}{T^{2}}\sum_{t=1}^{T}||\hat{\boldsymbol{\Sigma}}-\mathbf{r}_{t}\mathbf{r}_{t}^{T}||_{F}^{2}}{||\hat{\boldsymbol{\Sigma}}-\mathbf{T}||_{F}^{2}}\right) \]

Let’s now compare the estimation error measured in terms of squared error and PRIAL: \[ \begin{align} \textsf{SE} &= \|\hat{\boldsymbol{\Sigma}}^\textsf{sh} - \boldsymbol{\Sigma}\|_F^2\\ \textsf{PRIAL} &= \frac{\|\hat{\boldsymbol{\Sigma}} - \boldsymbol{\Sigma}\|_F^2 - \|\hat{\boldsymbol{\Sigma}}^\textsf{sh} - \boldsymbol{\Sigma}\|_F^2}{\|\hat{\boldsymbol{\Sigma}} - \boldsymbol{\Sigma}\|_F^2} \end{align} \]

Let’s try now the shrinkage estimator based on maximizing the Sharpe ratio: \[ \begin{array}{ll} \underset{\rho_{1}\geq0}{\textsf{maximize}} & \begin{array}{c} \frac{\hat{\boldsymbol{\mu}}^{T}(\hat{\boldsymbol{\Sigma}}^{\textsf{sh}})^{-1}\hat{\boldsymbol{\mu}}-\delta}{\sqrt{b\hat{\boldsymbol{\mu}}^{T}(\hat{\boldsymbol{\Sigma}}^{\textsf{sh}})^{-1}\hat{\boldsymbol{\Sigma}}(\hat{\boldsymbol{\Sigma}}^{\textsf{sh}})^{-1}\hat{\boldsymbol{\mu}}}}\end{array}\\ \textsf{subject to} & \begin{array}[t]{l} \hat{\boldsymbol{\Sigma}}^{\textsf{sh}}=\rho_{1}\mathbf{T}+\hat{\boldsymbol{\Sigma}}\\ \delta=D/(1-D)\\ D=\frac{1}{T}\textsf{Tr}(\hat{\boldsymbol{\Sigma}}(\hat{\boldsymbol{\Sigma}}^{\textsf{sh}})^{-1})\\ b=\frac{T}{\textsf{Tr}(\mathbf{W}(\mathbf{I}+\delta\mathbf{W})^{-2})} \end{array} \end{array} \] where \(\mathbf{W}=\mathbf{I}-\frac{1}{T}\mathbf{1}\mathbf{1}^{T}\).

We can now evaluate the performance:

We can now plot the performance measures as a function of the number of samples to get a more complete picture (we will use only the identity target):

# first generate all the data
set.seed(357)
N <- 100
T_max <- 300
mu <- runif(N)
U <- t(rmvnorm(n = round(0.7*N), sigma = 0.1*diag(N)))
Sigma <- U %*% t(U) + diag(N)
X <- rmvnorm(n = T_max, mean = mu, sigma = Sigma)

# now loop over subsets of the samples
SR_Sigma_vs_T <- SE_Sigma_vs_T <- NULL
T_sweep <- ceiling(seq(1.01*N, T_max, length.out = 20))
for (T in T_sweep) {
  X_ <- X[1:T, ]
  # sample estimates
  mu_sm <- colMeans(X_)
  Sigma_scm <- cov(X_)
  # methods: Ledoit-Wolf, max-SR
  T_I <- sum(diag(Sigma_scm))/N * diag(N)
  rho <- rho_LW(X_, Sigma_scm, T_I)
  Sigma_LW <- (1-rho)*Sigma_scm + rho*T_I
  rho <- rho_maxSR(mu_sm, Sigma_scm, T_I, T)
  Sigma_maxSR <- (1-rho)*Sigma_scm + rho*T_I
  # compute performance
  SE_Sigma_vs_T <- rbind(SE_Sigma_vs_T, c(norm(Sigma_scm - Sigma, "F")^2,
                                          norm(Sigma_LW - Sigma, "F")^2,
                                          norm(Sigma_maxSR - Sigma, "F")^2))
  SR_Sigma_vs_T <- rbind(SR_Sigma_vs_T, c(compute_SR(Sigma_scm, mu_sm, mu, Sigma),
                                          compute_SR(Sigma_LW, mu_sm, mu, Sigma),
                                          compute_SR(Sigma_maxSR, mu_sm, mu, Sigma)))
}
colnames(SR_Sigma_vs_T) <- colnames(SE_Sigma_vs_T) <- c("SCM", "LW", "max-SR")
rownames(SR_Sigma_vs_T) <- rownames(SE_Sigma_vs_T) <- paste("T =", T_sweep)
# compute PRIAL
PRIAL_Sigma_vs_T <- (SE_Sigma_vs_T[, 1] - SE_Sigma_vs_T[, -1])/SE_Sigma_vs_T[, 1]

# plots
matplot(T_sweep, SE_Sigma_vs_T, 
        main = "Squared error in estimation of Sigma", xlab = "T", ylab = "squared error",
        type = "b", pch = 20, col = rainbow(3))
legend("topright", inset=0.01, legend = colnames(SE_Sigma_vs_T), pch = 20, col = rainbow(3))

Shrinkage in \(\boldsymbol{\mu}\) and \(\boldsymbol{\Sigma}\) with market data

We start by loading some stock market data and dividing it into a training set and test set:

Let’s consider three estimations for \(\boldsymbol{\Sigma}\): the SCM, the shrinkage Ledoit-Wolf estimator, and the shrinkage max-SR estimator:

GMVP to compare shrinkage estimator for \(\boldsymbol{\Sigma}\)

We start with the Global Minimum Variance Portfolio (GMVP) since it does not require \(\boldsymbol{\mu}\) and we can focus on the estimation of \(\boldsymbol{\Sigma}\): \[ \mathbf{w}_\textsf{GMVP} = \frac{1}{\mathbf{1}^T\boldsymbol{\Sigma}^{-1}\mathbf{1}}\boldsymbol{\Sigma}^{-1}\mathbf{1} \]

We can see that in-sample the SCM performs the best in terms of standard deviation, obviously, while the others are similar or worse. However, out-of-sample the performance comparison is reversed and the shrinkage methods are better.

Let’s plot the wealth evolution over time:

One can see the clear superiority of the portfolio based on the shrinkage estimate whose \(\rho\) is computed using RMT in order to maximize the Sharpe ratio of the out-of-sample.

Max-Sharpe ratio portfolio to compare shrinkage estimators for \(\boldsymbol{\mu}\) and \(\boldsymbol{\Sigma}\)

Finally, let’s repeat a similar experiment but designing instead a maximum Sharpe ratio portfolio: \[ \begin{array}{ll} \underset{\mathbf{w}}{\textsf{minimize}} & \tilde{\mathbf{w}}^T\mathbf{\Sigma}\tilde{\mathbf{w}}\\ {\textsf{subject to}} & \tilde{\mathbf{w}}^T\boldsymbol{\mu} = 1\\ & \tilde{\mathbf{w}}\ge\mathbf{0} \end{array} \] and then \(\mathbf{w} = \tilde{\mathbf{w}}/(\mathbf{1}^T\tilde{\mathbf{w}})\). We will also consider the James-Stein shrinkage estimator for \(\boldsymbol{\mu}\):

library(CVXR)
# create function for the max-SR design
portolioMaxSharpeRatio <- function(mu, Sigma) {
  w_ <- Variable(nrow(Sigma))
  prob <- Problem(Minimize(quad_form(w_, Sigma)),
                  constraints = list(w_ >= 0, t(mu) %*% w_ == 1))
  result <- solve(prob)
  return(as.vector(result$getValue(w_)/sum(result$getValue(w_))))
}

# compute seven versions of max-SR portfolio
w_maxSR_scm <- portolioMaxSharpeRatio(mu_sm, Sigma_scm)
w_maxSR_LW <- portolioMaxSharpeRatio(mu_sm, Sigma_LW)
w_maxSR_maxSR <- portolioMaxSharpeRatio(mu_sm, Sigma_maxSR)
w_maxSR_JS_scm <- portolioMaxSharpeRatio(mu_JS, Sigma_scm)
w_maxSR_JS_LW <- portolioMaxSharpeRatio(mu_JS, Sigma_LW)
w_maxSR_JS_maxSR <- portolioMaxSharpeRatio(mu_JS, Sigma_maxSR)
w_JS_maxSR_JS_maxSR <- portolioMaxSharpeRatio(mu_JS, Sigma_JS_maxSR)
w_maxSR_all <- cbind(w_maxSR_scm, w_maxSR_LW, w_maxSR_maxSR,
                     w_maxSR_JS_scm, w_maxSR_JS_LW, w_maxSR_JS_maxSR, w_JS_maxSR_JS_maxSR)

# compute returns of the three portfolios
ret_maxSR_all <- xts(X_lin %*% w_maxSR_all, index(X_lin))
ret_maxSR_all_trn <- ret_maxSR_all[1:T_trn, ]
ret_maxSR_all_tst <- ret_maxSR_all[-c(1:T_trn), ]

# performance
table.AnnualizedReturns(ret_maxSR_all_trn)
#>                           w_maxSR_scm w_maxSR_LW w_maxSR_maxSR
#> Annualized Return              0.4534     0.4533        0.3032
#> Annualized Std Dev             0.1634     0.1632        0.1462
#> Annualized Sharpe (Rf=0%)      2.7741     2.7783        2.0742
#>                           w_maxSR_JS_scm w_maxSR_JS_LW w_maxSR_JS_maxSR
#> Annualized Return                 0.4317        0.4319           0.2817
#> Annualized Std Dev                0.1569        0.1568           0.1432
#> Annualized Sharpe (Rf=0%)         2.7508        2.7552           1.9679
#>                           w_JS_maxSR_JS_maxSR
#> Annualized Return                      0.2008
#> Annualized Std Dev                     0.1478
#> Annualized Sharpe (Rf=0%)              1.3592
table.AnnualizedReturns(ret_maxSR_all_tst)
#>                           w_maxSR_scm w_maxSR_LW w_maxSR_maxSR
#> Annualized Return              0.2170     0.2174        0.2555
#> Annualized Std Dev             0.0937     0.0941        0.0824
#> Annualized Sharpe (Rf=0%)      2.3160     2.3100        3.1010
#>                           w_maxSR_JS_scm w_maxSR_JS_LW w_maxSR_JS_maxSR
#> Annualized Return                 0.2197        0.2217           0.2483
#> Annualized Std Dev                0.0897        0.0901           0.0794
#> Annualized Sharpe (Rf=0%)         2.4504        2.4597           3.1255
#>                           w_JS_maxSR_JS_maxSR
#> Annualized Return                      0.2283
#> Annualized Std Dev                     0.0776
#> Annualized Sharpe (Rf=0%)              2.9426

Again, we can see totally different result in-sample and out-of-sample. In-sample, the portfolio based on the SCM has top performance, but this vanishes out-of-sample where the portfolio based on the max-SR shines in terms of Sharpe ratio (as expected!).

Let’s plot the wealth evolution over time:

The top performers by a huge difference are the ones in which \(\rho\) is obtained via the max-SR criterion. And among those three, the best is the one that fully uses the James-Stein’s estimator. Now let’s look at both the GMVP and max-SR portfolios together:

Indeed, using shrinkage clearly improves the performance and the maximum Sharpe ratio design improves over the GMVP which ignores \(\boldsymbol{\mu}\).

Black-Litterman model

The Black-Litterman model was created by Fisher Black and Robert Litterman in 1992 to resolve shortcomings of traditional Markovitz mean-variance asset allocation model. It allows to incorporate discretionary views on the expected returns to the estimation from historical data.

The log-returns are assumed normally distributed: \[ \mathbf{x}\sim\mathcal{N}(\boldsymbol{\mu}, \boldsymbol{\Sigma}), \] with the expected returns modeled as a random variable normally distributed \[ \boldsymbol{\mu}\sim\mathcal{N}(\boldsymbol{\pi}, \tau\boldsymbol{\Sigma}). \] We will consider two options for the market equilibrium \(\boldsymbol{\pi}\):

  • the sample mean: \(\boldsymbol{\pi} = \frac{1}{T}\sum_{t=1}^T \mathbf{x}_t\)
  • based on CAPM: \(\boldsymbol{\pi} = \delta \boldsymbol{\Sigma} \mathbf{w}_M\), where \(\delta = (E[R_{M,t}]-r_f)/\textsf{var}[R_{M,t}]\) is the risk aversion, \(\boldsymbol{\Sigma}\) is the covariance matrix of the returns, and \(\mathbf{w}_M\) is the normalized market capitalization.

The Black-Litterman model also introduces a mechanism to incorporate investor’s views into the input assumptions. The views are modeled as \[ \mathbf{P}\boldsymbol{\mu}\sim\mathcal{N}(\mathbf{v}, \boldsymbol{\Omega}). \] Employing Bayes formula one can finally find the posterior for the returns: \[ \mathbf{x} | \mathbf{v}, \boldsymbol{\Omega} \sim\mathcal{N}(\boldsymbol{\mu}_\textsf{BL}, \boldsymbol{\Sigma}_\textsf{BL}) \] where \[ \begin{align} \boldsymbol{\mu}_\textsf{BL} &= \boldsymbol{\pi} + \tau\boldsymbol{\Sigma}\mathbf{P}^T(\tau\mathbf{P}\boldsymbol{\Sigma}\mathbf{P}^T + \boldsymbol{\Omega})^{-1}(\mathbf{v} - \mathbf{P}\boldsymbol{\pi})\\ \boldsymbol{\Sigma}_\textsf{BL} &= (1+\tau)\boldsymbol{\Sigma} - \tau^2\boldsymbol{\Sigma}\mathbf{P}^T(\tau\mathbf{P}\boldsymbol{\Sigma}\mathbf{P}^T + \boldsymbol{\Omega})^{-1}\mathbf{P}\boldsymbol{\Sigma} \end{align} \]

Let’s start by loading the training and test sets:

Now let’s estimate the different expected returns \(\boldsymbol{\pi}\) with the training data:

Now we are ready to use Black-Litterman to compute the posterior mean and covariance of the returns:

# Black-Litterman
compute_BL <- function(pi, c_sweep) {
  tau <- 1/T_trn
  P <- rbind(c(0,0,0, 0,0, 0,0,1),
             c(0,1,0, 0,0, 0,0,0),
             c(0,0,0, 0,1, 0,0,0))
  v <- c(-0.00027, -0.00054, 0.00016)
  mu_BL <- NULL
  Sigma_BL <- list()
  for (c in c_sweep+1e-6) {
    Omega <- (1/c) * P %*% Sigma %*% t(P)
    mu_ <- pi + tau * Sigma %*% t(P) %*% solve(tau * P %*% Sigma %*% t(P) + Omega) %*% (v - P %*% pi)
    Sigma_ <- (1+tau)*Sigma - 
      (tau^2) * Sigma %*% t(P) %*% solve(tau * P %*% Sigma %*% t(P) + Omega) %*% P %*% Sigma
    mu_BL <- cbind(mu_BL, mu_)
    Sigma_BL <- c(Sigma_BL, list(Sigma_))
    }
  colnames(mu_BL) <- paste("c =", round(c_sweep))
  names(Sigma_BL) <- paste("c =", round(c_sweep))
  return(list(mu = mu_BL, Sigma = Sigma_BL))
}
c_sweep <- seq(0, T_trn/4, length.out = 10)
moments_sm_BL <- compute_BL(pi = mu_sm, c_sweep)
moments_capm_BL <- compute_BL(pi = mu_capm, c_sweep)
names(moments_sm_BL)
#> [1] "mu"    "Sigma"
print(moments_sm_BL$mu[, 1:5])
#>              c = 0        c = 29        c = 59        c = 88       c = 117
#> AAPL  0.0008904398  0.0008837597  0.0008774312  0.0008714271  0.0008657233
#> AMD  -0.0008874288 -0.0008780388 -0.0008691431 -0.0008607035 -0.0008526859
#> ADI   0.0005617339  0.0005565590  0.0005516565  0.0005470053  0.0005425867
#> AEZS -0.0012899591 -0.0012947978 -0.0012993818 -0.0013037307 -0.0013078622
#> A     0.0005742533  0.0005630573  0.0005524505  0.0005423877  0.0005328280
#> APD   0.0004674852  0.0004603498  0.0004535899  0.0004471767  0.0004410841
#> AA   -0.0002731132 -0.0002850071 -0.0002962750 -0.0003069651 -0.0003171207
#> CF    0.0011137992  0.0010763992  0.0010409676  0.0010073531  0.0009754193
names(moments_sm_BL$Sigma)
#>  [1] "c = 0"   "c = 29"  "c = 59"  "c = 88"  "c = 117" "c = 147" "c = 176"
#>  [8] "c = 205" "c = 235" "c = 264"
print(moments_sm_BL$Sigma[[1]])
#>              AAPL          AMD          ADI         AEZS            A
#> AAPL 3.044637e-04 0.0001949531 0.0001099570 6.153777e-05 0.0001438432
#> AMD  1.949531e-04 0.0010487974 0.0002788809 2.231771e-04 0.0003246761
#> ADI  1.099570e-04 0.0002788809 0.0002446124 1.479887e-04 0.0002107730
#> AEZS 6.153777e-05 0.0002231771 0.0001479887 3.335987e-03 0.0001922181
#> A    1.438432e-04 0.0003246761 0.0002107730 1.922181e-04 0.0004599589
#> APD  9.058100e-05 0.0001993283 0.0001257414 1.118831e-04 0.0001846213
#> AA   1.498914e-04 0.0003429236 0.0001971687 1.845851e-04 0.0002834053
#> CF   1.505383e-04 0.0002595812 0.0001560117 1.312042e-04 0.0002236725
#>               APD           AA           CF
#> AAPL 0.0000905810 0.0001498914 0.0001505383
#> AMD  0.0001993283 0.0003429236 0.0002595812
#> ADI  0.0001257414 0.0001971687 0.0001560117
#> AEZS 0.0001118831 0.0001845851 0.0001312042
#> A    0.0001846213 0.0002834053 0.0002236725
#> APD  0.0001978205 0.0001848751 0.0001532403
#> AA   0.0001848751 0.0004660323 0.0002628627
#> CF   0.0001532403 0.0002628627 0.0005782828

Finally, let’s compare the different estimations in terms of \(\boldsymbol{\mu}\) in the test data:

We can further design some portfolio for comparison:

library(PerformanceAnalytics)
library(CVXR)
portolioMarkowitz <- function(mu, Sigma, lmd = 0.5) {
  w <- Variable(nrow(Sigma))
  prob <- Problem(Maximize(t(mu) %*% w - lmd*quad_form(w, Sigma)),
                  constraints = list(w >= 0, sum(w) == 1))
  result <- solve(prob)
  return(as.vector(result$getValue(w)))
}

# portfolio based on BL with sample mean
w_sm_BL <- NULL
for (i in 1:length(c_sweep))
  w_sm_BL <- cbind(w_sm_BL, portolioMarkowitz(moments_sm_BL$mu[, i], moments_sm_BL$Sigma[[i]]))
colnames(w_sm_BL) <- colnames(moments_sm_BL$mu)
rownames(w_sm_BL) <- rownames(moments_sm_BL$mu)
w_sm_BL
#>             c = 0       c = 29       c = 59       c = 88      c = 117
#> AAPL 3.513616e-01 4.042001e-01 4.541870e-01 5.015850e-01 5.467752e-01
#> AMD  1.351994e-09 4.028070e-09 4.929919e-10 3.634594e-09 1.678314e-09
#> ADI  9.195681e-09 2.628387e-08 2.352472e-09 2.306694e-08 1.252940e-08
#> AEZS 1.175323e-09 3.435259e-09 3.882318e-10 3.044467e-09 1.421217e-09
#> A    7.912455e-09 2.182062e-08 1.805070e-09 1.846824e-08 1.004838e-08
#> APD  6.976509e-09 2.042599e-08 2.205070e-09 1.811402e-08 8.987464e-09
#> AA   2.073571e-09 6.107100e-09 7.467659e-10 5.399313e-09 2.482686e-09
#> CF   6.486384e-01 5.957999e-01 5.458130e-01 4.984149e-01 4.532247e-01
#>           c = 147      c = 176      c = 205      c = 235      c = 264
#> AAPL 5.896376e-01 6.304691e-01 6.694370e-01 7.066816e-01 7.421564e-01
#> AMD  1.779011e-09 4.874256e-10 1.179632e-09 2.657274e-09 1.608326e-09
#> ADI  1.603215e-08 4.754454e-09 1.010394e-08 2.320882e-08 1.470445e-08
#> AEZS 1.527581e-09 4.188994e-10 9.776744e-10 2.174912e-09 1.318776e-09
#> A    1.192012e-08 3.373305e-09 7.505394e-09 1.684865e-08 1.059850e-08
#> APD  1.028504e-08 2.891456e-09 6.765867e-09 1.545276e-08 9.384703e-09
#> AA   2.579968e-09 6.989491e-10 1.711054e-09 3.842002e-09 2.296500e-09
#> CF   4.103623e-01 3.695309e-01 3.305629e-01 2.933184e-01 2.578436e-01
chart.StackedBar(t(w_sm_BL), ylab = "w", space=0, border = NA,
                 main = "Portfolio (based on BL with sample mean) as a function of uncertainty")

# portfolio based on BL with CAMP
w_capm_BL <- NULL
for (i in 1:length(c_sweep))
  w_capm_BL <- cbind(w_capm_BL, portolioMarkowitz(moments_capm_BL$mu[, i], moments_capm_BL$Sigma[[i]]))
colnames(w_capm_BL) <- colnames(moments_capm_BL$mu)
rownames(w_capm_BL) <- rownames(moments_capm_BL$mu)
w_capm_BL
#>             c = 0       c = 29       c = 59       c = 88      c = 117
#> AAPL 5.238432e-01 5.323416e-01 5.367162e-01 5.383161e-01 5.388955e-01
#> AMD  6.741259e-02 4.163119e-02 1.724804e-02 3.662968e-07 2.658125e-07
#> ADI  1.792836e-02 3.110941e-02 3.879016e-02 4.103501e-02 3.411754e-02
#> AEZS 2.101271e-01 2.106285e-01 2.109787e-01 2.112547e-01 2.114950e-01
#> A    1.046323e-01 1.263009e-01 1.444388e-01 1.595406e-01 1.710667e-01
#> APD  4.605043e-07 4.426865e-06 6.284619e-03 1.272384e-02 2.017782e-02
#> AA   4.889685e-02 4.950629e-02 4.553904e-02 3.712926e-02 2.424671e-02
#> CF   2.715913e-02 8.477662e-03 4.403889e-06 1.657328e-07 4.563013e-07
#>           c = 147      c = 176      c = 205      c = 235      c = 264
#> AAPL 5.392697e-01 5.396670e-01 5.390785e-01 5.383159e-01 5.370888e-01
#> AMD  2.287161e-07 7.350301e-07 9.422779e-08 4.500727e-07 1.002073e-07
#> ADI  2.765112e-02 2.139788e-02 1.257194e-02 4.798770e-03 8.456774e-06
#> AEZS 2.117087e-01 2.119158e-01 2.120300e-01 2.121210e-01 2.121576e-01
#> A    1.822538e-01 1.927790e-01 1.999597e-01 2.064839e-01 2.119259e-01
#> APD  2.724840e-02 3.386259e-02 3.635360e-02 3.827302e-02 3.881765e-02
#> AA   1.186767e-02 3.758853e-04 5.910437e-06 6.232167e-06 1.231510e-06
#> CF   3.630449e-07 1.153990e-06 1.595844e-07 7.644872e-07 1.745917e-07
chart.StackedBar(t(w_capm_BL), ylab = "w", space=0, border = NA,
                 main = "Portfolio (based on BL with CAPM) as a function of uncertainty")

# performance
table.AnnualizedReturns(ret_sm_BL_trn)
#>                            c = 0 c = 29 c = 59 c = 88 c = 117 c = 147
#> Annualized Return         0.2421 0.2419 0.2412 0.2401  0.2387  0.2370
#> Annualized Std Dev        0.2966 0.2872 0.2792 0.2728  0.2676  0.2636
#> Annualized Sharpe (Rf=0%) 0.8163 0.8424 0.8637 0.8802  0.8919  0.8988
#>                           c = 176 c = 205 c = 235 c = 264
#> Annualized Return          0.2350  0.2329  0.2307  0.2283
#> Annualized Std Dev         0.2608  0.2589  0.2579  0.2577
#> Annualized Sharpe (Rf=0%)  0.9013  0.8997  0.8943  0.8858
table.AnnualizedReturns(ret_sm_BL_tst)
#>                            c = 0 c = 29 c = 59 c = 88 c = 117 c = 147
#> Annualized Return         0.0271 0.0425 0.0569 0.0705  0.0834  0.0954
#> Annualized Std Dev        0.2390 0.2313 0.2252 0.2205  0.2171  0.2148
#> Annualized Sharpe (Rf=0%) 0.1136 0.1838 0.2528 0.3198  0.3840  0.4443
#>                           c = 176 c = 205 c = 235 c = 264
#> Annualized Return          0.1068  0.1176  0.1279  0.1375
#> Annualized Std Dev         0.2135  0.2132  0.2136  0.2147
#> Annualized Sharpe (Rf=0%)  0.5004  0.5519  0.5986  0.6404
table.AnnualizedReturns(ret_capm_BL_trn)
#>                            c = 0 c = 29 c = 59 c = 88 c = 117 c = 147
#> Annualized Return         0.0131 0.0208 0.0300 0.0385  0.0411  0.0435
#> Annualized Std Dev        0.2905 0.2892 0.2879 0.2871  0.2873  0.2875
#> Annualized Sharpe (Rf=0%) 0.0449 0.0718 0.1043 0.1341  0.1429  0.1513
#>                           c = 176 c = 205 c = 235 c = 264
#> Annualized Return          0.0457  0.0455  0.0453  0.0451
#> Annualized Std Dev         0.2878  0.2882  0.2886  0.2890
#> Annualized Sharpe (Rf=0%)  0.1590  0.1580  0.1569  0.1560
table.AnnualizedReturns(ret_capm_BL_tst)
#>                             c = 0  c = 29  c = 59  c = 88 c = 117 c = 147
#> Annualized Return         -0.2860 -0.2810 -0.2764 -0.2732 -0.2722 -0.2713
#> Annualized Std Dev         0.3442  0.3442  0.3439  0.3436  0.3434  0.3432
#> Annualized Sharpe (Rf=0%) -0.8309 -0.8163 -0.8038 -0.7953 -0.7928 -0.7904
#>                           c = 176 c = 205 c = 235 c = 264
#> Annualized Return         -0.2704 -0.2709 -0.2714 -0.2717
#> Annualized Std Dev         0.3432  0.3434  0.3436  0.3437
#> Annualized Sharpe (Rf=0%) -0.7880 -0.7889 -0.7899 -0.7907

chart.CumReturns(ret_sm_BL_tst, 
            main = "Portfolios (based on BL with sample mean) as a function of uncertaintly of views", 
            wealth.index = TRUE, legend.loc = "topleft", colorset = rich10equal)