# basic finance
library(xts) # to manipulate time series of stock data
library(PerformanceAnalytics) # to compute performance measures
library(portfolioBacktest) # to conduct backtests
library(pob) # book package with financial data
# install with: devtools::install_github("dppalomar/pob")
library(fitHeavyTail) # to fit the multivariate skewed t distribution
library(highOrderPortfolios) # package for high-order portfolios
# plotting
library(ggplot2) # for nice plots
library(reshape2) # to reshape data
# optimization
library(CVXR)
library(quadprog)
library(nloptr)R Code for Portfolio Optimization
Chapter 9 – High-Order Portfolios
Daniel P. Palomar (2025). Portfolio Optimization: Theory and Application. Cambridge University Press.
Last update: March 14, 2025
Loading packages
The following packages are used in the examples:
Market data and benchmark portfolios
For illustrative purposes, we simply choose a few stocks from the S&P 500 index during 2015-2019.
data(SP500_2015to2020)
stock_prices <- SP500_2015to2020$stocks[, sample(ncol(SP500_2015to2020$stocks), 10)]
log(stock_prices) |>
fortify(melt = TRUE) |>
ggplot(aes(x = Index, y = Value, color = Series)) +
geom_line(linewidth = 1) +
scale_x_date(date_breaks = "6 months", date_labels = "%b %Y", date_minor_breaks = "1 week") +
labs(title = "Log-prices", x = "", y = "", color = "stocks")
As benchmark portfolios, we consider the following:
Global maximum return portfolio (GMRP): \[ \begin{array}{ll} \underset{\w}{\textm{maximize}} & \w^\T\bmu\\ \textm{subject to} & \bm{1}^\T\w=1, \quad \w\ge\bm{0}; \end{array} \]
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} \]
GMRP <- function(dataset, ...) {
N <- ncol(dataset$prices)
X <- diff(log(dataset$prices))[-1]
mu <- colMeans(X)
i_max <- which.max(mu)
w <- rep(0, N)
w[i_max] <- 1
return(w)
}
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)
}MVSK portfolio
The mean-variance-skewness-kurtosis (MVSK) portfolio is formulated as \[ \begin{array}{ll} \underset{\w}{\textm{minimize}} & - \lambda_1 \phi_1(\w) + \lambda_2 \phi_2(\w) - \lambda_3 \phi_3(\w) + \lambda_4 \phi_4(\w)\\ \textm{subject to} & \w \in \mathcal{W}, \end{array} \] where the hyper-parameters can be chosen, for example, as \(\lambda_1 = 1\), \(\lambda_2 = \frac{\gamma}{2}\), \(\lambda_3 = \frac{\gamma(\gamma + 1)}{6}\), and \(\lambda_4 = \frac{\gamma(\gamma + 1)(\gamma + 2)}{24}\).
We will consider three different ways to solve the problem, namely:
- simply by using an off-the-shelf nonlinear solver (e.g.,
nloptr); - implementing the SCA-Q-MVSK method described in (Palomar, 2025, Algorithm 9.2) (see (Zhou and Palomar, 2021) for more details); and
- conveniently using the MVSK solver implemented in the R package
highOrderPortfolios.
Recall that under the skewed \(t\) distribution (Palomar, 2025, Chapter 9):
- The moments are given by \[ \begin{aligned} \phi_1(\w) & = \w^\T\bmu + a_1 \w^\T\bm{\gamma}\\ \phi_2(\w) & = a_{21}\w^\T\bSigma\w + a_{22}(\w^\T\bm{\gamma})^2\\ \phi_3(\w) & = a_{31}(\w^\T\bm{\gamma})^3 + a_{32}(\w^\T\bm{\gamma})\w^\T\bSigma\w\\ \phi_4(\w) & = a_{41}(\w^\T\bm{\gamma})^4 + a_{42}(\w^\T\bm{\gamma})^2\w^\T\bSigma\w + a_{43}(\w^\T\bSigma\w)^2, \end{aligned} \] where \(a_1=\frac{\nu}{\nu-2}\), \(a_{21}=a_1\), \(a_{22}=\frac{2\nu^2}{(\nu-2)^2(\nu-4)}\), \(a_{31}=\frac{16\nu^3}{(\nu-2)^3(\nu-4)(\nu-6)}\), \(a_{32}=\frac{6\nu^2}{(\nu-2)^2(\nu-4)}\), \(a_{41}=\frac{(12\nu+120)\nu^4}{(\nu-2)^4(\nu-4)(\nu-6)(\nu-8)}\), \(a_{42}=\frac{6(2\nu+4)\nu^3}{(\nu-2)^3(\nu-4)(\nu-6)}\), and \(a_{43}=\frac{3\nu^2}{(\nu-2)(\nu-4)}\). Therefore, the code for the objective function \(f(\w) = - \lambda_1 \phi_1(\w) + \lambda_2 \phi_2(\w) - \lambda_3 \phi_3(\w) + \lambda_4 \phi_4(\w)\) is
compute_a <- function(nu) {
a21 <- a11 <- nu/(nu-2)
a22 <- (2 * nu ** 2)/((nu - 2) ** 2 * (nu - 4))
a31 <- 16 * (nu ** 3)/(((nu - 2) ** 3) * (nu - 4) * (nu - 6))
a32 <- 6 * (nu ** 2)/(((nu -2) ** 2) * (nu - 4))
a41 <- (12 * nu + 120) * (nu ** 4)/ (((nu - 2) ** 4) * (nu - 4) * (nu - 6) * (nu - 8))
a42 <- (2 * nu + 4) * (nu ** 3)/(((nu - 2) ** 3) * (nu - 4) * (nu - 6)) * 6
a43 <- (nu ** 2)/((nu - 2) * (nu - 4)) * 3
return(list("a11" = a11, "a21" = a21, "a22" = a22,
"a31" = a31, "a32" = a32,
"a41" = a41, "a42" = a42, "a43" = a43))
}
obj_value_skewt <- function(w, lambda, parameters) {
wT_mu <- t(w) %*% parameters$mu
wT_gamma <- t(w) %*% parameters$gamma
wT_Sigma_w <- sum((parameters$chol_Sigma %*% w) ** 2)
obj <- - lambda[1] * (wT_mu + parameters$a$a11 * wT_gamma) +
lambda[2] * (parameters$a$a21 * wT_Sigma_w + parameters$a$a22 * (wT_gamma**2) ) -
lambda[3] * (parameters$a$a31 * (wT_gamma**3) + parameters$a$a32 * wT_gamma * wT_Sigma_w) +
lambda[4] * (parameters$a$a41 * (wT_gamma**4) + parameters$a$a42 * wT_Sigma_w * (wT_gamma**2) + parameters$a$a43 * (wT_Sigma_w**2))
return(as.numeric(obj))
}- The gradients of the moments are \[ \begin{aligned} \nabla\phi_1(\w) &= \bmu + a_1 \bm{\gamma}\\ \nabla\phi_2(\w) &= 2a_{21}\bSigma\w + 2a_{22}(\w^\T\bm{\gamma})\bm{\gamma}\\ \nabla\phi_3(\w) &= 3a_{31}(\w^\T\bm{\gamma})^2\bm{\gamma} + a_{32}\left((\w^\T\bSigma\w)\bm{\gamma} + 2(\w^\T\bm{\gamma})\bSigma\w\right)\\ \nabla\phi_4(\w) &= 4a_{41}(\w^\T\bm{\gamma})^3\bm{\gamma}\\ &+ 2a_{42}\left((\w^\T\bm{\gamma})^2\bSigma\w + (\w^\T\bSigma\w)(\w^\T\bm{\gamma})\bm{\gamma}\right) + 4a_{43}(\w^\T\bSigma\w)\bSigma\w\\ \end{aligned} \] so the code for the gradient of the objective function is
grad_f_skewt <- function(w, lambda, parameters) {
wT_gamma <- as.numeric(t(w) %*% parameters$gamma)
wT_Sigma_w <- sum((parameters$chol_Sigma %*% w) ** 2)
grad <- - lambda[1] * (parameters$mu + parameters$a$a11 * parameters$gamma) +
lambda[2] * (2 * parameters$a$a21 * parameters$scatter %*% w +
parameters$a$a22 * 2 * wT_gamma * as.matrix(parameters$gamma)) -
lambda[3] * (parameters$a$a31 * 3 * (wT_gamma**2) * as.matrix(parameters$gamma) +
parameters$a$a32 * wT_Sigma_w * as.matrix(parameters$gamma) +
2 * wT_gamma * parameters$scatter %*% w) +
lambda[4] * (4 * parameters$a$a41 * (wT_gamma**3) *as.matrix(parameters$gamma) +
parameters$a$a42 * (2 * (wT_gamma**2) * parameters$scatter %*% w +
2 * wT_Sigma_w * wT_gamma * as.matrix(parameters$gamma)) +
4 * parameters$a$a43 * wT_Sigma_w * parameters$scatter %*% w)
return(grad)
}- The Hessians of the moments are \[ \begin{aligned} \nabla^2\phi_1(\w) &= \bm{0}\\ \nabla^2\phi_2(\w) &= 2a_{21}\bSigma + 2a_{22}\bm{\gamma}\bm{\gamma}^\T\\ \nabla^2\phi_3(\w) &= 6a_{31}(\w^\T\bm{\gamma})\bm{\gamma}\bm{\gamma}^\T + 2a_{32}\left(\bm{\gamma}\w^\T\bSigma + \bSigma\w\bm{\gamma}^\T + (\w^\T\bm{\gamma})\bSigma\right)\\ \nabla^2\phi_4(\w) &= 12a_{41}(\w^\T\bm{\gamma})^2\bm{\gamma}\bm{\gamma}^\T\\ &+ 2a_{42}\left(2(\w^\T\bm{\gamma})\bSigma\w\bm{\gamma}^\T + (\w^\T\bm{\gamma})^2\bSigma + 2(\w^\T\bm{\gamma})\bm{\gamma}\w^\T\bSigma + (\w^\T\bSigma\w)\bm{\gamma}\bm{\gamma}^\T\right)\\ &+ 4a_{43}\left(2\bSigma\w\w^\T\bSigma + (\w^\T\bSigma\w)\bSigma\right), \end{aligned} \] so the code for the Hessian of the objective function is
Hessian_f_skewt <- function(w, lambda, parameters) {
wT_gamma <- as.numeric(t(w) %*% parameters$gamma)
wT_Sigma_w <- sum((parameters$chol_Sigma %*% w) ** 2)
Sigma_w <- parameters$scatter %*% w
H_phi2 <- 2 * (parameters$a$a21 * parameters$scatter + parameters$a$a22 * parameters$gamma %*% t(parameters$gamma))
H_phi3 <- 6 * parameters$a$a31 * wT_gamma * parameters$gamma %*% t(parameters$gamma) + 2 * parameters$a$a32 * (parameters$gamma %*% t(Sigma_w) + Sigma_w %*% t(parameters$gamma) + wT_gamma * parameters$scatter)
H_phi4 <- 12 * parameters$a$a41 * (wT_gamma**2) * parameters$gamma %*% t(parameters$gamma) + 2 * parameters$a$a42 * (2 * wT_gamma * Sigma_w %*% t(parameters$gamma) + (wT_gamma**2) * parameters$scatter + 2 * wT_gamma * parameters$gamma %*% t(Sigma_w) + wT_Sigma_w * parameters$gamma %*% t(parameters$gamma)) + 4 * parameters$a$a43 * (wT_Sigma_w * parameters$scatter + 2 * Sigma_w %*% t(Sigma_w))
return(lambda[2] * H_phi2 - lambda[3] * H_phi3 + lambda[4] * H_phi4)
}1. Using off-the-shelf nonlinear solver
We can then write the code to solve the MVSK portfolio invoking the solver nloptr and using the package fitHeavyTail to fit the multivariate skewed \(t\) distribution with the function fit_mvst():
# hyper-parameters
xi <- 10
lambda <- c(0, xi/2, xi * (xi+1) / 6, xi * (xi + 1) * (xi + 2)/24)
MVSK_nloptr <- function(dataset, ...) {
N <- ncol(dataset$prices)
X <- diff(log(dataset$prices))[-1]
# fit skewed t distribution to data
options(nu_min = 8.5)
parameters <- fit_mvst(X)
parameters$a <- compute_a(parameters$nu)
parameters$chol_Sigma <- chol(parameters$scatter)
# define auxiliary function for equality constraints
eval_g_eq <- function(w, lambda, parameters) {
constr <- c(sum(w) - 1)
grad <- c(rep(1, N))
return(list("constraints" = constr,
"jacobian" = grad))
}
# call solver
res <- nloptr(x0 = pmax(design_GMVP(parameters$scatter), 1e-5), #rep(1/N, N), parameters$scatter
eval_f = obj_value_skewt,
eval_grad_f = grad_f_skewt,
lb = rep(0, N),
ub = rep(Inf, N),
eval_g_eq = eval_g_eq,
opts = list("algorithm" = "NLOPT_LD_SLSQP", maxeval = 2500, "xtol_rel" = 1e-4,
"local_opts" = list( "algorithm" = "NLOPT_LD_SLSQP",
"xtol_rel" = 1e-4,
"maxeval" = 2500)),
lambda = lambda,
parameters = parameters)
return(res$solution)
}2. Using the SCA-based method
The SCA-Q-MVSK method is described in (Palomar, 2025, Algorithm 9.2) (see (Zhou and Palomar, 2021) for more details). The corresponding R code is:
project_PSD <- function(Q, tau) {
N <- ncol(Q)
eigen_Q <- eigen(Q)
d <- eigen_Q$values
if (min(d) < 0) {
d[which(d < 0)] <- 0
Q <- eigen_Q$vectors %*% diag(d) %*% t(eigen_Q$vectors)
}
return(Q + tau * diag(N))
}
MVSK_Q_MVSK <- function(dataset, ...) {
N <- ncol(dataset$prices)
X <- diff(log(dataset$prices))[-1]
# fit skewed t distribution to data
options(nu_min = 8.5)
parameters <- fit_mvst(X)
parameters$a <- compute_a(parameters$nu)
parameters$chol_Sigma <- chol(parameters$scatter)
# main loop
wk <- pmax(design_GMVP(parameters$scatter), 1e-5) # initial point
iter <- 0
while(iter < 100) {
iter <- iter + 1
wk_old <- wk
Qk <- Hessian_f_skewt(wk, lambda, parameters)
qk <- grad_f_skewt(wk, lambda, parameters) - Qk %*% wk
Qk <- project_PSD(Qk, tau = 1e-8)
wk <- solve.QP(Dmat = Qk/norm(Qk,"2"), # for numerical stability we normalize with norm(Qk,"2")
dvec = -qk/norm(Qk,"2"),
Amat = t(rbind(rep(1, N), diag(N))),
bvec = c(1, rep(0, N)),
meq = 1)$solution
wk[which(wk < 0)] <- 0 # just in case
wk <- wk/sum(wk)
has_w_converged <- all(abs(wk - wk_old) <= .5 * 1e-8 * (abs(wk) + abs(wk_old)))
if (has_w_converged) break
}
return(wk)
}3. Using the R package highOrderPortfolios
The R package highOrderPortfolios conveniently implements algorithms to solve the MVSK portfolio and other variations, such as the MVSK tilting portfolio. See (Palomar, 2025, Chapter 9) and (Wang et al., 2023) for more details. The code is simply:
library(highOrderPortfolios)
MVSK_package <- function(dataset, ...) {
X <- diff(log(dataset$prices))[-1]
# fit skewed t distribution to data
X_skew_t_params <- estimate_skew_t(X)
# call solver
w <- design_MVSK_portfolio_via_skew_t(lambda, X_skew_t_params,
w_init = pmax(design_GMVP(X_skew_t_params$scatter), 1e-5),
method = "Q-MVSK")$w
return(w)
}Sanity check of the different MVSK portfolio methods
We first make sure that the three described ways of solving the MVSK portfolio design agree:
w_MVSK_nloptr <- MVSK_nloptr(dataset = list("prices" = stock_prices))
w_MVSK_Q_MVSK <- MVSK_Q_MVSK(dataset = list("prices" = stock_prices))
w_MVSK_package <- MVSK_package(dataset = list("prices" = stock_prices))
data.frame(
"stocks" = names(stock_prices),
"Nonlinear solver nloptr" = w_MVSK_nloptr,
"SCA-Q-MVSK" = w_MVSK_Q_MVSK,
"Package highOrderPortfolios" = w_MVSK_package,
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(fill = "Algorithms") +
ylab("weight") +
ggtitle("MVSK portfolio allocation for different algorithms")
Backtest
We are now ready to perform a backtest of the MVSK portfolio and the benchmarks. We will use the package portfolioBacktest.
data(SP500_2015to2020)
stock_prices <- SP500_2015to2020$stocks["2015::2019", sample(ncol(SP500_2015to2020$stocks), 20)]
bt <- portfolioBacktest(list("GMRP" = GMRP,
"GMVP" = GMVP,
"MVSK" = MVSK_Q_MVSK),
list("dataset1" = list("prices" = stock_prices)),
lookback = 252, optimize_every = 20, rebalance_every = 5)
listPortfoliosWithFailures(bt) # just in case there are errorsbacktestChartCumReturn(bt) +
scale_x_date(date_breaks = "6 months", date_labels = "%b %Y", date_minor_breaks = "1 week") +
ggtitle("Cumulative return")
backtestChartDrawdown(bt) +
scale_x_date(date_breaks = "6 months", date_labels = "%b %Y", date_minor_breaks = "1 week") +
ggtitle("Drawdown")
summary_bt <- backtestSummary(bt)
summaryTable(summary_bt, type = "DT",
measures = c("annual return", "annual volatility", "max drawdown",
"Sharpe ratio", "Sortino ratio", "Sterling ratio",
"VaR (0.95)", "CVaR (0.95)"))Multiple backtests
We now consider multiple randomized backtests (Palomar, 2025, Chapter 8) for a better characterization of the performance of the portfolios. Again, we will use the package portfolioBacktest.
We take a dataset of \(N=20\) stocks over the period 2015-2020 and generate 100 resamples each with \(N=8\) randomly selected stocks and a random period of two years. Then we perform a walk-forward backtest with a lookback window of 1 year, reoptimizing the portfolio every month.
stock_prices <- SP500_2015to2020$stocks[, sample(ncol(SP500_2015to2020$stocks), 20)]
# resample data
stock_prices_resampled <- financialDataResample(list("prices" = stock_prices),
num_datasets = 100, N_sample = 8, T_sample = 252*2)
# perform multiple backtest on the resampled data
bt <- portfolioBacktest(list("GMRP" = GMRP,
"GMVP" = GMVP,
"MVSK" = MVSK_Q_MVSK),
dataset_list = stock_prices_resampled,
lookback = 252, optimize_every = 20, rebalance_every = 5,
paral_datasets = 6)
listPortfoliosWithFailures(bt) # just in case there are errorssummary_bt <- backtestSummary(bt)
summaryTable(summary_bt, type = "simple",
measures = c("annual return", "annual volatility", "max drawdown",
"Sharpe ratio", "Sortino ratio", "Sterling ratio",
"VaR (0.95)", "CVaR (0.95)")) annual return annual volatility max drawdown Sharpe ratio Sortino ratio
GMRP 0.09 0.26 0.20 0.33 0.45
GMVP 0.13 0.13 0.12 1.05 1.48
MVSK 0.14 0.13 0.12 1.15 1.58
Sterling ratio VaR (0.95) CVaR (0.95)
GMRP 0.41 0.03 0.04
GMVP 1.18 0.01 0.02
MVSK 1.29 0.01 0.02summaryBarPlot(summary_bt, measures = c("max drawdown", "annual volatility"))
backtestBoxPlot(bt, "Sharpe ratio")
backtestBoxPlot(bt, "max drawdown")
References
Palomar, D. P. (2025). Portfolio optimization: Theory and application. Cambridge University Press.
Wang, X., Zhou, R., Ying, J., and Palomar, D. P. (2023). Efficient and scalable parametric high-order portfolios design via the skew-t distribution. IEEE Transactions on Signal Processing, 71, 3726–3740.
Zhou, R., and Palomar, D. P. (2021). Solving high-order portfolios via successive convex approximation algorithms. IEEE Transactions on Signal Processing, 69, 892–904.