Portfolio Optimization
Chapter 9: High-Order Portfolios

R code

Published

October 1, 2024

R code examples for Chapter 9 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:

# 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)

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:

  1. simply by using an off-the-shelf nonlinear solver (e.g., nloptr);
  2. implementing the SCA-Q-MVSK method described in (Palomar, 2024, Algorithm 9.2) (see (Zhou and Palomar, 2021) for more details); and
  3. conveniently using the MVSK solver implemented in the R package highOrderPortfolios.

Recall that under the skewed \(t\) distribution (Palomar, 2024, 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, 2024, 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, 2024, 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 errors
backtestChartCumReturn(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, 2024, 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 errors
summary_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.02
summaryBarPlot(summary_bt, measures = c("max drawdown", "annual volatility"))

backtestBoxPlot(bt, "Sharpe ratio")

backtestBoxPlot(bt, "max drawdown")

References

Palomar, D. P. (2024). 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.