Introduction

The purpose of this vignette is to demonstrate a sample of portfolio optimization problems that can be solved by using the ROI package. This vignette is based on joint work with Florian Schwendinger and Ronald Hochreiter which was presented at RFinance 2016, Chicago, USA, May 2016.

If previously no ROI version was installed, one should at least install the two plugins ROI.plugin.glpk, ROI.plugin.alabama and ROI.plugin.quadprog to run the examples below.

## NOTE: (Rglpk needs libglpk-dev on Debian/Ubuntu ("sudo apt-get install libglpk-dev"))

if (!require("ROI")) install.packages("ROI"); library("ROI")

if (!require("ROI.plugin.glpk")) install.packages("ROI.plugin.glpk"); library("ROI.plugin.glpk")

if (!require("ROI.plugin.quadprog")) install.packages("ROI.plugin.quadprog"); library("ROI.plugin.quadprog")

if (!require("ROI.plugin.alabama")) install.packages("ROI.plugin.alabama"); library("ROI.plugin.alabama")

Several R functions are created to implement the typical objectives and constraints used for portfolio optimization. All functions require a data.frame r_mat of returns. The mathematical formulation of the objectives and constraints is presented below. The default optimization in ROI is minimization. Hence for all maximization problems the argument maximum = TRUE in the OP() function must be specified.


Objectives implemented as R functions

Name Arguments Objectives Max/Min Type
reward_objective (r_mat) Portfolio return/reward max LP
markowitz_objective (r_mat) Variance of portfolio min QP
mad_objective (r_mat) Mean Absolute Deviation min LP
downside_var_objective (r_mat) Lower semi-variance of portfolio min QP
downside_mad_objective (r_mat) Lower semi-mean absolute deviation min LP
cvar_objective (r_mat, alpha, probs = NULL) Conditional Value-at-Risk min LP
minimax_young_objective (r_mat) Minimax portfolio max LP
quadratic_utility_objective (r_mat, lambda = 2) Quadratic Utility max LP
sharpe_objective (r_mat, rf = 0) Sharpe ratio max QP
omega_objective (r_mat, tau = 0) Omega max LP


Constraints implemented as R functions

Name Arguments Constraints Type
budget_constraint (r_mat, dir = "==", rhs = 1) Budget constraint,i.e., sum of the portfolio equals budget B (default B=rhs=1) LP
group_constraint (r_mat, index, coef.index = 1, dir = "==", rhs) Group constraint applied only the the index elements of \(x\) LP
turnover_constraint (r_mat, x0 = NULL, dir = "<=", rhs = 100) Constraint on the turnover from some initial weights \(x_0\) LP
cardinality_constraint (r_mat, dir = "<=", rhs = 100) Cardinality MILP
reward_constraint (r_mat, dir = ">=", rhs = 0) Target return constraint LP
markowitz_constraint (r_mat, dir = "<=", rhs = 1000) Constraint on variance of portfolio \(x^{\top}Qx \leq q^2\) QP
cvar_constraint (r_mat, alpha, probs = NULL, dir = "<=", rhs = 100) Constraints on the conditional Value-at-Risk LP

Objective functions

We consider a portfolio of assets with random returns. We denote the portfolio choice vector by \(x\) and by \(r\) the vector of random returns. \(N\) denotes the number of assets in the portfolio while \(S\) is the number of scenarios in the scenarios set.

Maximize expected return

\[ \max_{x \in \mathbb{R}^N} \hat \mu^\top x \] where \(\hat \mu\) is the vector of estimated mean asset returns.

In ROI we define a function which returns a linear objective is L_objective (note that the minus turns the maximization into a minimization problem):

reward_objective <- function(r_mat) {
    objective <- L_objective(colMeans(r_mat))
    list(objective = objective)
}

(see Example 1, Example 12, Example 13).

Minimize variance

The objective of the minimum of the variance portfolio is given by: \[\min_{x \in \mathbb{R}^N} x^\top Q x\] \[Q = \mathbb{C}ov(r)\]

markowitz_objective <- function(r_mat) {
    objective <- Q_objective(Q = 2 * cov(r_mat), 
                             L = rep(0, NCOL(r_mat)))
    list(objective = objective)
}

(see Example 2, Example 3).

Minimize mean absolute deviation

The minimization of the mean absolute deviation can be written as a linear problem. \[ \begin{eqnarray*} \min_{x\in \mathbb{R}^N, y\in \mathbb{R}^S, z\in \mathbb{R}^S}&& \frac{1}{S} \sum_{s = 1}^S (y_s + z_s)\\ s.t.&&y_s - z_s = (r_{s} - \hat \mu)^\top x,\\ && y_s \geq 0, z_s \geq 0. \end{eqnarray*} \] Note that the bounds \(y_s\geq 0\) and \(z_s\geq 0\) must not be additionally specified as these are the default in ROI.

mad_objective <- function(r_mat){
    x.names <- colnames(r_mat)
    N <- NCOL(r_mat)
    S <- nrow(r_mat)
    mu <- colMeans(r_mat)
    Amat <-  cbind(sweep(as.matrix(r_mat), 2,  mu), - diag(S), diag(S))
    var.names <- c(x.names, 
                   paste0("y_mad_aux", seq_len(S)),
                   paste0("z_mad_aux", seq_len(S)))
    
    constraint <-  L_constraint(L = Amat, dir = rep("==", S), 
                                rhs = rep(0, S), 
                                names = var.names)

    objective <- L_objective(L = c(rep(0, N), rep(1/S, 2 * S)))
    
    list(objective = objective, constraint = constraint)
}

(See Example 4)

Minimize lower semi-variance

A downside risk measure is the lower semi-variance, which is the expected squared deviation from the mean, calculated over those points that are no greater than the mean. Only returns that are below the mean contribute to the portfolio risk. Let \(r_s\) denote the vector of random returns of the \(N\) assets in scenario \(s\). The problem can be formulated as (using the auxiliary variable \(z\in \mathbb{R}^S\)): \[\min_{x\in \mathbb{R}^N} \frac{1}{S} \sum_{s = 1}^S \mathrm{max}(\hat\mu^\top x - r_s^\top x, 0)^2.\] The problem can be reformulated as a QP: \[ \begin{eqnarray*} \min_{x\in \mathbb{R}^N, z\in \mathbb{R}^S}&& \frac{1}{S} \sum_{s = 1}^S z_s^2\\ s.t.&&z_s \geq (\hat\mu - r_s)^\top x)\\ && z_s \geq 0. \end{eqnarray*} \]

downside_var_objective <- function(r_mat, x.names = NULL) {
    x.names <- colnames(r_mat)
    N <- NCOL(r_mat)
    S <- NROW(r_mat)
    mu <- colMeans(r_mat)
    Amat <- cbind(sweep(as.matrix(r_mat), 2, mu), diag(S))
    var.names <- c(x.names, paste0("z_dvar_aux", seq_len(S)))
    
    constraint <- L_constraint(L = Amat, dir = rep(">=", S), 
                               rhs = rep(0, S), 
                               names = var.names)
    objective <- Q_objective(
             Q = 2 * diag(c(rep(1e-05, N), rep(1/S, S))),
             L = rep(0, N + S))
    
    list(objective = objective, constraint = constraint)
}

(see Example 5).

Minimize lower semi-absolute deviation

Let \(r_s\) denote the vector of random returns of the \(N\) assets in scenario \(s\). The downside risk measure given by the expected absolute deviation from the mean, calculated over those points that are no greater than the mean is given by: \[\frac{1}{S} \sum_{s = 1}^S \vert \mathrm{max}(\hat\mu^\top x - r_s^\top x, 0) \vert\] The minimization of this risk measure can be reformulated as a linear problem: \[ \begin{eqnarray*} \min_{x\in \mathbb{R}^N, z\in \mathbb{R}^S}&& \frac{1}{S} \sum_{s = 1}^S z_s,\\ s.t.&& z_s \geq (\hat\mu - r_s)^\top x\\ && z_s \geq 0. \end{eqnarray*} \]

downside_mad_objective <- function(r_mat){
    x.names <- colnames(r_mat)
    N <- NCOL(r_mat)
    S <- NROW(r_mat)
    mu <- colMeans(r_mat)
    Amat <- cbind(sweep(as.matrix(r_mat), 2, mu), diag(S))
    var.names <- c(x.names, paste0("z_dmad_aux", seq_len(S)))
    
    constraint <- L_constraint(L = Amat, dir = rep(">=", S),
                               rhs = rep(0, S), 
                               names = var.names)
    
    objective  <- L_objective(L = c(rep(0, N), rep(1/S,  S)))
    
    list(objective = objective, constraint = constraint)    
}

(see Example 6).

Minimize conditional value at risk

The conditional value at risk/expected shortfall is the average of the losses that exceed the \(\alpha\)-quantile of the loss distribution, also called \(\alpha\) Value-at-Risk (VaR). Defining the loss distribution as minus of the return of the portfolio, the problem can be formulated as follows: \[ \begin{eqnarray*} \min_{x\in \mathbb{R}^N, z\in \mathbb{R}^S, \gamma\in \mathbb{R}}& \gamma + \frac{1}{(1-\alpha)} \sum_{s = 1}^S p_s z_s\\ s.t.\ &z_s \geq - r_s^\top x - \gamma,\\ &z_s \geq 0, \end{eqnarray*} \] where VaR is a minimizer of the function over \(\gamma\) and \(p=(p_1, \dots, p_S)\) is a vector of probabilities associated with the scenarios.

cvar_objective <- function(r_mat, alpha, probs = NULL) {
  x.names <- colnames(r_mat)
  N <- NCOL(r_mat)
  S <- NROW(r_mat)
  mu <- colMeans(r_mat)
  if (is.null(probs)) probs <- rep(1/S, S)
  if (alpha < 0.5) alpha <- 1 - alpha
  
  Amat <- cbind(as.matrix(r_mat),  diag(S), 1)
  var.names <- c(x.names, paste0("z_cvar_aux", seq_len(S)), "gamma")
  
  ## set bounds for gamma (-Inf, Inf) 
  bnds <- ROI::V_bound(li = c(N + S + 1), lb = c( -Inf),
                       ui = c(N + S + 1), ub = c(  Inf))
  
  constraint <- L_constraint(L = Amat, dir = rep(">=", S), 
                            rhs = rep(0, S), 
                            names = var.names)

  objective <- L_objective(c(rep(0, N), probs/(1 - alpha), 1))

  list(objective = objective, constraint = constraint, bounds = bnds)
}

Note that if no value is specified for \(p_1, \cdots, p_S\) the default is equal weights.

(see Example 7, Example 14).

Minimax Portfolio

The optimal portfolio is defined as that one that minimizes the maximum loss over all past historical periods, subject to a restriction on the minimum acceptable average return across all observed periods of time. The minimax portfolio maximizes the minimum gain and can be seen as a limiting case of CVaR for \(\alpha \rightarrow 1\). Let \(M_p=\min_s r_s^\top x\) denote the minimum gain of the portfolio. \[ \begin{eqnarray*} \max_{x \in \mathbb{R}^N, M_p\in \mathbb{R}}&& M_p\\ s.t.&% r_s^\top x - M_p \geq 0. \end{eqnarray*} \] Young, M.R., 1998. A minimax portfolio selection rule with linear programming solution. Management science, 44(5), pp.673-683.

minimax_young_objective <- function(r_mat){
    x.names <- colnames(r_mat)
    N <- NCOL(r_mat)
    S <- NROW(r_mat)
    Amat <- cbind(as.matrix(r_mat), -1)
    
    bnds <- ROI::V_bound(li = c(N + 1), lb = c( -Inf),
                          ui = c(N + 1), ub = c(  Inf))
    constraint <- L_constraint(L = Amat, dir=rep(">=", S), rhs=rep(0, S), 
                                      names=c(x.names,  "mp_aux"))
  
    objective <- L_objective(c(rep(0, N), 1))
  
    list(objective = objective, constraint = constraint, bounds = bnds)     
}

(see Example 8).

Maximize quadratic utility

Let \(\lambda\) denote the risk aversion parameter. Typical risk aversion parameters lie between 2 and 4. Then the quadratic utility is given by: \[\begin{eqnarray*} \max_{x\in\mathbb{R}^N}&& \hat \mu^\top x - \frac{\lambda}{2} x^\top Q x \end{eqnarray*}\]
quadratic_utility_objective <- function(r_mat, lambda = 2){
    objective <- Q_objective(Q = - lambda * cov(r_mat), 
                             L = colMeans(r_mat))
    list(objective = objective)
}

(see Example 9).

Maximize Sharpe ratio

Let \(r_f\) denote the risk free rate (when not specified is set by default to zero). The Sharpe ratio of a random return of an asset \(\tilde r\) is defined as: \[ \frac{\mathbb{E}(\tilde r)- r_f}{\sqrt{\mathbb{V}ar(\tilde r)}} \] The budget normalization constraint is applied and the problem can be formulated as a QP: \[ \begin{eqnarray*} \max_{y\in \mathbb{R}^N, \kappa\in \mathbb{R}}& - y^\top Q y\\ s.t.&(\hat\mu - r_f \mathit{1})^\top y = \text{1},\\ & \mathit{1}^\top y = \kappa,\\ & \kappa > 0. \end{eqnarray*} \] The optimal solution is given by \(x^* = y^* / \kappa^*\).

sharpe_objective <- function(r_mat, rf = 0){
    N <- NCOL(r_mat)
    S <- NROW(r_mat)
    mu <- colMeans(r_mat)
    Amat <- rbind(c(mu - rf, 0),
                  c(rep(0, N), 1),
                  c(rep(1, N), -1))
    var.names <- c(paste0("y_sharpe_aux", seq_len(N)), "kappa_sharpe")
    
    constraint <- L_constraint(L = Amat, dir = c("==", ">", "=="),
                               rhs = c(1, 0, 0), names = var.names)
    
    mat <- matrix(0, ncol = N + 1, nrow = N + 1)
    mat[1:N, 1:N] <- 2 * cov(r_mat)
    mat[N + 1, N + 1] <-  1e-04

    objective <- Q_objective(Q = - mat, L = c(rep(0, N), 0))
  
    list(objective = objective, constraint = constraint)     
}

(see Example 10).

Maximize Omega

The Omega Ratio is a risk-return performance measure of an investment asset, portfolio, or strategy. The Omega Ratio, introduced in 2002 by Keating and Shadwick, is defined as the probability weighted ratio of gains versus losses for some threshold return target \(\tau\). \[\Omega(\tilde r) = \frac{\int_\tau^{+\infty}\left(1 - F(r)\right)\mathrm{d}r}{\int_{-\infty}^\tau F(r)\mathrm{d}r} \] where \(F(r)\) denotes the cumulative distribution function of \(\tilde r\). The maximization of \(\Omega(\tilde r)\) can be formulated as an linear problem (note that the budget normalization constraint and the target return constraint are already part of the model formulation). \[\begin{eqnarray*} \max_{y\in \mathbb{R}^N, u\in \mathbb{R}^S, z\in \mathbb{R}}& \hat \mu^\top y - \tau z\\ s.t.\ & u_s \geq \tau - r_s^\top y , \\ & u_i\geq 0,\\ & \mathit{1}^\top u = 1,\\ & \mathit{1}^\top y = z,\\ & \hat \mu^\top y \geq \tau z,\\ & z \geq 0 . \end{eqnarray*}\]

The Omega optimal portfolio is given by \(x^* = y^* /z^*\), where \(z\) is the homogenizing variable. If not explicitly specified, the default for the threshold is \(\tau = 0\).

omega_objective <- function(r_mat, tau = 0){
  ## variables y_1, ... y_N, u_1, .., u_S, z
  N <- NCOL(r_mat)
  S <- NROW(r_mat)
  mu <- colMeans(r_mat)
  
  Amat <- rbind(cbind(as.matrix(r_mat), diag(S), 0),# u_s >= tau - r_s'y
                c(rep(0, N), rep(1, S),  0), # sum(u) = 1
                c(rep(1, N), rep(0, S), -1), # sum(y) = z
                c(mu,        rep(0, S), - tau),
                c(rep(0, N), rep(0, S),  1)) # mu'y  >= tau * z
  var.names <- c(paste0("y_omega_aux", seq_len(N)),
                 paste0("u_omega_aux", seq_len(S)), 
                 "z_omega")
  constraint <-  L_constraint(L = Amat,  
                              dir = c(rep(">=", S), "==", "==", ">=", ">"),
                              rhs = c(rep(tau, S), 1, 0, 0, 1e-05), 
                              names=var.names)
  objective  <- L_objective(L = c(mu, rep(0, S), -tau))
  
  list(objective = objective, constraint = constraint)    
  ## x* <- y*/z*
}

(see Example 11).

Constraints

Box constraints

Can be set by the argument bounds in OP(). Note that the default in ROI is \(0 \leq x_i \leq \infty\).

Budget constraint

\[\sum_{i=1}^N x_i = B \qquad \sum_{i=1}^N x_i \leq B_u \qquad \sum_{i=1}^N x_i \geq B_l.\] Budget normalization constraint is obtained for B = 1 (default).

budget_constraint <- function(r_mat, dir = "==", rhs = 1) {
    x.names <- colnames(r_mat)
    L_constraint(L = rep(1, NCOL(r_mat)), 
                 dir = dir,  rhs = rhs, names = x.names)
}

Group constraints

Group constraints are linear constraints which apply to only some elements of the portfolio choice vector \(x\): \[ ax_i + bx_j = c \] where the \(=\) can be replaced by inequalities.

group_constraint <- function(r_mat, index, coef.index = 1, dir = "==", rhs) {
  ## index = (i, j)
  ## coef.index = c(a,b)
  ## rhs = c
  x.names <- colnames(r_mat)
  N <- NCOL(r_mat)
  L <- rep(0, N)
  L[index] <- coef.index
  L_constraint(L = L, dir = dir, rhs = rhs, names = x.names)
}

(see Example 1, Example 12).

Turnover constraint

For a target turnover \(L\) and some initial weights \(x_0\), the turnover contraint is \[\sum_{i=1}^N \vert x_i - x_{i0}\vert \leq L.\] This can be reformulated in terms of auxiliary variables \(y^+, y^- \in \mathbb{R}^{+N}\) such that \[\begin{eqnarray*} &&y^+_i- y_i^- = x_i - x_{0i}\\ &&\sum_{i=1}^N (y_i^+ + y_i^-)\leq L\\ &&y^+ \geq 0, \, y^-\geq 0 \end{eqnarray*}\]

If no initial weights are specified, the default is equal weights. Note: does not work out of the box with Omega and Sharpe objectives.

turnover_constraint <- function(r_mat, x0 = NULL, dir = "<=", rhs = 100) {
    x.names <- colnames(r_mat)
    N <- NCOL(r_mat)
    S <- NROW(r_mat)
    if (is.null(x0)) x0 <- rep(1/N, N)
    Amat <- cbind(diag(N), - diag(N), diag(N))
    var.names <- c(x.names,  
                   paste0("y_plus_aux", seq_len(N)), 
                   paste0("y_minus_aux", seq_len(N)))
    
    rbind(L_constraint(L = Amat, dir = rep("==", N), rhs = x0, 
                       names = var.names),
          L_constraint(c(rep(0, N), rep(1, N), rep(1, N) ), 
                       dir = dir, rhs = rhs, names = var.names))
}

(see Example 12).

Cardinality constraint

The cardinality constraint assigns bounds \(P_{min}\) and \(P_{max}\) on the number of assets in a portfolio. Binary auxiliary variables \(z\in \mathbb{R}^N\) are introduced such that \(P_{min} \leq \sum_{i=1}^N z_i\leq P_{max}\) and \(x_i - z_i \leq 0\).

cardinality_constraint <- function(r_mat, dir = "<=", rhs = 100) {
  x.names <- colnames(r_mat)
  N <- NCOL(r_mat)
  Amat <- cbind(diag(N), -diag(N))
  var.names <- c(x.names, paste0("z_card_aux", seq_len(N)))
  cat("Variable types for z_card_aux must be set to binary.\n")
  rbind(L_constraint(L = Amat, dir = rep("<=", N), 
                     rhs = rep(0, N), names = var.names), 
        L_constraint(L = c(rep(0, N), rep(1, N)), dir = dir,  
                     rhs = rhs, names = var.names))
}

(see Example 14).

Target return constraint

Also typical objective functions can be used as constraints. For example, the expected return must not fall below a value \(\tau\): \[\hat\mu^\top x \geq \tau\] Note: does not work with Sharpe or Omega (it is included in both Sharpe and Omega maximization problems).

reward_constraint <- function(r_mat, dir = ">=", rhs = 0) {
  x.names <- colnames(r_mat)
  L_constraint(L = colMeans(r_mat), dir = dir,  
               rhs = rhs, names = x.names)
}

(see Example 3, Example 6).

Risk constraints – Variance

Typically constraints of the form \(x^\top Q x \leq q^2\). Note however that this is a quadratic constraint and needs appropriate solvers (e.g., in ROI.plugin.alabama). Note: does not work with Omega.

markowitz_constraint <- function(r_mat, dir = "<", rhs = 1000) {
    x.names <- colnames(r_mat)
    N <- NCOL(r_mat)
    Q_constraint(Q = 2 * cov(r_mat), L = rep(0, N), 
                 dir = dir, rhs = rhs, names=x.names)
}

(see Example 13).

Risk constraints – CVaR

Typically CVaR\((x, \alpha) \leq q\). Note: does not work with Omega and Sharpe.

cvar_constraint <- function(r_mat, alpha, probs = NULL, dir = "<=", rhs = 100) {
  x.names <- colnames(r_mat)
  N <- NCOL(r_mat)
  S <- NROW(r_mat)
  if (alpha < 0.5) alpha <- 1 - alpha
  if (is.null(probs)) probs <- rep(1/S, S)
  
  Amat <- cbind(as.matrix(r_mat), diag(S), 1)
  var.names <- c(x.names, paste0("z_cvar_aux", seq_len(S)), "gamma")
  # set bounds for gama
  bnds <- ROI::V_bound(li = c(N + S + 1), lb = c( -Inf),
                       ui = c(N + S + 1), ub = c(  Inf))
  rbind(L_constraint(L = Amat, dir = rep(">=", S), 
                     rhs=rep(0, S), names = var.names), 
        L_constraint(c(rep(0, N), rep(1 / ((1 - alpha) * S), S), 1), 
                     dir = dir, rhs = rhs, names = var.names))
}    

(see Example 12).

Data

For illustration purposes daily log returns for the past 180 days of the 30 companies included in the 2018 Dow Jones Industrial Average will be used.

if (!require("quantmod")) install.packages("quantmod"); library("quantmod")
Tickers     <- c("MMM", "AXP","AAPL", "BA","CAT","CVX","CSCO", "KO", "DWDP", "XOM",
                 "GS","HD", "IBM", "INTC", "JNJ", "JPM","MCD","MRK","MSFT",
                 "NKE", "PFE","PG",
                 "TRV", "UNH","UTX", "VZ", "V", "WMT", "WBA", "DIS")

cached_file <- "cache/cached_portfolio_data.rda"
if ( file.exists(cached_file) ) {
  load(cached_file)
} else {
  download_tickers <- function(x) {
    download_ticker_error_handler <- function(e) {
      warning(paste(x, 'not found'))
      return(NA)
    }
    tryCatch(getSymbols(x, from = Sys.Date() - 180, auto.assign = FALSE),
             error = download_ticker_error_handler)
  }
  Historicals <- lapply(Tickers, download_tickers)

  prices <- do.call("cbind",lapply(Historicals, function(x) x[, 6]))

  djia2018 <- as.data.frame((log(lag(prices)) - log(prices))[-1, ])
  save(Historicals, prices, djia2018, file = cached_file)
}

We consider a portfolio of assets with random returns. We denote the portfolio choice vector by \(x\) and by \(r\) the vector random returns. \(N\) denotes the number of assets in the portfolio while \(S\) is the number of scenarios in the scenarios set.

Examples

Example 1: Maximize expected return subject to budget normalization and group constraints

The following optimization problem: \[ \begin{eqnarray*} \max_{x \in \mathbb{r}^N}&& \hat \mu^\top x\\ \sum_{x=1}^N x_i &=& 1\\ x_i &\geq& 0\\ x_3 + x_{17} &\leq& 0.5 \end{eqnarray*} \] can be set up easily in ROI:

lp  <- OP(objective  =  reward_objective(djia2018)$objective,
          constraints = rbind(budget_constraint(djia2018),
                              group_constraint(djia2018, index = c(3, 17), dir = "==", rhs = 0.5)),
          maximum = T)

(Note that the bounds can be omitted as this is the default in ROI). To perfom the optimization, the ROI_solve() function is called:

(sol <- ROI_solve(lp, solver = "glpk"))
## Optimal solution found.
## The objective value is: 1.516122e-03
solution(sol)
##  MMM.Adjusted  AXP.Adjusted AAPL.Adjusted   BA.Adjusted  CAT.Adjusted 
##           0.0           0.0           0.5           0.0           0.0 
##  CVX.Adjusted CSCO.Adjusted   KO.Adjusted DWDP.Adjusted  XOM.Adjusted 
##           0.0           0.0           0.0           0.0           0.0 
##   GS.Adjusted   HD.Adjusted  IBM.Adjusted INTC.Adjusted  JNJ.Adjusted 
##           0.5           0.0           0.0           0.0           0.0 
##  JPM.Adjusted  MCD.Adjusted  MRK.Adjusted MSFT.Adjusted  NKE.Adjusted 
##           0.0           0.0           0.0           0.0           0.0 
##  PFE.Adjusted   PG.Adjusted  TRV.Adjusted  UNH.Adjusted  UTX.Adjusted 
##           0.0           0.0           0.0           0.0           0.0 
##   VZ.Adjusted    V.Adjusted  WMT.Adjusted  WBA.Adjusted  DIS.Adjusted 
##           0.0           0.0           0.0           0.0           0.0

Example 2: Minimum variance portfolio

lp  <- OP(objective  =  markowitz_objective(djia2018)$objective,
          constraints = rbind(budget_constraint(djia2018)))

(Note that the bounds can be omitted as this is the default in ROI). To perfom the optimization, the ROI_solve() function is called:

(sol <- ROI_solve(lp, solver = "quadprog"))
## Optimal solution found.
## The objective value is: 4.345243e-05
round(solution(sol), 3)
##  MMM.Adjusted  AXP.Adjusted AAPL.Adjusted   BA.Adjusted  CAT.Adjusted 
##         0.000         0.000         0.000         0.000         0.000 
##  CVX.Adjusted CSCO.Adjusted   KO.Adjusted DWDP.Adjusted  XOM.Adjusted 
##         0.051         0.000         0.212         0.000         0.069 
##   GS.Adjusted   HD.Adjusted  IBM.Adjusted INTC.Adjusted  JNJ.Adjusted 
##         0.000         0.068         0.000         0.000         0.076 
##  JPM.Adjusted  MCD.Adjusted  MRK.Adjusted MSFT.Adjusted  NKE.Adjusted 
##         0.040         0.067         0.000         0.000         0.000 
##  PFE.Adjusted   PG.Adjusted  TRV.Adjusted  UNH.Adjusted  UTX.Adjusted 
##         0.000         0.079         0.000         0.000         0.048 
##   VZ.Adjusted    V.Adjusted  WMT.Adjusted  WBA.Adjusted  DIS.Adjusted 
##         0.142         0.013         0.012         0.018         0.107

Example 3: Minimum variance portfolio with budget normalization and target return constraint

lp  <- OP(objective  =  markowitz_objective(djia2018)$objective,
          constraints = rbind(budget_constraint(djia2018),
                              reward_constraint(djia2018, rhs = 0.001)))

To perfom the optimization, the ROI_solve() function is called:

(sol <- ROI_solve(lp, solver = "quadprog"))
## Optimal solution found.
## The objective value is: 8.484419e-05
round(solution(sol), 3)
##  MMM.Adjusted  AXP.Adjusted AAPL.Adjusted   BA.Adjusted  CAT.Adjusted 
##         0.000         0.000         0.000         0.000         0.000 
##  CVX.Adjusted CSCO.Adjusted   KO.Adjusted DWDP.Adjusted  XOM.Adjusted 
##         0.102         0.000         0.135         0.071         0.000 
##   GS.Adjusted   HD.Adjusted  IBM.Adjusted INTC.Adjusted  JNJ.Adjusted 
##         0.225         0.219         0.157         0.000         0.000 
##  JPM.Adjusted  MCD.Adjusted  MRK.Adjusted MSFT.Adjusted  NKE.Adjusted 
##         0.000         0.020         0.000         0.000         0.000 
##  PFE.Adjusted   PG.Adjusted  TRV.Adjusted  UNH.Adjusted  UTX.Adjusted 
##         0.000         0.000         0.071         0.000         0.000 
##   VZ.Adjusted    V.Adjusted  WMT.Adjusted  WBA.Adjusted  DIS.Adjusted 
##         0.000         0.000         0.000         0.000         0.000

Example 4: Minimize mean absolute deviation portfolio

The minimum MAD portfolio is given by:

lp  <- OP(objective  =  mad_objective(djia2018)$objective,
          constraints = rbind(mad_objective(djia2018)$constraint,
                              budget_constraint(djia2018),
                              use.names = TRUE))

(Note that the bounds can be omitted as this is the default in ROI). To perfom the optimization, the ROI_solve() function is called:

(sol <- ROI_solve(lp, solver = "glpk"))
## Optimal solution found.
## The objective value is: 4.587062e-03

The optimal weights are:

round(solution(sol)[seq_len(NCOL(djia2018))], 3)
##  MMM.Adjusted  AXP.Adjusted AAPL.Adjusted   BA.Adjusted  CAT.Adjusted 
##         0.000         0.026         0.027         0.000         0.000 
##  CVX.Adjusted CSCO.Adjusted   KO.Adjusted DWDP.Adjusted  XOM.Adjusted 
##         0.053         0.000         0.185         0.000         0.033 
##   GS.Adjusted   HD.Adjusted  IBM.Adjusted INTC.Adjusted  JNJ.Adjusted 
##         0.000         0.132         0.000         0.000         0.000 
##  JPM.Adjusted  MCD.Adjusted  MRK.Adjusted MSFT.Adjusted  NKE.Adjusted 
##         0.039         0.054         0.033         0.000         0.000 
##  PFE.Adjusted   PG.Adjusted  TRV.Adjusted  UNH.Adjusted  UTX.Adjusted 
##         0.018         0.113         0.000         0.000         0.000 
##   VZ.Adjusted    V.Adjusted  WMT.Adjusted  WBA.Adjusted  DIS.Adjusted 
##         0.088         0.000         0.032         0.041         0.126

Example 5: The minimum lower semi-variance portfolio

The minimum lower semi-variance portfolio is given by:

tmp <- downside_var_objective(djia2018)
lp  <- OP(objective  = tmp$objective,
          constraints = rbind(tmp$constraint,
                              budget_constraint(djia2018), 
                              use.names = TRUE))

To perfom the optimization, the ROI_solve() function is called:

(sol <- ROI_solve(lp, solver = "quadprog"))
## Optimal solution found.
## The objective value is: 1.804649e-05

The optimal weights are:

round(solution(sol)[seq_len(NCOL(djia2018))], 3)
##  MMM.Adjusted  AXP.Adjusted AAPL.Adjusted   BA.Adjusted  CAT.Adjusted 
##         0.000         0.000         0.000         0.000         0.000 
##  CVX.Adjusted CSCO.Adjusted   KO.Adjusted DWDP.Adjusted  XOM.Adjusted 
##         0.044         0.000         0.100         0.000         0.078 
##   GS.Adjusted   HD.Adjusted  IBM.Adjusted INTC.Adjusted  JNJ.Adjusted 
##         0.000         0.046         0.019         0.000         0.013 
##  JPM.Adjusted  MCD.Adjusted  MRK.Adjusted MSFT.Adjusted  NKE.Adjusted 
##         0.086         0.064         0.006         0.000         0.014 
##  PFE.Adjusted   PG.Adjusted  TRV.Adjusted  UNH.Adjusted  UTX.Adjusted 
##         0.079         0.048         0.000         0.000         0.074 
##   VZ.Adjusted    V.Adjusted  WMT.Adjusted  WBA.Adjusted  DIS.Adjusted 
##         0.144         0.000         0.000         0.031         0.153

Example 6: Minimum lower semi-absolute deviation portfolio with target return constraint

The portfolio with minimum lower semi-absolute deviation and minimum 0.001 target return is given by:

tmp <- downside_mad_objective(djia2018)
lp  <- OP(objective  = tmp$objective,
          constraints = rbind(tmp$constraint,
                              budget_constraint(djia2018),
                              reward_constraint(djia2018, rhs = 0.001),
                              use.names = TRUE))

To perfom the optimization, the ROI_solve() function is called:

(sol <- ROI_solve(lp, solver = "glpk"))
## Optimal solution found.
## The objective value is: 3.290369e-03
round(solution(sol)[seq_len(NCOL(djia2018))], 3)
##  MMM.Adjusted  AXP.Adjusted AAPL.Adjusted   BA.Adjusted  CAT.Adjusted 
##         0.000         0.000         0.000         0.000         0.000 
##  CVX.Adjusted CSCO.Adjusted   KO.Adjusted DWDP.Adjusted  XOM.Adjusted 
##         0.068         0.000         0.093         0.006         0.029 
##   GS.Adjusted   HD.Adjusted  IBM.Adjusted INTC.Adjusted  JNJ.Adjusted 
##         0.283         0.283         0.130         0.000         0.000 
##  JPM.Adjusted  MCD.Adjusted  MRK.Adjusted MSFT.Adjusted  NKE.Adjusted 
##         0.000         0.040         0.000         0.000         0.000 
##  PFE.Adjusted   PG.Adjusted  TRV.Adjusted  UNH.Adjusted  UTX.Adjusted 
##         0.000         0.000         0.047         0.000         0.000 
##   VZ.Adjusted    V.Adjusted  WMT.Adjusted  WBA.Adjusted  DIS.Adjusted 
##         0.021         0.000         0.000         0.000         0.000

Example 7: Minimum 95% & 99% conditional value at risk:

The portfolio with minimum 95% conditional value at risk can be obtained as follows:

tmp <- cvar_objective(djia2018, alpha = 0.95)
lp  <- OP(objective  =  tmp$objective,
          constraints = rbind(tmp$constraint,
                              budget_constraint(djia2018),
                              use.names = TRUE),
          bounds = tmp$bounds)

To perfom the optimization, the ROI_solve() function is called:

(sol <- ROI_solve(lp, solver = "glpk"))
## Optimal solution found.
## The objective value is: 1.119712e-02
round(solution(sol)[seq_len(NCOL(djia2018))], 3)
##  MMM.Adjusted  AXP.Adjusted AAPL.Adjusted   BA.Adjusted  CAT.Adjusted 
##         0.000         0.000         0.000         0.000         0.000 
##  CVX.Adjusted CSCO.Adjusted   KO.Adjusted DWDP.Adjusted  XOM.Adjusted 
##         0.000         0.000         0.054         0.000         0.163 
##   GS.Adjusted   HD.Adjusted  IBM.Adjusted INTC.Adjusted  JNJ.Adjusted 
##         0.000         0.030         0.155         0.000         0.000 
##  JPM.Adjusted  MCD.Adjusted  MRK.Adjusted MSFT.Adjusted  NKE.Adjusted 
##         0.006         0.000         0.000         0.000         0.000 
##  PFE.Adjusted   PG.Adjusted  TRV.Adjusted  UNH.Adjusted  UTX.Adjusted 
##         0.046         0.070         0.000         0.000         0.074 
##   VZ.Adjusted    V.Adjusted  WMT.Adjusted  WBA.Adjusted  DIS.Adjusted 
##         0.211         0.000         0.000         0.018         0.172

The solution for 99% CVaR is given by:

tmp99 <- cvar_objective(djia2018, alpha = 0.99)
lp99  <- OP(objective  =  tmp99$objective,
          constraints = rbind(tmp99$constraint,
                              budget_constraint(djia2018),
                              use.names = TRUE),
          bounds = tmp$bounds)
sol99 <- ROI_solve(lp99, solver = "glpk")
round(solution(sol99)[seq_len(NCOL(djia2018))], 3)
##  MMM.Adjusted  AXP.Adjusted AAPL.Adjusted   BA.Adjusted  CAT.Adjusted 
##         0.000         0.000         0.000         0.000         0.000 
##  CVX.Adjusted CSCO.Adjusted   KO.Adjusted DWDP.Adjusted  XOM.Adjusted 
##         0.035         0.000         0.005         0.000         0.217 
##   GS.Adjusted   HD.Adjusted  IBM.Adjusted INTC.Adjusted  JNJ.Adjusted 
##         0.000         0.069         0.133         0.000         0.000 
##  JPM.Adjusted  MCD.Adjusted  MRK.Adjusted MSFT.Adjusted  NKE.Adjusted 
##         0.016         0.000         0.000         0.000         0.000 
##  PFE.Adjusted   PG.Adjusted  TRV.Adjusted  UNH.Adjusted  UTX.Adjusted 
##         0.071         0.100         0.000         0.000         0.000 
##   VZ.Adjusted    V.Adjusted  WMT.Adjusted  WBA.Adjusted  DIS.Adjusted 
##         0.242         0.000         0.000         0.007         0.105

The 95% or the 99% VaR are a by-product of the optimization problem and can be extracted from the solution

VaR <- c("0.95" = sol$solution["gamma"], 
         "0.99" = sol99$solution["gamma"]) 
VaR
## 0.95.gamma 0.99.gamma 
## 0.01064806 0.01169682

Example 8: The minimax portfolio

tmp <- minimax_young_objective(djia2018)
lp  <- OP(objective  =  tmp$objective,
          constraints = rbind(tmp$constraint,
                              budget_constraint(djia2018),
                              use.names = TRUE),
          bounds =  tmp$bounds,
          maximum =  T)

To perfom the optimization, the ROI_solve() function is called:

(sol <- ROI_solve(lp, solver = "glpk"))
## Optimal solution found.
## The objective value is: -1.169682e-02
round(solution(sol)[seq_len(NCOL(djia2018))], 3)
##  MMM.Adjusted  AXP.Adjusted AAPL.Adjusted   BA.Adjusted  CAT.Adjusted 
##         0.000         0.000         0.000         0.000         0.000 
##  CVX.Adjusted CSCO.Adjusted   KO.Adjusted DWDP.Adjusted  XOM.Adjusted 
##         0.035         0.000         0.005         0.000         0.217 
##   GS.Adjusted   HD.Adjusted  IBM.Adjusted INTC.Adjusted  JNJ.Adjusted 
##         0.000         0.069         0.133         0.000         0.000 
##  JPM.Adjusted  MCD.Adjusted  MRK.Adjusted MSFT.Adjusted  NKE.Adjusted 
##         0.016         0.000         0.000         0.000         0.000 
##  PFE.Adjusted   PG.Adjusted  TRV.Adjusted  UNH.Adjusted  UTX.Adjusted 
##         0.071         0.100         0.000         0.000         0.000 
##   VZ.Adjusted    V.Adjusted  WMT.Adjusted  WBA.Adjusted  DIS.Adjusted 
##         0.242         0.000         0.000         0.007         0.105

Example 9: Maximize quadratic utility where short selling is allowed but the weights should not be less than -1.

tmp <- quadratic_utility_objective(djia2018)
lp  <- OP(objective  =  tmp$objective,
          constraints = budget_constraint(djia2018),
          bounds = V_bound(li = seq_len(NCOL(djia2018)), 
                           lb = rep(-1, NCOL(djia2018))),
          maximum = T)
(sol <- ROI_solve(lp, solver = "quadprog"))
## Optimal solution found.
## The objective value is: 3.018209e-02
round(solution(sol)[seq_len(NCOL(djia2018))], 3)
##  MMM.Adjusted  AXP.Adjusted AAPL.Adjusted   BA.Adjusted  CAT.Adjusted 
##        -1.000        -1.000        -0.173        -0.221        -1.000 
##  CVX.Adjusted CSCO.Adjusted   KO.Adjusted DWDP.Adjusted  XOM.Adjusted 
##         1.700        -1.000        -1.000         3.359        -1.000 
##   GS.Adjusted   HD.Adjusted  IBM.Adjusted INTC.Adjusted  JNJ.Adjusted 
##         6.592         3.561         3.906         1.708        -1.000 
##  JPM.Adjusted  MCD.Adjusted  MRK.Adjusted MSFT.Adjusted  NKE.Adjusted 
##        -1.000        -1.000        -1.000        -1.000        -1.000 
##  PFE.Adjusted   PG.Adjusted  TRV.Adjusted  UNH.Adjusted  UTX.Adjusted 
##        -1.000        -1.000         1.567        -1.000        -1.000 
##   VZ.Adjusted    V.Adjusted  WMT.Adjusted  WBA.Adjusted  DIS.Adjusted 
##        -1.000        -1.000        -1.000        -1.000        -1.000

Example 10: Maximize Sharpe ratio with shortselling

N <- NCOL(djia2018)
tmp <-sharpe_objective(djia2018)

lp <- OP(maximum = T)
objective(lp) <- tmp$objective

When imposing constraints on the portfolio choice vector \(x\) using the Sharpe objective, one must keep in mind that \(x^* = y^* / \kappa^*\), where \(y\) and \(\kappa\) are the values optimized. Hence, the constraint \(x_i \geq -1\) must be transformed to: \[ y_i/\kappa \geq -1 \Rightarrow y_i \geq -\kappa \Rightarrow y_i + \kappa \geq 0 \]:

mat <- cbind(diag(N),  1)
shortsell_constraint <- L_constraint(mat, dir = rep(">=", N),
                                     rhs = rep(0, N))

constraints(lp) <- rbind(tmp$constraint,
                         shortsell_constraint)

Moreover, the lower bounds for the auxiliary variables must be set to -Inf.

bounds(lp) <- V_bound(li = seq_len(N +1), 
                      lb = c(rep(-Inf, N), 0))
(sol <- ROI_solve(lp, solver = "quadprog"))
## Optimal solution found.
## The objective value is: -6.311826e+00

The optimal solution is given by \(x^* = y^* / \kappa^*\):

sol_sharpe <- solution(sol)
x_opt <- round(sol_sharpe[1:30]/sol_sharpe["kappa_sharpe"], 3)
names(x_opt) <-colnames(djia2018)
x_opt
##  MMM.Adjusted  AXP.Adjusted AAPL.Adjusted   BA.Adjusted  CAT.Adjusted 
##        -0.970        -1.000         0.831        -0.037        -0.519 
##  CVX.Adjusted CSCO.Adjusted   KO.Adjusted DWDP.Adjusted  XOM.Adjusted 
##         1.352        -1.000         2.387         1.414        -0.399 
##   GS.Adjusted   HD.Adjusted  IBM.Adjusted INTC.Adjusted  JNJ.Adjusted 
##         2.925         1.828         1.861         0.993        -1.000 
##  JPM.Adjusted  MCD.Adjusted  MRK.Adjusted MSFT.Adjusted  NKE.Adjusted 
##        -1.000        -0.243        -1.000        -0.560        -0.364 
##  PFE.Adjusted   PG.Adjusted  TRV.Adjusted  UNH.Adjusted  UTX.Adjusted 
##        -1.000        -0.916         1.861        -1.000        -0.303 
##   VZ.Adjusted    V.Adjusted  WMT.Adjusted  WBA.Adjusted  DIS.Adjusted 
##        -1.000        -0.622        -0.185        -0.617        -0.715

Example 11: Maximize Omega

tmp <- omega_objective(djia2018)
lp  <- OP(objective  =  tmp$objective,
          constraints = tmp$constraint,
          maximum = T)
(sol <- ROI_solve(lp, solver = "glpk"))
## Optimal solution found.
## The objective value is: 4.407022e-03

The optimal solution is given by \(x^* = y^* / z^*\).

sol_omega <- solution(sol)
x_opt <- round(sol_omega[1:30]/sol_omega["z_omega"], 3)
names(x_opt) <-colnames(djia2018)
x_opt
##  MMM.Adjusted  AXP.Adjusted AAPL.Adjusted   BA.Adjusted  CAT.Adjusted 
##         0.000         0.000         0.000         0.000         0.000 
##  CVX.Adjusted CSCO.Adjusted   KO.Adjusted DWDP.Adjusted  XOM.Adjusted 
##         0.000         0.000         0.000         0.175         0.000 
##   GS.Adjusted   HD.Adjusted  IBM.Adjusted INTC.Adjusted  JNJ.Adjusted 
##         0.511         0.166         0.149         0.000         0.000 
##  JPM.Adjusted  MCD.Adjusted  MRK.Adjusted MSFT.Adjusted  NKE.Adjusted 
##         0.000         0.000         0.000         0.000         0.000 
##  PFE.Adjusted   PG.Adjusted  TRV.Adjusted  UNH.Adjusted  UTX.Adjusted 
##         0.000         0.000         0.000         0.000         0.000 
##   VZ.Adjusted    V.Adjusted  WMT.Adjusted  WBA.Adjusted  DIS.Adjusted 
##         0.000         0.000         0.000         0.000         0.000

Example 12: Maximize expected returns wrt to cvar, turnover, short selling and group constraints

\[\begin{eqnarray*} \max_{x\in\mathbb{R}^N}& \hat\mu^\top x \\ s.t.& \mathrm{CVaR}(x, 0.95) \leq 0.02 \\ & x_i \geq -1\\ & x_2 + x_{10} + x_{20}\leq 0.5 \\ & \sum_{i=1}^N \vert x_i - x_{0i} \vert \leq 0.5 \end{eqnarray*}\]

Given that the constraints have more (auxiliary) variables than the objective, first the constraints are defined:

lp <- OP(maximum = T)

constraints(lp) <- rbind(
  budget_constraint(djia2018),
  group_constraint(djia2018, index = c(2, 10, 20), dir = "<=", rhs = 0.5),
  turnover_constraint(djia2018, dir = "<=", rhs = 0.5),
  cvar_constraint(djia2018, alpha = 0.95, rhs = 0.02),
  use.names = TRUE)

Next, the vector of the linear objective is filled with zeros to match the dimension of the constraints:

obj <- c(terms(reward_objective(djia2018)$objective)$L)
objective(lp) <- c(obj, double(ncol(constraints(lp)) - length(obj)))
(sol <- ROI_solve(lp, solver = "glpk"))
## Optimal solution found.
## The objective value is: 8.153523e-04
solution(sol)[1:30]
##  MMM.Adjusted  AXP.Adjusted AAPL.Adjusted   BA.Adjusted  CAT.Adjusted 
##    0.03333333    0.03333333    0.03333333    0.03333333    0.03333333 
##  CVX.Adjusted CSCO.Adjusted   KO.Adjusted DWDP.Adjusted  XOM.Adjusted 
##    0.03333333    0.03333333    0.03333333    0.07703204    0.03333333 
##   GS.Adjusted   HD.Adjusted  IBM.Adjusted INTC.Adjusted  JNJ.Adjusted 
##    0.23963463    0.03333333    0.03333333    0.03333333    0.00000000 
##  JPM.Adjusted  MCD.Adjusted  MRK.Adjusted MSFT.Adjusted  NKE.Adjusted 
##    0.03333333    0.03333333    0.00000000    0.03333333    0.03333333 
##  PFE.Adjusted   PG.Adjusted  TRV.Adjusted  UNH.Adjusted  UTX.Adjusted 
##    0.00000000    0.00000000    0.03333333    0.01666667    0.03333333 
##   VZ.Adjusted    V.Adjusted  WMT.Adjusted  WBA.Adjusted  DIS.Adjusted 
##    0.00000000    0.03333333    0.00000000    0.00000000    0.03333333

Example 13: Maximize expected return wrt variance constraint

\[\begin{eqnarray*} \max_{x\in\mathbb{R}^N}&& \hat\mu^\top x \\ s.t.&& x^{\top}Qx \leq 0.1 \\ &&\sum_{i=1}^N x_i =1\\ && x_i \geq 0 \end{eqnarray*}\]
p <- OP(maximum = T)

constraints(p) <- rbind(
  markowitz_constraint(djia2018, dir = "<=", rhs = 0.5^2),
  budget_constraint(djia2018),
  use.names = TRUE)

objective(p) <- reward_objective(djia2018)$objective 

When using the alabama solver starting values need to be chosen. A good strategy is to use several starting values and compare the results.

mstart <- lapply(1:10, function(x) runif(length(objective(p))))
solus <- lapply(mstart, function(s) ROI_solve(p, solver = "alabama", start = s))
best_solution <- which.max(sapply(solus, solution, type = "objval"))  
round(solution(solus[[best_solution]]), 3)
##  MMM.Adjusted  AXP.Adjusted AAPL.Adjusted   BA.Adjusted  CAT.Adjusted 
##         0.000         0.000         0.000         0.000         0.000 
##  CVX.Adjusted CSCO.Adjusted   KO.Adjusted DWDP.Adjusted  XOM.Adjusted 
##         0.000         0.000         0.000         0.000         0.000 
##   GS.Adjusted   HD.Adjusted  IBM.Adjusted INTC.Adjusted  JNJ.Adjusted 
##         0.999         0.000         0.000         0.000         0.000 
##  JPM.Adjusted  MCD.Adjusted  MRK.Adjusted MSFT.Adjusted  NKE.Adjusted 
##         0.000         0.000         0.000         0.000         0.000 
##  PFE.Adjusted   PG.Adjusted  TRV.Adjusted  UNH.Adjusted  UTX.Adjusted 
##         0.000         0.000         0.000         0.000         0.000 
##   VZ.Adjusted    V.Adjusted  WMT.Adjusted  WBA.Adjusted  DIS.Adjusted 
##         0.000         0.000         0.000         0.000         0.000

Example 14: Minimize 99% CVaR wrt cardinality constraints

\[\begin{eqnarray*} \min_{x\in \mathbb{R}^N, z\in \mathbb{R}^S, \gamma\in \mathbb{R}}& &\gamma + \frac{1}{(1-0.99)} \sum_{s = 1}^S p_s z_s\\ s.t.&&z_s \geq - r_s^\top x - \gamma,\\ &&z_s \geq 0,\\ && \sum_{i=1}^N u_i\leq 6,\\ && x_i - u_i \leq 0. \end{eqnarray*}\]

where \(u\in \mathbb{R}^N\)are binary auxiliary variables.

tmp <- cvar_objective(djia2018, 0.99)
lp <- OP()

constraints(lp) <- rbind(
  tmp$constraint,
  budget_constraint(djia2018),
  cardinality_constraint(djia2018, dir = "<=", rhs = 6),
  use.names = TRUE)
## Variable types for z_card_aux must be set to binary.
obj <-c((tmp$objective)$L)
objective(lp) <- c(obj, double(NCOL(constraints(lp)) - length(obj)))

types(lp) <- rep("C",  NCOL(constraints(lp) ))
types(lp)[grep("z_card_aux", constraints(lp)$names)] <- "B"

(sol <- ROI_solve(lp, solver = "glpk"))
## Optimal solution found.
## The objective value is: 1.225243e-02
round(solution(sol)[1:30], 3)
##  MMM.Adjusted  AXP.Adjusted AAPL.Adjusted   BA.Adjusted  CAT.Adjusted 
##         0.000         0.000         0.000         0.000         0.000 
##  CVX.Adjusted CSCO.Adjusted   KO.Adjusted DWDP.Adjusted  XOM.Adjusted 
##         0.000         0.000         0.000         0.000         0.284 
##   GS.Adjusted   HD.Adjusted  IBM.Adjusted INTC.Adjusted  JNJ.Adjusted 
##         0.000         0.106         0.154         0.000         0.000 
##  JPM.Adjusted  MCD.Adjusted  MRK.Adjusted MSFT.Adjusted  NKE.Adjusted 
##         0.000         0.000         0.000         0.000         0.000 
##  PFE.Adjusted   PG.Adjusted  TRV.Adjusted  UNH.Adjusted  UTX.Adjusted 
##         0.000         0.123         0.000         0.000         0.000 
##   VZ.Adjusted    V.Adjusted  WMT.Adjusted  WBA.Adjusted  DIS.Adjusted 
##         0.227         0.000         0.000         0.000         0.107