# basic finance
library(xts) # to manipulate time series of stock data
library(quantmod) # to download stock data
library(PerformanceAnalytics) # to compute performance measures
library(portfolioBacktest) # to conduct backtests
library(pob) # book package with financial data
# plotting
library(ggplot2) # for nice plots
library(reshape2) # to reshape data
# optimization
library(CVXR)
library(nloptr)Portfolio Optimization
Chapter 7: Modern Portfolio Theory
R code
R code examples for Chapter 7 of the book:
Daniel P. Palomar (2024). Portfolio Optimization: Theory and Application. Cambridge University Press.
Loading packages
The following packages are used in the examples:
Mean-variance portfolio (MVP)
Vanilla MVP
We explore the mean-variance portfolio (MVP), \[ \begin{array}{ll} \underset{\w}{\textm{maximize}} & \w^\T\bmu - \frac{\lambda}{2}\w^\T\bSigma\w\\ \textm{subject to} & \bm{1}^\T\w=1, \quad \w\ge\bm{0}, \end{array} \] with different values of the hyper-parameter \(\lambda\).
data(SP500_2015to2020)
stock_prices <- SP500_2015to2020$stocks["2019::", c("AAPL", "AMZN", "AMD", "GM", "GOOGL", "MGM", "MSFT", "QCOM", "TSCO", "UPS")]
T <- nrow(stock_prices)
T_trn <- round(0.70*T)
#
# Define portfolios
#
EWP <- function(data, ...) {
N <- ncol(data$prices)
return(rep(1/N, N))
}
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))
return(w)
}
MVP_lmd1 <- function(dataset, ...) {
N <- ncol(dataset$prices)
X <- diff(log(dataset$prices))[-1]
mu <- colMeans(X)
Sigma <- cov(X)
w <- design_MVP(mu, Sigma, lambda = 1)
return(w)
}
MVP_lmd4 <- function(dataset, ...) {
N <- ncol(dataset$prices)
X <- diff(log(dataset$prices))[-1]
mu <- colMeans(X)
Sigma <- cov(X)
w <- design_MVP(mu, Sigma, lambda = 4)
return(w)
}
MVP_lmd16 <- function(dataset, ...) {
N <- ncol(dataset$prices)
X <- diff(log(dataset$prices))[-1]
mu <- colMeans(X)
Sigma <- cov(X)
w <- design_MVP(mu, Sigma, lambda = 16)
return(w)
}
MVP_lmd64 <- function(dataset, ...) {
N <- ncol(dataset$prices)
X <- diff(log(dataset$prices))[-1]
mu <- colMeans(X)
Sigma <- cov(X)
w <- design_MVP(mu, Sigma, lambda = 64)
return(w)
}
# single backtest
bt <- portfolioBacktest(portfolio_funs = list("1/N" = EWP,
"MVP (lmd = 1)" = MVP_lmd1,
"MVP (lmd = 4)" = MVP_lmd4,
"MVP (lmd = 16)" = MVP_lmd16,
"MVP (lmd = 64)" = MVP_lmd64),
dataset_list = list("dataset1" = list("prices" = stock_prices)),
lookback = T_trn, optimize_every = 30, rebalance_every = 1)Backtesting 5 portfolios over 1 datasets (periodicity = daily data)...# barplot
data.frame(
"stocks" = names(stock_prices),
"1/N" = as.numeric(bt$`1/N`$dataset1$w_optimized[1, ]),
"MVP (lmd = 1)" = as.numeric(bt$`MVP (lmd = 1)`$dataset1$w_optimized[1, ]),
"MVP (lmd = 4)" = as.numeric(bt$`MVP (lmd = 4)`$dataset1$w_optimized[1, ]),
"MVP (lmd = 16)" = as.numeric(bt$`MVP (lmd = 16)`$dataset1$w_optimized[1, ]),
"MVP (lmd = 64)" = as.numeric(bt$`MVP (lmd = 64)`$dataset1$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) +
labs(title = "Capital weight allocation", y = "weight", fill = "portfolios")
backtestChartCumReturn(bt) +
scale_x_date(date_breaks = "1 month", date_labels = "%b %Y", date_minor_breaks = "1 day") +
ggtitle("Cumulative return (out of sample)")
backtestChartDrawdown(bt) +
scale_x_date(date_breaks = "1 month", date_labels = "%b %Y", date_minor_breaks = "1 day") +
ggtitle("Drawdown (out of sample)")
summary_bt <- backtestSummary(bt)
summaryTable(summary_bt, type = "DT",
measures = c("Sharpe ratio", "annual return", "annual volatility", "max drawdown"))MVP with diversification heuristics
We now consider the MVP with two diversification heuristics:
- upper bound: \(\w\leq0.25\times\bm{1}\); and
- diversification constraint: \(\|\w\|_2^2 \leq 2/N\).
data(SP500_2015to2020)
stock_prices <- SP500_2015to2020$stocks["2019::", c("AAPL", "AMZN", "AMD", "GM", "GOOGL", "MGM", "MSFT", "QCOM", "TSCO", "UPS")]
T <- nrow(stock_prices)
T_trn <- round(0.70*T)
#
# Define portfolios
#
design_GMVP <- function(Sigma) {
w <- Variable(nrow(Sigma))
prob <- Problem(Minimize(quad_form(w, Sigma)),
constraints = list(w >= 0, sum(w) == 1))
result <- solve(prob)
w <- as.vector(result$getValue(w))
return(w)
}
GMVP <- function(dataset, ...) {
N <- ncol(dataset$prices)
X <- diff(log(dataset$prices))[-1]
Sigma <- cov(X)
w <- design_GMVP(Sigma)
return(w)
}
MVP <- function(dataset, ...) {
N <- ncol(dataset$prices)
X <- diff(log(dataset$prices))[-1]
mu <- colMeans(X)
Sigma <- cov(X)
w <- design_MVP(mu, Sigma)
return(w)
}
MVP_ub <- function(dataset, ...) {
N <- ncol(dataset$prices)
X <- diff(log(dataset$prices))[-1]
mu <- colMeans(X)
Sigma <- cov(X)
w <- design_MVP(mu, Sigma, ub = 0.25)
return(w)
}
design_MVP_diversification <- function(mu, Sigma, lmd = 1, D) {
w <- Variable(nrow(Sigma))
prob <- Problem(Maximize(t(mu) %*% w - (lmd/2)*quad_form(w, Sigma)),
constraints = list(w >= 0, sum(w) == 1, sum(w^2) <= D))
result <- solve(prob)
w <- as.vector(result$getValue(w))
return(w)
}
MVP_diversification <- function(dataset, ...) {
N <- ncol(dataset$prices)
X <- diff(log(dataset$prices))[-1]
mu <- colMeans(X)
Sigma <- cov(X)
w <- design_MVP_diversification(mu, Sigma, D = 2/N)
return(w)
}
# single backtest
bt <- portfolioBacktest(portfolio_funs = list("1/N" = EWP,
"GMVP" = GMVP,
"MVP" = MVP,
"MVP with upper bound" = MVP_ub,
"MVP with diversific. const." = MVP_diversification),
dataset_list = list("dataset1" = list("prices" = stock_prices)),
lookback = T_trn, optimize_every = 30, rebalance_every = 1)Backtesting 5 portfolios over 1 datasets (periodicity = daily data)...# barplot
data.frame(
"stocks" = names(stock_prices),
"1/N" = as.numeric(bt$`1/N`$dataset1$w_optimized[1, ]),
"GMVP" = as.numeric(bt$`GMVP`$dataset1$w_optimized[1, ]),
"MVP" = as.numeric(bt$`MVP`$dataset1$w_optimized[1, ]),
"MVP with upper bound" = as.numeric(bt$`MVP with upper bound`$dataset1$w_optimized[1, ]),
"MVP with diversific. const." = as.numeric(bt$`MVP with diversific. const.`$dataset1$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) +
labs(title = "Capital weight allocation", y = "weight", fill = "portfolios")
backtestChartCumReturn(bt) +
scale_x_date(date_breaks = "1 month", date_labels = "%b %Y", date_minor_breaks = "1 day") +
ggtitle("Cumulative return (out of sample)")
backtestChartDrawdown(bt) +
scale_x_date(date_breaks = "1 month", date_labels = "%b %Y", date_minor_breaks = "1 day") +
ggtitle("Drawdown (out of sample)")
summary_bt <- backtestSummary(bt)
summaryTable(summary_bt, type = "DT",
measures = c("Sharpe ratio", "annual return", "annual volatility", "max drawdown"))Maximum Sharpe ratio portfolio (MSRP)
Recall the maximum Sharpe ratio portfolio (MSRP) formulation: \[ \begin{array}{ll} \underset{\w}{\textm{maximize}} & \dfrac{\w^\T\bmu - r_\textm{f}}{\sqrt{\w^\T\bSigma\w}}\\ \textm{subject to} & \begin{array}{l} \bm{1}^\T\w=1, \quad \w\ge\bm{0}.\end{array} \end{array} \]
We will solve it numerically via different methods: - general-purpose nonlinear solver - bisection method - Dinkelbach method - Schaible transform
First, let’s obtain the mean vector \(\bmu\) and covariance matrix \(\bSigma\) from the market data:
data(SP500_2015to2020)
stock_prices <- SP500_2015to2020$stocks["2019::", c("AAPL", "AMZN", "AMD", "GM", "GOOGL", "MGM", "MSFT", "QCOM", "TSCO", "UPS")]
X <- diff(log(stock_prices))[-1]
N <- ncol(X)
mu <- colMeans(X)
Sigma <- cov(X)MSRP via general-purpose nonlinear solver
We will solve this problem with the general-purpose nonlinear solver nloptr in R:
library(nloptr)
# define the nonconvex objective function
fn_SR <- function(w) {
return(as.numeric(t(w) %*% mu / sqrt(t(w) %*% Sigma %*% w)))
}
neg_fn_SR <- function(w) return(-fn_SR(w))
# initial point
w0 <- rep(1/N, N)
res <- nloptr::slsqp(w0, neg_fn_SR,
lower = rep(0, N), # w >= 0
heq = function(w) return(sum(w) - 1)) # sum(w) = 1
w_nonlinear_solver <- res$par
res$par
[1] 5.610928e-01 1.563860e-01 1.884386e-01 4.271107e-17 0.000000e+00
[6] 6.552656e-17 8.326145e-17 0.000000e+00 7.884021e-02 1.524237e-02
$value
[1] -0.1071244
$iter
[1] 28
$convergence
[1] 4
$message
[1] "NLOPT_XTOL_REACHED: Optimization stopped because xtol_rel or xtol_abs (above) was reached."MSRP via bisection method
We are going to solve the nonconvex problem \[ \begin{array}{ll} \underset{\w,t}{\textm{maximize}} & t\\ \textm{subject to} & t \leq \dfrac{\w^\T\bmu}{\sqrt{\w^\T\bSigma\w}}\\ & \bm{1}^\T\w=1, \quad \w\geq\bm{0}. \end{array} \] via bisection on \(t\) with the following (convex) SOCP problem for a given \(t\): \[ \begin{array}{ll} \underset{\;}{\textm{find}} & \w\\ \textm{subject to} & t \left\Vert \bSigma^{1/2}\w\right\Vert_{2}\leq\w^\T\bmu\\ & \bm{1}^\T\w=1, \quad \w\geq\bm{0}. \end{array} \]
# define the inner solver based on an SOCP solver
# (we will simply use CVXR for convenience, see: https://cvxr.rbind.io/cvxr_functions/)
library(CVXR)
# square-root of matrix Sigma
Sigma_12 <- chol(Sigma)
max(abs(t(Sigma_12) %*% Sigma_12 - Sigma)) # sanity check[1] 1.084202e-19# create function for MVP
SOCP_bisection <- function(t) {
w <- Variable(nrow(Sigma))
prob <- Problem(Maximize(0),
constraints = list(t*cvxr_norm(Sigma_12 %*% w, 2) <= t(mu) %*% w,
sum(w) == 1,
w >= 0))
result <- solve(prob)
return(list("status" = result$status, "w" = as.vector(result$getValue(w))))
}
# now run the bisection algorithm
t_lb <- 0 # for sure the problem is feasible in this case
t_ub <- 10 # a tighter upper bound could be chosen (10 is a conservative choice)
while(t_ub - t_lb > 1e-6) {
t <- (t_ub + t_lb)/2 # midpoint
if(SOCP_bisection(t)$status == "infeasible")
t_ub <- t
else
t_lb <- t
}
w_bisection <- SOCP_bisection(t_lb)$w
# comparison between two solutions
round(cbind(w_nonlinear_solver, w_bisection), digits = 3) w_nonlinear_solver w_bisection
[1,] 0.561 0.560
[2,] 0.156 0.156
[3,] 0.188 0.188
[4,] 0.000 0.000
[5,] 0.000 0.000
[6,] 0.000 0.000
[7,] 0.000 0.000
[8,] 0.000 0.000
[9,] 0.079 0.079
[10,] 0.015 0.017# Sharpe ratio of two solutions
c("nonlinear_solver" = fn_SR(w_nonlinear_solver),
"bisection" = fn_SR(w_bisection))nonlinear_solver bisection
0.1071244 0.1071230 MSRP via Dinkelbach method
We are going to solve the nonconvex problem \[ \begin{array}{ll} \underset{\w}{\textm{maximize}} & \dfrac{\w^\T\bmu}{\sqrt{\w^\T\bSigma\w}}\\ \textm{subject to} & \begin{array}{l} \bm{1}^\T\w=1, \quad \w\ge\bm{0}.\end{array} \end{array} \] by iteratively solving the (convex) SOCP problem for a given \(y^{(k)}\): \[ \begin{array}{ll} \underset{\w}{\textm{maximize}} & \w^\T\bmu - y^{(k)} \left\Vert \bSigma^{1/2}\w\right\Vert_{2}\\ \textm{subject to} & \begin{array}{l} \bm{1}^\T\w=1, \quad \w\ge\bm{0}.\end{array} \end{array} \] where \[y^{(k)} = \frac{\mathbf{w}^{(k)\T}\bmu}{\sqrt{\w^{(k)\T}\bSigma\w^{(k)}}}.\]
# define the inner solver based on an SOCP solver
# (we will simply use CVXR for convenience, see: https://cvxr.rbind.io/cvxr_functions/)
library(CVXR)
# square-root of matrix Sigma
Sigma_12 <- chol(Sigma)
max(abs(t(Sigma_12) %*% Sigma_12 - Sigma)) # sanity check[1] 1.084202e-19# create function for MVP
SOCP_Dinkelbach <- function(y) {
w <- Variable(nrow(Sigma))
prob <- Problem(Maximize(t(mu) %*% w - y*cvxr_norm(Sigma_12 %*% w, 2)),
constraints = list(sum(w) == 1,
w >= 0))
result <- CVXR::solve(prob)
return(as.vector(result$getValue(w)))
}
# initial point (has to satisfy t(w_k) %*% mu>=0)
i_max <- which.max(mu)
w_k <- rep(0, N)
w_k[i_max] <- 1
# now the iterative Dinkelbach algorithm
k <- 1
while(k == 1 || max(abs(w_k - w_prev)) > 1e-6) {
w_prev <- w_k
y_k <- as.numeric(t(w_k) %*% mu / sqrt(t(w_k) %*% Sigma %*% w_k))
w_k <- SOCP_Dinkelbach(y_k)
k <- k + 1
}
w_Dinkelbach <- w_k
cat("Number of iterarions:", k-1)Number of iterarions: 5# comparison among three solutions
round(cbind(w_nonlinear_solver, w_bisection, w_Dinkelbach), digits = 3) w_nonlinear_solver w_bisection w_Dinkelbach
[1,] 0.561 0.560 0.561
[2,] 0.156 0.156 0.156
[3,] 0.188 0.188 0.188
[4,] 0.000 0.000 0.000
[5,] 0.000 0.000 0.000
[6,] 0.000 0.000 0.000
[7,] 0.000 0.000 0.000
[8,] 0.000 0.000 0.000
[9,] 0.079 0.079 0.079
[10,] 0.015 0.017 0.015# Sharpe ratio of three solutions
c("nonlinear_solver" = fn_SR(w_nonlinear_solver),
"bisection" = fn_SR(w_bisection),
"Dinkelbach" = fn_SR(w_Dinkelbach))nonlinear_solver bisection Dinkelbach
0.1071244 0.1071230 0.1071244 MSRP via Schaible transform
The maximum Sharpe ratio portfolio (MSRP) is the nonconvex problem \[ \begin{array}{ll} \underset{\w}{\textm{maximize}} & \dfrac{\w^\T\bmu}{\sqrt{\w^\T\bSigma\w}}\\ \textm{subject to} & \begin{array}{l} \bm{1}^\T\w=1, \quad \w\ge\bm{0}.\end{array} \end{array} \] that can be rewritten in convex form as \[ \begin{array}{ll} \underset{\w}{\textm{minimize}} & \tilde{\w}^\T\bSigma\tilde{\w}\\ \textm{subject to} & \tilde{\w}^\T\bmu = 1\\ & \tilde{\w}\ge\bm{0} \end{array} \] and then \(\w = \tilde{\w}/(\bm{1}^\T\tilde{\w})\).
This is a quadratic problem (QP) and we can conveniently use CVXR (although one is advised to use a specific QP solver like quadprog for speed and stability):
# create function for MSRP
MSRP <- 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)
w <- as.vector(result$getValue(w_)/sum(result$getValue(w_)))
names(w) <- colnames(Sigma)
return(w)
}
# this function can now be used as
w_MSRP <- MSRP(mu, Sigma)
# comparison among solutions
round(cbind(w_nonlinear_solver, w_bisection, w_Dinkelbach, w_MSRP), digits = 3) w_nonlinear_solver w_bisection w_Dinkelbach w_MSRP
AAPL 0.561 0.560 0.561 0.561
AMZN 0.156 0.156 0.156 0.156
AMD 0.188 0.188 0.188 0.188
GM 0.000 0.000 0.000 0.000
GOOGL 0.000 0.000 0.000 0.000
MGM 0.000 0.000 0.000 0.000
MSFT 0.000 0.000 0.000 0.000
QCOM 0.000 0.000 0.000 0.000
TSCO 0.079 0.079 0.079 0.079
UPS 0.015 0.017 0.015 0.015# Sharpe ratio of different solutions
c("nonlinear_solver" = fn_SR(w_nonlinear_solver),
"bisection" = fn_SR(w_bisection),
"Dinkelbach" = fn_SR(w_Dinkelbach),
"Schaible" = fn_SR(w_MSRP))nonlinear_solver bisection Dinkelbach Schaible
0.1071244 0.1071230 0.1071244 0.1071244 All methods produce the optimal solution.
Numerical experiments
We now compare the following portfolios:
global minimum variance portfolio (GMVP): \[ \begin{array}{ll} \underset{\w}{\textm{minimize}} & \w^\T\bSigma\w\\ \textm{subject to} & \bm{1}^\T\w=1, \quad \w\ge\bm{0}; \end{array} \]
mean-variance portfolio (MVP): \[ \begin{array}{ll} \underset{\w}{\textm{maximize}} & \w^\T\bmu - \frac{\lambda}{2}\w^\T\bSigma\w\\ \textm{subject to} & \bm{1}^\T\w=1, \quad \w\ge\bm{0}; \end{array} \]
maximum Sharpe ratio portfolio (MSRP): \[ \begin{array}{ll} \underset{\w}{\textm{maximize}} & \dfrac{\w^\T\bmu - r_\textm{f}}{\sqrt{\w^\T\bSigma\w}}\\ \textm{subject to} & \begin{array}{l} \bm{1}^\T\w=1, \quad \w\ge\bm{0}.\end{array} \end{array} \]
data(SP500_2015to2020)
stock_prices <- SP500_2015to2020$stocks["2019::", c("AAPL", "AMZN", "AMD", "GM", "GOOGL", "MGM", "MSFT", "QCOM", "TSCO", "UPS")]
T <- nrow(stock_prices)
T_trn <- round(0.70*T)
#
# Define portfolios
#
GMVP <- function(dataset, ...) {
N <- ncol(dataset$prices)
X <- diff(log(dataset$prices))[-1]
Sigma <- cov(X)
w <- design_GMVP(Sigma)
return(w)
}
MVP <- function(dataset, ...) {
N <- ncol(dataset$prices)
X <- diff(log(dataset$prices))[-1]
mu <- colMeans(X)
Sigma <- cov(X)
w <- design_MVP(mu, Sigma, lambda = 4)
return(w)
}
design_MSRP <- 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)
w <- as.vector(result$getValue(w_)/sum(result$getValue(w_)))
names(w) <- colnames(Sigma)
return(w)
}
MSRP <- function(dataset, ...) {
N <- ncol(dataset$prices)
X <- diff(log(dataset$prices))[-1]
mu <- colMeans(X)
Sigma <- cov(X)
w <- design_MSRP(mu, Sigma)
return(w)
}
# single backtest
bt <- portfolioBacktest(portfolio_funs = list("GMVP" = GMVP,
"MVP" = MVP,
"MSRP" = MSRP),
dataset_list = list("dataset1" = list("prices" = stock_prices)),
lookback = T_trn, optimize_every = 30, rebalance_every = 1)Backtesting 3 portfolios over 1 datasets (periodicity = daily data)...# barplot
data.frame(
"stocks" = names(stock_prices),
"GMVP" = as.numeric(bt$`GMVP`$dataset1$w_optimized[1, ]),
"MVP" = as.numeric(bt$`MVP`$dataset1$w_optimized[1, ]),
"MSRP" = as.numeric(bt$`MSRP`$dataset1$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) +
labs(title = "Capital weight allocation", y = "weight", fill = "portfolios")
backtestChartCumReturn(bt) +
scale_x_date(date_breaks = "1 month", date_labels = "%b %Y", date_minor_breaks = "1 day") +
ggtitle("Cumulative return (out of sample)")
backtestChartDrawdown(bt) +
scale_x_date(date_breaks = "1 month", date_labels = "%b %Y", date_minor_breaks = "1 day") +
ggtitle("Drawdown (out of sample)")
summary_bt <- backtestSummary(bt)
summaryTable(summary_bt, type = "DT",
measures = c("Sharpe ratio", "annual return", "annual volatility", "max drawdown"))Universal algorithm for portfolio optimization
MSRP
We consider two methods for the resolution of the MSRP:
Via the Schaible transform, i.e., solving \[ \begin{array}{ll} \underset{\bm{y}}{\textm{minimize}} & \bm{y}^\T\bSigma\bm{y}\\ \textm{subject to} & \bm{y}^\T\left(\bmu - r_\textm{f}\bm{1}\right) = 1\\ & \bm{y}\ge\bm{0}, \end{array} \] with \(\w = \bm{y}/\left(\bm{1}^\T\bm{y}\right)\).
Via the universal iterative SQP-MVP algorithm, i.e., we iteratively obtain \(\w^{k+1}\), for \(k=0,1,\dots,\), by solving \[ \begin{array}{ll} \underset{\w}{\textm{maximize}} & \w^\T\bmu - \dfrac{\lambda^k}{2}\w^\T\bSigma\w\\ \textm{subject to} & \bm{1}^\T\w=1, \quad \w\ge\bm{0}, \end{array} \] with \[ \lambda^k = \dfrac{(\w^k)^\T\bmu - r_\textm{f}}{(\w^k)^\T\bSigma\w^k}. \]
data(SP500_2015to2020)
stock_prices <- SP500_2015to2020$stocks["2019::", c("AAPL", "AMZN", "AMD", "GM", "GOOGL", "MGM", "MSFT", "QCOM", "TSCO", "UPS")]
X <- diff(log(stock_prices))[-1]
T <- nrow(X)
N <- ncol(X)
mu <- colMeans(X)
Sigma <- cov(X)
# using Schaible transform
w_ <- Variable(N)
prob <- Problem(Minimize(quad_form(w_, Sigma)),
constraints = list(w_ >= 0, t(mu) %*% w_ == 1))
result <- solve(prob)
w_cvx <- as.vector(result$getValue(w_)/sum(result$getValue(w_)))
obj_cvx <- as.numeric(t(w_cvx) %*% mu / sqrt(t(w_cvx) %*% Sigma %*% w_cvx))
# using the SQP-MVP algorithm
w <- rep(1/N, N)
obj_sqp <- c(t(w) %*% mu / sqrt(t(w) %*% Sigma %*% w))
k <- 0
for (k in 0:20) {
lmd_k <- as.numeric(t(w) %*% mu / (t(w) %*% Sigma %*% w))
w_prev <- w
w <- Variable(N)
prob <- Problem(Maximize(t(w) %*% mu - (lmd_k/2) * quad_form(w, Sigma)),
constraints = list(sum(w) == 1, w >= 0))
result <- solve(prob)
w <- as.vector(result$getValue(w))
obj_sqp <- c(obj_sqp, t(w) %*% mu / sqrt(t(w) %*% Sigma %*% w))
k <- k + 1
if (max(abs(w - w_prev)) < 1e-4)
break
}
# plot
data.frame(
"k" = c(0:k),
"SQP iterative method" = obj_sqp,
check.names = FALSE) |>
melt(id.vars = "k") |>
ggplot(aes(x = k, y = value, color = variable, fill = variable)) +
geom_hline(yintercept = obj_cvx, size = 0.8) +
geom_line(size = 0.8, show.legend = FALSE) +
geom_point(size = 1.5, show.legend = FALSE) +
ggtitle("Convergence") + ylab("objective value") + xlab("iteration")
Mean-volatility portfolio
We now consider two methods for the resolution of the mean-volatility portfolio:
Via an SOCP solver, i.e., solving \[ \begin{array}{ll} \underset{\w}{\textm{maximize}} & \w^\T\bmu - \kappa\sqrt{\w^\T\bSigma\w}\\ \textm{subject to} & \bm{1}^\T\w=1, \quad \w\ge\bm{0}, \end{array} \]
Via the universal iterative SQP-MVP algorithm, i.e., we iteratively obtain \(\w^{k+1}\), for \(k=0,1,\dots,\), by solving \[ \begin{array}{ll} \underset{\w}{\textm{maximize}} & \w^\T\bmu - \dfrac{\lambda^k}{2}\w^\T\bSigma\w\\ \textm{subject to} & \bm{1}^\T\w=1, \quad \w\ge\bm{0}, \end{array} \] with \[ \lambda^k = \kappa/\sqrt{(\w^k)^\T\bSigma\w^k}. \]
data(SP500_2015to2020)
stock_prices <- SP500_2015to2020$stocks["2019::", c("AAPL", "AMZN", "AMD", "GM", "GOOGL", "MGM", "MSFT", "QCOM", "TSCO", "UPS")]
X <- diff(log(stock_prices))[-1]
T <- nrow(X)
N <- ncol(X)
mu <- colMeans(X)
Sigma <- cov(X)
kappa <- 1.0
# using CVX
Sigma_12 <- chol(Sigma)
#max(abs(t(Sigma_12) %*% Sigma_12 - Sigma)) # sanity check
w <- Variable(N)
prob <- Problem(Maximize(t(w) %*% mu - kappa*cvxr_norm(Sigma_12 %*% w, 2)),
constraints = list(sum(w) == 1, w >= 0))
result <- solve(prob) #result$status
w_cvx <- as.vector(result$getValue(w))
obj_cvx <- as.numeric(t(w_cvx) %*% mu - kappa*sqrt(t(w_cvx) %*% Sigma %*% w_cvx)) # same as: result$value
# using the SQP-MVP algorithm
w <- rep(1/N, N)
obj_sqp <- c(t(w) %*% mu - sqrt(t(w) %*% Sigma %*% w))
k <- 0
for (k in 0:20) {
lmd_k <- kappa/sqrt(t(w) %*% Sigma %*% w)
w_prev <- w
w <- Variable(N)
prob <- Problem(Maximize(t(w) %*% mu - (lmd_k/2) * quad_form(w, Sigma)),
constraints = list(sum(w) == 1, w >= 0))
result <- solve(prob)
w <- as.vector(result$getValue(w))
obj_sqp <- c(obj_sqp, t(w) %*% mu - kappa*sqrt(t(w) %*% Sigma %*% w))
k <- k + 1
if (max(abs(w - w_prev)) < 1e-4)
break
}
# plot
data.frame(
"k" = c(0:k),
#"SOCP solver" = rep(obj_cvx, k+1),
"SQP iterative method" = obj_sqp,
check.names = FALSE) |>
melt(id.vars = "k") |>
ggplot(aes(x = k, y = value, color = variable, fill = variable)) +
geom_hline(yintercept = obj_cvx, size = 0.8) +
geom_line(size = 0.8, show.legend = FALSE) +
geom_point(size = 1.5, show.legend = FALSE) +
ggtitle("Convergence") + ylab("objective value") + xlab("iteration")