# 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")
Naive versus robust portfolios
We now compute the empirical distribution of the achieved mean–variance objective, as well as the Sharpe ratio, calculated over 1,000 Monte Carlo noisy observations. We can see that the robust designs avoid the worst-case realizations (left tail in the distributions) at the expense of not achieving the best-case realizations (right tail); thus, they are more stable and robust as expected.
Code
library(pob)
library(mvtnorm)
library(ggridges)
MVP_noisy <- function(dataset, ...) {
X <- diff(log(dataset$prices))[-1]
X_noisy <- rmvnorm(n = nrow(X), mean = colMeans(X), sigma = cov(X))
design_MVP(mu = colMeans(X_noisy), Sigma = cov(X_noisy), lambda = 1)
}
MVP_noisy_robust_mu_box <- function(dataset, ...) {
X <- diff(log(dataset$prices))[-1]
X_noisy <- rmvnorm(n = nrow(X), mean = colMeans(X), sigma = cov(X))
design_MVP_robust_mu_box(mu_hat = colMeans(X_noisy), delta = 0.2*sqrt(diag(cov(X_noisy))),
Sigma_hat = cov(X_noisy),
lambda = 1)
}
MVP_noisy_robust_mu_ellipsoid <- function(dataset, ...) {
X <- diff(log(dataset$prices))[-1]
X_noisy <- rmvnorm(n = nrow(X), mean = colMeans(X), sigma = cov(X))
design_MVP_robust_mu_ellipsoid(mu_hat = colMeans(X_noisy), kappa_mu = 0.2,
Sigma_hat = cov(X_noisy),
lambda = 1)
}
# parameters
num_realiz_histogram <- 1000
N <- 50
lookback_months <- 12
set.seed(42)
stock_prices <- SP500_2015to2020$stocks["2019::", sample(471, N)]
# true mu and Sigma in the first training set
X_ <- diff(log(stock_prices[1:(lookback_months*21), ]))[-1]
mu_true <- colMeans(X_)
Sigma_true <- cov(X_)
w_clairvoyant <- design_MVP(mu = mu_true, Sigma = Sigma_true, lambda = 1)
# backtests for naive MVP
set.seed(42)
bt_naive <- portfolioBacktest(rep(list(MVP_noisy), num_realiz_histogram),
list(list("prices" = stock_prices)), price_name = "prices",
lookback = lookback_months*21, optimize_every = 10000, rebalance_every = 10000,
paral_portfolios = 6, show_progress_bar = TRUE,
return_returns = FALSE)Backtesting 1000 portfolios over 1 datasets (periodicity = daily data)...Code
listPortfoliosWithFailures(bt_naive)
# backtests for robust MVP under box
set.seed(42)
bt_robust_box <- portfolioBacktest(rep(list(MVP_noisy_robust_mu_box), num_realiz_histogram),
list(list("prices" = stock_prices)), price_name = "prices",
lookback = lookback_months*21, optimize_every = 10000, rebalance_every = 10000,
paral_portfolios = 6, show_progress_bar = TRUE,
return_returns = FALSE)Backtesting 1000 portfolios over 1 datasets (periodicity = daily data)...Code
listPortfoliosWithFailures(bt_robust_box)
# backtests for robust MVP under ellipsoid
set.seed(42)
bt_robust_ellipsoid <- portfolioBacktest(rep(list(MVP_noisy_robust_mu_ellipsoid), num_realiz_histogram),
list(list("prices" = stock_prices)), price_name = "prices",
lookback = lookback_months*21, optimize_every = 10000, rebalance_every = 10000,
paral_portfolios = 6, show_progress_bar = TRUE,
return_returns = FALSE)Backtesting 1000 portfolios over 1 datasets (periodicity = daily data)...Code
listPortfoliosWithFailures(bt_robust_ellipsoid)
# evaluation
fun_mean_variance <- function(w) as.numeric(t(w) %*% mu_true - (1/2) * t(w) %*% Sigma_true %*% w)
fun_Sharpe_ratio <- function(w) as.numeric(t(w) %*% mu_true / sqrt(t(w) %*% Sigma_true %*% w))
mean_variance_clairvoyant <- fun_mean_variance(w_clairvoyant)
mean_variance_naive <- sapply(bt_naive, function(bt_realiz) fun_mean_variance(as.vector(bt_realiz$data1$w_optimized[1, ])))
mean_variance_robust_ellipsoid <- sapply(bt_robust_ellipsoid, function(bt_realiz) fun_mean_variance(as.vector(bt_realiz$data1$w_optimized[1, ])))
mean_variance_robust_box <- sapply(bt_robust_box, function(bt_realiz) fun_mean_variance(as.vector(bt_realiz$data1$w_optimized[1, ])))
SR_clairvoyant <- fun_Sharpe_ratio(w_clairvoyant)
SR_naive <- sapply(bt_naive, function(bt_realiz) fun_Sharpe_ratio(as.vector(bt_realiz$data1$w_optimized[1, ])))
SR_robust_ellipsoid <- sapply(bt_robust_ellipsoid, function(bt_realiz) fun_Sharpe_ratio(as.vector(bt_realiz$data1$w_optimized[1, ])))
SR_robust_box <- sapply(bt_robust_box, function(bt_realiz) fun_Sharpe_ratio(as.vector(bt_realiz$data1$w_optimized[1, ])))
df <- data.frame("method" = "Naive MVP",
"realization" = 1:num_realiz_histogram,
"mean-variance" = mean_variance_naive,
"Sharpe ratio" = SR_naive,
check.names = FALSE)
df <- rbind(df, data.frame("method" = "Robust MVP (box uncertainty set)",
"realization" = 1:num_realiz_histogram,
"mean-variance" = mean_variance_robust_box,
"Sharpe ratio" = SR_robust_box,
check.names = FALSE))
df <- rbind(df, data.frame("method" = "Robust MVP (ellipsoidal uncertainty set)",
"realization" = 1:num_realiz_histogram,
"mean-variance" = mean_variance_robust_ellipsoid,
"Sharpe ratio" = SR_robust_ellipsoid,
check.names = FALSE))
df$method <- factor(df$method, levels = unique(df$method))Code
ggplot(df, aes(x = `mean-variance`, y = method, fill = method)) +
geom_density_ridges(show.legend = FALSE) +
scale_y_discrete(limits = rev(levels(df$method)), expand = c(0, 0.2, 0, 1.6)) +
theme(axis.text.y = element_text(colour = "black", size = 10)) +
xlim(c(-0.001, NA)) +
labs(title = "Empirical distribution of mean-variance value", y = NULL)
Code
ggplot(df, aes(x = `Sharpe ratio`, y = method, fill = method)) +
geom_density_ridges(show.legend = FALSE) +
scale_y_discrete(limits = rev(levels(df$method)), expand = c(0, 0.2, 0, 1.6)) +
theme(axis.text.y = element_text(colour = "black", size = 10)) +
xlim(c(-0.05, NA)) +
labs(title = "Empirical distribution of Sharpe ratio", y = NULL)
Backtests
Backtest of naive versus robust mean–variance portfolios:
Code
library(pob)
MVP_naive <- function(dataset, ...) {
X <- diff(log(dataset$prices))[-1]
design_MVP(mu = colMeans(X), Sigma = cov(X), lambda = 1)
}
MVP_robust_mu_box <- function(dataset, ...) {
X <- diff(log(dataset$prices))[-1]
Sigma <- cov(X)
# robust design
design_MVP_robust_mu_box(mu_hat = colMeans(X), delta = 0.2*sqrt(diag(Sigma)),
Sigma_hat = Sigma,
lambda = 1)
}
MVP_robust_mu_ellipsoid <- function(dataset, ...) {
X <- diff(log(dataset$prices))[-1]
design_MVP_robust_mu_ellipsoid(mu_hat = colMeans(X), kappa_mu = 0.2,
Sigma_hat = cov(X),
lambda = 1)
}
# single backtest
N <- 50
lookback_months <- 6
set.seed(42)
stock_prices <- SP500_2015to2020$stocks["2017::", sample(471, N)]
bt <- portfolioBacktest(list("Naive MVP" = MVP_naive,
"Robust MVP (box uncertainty set)" = MVP_robust_mu_box,
"Robust MVP (ellipsoidal uncertainty set)" = MVP_robust_mu_ellipsoid),
list(list("prices" = stock_prices)), price_name = "prices",
lookback = lookback_months*21, optimize_every = 21, rebalance_every = 1,
return_portfolio = FALSE)Backtesting 3 portfolios over 1 datasets (periodicity = daily data)...Code
listPortfoliosWithFailures(bt)
bt <- lapply(bt, function(pf) {pf$data1$X_lin <- NULL; return(pf)}) # reduce object size Code
backtestChartCumReturn(bt) +
scale_x_date(date_breaks = "6 months", date_labels = "%b %Y", date_minor_breaks = "1 week") +
theme(legend.title = element_blank()) +
ggtitle("Cumulative P&L")
Code
backtestChartDrawdown(bt) +
scale_x_date(date_breaks = "6 months", date_labels = "%b %Y", date_minor_breaks = "1 week") +
theme(legend.title = element_blank()) +
ggtitle("Drawdown")
To be more exhaustive, we now run multiple backtests of naive versus robust mean–variance portfolios:
Code
library(pob)
# multiple backtests
N <- 50
lookback_months <- 6
set.seed(42)
stock_prices_resampled <- financialDataResample(list("prices" = SP500_2015to2020$stocks),
num_datasets = 200, N_sample = N, T_sample = 12*21)200 datasets resampled (with N = 50 instruments and length T = 252) from the original data between 2015-01-05 and 2020-09-22.Code
bt <- portfolioBacktest(list("Naive MVP" = MVP_naive,
"Robust MVP (box uncertainty set)" = MVP_robust_mu_box,
"Robust MVP (ellipsoidal uncertainty set)" = MVP_robust_mu_ellipsoid),
dataset_list = stock_prices_resampled,
lookback = lookback_months*21, optimize_every = 21, rebalance_every = 1,
return_portfolio = FALSE, return_returns = FALSE,
paral_datasets = 6, show_progress_bar = TRUE)Backtesting 3 portfolios over 200 datasets (periodicity = daily data)... Backtesting function "Naive MVP " (1/3) Backtesting function "Robust MVP (box uncertainty set)" (2/3) Backtesting function "Robust MVP (ellipsoidal uncertainty set)" (3/3)Code
listPortfoliosWithFailures(bt)Code
backtestBoxPlot(bt, "Sharpe ratio")
Code
backtestBoxPlot(bt, "max drawdown")
Portfolio resampling
Resampling methods
Estimating a parameter \(\bm{\theta}\) with \(\hat{\bm{\theta}}\) is insufficient without knowing the estimate’s accuracy. Confidence intervals traditionally provide this accuracy but rely heavily on theoretical mathematics. Resampling methods, like cross-validation and the bootstrap, use computer-based techniques to assess statistical accuracy without complex formulas (Efron and Tibshirani, 1993).
The bootstrap, introduced by Efron in 1979, mimics the original sampling process by resampling the observed data with replacement to create multiple bootstrap samples. Each sample leads to different realizations of the estimate, \(\hat{\bm{\theta}}^{*(b)}\), from which accuracy measures (bias, variance, confidence intervals) are empirically derived. This method allows for the assessment of estimator accuracy without needing additional data, with the statistical behavior of \(\hat{\bm{\theta}}^{*(b)}\) faithfully representing that of \(\hat{\bm{\theta}}\) compared to the true parameter \(\bm{\theta}\), especially as the number of resamples \(B\) increases (Efron and Tibshirani, 1993).
To be more precise, the original \(n\) observations \(\bm{x}_1,\dots,\bm{x}_n\) are sampled \(n\) times with replacement (some samples will be selected multiple times while others will not be used). This procedure is repeated \(B\) times to obtain the bootstrap samples, \[ \big(\bm{x}_1,\dots,\bm{x}_n\big) \rightarrow \left(\bm{x}_1^{*(b)},\dots,\bm{x}_n^{*(b)}\right), \quad b=1,\dots,B, \] each of which leads to a different realization of the estimation (bootstrap replicates), \[ \hat{\bm{\theta}}^{*(b)} = f\left(\bm{x}_1^{*(b)},\dots,\bm{x}_n^{*(b)}\right), \quad b=1,\dots,B, \] from which measures of accuracy of the estimator (bias, variance, confidence intervals, etc.) can then be empirically obtained.
The key theoretical result is that the statistical behavior of the random resampled estimates \(\hat{\bm{\theta}}^{*(b)}\) compared to \(\hat{\bm{\theta}}\) (taken as the true parameter) faithfully represent the statistics of the random estimates \(\hat{\bm{\theta}}\) compared to the true (unknown) parameter \(\bm{\theta}\). More exactly, the estimations of accuracy are asymptotically consistent as \(B\rightarrow\infty\) (under some technical conditions) (Efron and Tibshirani, 1993). This is a rather surprising result that allows the empirical assessment of the accuracy of the estimator without having access to the true parameter \(\bm{\theta}\).
We now illustrate the magic of the bootstrap to estimate the accuracy of the sample mean estimator (from \(n=100\) observations). In this case, the empirical distribution of the bias of the estimator is computed via \(B=1,000\) bootstraps, producing an accurate histogram compared to the true distribution. In practice, confidence intervals may suffice to assess the accuracy and fewer bootstraps may be used (even less bootstraps are necessary to compute the standard deviation of the bias).
Code
library(ggplot2)
B <- 1000 # number of boostraps
n <- 100 # number of observations for the estimation
mu_true <- 0
sd_true <- 1
set.seed(42)
# nominal sample
x_nominal <- rnorm(n = n, mean = mu_true, sd = sd_true)
mu_nominal <- mean(x_nominal)
sd_nominal <- sd(x_nominal)
df <- data.frame()
# nonparametric bootstrap
mu_hat_resampled <- c()
for (i in 1:B) {
x_resampled <- sample(x_nominal, size = n, replace = TRUE)
mu_hat_resampled <- c(mu_hat_resampled, mean(x_resampled))
}
df <- rbind(df, data.frame("bias" = mu_hat_resampled - mu_nominal,
"distribution" = "nonparametric bootstrap"))
# plot
ggplot(df[df$distribution != "true", ], aes(x = bias, fill = distribution)) +
geom_histogram(aes(y = ..density..), position = "identity", alpha = 0.7, color = "gray31", show.legend = FALSE) +
stat_function(fun = dnorm, args = list(mean = mu_true, sd = sd_true/sqrt(n)), color = "black", linewidth = 1.2) +
annotate("text", x = 0.13, y = 3.58, label = "true distribution") +
geom_segment(aes(x = 0.1, y = 3.45, xend = 0.088, yend = 2.8)) +
xlim(-0.4, 0.4) +
theme(axis.text.y = element_blank()) +
theme(legend.title = element_blank()) +
ggtitle("Empirical distribution of the mean estimation bias")
Portfolio resampling
In portfolio design, an optimization problem is formulated based on \(T\) observations of the assets’ returns \(\bm{x}_1,\dots,\bm{x}_T\) whose solution is supposedly an optimal portfolio \(\w\). As previous discussed, this solution is very sensitive to the inherent noise in the observed data or in the noise in the estimated parameters used in the portfolio formulation, such as the mean vector \(\hat{\bmu}\) and covariance matrix \(\hat{\bSigma}\).
Fortunately, we can capitalize on the results from the past half century in statistics; in particular, we can use resampling techniques, such as the bootstrap and bagging, to improve the portfolio design.
Portfolio bagging: The technique of aggregating portfolios was considered in 1998 via a bagging procedure (Michaud and Michaud, 1998, 2007, 2008; Scherer, 2002):
- Resample the original data \((\bm{x}_1,\dots,\bm{x}_T)\) \(B\) times via the bootstrap method and estimate \(B\) different versions of the mean vector and covariance matrix: \(\hat{\bmu}^{*(b)}\) and \(\hat{\bSigma}^{*(b)}\) for \(b=1,\dots,B\);
- Solve the optimal portfolio \(\w^{*(b)}\) for each bootstrap sample; and
- Average the portfolios via bagging: \[ \w^{\textm{bag}} = \frac{1}{B}\sum_{b=1}^B \w^{*(b)}. \]
Numerical experiments
Sensitivity
First, let’s write the code for some resampled portfolios, namely, bagged portfolio, subsample resampled (SSR) portfolio, and their combination SSR bagged:
Code
bootstrap_portfolio <- function(portfolio_fun, X,
num_bootstraps = 20, same_seed = TRUE, SSR_exp = 0.7,
type_bootstrap = c("vanilla", "subset resampling", "vanilla bootstrap + subset resampling"),
type_aggregation = c("mean", "gmedian"), force_sum_abs_w_1 = FALSE, cardinality = NULL
) {
type_bootstrap <- match.arg(type_bootstrap)
X <- as.matrix(X)
if (is.null(colnames(X)))
colnames(X) <- paste0("asset #", 1:ncol(X))
# generate bootstrap samples
X_sample_list <- switch(type_bootstrap,
"vanilla" = {
replicate(num_bootstraps, X[sample(nrow(X), replace = TRUE), ], simplify = FALSE)
},
"subset resampling" = {
replicate(num_bootstraps, X[, sample(ncol(X), ceiling(ncol(X)^SSR_exp))], simplify = FALSE)
},
"vanilla bootstrap + subset resampling" = {
replicate(num_bootstraps, X[sample(nrow(X), replace = TRUE), sample(ncol(X), ceiling(ncol(X)^SSR_exp))], simplify = FALSE)
},
stop("Bootstrap method unknown")
)
# design portfolios
w_list <- lapply(X_sample_list, portfolio_fun)
# include zeros in case of subspace resampling
if (grepl("subset resampling", type_bootstrap)) {
w_list <- lapply(w_list, function(w) {
if (is.null(names(w)))
stop("w does not have names in subset resampling")
w_full <- rep(0, ncol(X))
names(w_full) <- colnames(X)
w_full[names(w)] <- w
return(w_full)
})
}
# aggregate portfolios
w_matrix <- do.call(rbind, w_list)
w <- colMeans(w_matrix)
return(w)
}
design_MVP_bagged <- function(X, B, lambda = 1) {
bootstrap_portfolio(
portfolio_fun = function(X_sample) design_MVP(mu = colMeans(X_sample), Sigma = cov(X_sample), lambda = lambda),
X = X,
num_bootstraps = B,
type_bootstrap = "vanilla",
type_aggregation = "mean",
cardinality = ncol(X)
)
}
design_MVP_SSR <- function(X, B, lambda = 1) {
bootstrap_portfolio(
portfolio_fun = function(X_sample) design_MVP(mu = colMeans(X_sample), Sigma = cov(X_sample), lambda = lambda),
X = X,
num_bootstraps = B,
type_bootstrap = "subset resampling",
type_aggregation = "mean",
cardinality = ncol(X)
)
}
design_MVP_SSR_bagged <- function(X, B, lambda = 1) {
bootstrap_portfolio(
portfolio_fun = function(X_sample) design_MVP(mu = colMeans(X_sample), Sigma = cov(X_sample), lambda = lambda),
X = X,
num_bootstraps = B,
type_bootstrap = "vanilla bootstrap + subset resampling",
type_aggregation = "mean",
cardinality = ncol(X)
)
}We can now observe the improved sensitivity of the bagged mean–variance portfolios:
Code
library(pob)
library(mvtnorm)
stock_prices <- SP500_2015to2020$stocks["2019::", c("AAPL", "AMZN", "AMD", "GM", "GOOGL", "MGM", "MSFT", "QCOM", "TSCO", "UPS")]
MVP_bagged_noisy <- function(dataset, ...) {
X <- diff(log(dataset$prices))[-1]
X_noisy <- rmvnorm(n = nrow(X), mean = colMeans(X), sigma = cov(X))
design_MVP_bagged(X_noisy, B = 20, lambda = 1)
}
# backtests
set.seed(42)
bt <- portfolioBacktest(list("#1" = MVP_bagged_noisy,
"#2" = MVP_bagged_noisy,
"#3" = MVP_bagged_noisy,
"#4" = MVP_bagged_noisy,
"#5" = MVP_bagged_noisy,
"#6" = MVP_bagged_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)
# plot
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")
Naive versus resampled portfolios
We now compute the empirical distribution of the achieved mean–variance objective, as well as the Sharpe ratio, calculated over 1,000 Monte Carlo noisy observations. We can observe that resampled portfolios are more stable and do not suffer from extreme bad realizations (unlike the naive portfolio). However, the naive portfolio can be superior for some realizations.
Code
library(pob)
library(mvtnorm)
library(ggridges)
MVP_noisy <- function(dataset, ...) {
X <- diff(log(dataset$prices))[-1]
X_noisy <- rmvnorm(n = nrow(X), mean = colMeans(X), sigma = cov(X))
design_MVP(mu = colMeans(X_noisy), Sigma = cov(X_noisy), lambda = 1)
}
MVP_bagged_noisy <- function(dataset, ...) {
X <- diff(log(dataset$prices))[-1]
X_noisy <- rmvnorm(n = nrow(X), mean = colMeans(X), sigma = cov(X))
design_MVP_bagged(X_noisy, B = 20, lambda = 1)
}
MVP_SSR_noisy <- function(dataset, ...) {
X <- diff(log(dataset$prices))[-1]
X_noisy <- rmvnorm(n = nrow(X), mean = colMeans(X), sigma = cov(X))
design_MVP_SSR(X_noisy, B = 20, lambda = 1)
}
MVP_SSR_bagged_noisy <- function(dataset, ...) {
X <- diff(log(dataset$prices))[-1]
X_noisy <- rmvnorm(n = nrow(X), mean = colMeans(X), sigma = cov(X))
design_MVP_SSR_bagged(X_noisy, B = 20, lambda = 1)
}
# parameters
num_realiz_histogram <- 1000
N <- 50
lookback_months <- 12
set.seed(42)
stock_prices <- SP500_2015to2020$stocks["2019::", sample(471, N)]
# true mu and Sigma in the first training set
X_ <- diff(log(stock_prices[1:(lookback_months*21), ]))[-1]
mu_true <- colMeans(X_)
Sigma_true <- cov(X_)
w_clairvoyant <- design_MVP(mu = mu_true, Sigma = Sigma_true, lambda = 1)
# backtests for naive MVP
set.seed(42)
bt_naive <- portfolioBacktest(rep(list(MVP_noisy), num_realiz_histogram),
list(list("prices" = stock_prices)), price_name = "prices",
lookback = lookback_months*21, optimize_every = 10000, rebalance_every = 10000,
paral_portfolios = 6, show_progress_bar = TRUE,
return_returns = FALSE)Backtesting 1000 portfolios over 1 datasets (periodicity = daily data)...Code
listPortfoliosWithFailures(bt_naive)
# backtests for bagged MVP
set.seed(42)
bt_bagged <- portfolioBacktest(rep(list(MVP_bagged_noisy), num_realiz_histogram),
list(list("prices" = stock_prices)), price_name = "prices",
lookback = lookback_months*21, optimize_every = 10000, rebalance_every = 10000,
paral_portfolios = 6, show_progress_bar = TRUE,
return_returns = FALSE)Backtesting 1000 portfolios over 1 datasets (periodicity = daily data)...Code
listPortfoliosWithFailures(bt_bagged)
# backtests for SSR MVP
set.seed(42)
bt_SSR <- portfolioBacktest(rep(list(MVP_SSR_noisy), num_realiz_histogram),
list(list("prices" = stock_prices)), price_name = "prices",
lookback = lookback_months*21, optimize_every = 10000, rebalance_every = 10000,
paral_portfolios = 6, show_progress_bar = TRUE,
return_returns = FALSE)Backtesting 1000 portfolios over 1 datasets (periodicity = daily data)...Code
listPortfoliosWithFailures(bt_SSR)
# backtests for SSR bagged MVP
set.seed(42)
bt_SSR_bagged <- portfolioBacktest(rep(list(MVP_SSR_bagged_noisy), num_realiz_histogram),
list(list("prices" = stock_prices)), price_name = "prices",
lookback = lookback_months*21, optimize_every = 10000, rebalance_every = 10000,
paral_portfolios = 6, show_progress_bar = TRUE,
return_returns = FALSE)Backtesting 1000 portfolios over 1 datasets (periodicity = daily data)...Code
listPortfoliosWithFailures(bt_SSR_bagged)
# evaluation
fun_mean_variance <- function(w) as.numeric(t(w) %*% mu_true - (1/2) * t(w) %*% Sigma_true %*% w)
fun_Sharpe_ratio <- function(w) as.numeric(t(w) %*% mu_true / sqrt(t(w) %*% Sigma_true %*% w))
mean_variance_clairvoyant <- fun_mean_variance(w_clairvoyant)
mean_variance_naive <- sapply(bt_naive, function(bt_realiz) fun_mean_variance(as.vector(bt_realiz$data1$w_optimized[1, ])))
mean_variance_bagged <- sapply(bt_bagged, function(bt_realiz) fun_mean_variance(as.vector(bt_realiz$data1$w_optimized[1, ])))
mean_variance_SSR <- sapply(bt_SSR, function(bt_realiz) fun_mean_variance(as.vector(bt_realiz$data1$w_optimized[1, ])))
mean_variance_SSR_bagged <- sapply(bt_SSR_bagged, function(bt_realiz) fun_mean_variance(as.vector(bt_realiz$data1$w_optimized[1, ])))
SR_clairvoyant <- fun_Sharpe_ratio(w_clairvoyant)
SR_naive <- sapply(bt_naive, function(bt_realiz) fun_Sharpe_ratio(as.vector(bt_realiz$data1$w_optimized[1, ])))
SR_bagged <- sapply(bt_bagged, function(bt_realiz) fun_Sharpe_ratio(as.vector(bt_realiz$data1$w_optimized[1, ])))
SR_SSR <- sapply(bt_SSR, function(bt_realiz) fun_Sharpe_ratio(as.vector(bt_realiz$data1$w_optimized[1, ])))
SR_SSR_bagged <- sapply(bt_SSR_bagged, function(bt_realiz) fun_Sharpe_ratio(as.vector(bt_realiz$data1$w_optimized[1, ])))
df <- data.frame("method" = "Naive MVP",
"realization" = 1:num_realiz_histogram,
"mean-variance" = mean_variance_naive,
"Sharpe ratio" = SR_naive,
check.names = FALSE)
df <- rbind(df, data.frame("method" = "Bagged MVP",
"realization" = 1:num_realiz_histogram,
"mean-variance" = mean_variance_bagged,
"Sharpe ratio" = SR_bagged,
check.names = FALSE))
df <- rbind(df, data.frame("method" = "Subset resampled MVP",
"realization" = 1:num_realiz_histogram,
"mean-variance" = mean_variance_SSR,
"Sharpe ratio" = SR_SSR,
check.names = FALSE))
df <- rbind(df, data.frame("method" = "Subset bagged MVP",
"realization" = 1:num_realiz_histogram,
"mean-variance" = mean_variance_SSR_bagged,
"Sharpe ratio" = SR_SSR_bagged,
check.names = FALSE))
df$method <- factor(df$method, levels = unique(df$method))Code
ggplot(df, aes(x = `mean-variance`, y = method, fill = method)) +
geom_density_ridges(show.legend = FALSE) +
scale_y_discrete(limits = rev(levels(df$method)), expand = c(0, 0.2, 0, 1.6)) +
theme(axis.text.y = element_text(colour = "black", size = 10)) +
xlim(c(-0.001, NA)) +
labs(title = "Empirical distribution of mean-variance value", y = NULL)
Code
ggplot(df, aes(x = `Sharpe ratio`, y = method, fill = method)) +
geom_density_ridges(show.legend = FALSE) +
scale_y_discrete(limits = rev(levels(df$method)), expand = c(0, 0.2, 0, 1.6)) +
theme(axis.text.y = element_text(colour = "black", size = 10)) +
xlim(c(-0.05, NA)) +
labs(title = "Empirical distribution of Sharpe ratio", y = NULL)
Backtests
Backtest of naive versus resampled mean–variance portfolios:
Code
library(pob)
MVP_naive <- function(dataset, ...) {
X <- diff(log(dataset$prices))[-1]
design_MVP(mu = colMeans(X), Sigma = cov(X), lambda = 1)
}
MVP_bagged <- function(dataset, ...) {
X <- diff(log(dataset$prices))[-1]
design_MVP_bagged(X, B = 20, lambda = 1)
}
MVP_SSR <- function(dataset, ...) {
X <- diff(log(dataset$prices))[-1]
design_MVP_SSR(X, B = 20, lambda = 1)
}
MVP_SSR_bagged <- function(dataset, ...) {
X <- diff(log(dataset$prices))[-1]
design_MVP_SSR_bagged(X, B = 20, lambda = 1)
}
# single backtest
N <- 50
lookback_months <- 6
set.seed(42)
stock_prices <- SP500_2015to2020$stocks["2017::", sample(471, N)]
bt <- portfolioBacktest(list("Naive MVP" = MVP_naive,
"Bagged MVP" = MVP_bagged,
"Subset resampled MVP" = MVP_SSR,
"Subset bagged MVP" = MVP_SSR_bagged),
list(list("prices" = stock_prices)), price_name = "prices",
lookback = lookback_months*21, optimize_every = 21, rebalance_every = 1,
return_portfolio = FALSE)Backtesting 4 portfolios over 1 datasets (periodicity = daily data)...Code
listPortfoliosWithFailures(bt)
bt <- lapply(bt, function(pf) {pf$data1$X_lin <- NULL; return(pf)}) # reduce object sizeCode
backtestChartCumReturn(bt) +
scale_x_date(date_breaks = "6 months", date_labels = "%b %Y", date_minor_breaks = "1 week") +
theme(legend.title = element_blank()) +
ggtitle("Cumulative P&L")
Code
backtestChartDrawdown(bt) +
scale_x_date(date_breaks = "6 months", date_labels = "%b %Y", date_minor_breaks = "1 week") +
theme(legend.title = element_blank()) +
ggtitle("Drawdown")
To be more exhaustive, we now run multiple backtests of naive versus resampled mean–variance portfolios:
Code
library(pob)
# multiple backtests
N <- 50
lookback_months <- 6
set.seed(42)
stock_prices_resampled <- financialDataResample(list("prices" = SP500_2015to2020$stocks),
num_datasets = 200, N_sample = N, T_sample = 12*21)200 datasets resampled (with N = 50 instruments and length T = 252) from the original data between 2015-01-05 and 2020-09-22.Code
bt <- portfolioBacktest(list("Naive MVP" = MVP_naive,
"Bagged MVP" = MVP_bagged,
"Subset resampled MVP" = MVP_SSR,
"Subset bagged MVP" = MVP_SSR_bagged),
dataset_list = stock_prices_resampled,
lookback = lookback_months*21, optimize_every = 21, rebalance_every = 1,
return_portfolio = FALSE, return_returns = FALSE,
paral_datasets = 6, show_progress_bar = TRUE)Backtesting 4 portfolios over 200 datasets (periodicity = daily data)... Backtesting function "Naive MVP " (1/4) Backtesting function "Bagged MVP " (2/4) Backtesting function "Subset resampled MVP" (3/4) Backtesting function "Subset bagged MVP" (4/4)Code
listPortfoliosWithFailures(bt)Code
backtestBoxPlot(bt, "Sharpe ratio")
Code
backtestBoxPlot(bt, "max drawdown")
References
Efron, B., and Tibshirani, R. J. (1993). An introduction to the bootstrap. Springer.
Michaud, R. O., and Michaud, R. O. (1998). Efficient asset management: A practical guide to stock portfolio optimization and asset allocation. Harvard Business School Press.
Michaud, R. O., and Michaud, R. O. (2007). Estimation error and portfolio optimization: A resampling solution. Journal of Investment and Management, 6(1), 8–28.
Michaud, R. O., and Michaud, R. O. (2008). Efficient asset management: A practical guide to stock portfolio optimization and asset allocation. Oxford University Press.
Scherer, B. (2002). Portfolio resampling: Review and critique. Financial Analysts Journal, 58(6), 98–109.