# basic finance
library(xts) # to manipulate time series of stock data
library(portfolioBacktest) # to conduct backtests
library(pob) # book package with financial data
library(quantmod) # to download market data from the Internet
library(mvtnorm) # to generate synthetic data
# plotting
library(ggplot2) # for nice plots
library(reshape2) # to reshape data
library(ggridges) # to plot the geom_density_ridges
# optimization
#library(quadprog)
library(CVXR)Portfolio Optimization
Chapter 14: Robust Portfolios
R code
R code examples for Chapter 14 of the book:
Daniel P. Palomar (2024). Portfolio Optimization: Theory and Application. Cambridge University Press.
Loading packages
First load the packages used in the examples:
Introduction
Markowitz’s mean-variance portfolio optimization balances expected return \(\w^\T\bmu\) against risk, measured by variance \(\w^\T\bSigma\w\): \[ \begin{array}{ll} \underset{\w}{\text{maximize}} & \w^\T\bmu - \frac{\lambda}{2}\w^\T\bSigma\w\\ \text{subject to} & \w \in \mathcal{W}, \end{array} \] where \(\bmu\) and \(\bSigma\) are parameters, \(\lambda\) controls risk aversion, and \(\mathcal{W}\) represents constraints like \(\{\w \mid \bm{1}^\T\w=1, \w\ge\bm{0} \}\).
In practice, \(\bmu\) and \(\bSigma\) are estimated from historical data, leading to estimation errors due to limited data and non-stationarity. This noise, especially in \(\hat{\bmu}\), results in erratic portfolio designs, a problem known as the “Markowitz optimization enigma.” Michaud highlighted that portfolio optimization often maximizes estimation errors, rendering “optimal” portfolios financially meaningless.
Next figure shows the sensitivity of the mean-variance portfolio to estimation errors, demonstrating its impracticality due to erratic behavior.
Code
library(pob)
library(CVXR)
library(mvtnorm)
data(SP500_2015to2020)
stock_prices <- SP500_2015to2020$stocks["2019::", c("AAPL", "AMZN", "AMD", "GM", "GOOGL", "MGM", "MSFT", "QCOM", "TSCO", "UPS")]
design_MVP <- function(mu, Sigma, lambda = 1, ub = Inf) {
w <- Variable(nrow(Sigma))
prob <- Problem(Maximize(t(mu) %*% w - (lambda/2)*quad_form(w, Sigma)),
constraints = list(w >= 0, sum(w) == 1, w <= ub))
result <- solve(prob)
w <- as.vector(result$getValue(w))
names(w) <- colnames(Sigma)
return(w)
}
MVP_noisy <- function(dataset, ...) {
X <- diff(log(dataset$prices))[-1]
mu <- colMeans(X)
Sigma <- cov(X)
# noisy version
X_noisy <- rmvnorm(n = nrow(X), mean = mu, sigma = Sigma)
mu_noisy <- colMeans(X_noisy)
Sigma_noisy <- cov(X_noisy)
design_MVP(mu = mu_noisy, Sigma = Sigma_noisy, lambda = 1)
}
# backtests
set.seed(42)
bt <- portfolioBacktest(list("#1" = MVP_noisy,
"#2" = MVP_noisy,
"#3" = MVP_noisy,
"#4" = MVP_noisy,
"#5" = MVP_noisy,
"#6" = MVP_noisy),
list(list("prices" = stock_prices)), price_name = "prices",
lookback = 6*21, optimize_every = 10000, rebalance_every = 10000)Backtesting 6 portfolios over 1 datasets (periodicity = daily data)...Code
listPortfoliosWithFailures(bt)
data.frame(
"stocks" = names(stock_prices),
"#1" = as.numeric(bt$`#1`$data1$w_optimized[1, ]),
"#2" = as.numeric(bt$`#2`$data1$w_optimized[1, ]),
"#3" = as.numeric(bt$`#3`$data1$w_optimized[1, ]),
"#4" = as.numeric(bt$`#4`$data1$w_optimized[1, ]),
"#5" = as.numeric(bt$`#5`$data1$w_optimized[1, ]),
"#6" = as.numeric(bt$`#6`$data1$w_optimized[1, ]),
check.names = FALSE) |>
melt(id.vars = "stocks") |>
ggplot(aes(x = stocks, y = value, fill = variable)) +
geom_bar(stat = "identity", position = "dodge", color = "black", width = 0.8) +
theme(legend.title = element_text()) + labs(fill = "realization") +
labs(title = "Portfolio for different estimation errors", y = "weight")
Robust worst-case portfolio optimization
To study the uncertainty in \(\bmu\), we consider the global maximum return portfolio (GMRP) for an estimated mean vector \(\hat{\bmu}\): \[ \begin{array}{ll} \underset{\w}{\textm{maximize}} & \w^\T\hat{\bmu}\\ \textm{subject to} & \bm{1}^\T\w=1, \quad \w\ge\bm{0}, \end{array} \] which is is known to be terribly sensitive to estimation errors.
We will instead use the worst-case formulation: \[ \begin{array}{ll} \underset{\w}{\textm{maximize}} & \textm{inf}_{\bmu \in \mathcal{U}_{\bmu}} \; \w^\T\bmu\\ \textm{subject to} & \bm{1}^\T\w=1, \quad \w\ge\bm{0}, \end{array} \] where typical choices for the uncertainty region \(\mathcal{U}_{\bmu}\) are the ellipsoidal and box sets considered next.
Ellipsoidal uncertainty set
The ellipsoidal uncertainty set for \(\bmu\) is \[ \mathcal{U}_{\bmu} = \left\{\bmu = \hat{\bmu} + \kappa \bm{S}^{1/2}\bm{u} \mid \|\bm{u}\|_2 \le 1\right\}, \] where \(\bm{S}^{1/2}\) is the symmetric square-root matrix of the shape \(\bm{S}\), e.g., \(\bm{S}=(1/T)\bSigma\), and \(\kappa\) determines the size of the ellipsoid.
The corresponding worst-case value of \(\w^\T\bmu\) is \[ \begin{aligned} \underset{\bmu \in \mathcal{U}_{\bmu}}{\text{inf}} \; \w^\T\bmu &= \w^\T\hat{\bmu} - \kappa\; \|\bm{S}^{1/2}\w\|_2. \end{aligned} \]
Box uncertainty set
The box uncertainty set for \(\bmu\) is \[ \mathcal{U}_{\bmu} = \left\{\bmu \mid -\bm{\delta} \le \bmu - \hat{\bmu} \le \bm{\delta}\right\}, \] where \(\bm{\delta}\) is the half-width of the box in all dimensions.
The corresponding worst-case value of \(\w^\T\bmu\) is \[ \begin{aligned} \underset{\bmu \in \mathcal{U}_{\bmu}}{\text{inf}} \; \w^\T\bmu &= \w^\T\hat{\bmu} - |\w|^\T\bm{\delta}. \end{aligned} \]
Numerical experiments
Sensitivity
First, let’s write the code for the robust mean–variance formulations with an ellipsoidal and box uncertainty sets for the mean vector:
Code
library(CVXR)
design_MVP_robust_mu_box <- function(mu_hat, delta,
Sigma_hat,
lambda = 1) {
w <- Variable(length(mu_hat))
prob <- Problem(Maximize(t(w) %*% mu_hat - t(abs(w)) %*% delta - (lambda/2)*quad_form(w, Sigma_hat)),
constraints = list(w >= 0, sum(w) == 1))
result <- solve(prob)
return(as.vector(result$getValue(w)))
}
design_MVP_robust_mu_ellipsoid <- function(mu_hat, S_mu = Sigma_hat, kappa_mu = 0.1,
Sigma_hat,
lambda = 1) {
S12 <- chol(S_mu) # t(S12) %*% S12 = Sigma
w <- Variable(length(mu_hat))
prob <- Problem(Maximize(t(w) %*% mu_hat - kappa_mu*norm2(S12 %*% w)
- (lambda/2) * (quad_form(w, Sigma_hat))),
constraints = list(w >= 0, sum(w) == 1))
result <- solve(prob)
return(as.vector(result$getValue(w)))
}We can now observe the improved sensitivity of the robust mean–variance portfolio under an ellipsoidal uncertainty set for \(\bm{\mu}\):
Code
library(pob)
library(mvtnorm)
data(SP500_2015to2020)
stock_prices <- SP500_2015to2020$stocks["2019::", c("AAPL", "AMZN", "AMD", "GM", "GOOGL", "MGM", "MSFT", "QCOM", "TSCO", "UPS")]
MVP_noisy_robust_mu_ellipsoid <- function(dataset, ...) {
X <- diff(log(dataset$prices))[-1]
mu <- colMeans(X)
Sigma <- cov(X)
# noisy version
X_noisy <- rmvnorm(n = nrow(X), mean = mu, sigma = Sigma)
mu_noisy <- colMeans(X_noisy)
Sigma_noisy <- cov(X_noisy)
# robust design
design_MVP_robust_mu_ellipsoid(mu_hat = mu_noisy, kappa_mu = 0.5,
Sigma_hat = Sigma_noisy,
lambda = 1)
}
# backtests
set.seed(42)
bt <- portfolioBacktest(list("#1" = MVP_noisy_robust_mu_ellipsoid,
"#2" = MVP_noisy_robust_mu_ellipsoid,
"#3" = MVP_noisy_robust_mu_ellipsoid,
"#4" = MVP_noisy_robust_mu_ellipsoid,
"#5" = MVP_noisy_robust_mu_ellipsoid,
"#6" = MVP_noisy_robust_mu_ellipsoid),
list(list("prices" = stock_prices)), price_name = "prices",
lookback = 6*21, optimize_every = 10000, rebalance_every = 10000)Backtesting 6 portfolios over 1 datasets (periodicity = daily data)...Code
listPortfoliosWithFailures(bt)
data.frame(
"stocks" = names(stock_prices),
"#1" = as.numeric(bt$`#1`$data1$w_optimized[1, ]),
"#2" = as.numeric(bt$`#2`$data1$w_optimized[1, ]),
"#3" = as.numeric(bt$`#3`$data1$w_optimized[1, ]),
"#4" = as.numeric(bt$`#4`$data1$w_optimized[1, ]),
"#5" = as.numeric(bt$`#5`$data1$w_optimized[1, ]),
"#6" = as.numeric(bt$`#6`$data1$w_optimized[1, ]),
check.names = FALSE) |>
melt(id.vars = "stocks") |>
ggplot(aes(x = stocks, y = value, fill = variable)) +
geom_bar(stat = "identity", position = "dodge", color = "black", width = 0.8) +
theme(legend.title = element_text()) + labs(fill = "realization") +
labs(title = "Portfolio for different estimation errors", y = "weight")