We will load some stock market data and divide it into a training part (for the estimation of the expected return \(\boldsymbol{\mu}\) and covariance matrix \(\boldsymbol{\Sigma}\), and subsequent portfolio design) and a test part (for the out-of-sample performance evaluation).

In particular, we will start by loading some stock data from three different sectors:

• sector technology: AAPL, AMD, ADI.
• sector healthcare: ABBV, AEZS, A.
• sector basic materials: APD, AA, CF.

library(xts)  # to manipulate time series of stock data
library(quantmod)  # to download stock data
library(PerformanceAnalytics)  # to compute performance measures

# download data from YahooFinance
stock_namelist <- c("AAPL", "AMD", "ADI",  "ABBV", "AEZS", "A",  "APD", "AA","CF")
prices <- xts()
for (i in 1:length(stock_namelist)) {
  tmp <- Ad(getSymbols(stock_namelist[i], from = "2013-01-01", to = "2018-12-31", auto.assign = FALSE))
  tmp <- na.approx(tmp, na.rm = FALSE)  # interpolate NAs
  prices <- cbind(prices, tmp)
}
colnames(prices) <- stock_namelist
indexClass(prices) <- "Date"

str(prices)

R>> An 'xts' object on 2013-01-02/2016-12-30 containing:
R>>   Data: num [1:1008, 1:9] 55.5 54.8 53.2 52.9 53.1 ...
R>>  - attr(*, "dimnames")=List of 2
R>>   ..$ : NULL
R>>   ..$ : chr [1:9] "AAPL" "AMD" "ADI" "ABBV" ...
R>>   Indexed by objects of class: [Date] TZ: UTC
R>>   xts Attributes:  
R>>  NULL

head(prices)

R>>                AAPL  AMD      ADI     ABBV AEZS        A      APD       AA       CF
R>> 2013-01-02 55.47129 2.53 37.50844 27.60817  253 27.97831 66.43760 20.62187 29.44429
R>> 2013-01-03 54.77112 2.49 36.90318 27.38020  254 28.07852 66.20545 20.80537 29.30506
R>> 2013-01-04 53.24548 2.59 36.24679 27.03431  257 28.63300 67.09529 21.24121 29.96146
R>> 2013-01-07 52.93227 2.67 36.35761 27.08934  259 28.42591 67.03341 20.87419 29.84922
R>> 2013-01-08 53.07474 2.67 35.98254 26.49976  255 28.19877 67.15720 20.87419 29.41162
R>> 2013-01-09 52.24524 2.63 35.88875 26.64913  258 28.96036 68.06254 20.82831 30.44025

tail(prices)

R>>                AAPL   AMD      ADI     ABBV AEZS        A      APD    AA       CF
R>> 2016-12-22 111.8445 11.60 70.07727 55.82264 4.10 44.91877 136.0360 29.75 27.82244
R>> 2016-12-23 112.0657 11.58 70.44936 56.43826 4.10 45.14350 136.3992 29.71 28.37686
R>> 2016-12-27 112.7774 12.07 70.89778 56.58311 4.05 45.44642 137.2373 29.65 29.52265
R>> 2016-12-28 112.2965 11.55 70.18222 56.37488 3.55 44.67448 135.1328 29.43 29.22696
R>> 2016-12-29 112.2676 11.59 70.20129 56.79134 3.60 44.72544 135.2359 28.89 29.47645
R>> 2016-12-30 111.3924 11.34 69.28539 56.69175 3.60 44.64704 134.7206 28.08 29.08836

# compute data matrix with log-returns
X <- diff(log(prices))[-1]
N <- ncol(X)  # number of stocks
T <- nrow(X)  # number of days

We can take a look at the prices of the stocks:

plot(prices/rep(prices[1, ], each = nrow(prices)), col = rainbow10equal, legend.loc = "topleft",
     main = "Normalized prices")

We can now use obtain the sample estimates from the returns \(\mathbf{x}_t\) (i.e., sample means and sample covariance matrix) as \[ \begin{align} \hat{\boldsymbol{\mu}} & = \frac{1}{T}\sum_{t=1}^T \mathbf{x}_t\\ \hat{\boldsymbol{\Sigma}} & = \frac{1}{T-1}\sum_{t=1}^T (\mathbf{x}_t - \hat{\boldsymbol{\mu}})(\mathbf{x}_t - \hat{\boldsymbol{\mu}})^T \end{align} \]

mu <- colMeans(X)
Sigma <- cov(X)

Now that we have some realistic parameters \(\boldsymbol{\mu}\) and \(\boldsymbol{\Sigma}\), we will design a portfolio with some small noise in the parameters to observe the effect.

Let’s start with the global maximum return portfolio (GMRP) formulated as \[ \begin{array}{ll} \underset{\mathbf{w}}{\textsf{maximize}} & \boldsymbol{\mu}^T\mathbf{w}\\ \textsf{subject to} & \mathbf{1}^T\mathbf{w} = 1\\ & \mathbf{w}\ge\mathbf{0}. \end{array} \] Of course the optimal solution is to allocate all the budget on the single stock with largest expected return, not very diversified. We can still implement it using a solver:

library(CVXR)

portfolioMaxReturn <- function(mu) {
  w <- Variable(length(mu))
  prob <- Problem(Maximize(t(w) %*% mu), 
                  constraints = list(w >= 0, sum(w) == 1))
  result <- solve(prob)
  return(as.vector(result$getValue(w)))
}

w_GMRP <- portfolioMaxReturn(mu)
names(w_GMRP) <- colnames(X)
w_all_GMRP <- cbind(w_GMRP)

# plot allocation
barplot(t(w_all_GMRP), col = rainbow8equal[1], legend = colnames(w_all_GMRP), beside = TRUE,
        main = "Global maximum return portfolio allocation", xlab = "stocks", ylab = "dollars")

As expected, all the budget is allocated to a single stock.

Let’s now introduce some estimation error in \(\boldsymbol{\mu}\) to see the effect:

library(mvtnorm)

# generate a few noisy estimations of mu
set.seed(357)
for (i in 1:6) {
  mu_noisy <- as.vector(rmvnorm(n = 1, mean = mu, sigma = (1/T)*Sigma))
  w_GMRP_noisy <- portfolioMaxReturn(mu_noisy)
  w_all_GMRP <- cbind(w_all_GMRP, w_GMRP_noisy)
}

# plot to compare the allocations
barplot(t(w_all_GMRP), col = rainbow8equal[1:7], legend = colnames(w_all_GMRP), beside = TRUE,
        main = "Global maximum return portfolio allocation", xlab = "stocks", ylab = "dollars")

We can see that the optimal allocation changes totally from realization to realization, which is highly undesirable (in fact, it is totally unacceptable).

Let’s continue with the GMVP formulated as \[ \begin{array}{ll} \underset{\mathbf{w}}{\textsf{minimize}} & \mathbf{w}^T\mathbf{\Sigma}\mathbf{w}\\ {\textsf{subject to}} & \mathbf{1}^T\mathbf{w} = 1\\ & \mathbf{w}\ge\mathbf{0} \end{array} \]

portfolioGMVP <- function(Sigma) {
  w <- Variable(nrow(Sigma))
  prob <- Problem(Minimize(quad_form(w, Sigma)), 
                  constraints = list(w >= 0, sum(w) == 1))
  result <- solve(prob)
  return(as.vector(result$getValue(w)))
}

w_GMVP <- portfolioGMVP(Sigma)
names(w_GMVP) <- colnames(X)
w_all_GMVP <- cbind(w_GMVP)

# plot allocation
barplot(t(w_all_GMVP), col = rainbow8equal[1], legend = colnames(w_all_GMVP), beside = TRUE,
        main = "Global minimum variance portfolio allocation", xlab = "stocks", ylab = "dollars")

The GMVP is totally different from the maximum return portfolio. It’s quite the opposite: the budget is allocated to most of the stocks to diversity the risk.

Let’s now introduce some estimation error in \(\boldsymbol{\Sigma}\) to see the effect:

library(PerformanceAnalytics)
# generate a few noisy estimations of Sigma
set.seed(357)
for (i in 1:6) {
  X_noisy <- rmvnorm(n = T, mean = rep(0, N), sigma = Sigma)
  Sigma_noisy <- cov(X_noisy)
  
  w_GMVP_noisy <- portfolioGMVP(Sigma_noisy)
  w_all_GMVP <- cbind(w_all_GMVP, w_GMVP_noisy)
}

# plot to compare the allocations
barplot(t(w_all_GMVP), col = rainbow8equal[1:7], legend = colnames(w_all_GMVP), beside = TRUE,
        main = "Global minimum variance portfolio allocation", xlab = "stocks", ylab = "dollars")

chart.StackedBar(t(w_all_GMVP), main = "GMVP allocation", ylab = "w", space = 0, border = NA)

We can see that the optimal allocation still changes from realization to realization, although in this case is not that extreme.

Finally, let’s take a look at Markowitz’s mean-variance portfolio formulated as \[ \begin{array}{ll} \underset{\mathbf{w}}{\textsf{maximize}} & \boldsymbol{\mu}^T\mathbf{w} -\lambda\mathbf{w}^T\mathbf{\Sigma}\mathbf{w}\\ {\textsf{subject to}} & \mathbf{1}^T\mathbf{w} = 1\\ & \mathbf{w}\ge\mathbf{0} \end{array} \]

portfolioMarkowitz <- 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)))
}

w_Markowitz <- portfolioMarkowitz(mu, Sigma)
names(w_Markowitz) <- colnames(X)
w_all_Markowitz <- cbind(w_Markowitz)

# plot allocation
barplot(t(w_all_Markowitz), col = rainbow8equal[1], legend = colnames(w_all_Markowitz), beside = TRUE,
        main = "Markowitz portfolio allocation", xlab = "stocks", ylab = "dollars")

The mean-variance portfolio seems to be not very diversified. This is one of the reasons why practitioners do not use it. The other reason is the sensitivity or lack of robustness to the estimation of the parameters as we will explore next.

# generate a few noisy estimations of mu and Sigma
set.seed(357)
for (i in 1:6) {
  X_noisy <- rmvnorm(n = T, mean = mu, sigma = Sigma)
  mu_noisy <- colMeans(X_noisy)
  Sigma_noisy <- cov(X_noisy)
  
  w_Markowitz_noisy <- portfolioMarkowitz(mu_noisy, Sigma_noisy)
  w_all_Markowitz <- cbind(w_all_Markowitz, w_Markowitz_noisy)
}

# plot to compare the allocations
barplot(t(w_all_Markowitz), col = rainbow8equal[1:7], legend = colnames(w_all_Markowitz), beside = TRUE,
        main = "Markowitz portfolio allocation", xlab = "stocks", ylab = "dollars")

chart.StackedBar(t(w_all_Markowitz), main = "Markowitz portfolio allocation", ylab = "w", space = 0, border = NA)

Again, the allocations are totally different from realization to realization. Totally unacceptable. Who is going to trust these allocations to invest their money?

The robust formulation of the global maximum return portfolio (GMRP) reads \[ \begin{array}{ll} \underset{\mathbf{w}}{\textsf{maximize}} & \underset{\boldsymbol{\mu}\in\mathcal{U}_{\boldsymbol{\mu}}}{\min} \mathbf{w}^T\boldsymbol{\mu}\\ {\textsf{subject to}} & \mathbf{1}^T\mathbf{w} = 1\\ & \mathbf{w}\ge\mathbf{0}. \end{array} \]

• For an elliptical uncertainty set \[ \mathcal{U}_\boldsymbol{\mu}^e = \{\boldsymbol{\mu}=\hat{\boldsymbol{\mu}}+\kappa\mathbf{S}^{1/2}\mathbf{u}\ | \ \left\Vert \mathbf{u}\right\Vert_{2}\leq 1\}, \] where typically \(\mathbf{S}=\hat{\mathbf{\Sigma}}\), the robust formulation becomes \[ \begin{array}{ll} \underset{\mathbf{w}}{\textsf{maximize}} & \mathbf{w}^T\hat{\boldsymbol{\mu}} - \kappa\left\Vert \mathbf{S}^{1/2}\mathbf{w}\right\Vert _{2}\\ {\textsf{subject to}} & \mathbf{1}^T\mathbf{w} = 1\\ & \mathbf{w}\ge\mathbf{0}. \end{array} \]

• For a box uncertainty set \[ \mathcal{U}_\boldsymbol{\mu}^b = \{\boldsymbol{\mu}\ | -\boldsymbol{\delta}\leq\boldsymbol{\mu}-\hat{\boldsymbol{\mu}}\leq\boldsymbol{\delta}\}, \] where typically \(\delta_i = \kappa\hat{\sigma}_i\), the robust formulation becomes \[ \begin{array}{ll} \underset{\mathbf{w}}{\textsf{maximize}} & \mathbf{w}^T\hat{\boldsymbol{\mu}} - |\mathbf{w}|^T\boldsymbol{\delta}\\ {\textsf{subject to}} & \mathbf{1}^T\mathbf{w} = 1\\ & \mathbf{w}\ge\mathbf{0}. \end{array} \]

Let’s start implementing the robust solution for the box constraint:

portfolioMaxReturnRobustBox <- function(mu_hat, delta) {
  w <- Variable(length(mu_hat))
  prob <- Problem(Maximize( t(w) %*% mu_hat - t(abs(w)) %*% delta ), 
                  constraints = list(w >= 0, sum(w) == 1))
  result <- solve(prob)
  return(as.vector(result$getValue(w)))
}

# clairvoyant solution
w_GMRP <- portfolioMaxReturn(mu)
names(w_GMRP) <- colnames(X)
w_all_GMRP_robust_box <- cbind(w_GMRP)

# multiple robust solutions
kappa <- 0.1
delta <- kappa*sqrt(diag(Sigma_noisy))
set.seed(357)
for (i in 1:6) {
  X_noisy <- rmvnorm(n = T, mean = mu, sigma = Sigma)
  mu_noisy <- colMeans(X_noisy)
  Sigma_noisy <- cov(X_noisy)
  
  w_GMRP_robust_box_noisy <- portfolioMaxReturnRobustBox(mu_noisy, delta)
  w_all_GMRP_robust_box <- cbind(w_all_GMRP_robust_box, w_GMRP_robust_box_noisy)
}

# plot to compare the allocations
barplot(t(w_all_GMRP_robust_box), col = rainbow8equal[1:7], legend = colnames(w_all_GMRP_robust_box), 
        beside = TRUE, args.legend = list(bg = "white"),
        main = "Robust (box) global maximum return portfolio allocation", xlab = "stocks", ylab = "dollars")

Not much improvement with respect to the naive case. The sensitivity is still there. In fact, for the particular case of \(\boldsymbol{\delta}=\mathbf{1}\), it follows that \(|\mathbf{w}|^T\boldsymbol{\delta}=\mathbf{w}^T\mathbf{1}=1\)!

Let’s try the robust solution for the elliptical constraint:

portfolioMaxReturnRobustEllipsoid <- function(mu_hat, S, kappa = 0.1) {
  S12 <- chol(S)  # t(S12) %*% S12 = Sigma
  w <- Variable(length(mu_hat))
  prob <- Problem(Maximize( t(w) %*% mu_hat - kappa*norm2(S12 %*% w) ), 
                  constraints = list(w >= 0, sum(w) == 1))
  result <- solve(prob)
  return(as.vector(result$getValue(w)))
}

# clairvoyant solution
w_GMRP <- portfolioMaxReturn(mu)
names(w_GMRP) <- colnames(X)
w_all_GMRP_robust_ellipsoid <- cbind(w_GMRP)

# multiple robust solutions
kappa <- 0.2
set.seed(357)
for (i in 1:6) {
  X_noisy <- rmvnorm(n = T, mean = mu, sigma = Sigma)
  mu_noisy <- colMeans(X_noisy)
  Sigma_noisy <- cov(X_noisy)
  
  w_GMRP_robust_ellipsoid_noisy <- portfolioMaxReturnRobustEllipsoid(mu_noisy, Sigma_noisy, kappa)
  w_all_GMRP_robust_ellipsoid <- cbind(w_all_GMRP_robust_ellipsoid, w_GMRP_robust_ellipsoid_noisy)
}

# plot to compare the allocations
barplot(t(w_all_GMRP_robust_ellipsoid), col = rainbow8equal[1:7], legend = colnames(w_all_GMRP_robust_ellipsoid), 
        beside = TRUE, args.legend = list(bg = "white"),
        main = "Robust (ellipsoid) global maximum return portfolio allocation", xlab = "stocks", ylab = "dollars")

chart.StackedBar(t(w_all_GMRP_robust_ellipsoid), 
                 main = "Robust (ellipsoid) global maximum return portfolio allocation", 
                 ylab = "w", space = 0, border = NA)

The elliptical uncertainty set seems to work better for the maximum return portfolio. It definitely provides more diversity compared to the naive solution, as well as being more stable (less sensitive). However, it is still too sensitive.

Summarizing, for the robust global maximum return portfolio, it seems that the ellipsoid uncertainty in \(\boldsymbol{\mu}\) is more stable (and diverse) than the box uncertainly (which is totally unstable). Nevertheless, the allocation is still too sensitive owing to the fact that the global maximum return portfolio is probably not a good formulation as it ignores the covariance matrix.

The robust formulation of the GMVP reads \[ \begin{array}{ll} \underset{\mathbf{w}}{\textsf{minimize}} & \underset{\boldsymbol{\Sigma}\in\mathcal{U}_{\boldsymbol{\Sigma}}}{\max} \mathbf{w}^T\mathbf{\Sigma}\mathbf{w}\\ {\textsf{subject to}} & \mathbf{1}^T\mathbf{w} = 1\\ & \mathbf{w}\ge\mathbf{0}. \end{array} \]

• For a spherical uncertainty set on the data matrix \(\mathbf{X}\) (assuming we use the sample covariance matrix \(\hat{\boldsymbol{\Sigma}}=\frac{1}{T-1}\mathbf{X}^T\mathbf{X}\)) \[ \mathcal{U}_\mathbf{X}^s = \{\mathbf{X}\ | \ \|\mathbf{X}-\hat{\mathbf{X}}\|_F \le \delta\}, \] the robust formulation becomes \[ \begin{array}{ll} \underset{\mathbf{w}}{\textsf{minimize}} & \|\hat{\mathbf{X}}\mathbf{w}\|_2 + \delta \|\mathbf{w}\|_2\\ {\textsf{subject to}} & \mathbf{1}^T\mathbf{w} = 1\\ & \mathbf{w}\ge\mathbf{0}. \end{array} \]

• For a box uncertainty set on the covariance matrix \[ \mathcal{U}_\boldsymbol{\Sigma}^b = \{\boldsymbol{\Sigma}\ | \ \underline{\boldsymbol{\Sigma}} \le \boldsymbol{\Sigma} \le \overline{\boldsymbol{\Sigma}},\ \boldsymbol{\Sigma} \succeq \mathbf{0}\}, \] where a choice for the bounds is \(\underline{\boldsymbol{\Sigma}}=\delta^{-1}\hat{\boldsymbol{\Sigma}}\) and \(\overline{\boldsymbol{\Sigma}}=\delta\hat{\boldsymbol{\Sigma}}\) with \(\delta>1\), the robust formulation becomes

\[ \begin{array}{ll} \underset{\mathbf{w},\underline{\boldsymbol{\Lambda}},\overline{\boldsymbol{\Lambda}}}{\textsf{minimize}} & \textsf{Tr}(\overline{\boldsymbol{\Lambda}}\ \overline{\boldsymbol{\Sigma}}) - \textsf{Tr}(\underline{\boldsymbol{\Lambda}}\ \underline{\boldsymbol{\Sigma}})\\ \textsf{subject to} & \left[\begin{array}{cc}\overline{\boldsymbol{\Lambda}}-\underline{\boldsymbol{\Lambda}} & \mathbf{w}\\ \mathbf{w}^{T} & 1\end{array} \right]\succeq\mathbf{0}\\ & \overline{\boldsymbol{\Lambda}}\ge\mathbf{0},\ \underline{\boldsymbol{\Lambda}}\ge\mathbf{0}\\ & \mathbf{1}^T\mathbf{w} = 1\\ & \mathbf{w}\ge\mathbf{0}. \end{array} \]

• For an elliptical uncertainty set on the covariance matrix \[ \mathcal{U}_{\boldsymbol{\Sigma}}^{e} = \{\boldsymbol{\Sigma}\ |\ (\mathsf{vec}(\boldsymbol{\Sigma})-\mathsf{vec}(\hat{\boldsymbol{\Sigma}}))^{T}\mathbf{S}_{\boldsymbol{\Sigma}}^{-1}(\mathsf{vec}(\boldsymbol{\Sigma})-\mathsf{vec}(\hat{\boldsymbol{\Sigma}}))\leq\delta^{2},\ \boldsymbol{\Sigma}\succeq\mathbf{0}\}, \] where one simple choice for \(\mathbf{S}_{\boldsymbol{\Sigma}}\) is the identity matrix, the robust formulation becomes \[ \begin{array}{ll} \underset{\mathbf{w},\mathbf{X},\mathbf{Z}}{\textsf{minimize}} & \textsf{Tr}\left(\hat{\boldsymbol{\Sigma}}\left(\mathbf{X}+\mathbf{Z}\right)\right) + \delta\left\Vert \mathbf{S}_{\boldsymbol{\Sigma}}^{1/2}\left(\textsf{vec}(\mathbf{X})+\textsf{vec}(\mathbf{Z})\right)\right\Vert _{2}\\ \textsf{subject to} & \mathbf{w}^{T}\mathbf{1}=1,\quad \mathbf{w}\ge\mathbf{0}\\ & \left[\begin{array}{cc}\mathbf{X} & \mathbf{w}\\ \mathbf{w}^{T} & 1\end{array} \right]\succeq\mathbf{0}\\ & \mathbf{Z}\succeq\mathbf{0}. \end{array} \]

Let’s start implementing the robust formulation for the box constraint on \(\boldsymbol{\Sigma}\):

# define function (check: https://cvxr.rbind.io/cvxr_functions/)
portfolioGMVPRobustBox <- function(Sigma_lb, Sigma_ub) {
  N <- nrow(Sigma_lb)
  w         <- Variable(N)
  Lambda_ub <- Variable(N, N)
  Lambda_lb <- Variable(N, N)
  Tmp       <- Variable(N+1, N+1, PSD = TRUE)
  prob <- Problem(Minimize(matrix_trace(Lambda_ub %*% Sigma_ub) - matrix_trace(Lambda_lb %*% Sigma_lb)),
                  constraints = list(w >= 0, sum(w) == 1,
                                     Lambda_ub >= 0, Lambda_lb >= 0,
                                     bmat(list(list(Lambda_ub - Lambda_lb, w),
                                               list(t(w),                  1))) == Tmp))
  result <- solve(prob)
  return(as.vector(result$getValue(w)))
}

# clairvoyant solution
w_GMVP <- portfolioGMVP(Sigma)
names(w_GMVP) <- colnames(X)
w_all_GMVP_robust_box <- cbind(w_GMVP)

# multiple robust solutions
delta <- 1.1
set.seed(357)
for (i in 1:6) {
  X_noisy <- rmvnorm(n = T, mean = mu, sigma = Sigma)
  mu_noisy <- colMeans(X_noisy)
  Sigma_noisy <- cov(X_noisy)
  Sigma_lb <- ifelse(Sigma_noisy >= 0,  (1/delta)*Sigma_noisy, delta*Sigma_noisy)
  Sigma_ub <- ifelse(Sigma_noisy >= 0,  delta*Sigma_noisy, (1/delta)*Sigma_noisy)
  w_GMVP_robust_box_noisy <- portfolioGMVPRobustBox(Sigma_lb, Sigma_ub)
  w_all_GMVP_robust_box <- cbind(w_all_GMVP_robust_box, w_GMVP_robust_box_noisy)
}

# plot to compare the allocations
barplot(t(w_all_GMVP_robust_box), col = rainbow8equal[1:7], legend = colnames(w_all_GMVP_robust_box), 
        beside = TRUE, args.legend = list(bg = "white", x = "topleft"),
        main = "Robust (box) GMVP allocation", xlab = "stocks", ylab = "dollars")

chart.StackedBar(t(w_all_GMVP), main = "Naive GMVP allocation", ylab = "w", space = 0, border = NA)

chart.StackedBar(t(w_all_GMVP_robust_box), 
                 main = "Robust (box) GMVP allocation", ylab = "w", space = 0, border = NA)

The naive GMVP was already kind of stable. Nevertheless, the robust GMVP becomes even more stable (less sensitive to the estimation errors). Better results could be obtained by choosing more appropriate upper and lower bounds on \(\boldsymbol{\Sigma}\).

Let’s consider now the robust formulation for the spherical uncertainty region on the data matrix \(\mathbf{X}\):

# define function
portfolioGMVPRobustSphereX <- function(X_hat, delta) {
  N <- ncol(X_hat)
  w <- Variable(N)
  X_ <- scale(X_hat, center = TRUE, scale = FALSE)  # demean
  prob <- Problem(Minimize(norm2(X_ %*% w) + delta*norm2(w)),
                  constraints = list(w >= 0, sum(w) == 1))
  result <- solve(prob)
  return(as.vector(result$getValue(w)))
}

# clairvoyant solution
w_GMVP <- portfolioGMVP(Sigma)
names(w_GMVP) <- colnames(X)
w_all_GMVP_robust_sphereX <- cbind(w_GMVP)

# multiple robust solutions
delta <- 0.1
set.seed(357)
for (i in 1:6) {
  X_noisy <- rmvnorm(n = T, mean = mu, sigma = Sigma)

  w_GMVP_robust_sphereX_noisy <- portfolioGMVPRobustSphereX(X_noisy, delta)
  w_all_GMVP_robust_sphereX <- cbind(w_all_GMVP_robust_sphereX, w_GMVP_robust_sphereX_noisy)
}

# plot to compare the allocations
barplot(t(w_all_GMVP_robust_sphereX), col = rainbow8equal[1:7], legend = colnames(w_all_GMVP_robust_sphereX), 
        beside = TRUE, args.legend = list(bg = "white", x = "topleft"),
        main = "Robust (sphere in X) GMVP allocation", xlab = "stocks", ylab = "dollars")

chart.StackedBar(t(w_all_GMVP), main = "Naive GMVP allocation", ylab = "w", space = 0, border = NA)

chart.StackedBar(t(w_all_GMVP_robust_sphereX), 
                 main = "Robust (sphere in X) GMVP allocation", ylab = "w", space = 0, border = NA)