knitr::opts_chunk$set(
  echo = TRUE,
  collapse = TRUE,
  comment = "#>"
)
Sys.setlocale("LC_TIME", "English")

\hypersetup{linkcolor=black} \tableofcontents \newpage

1 Introduction

CVXR is an R package that provides an object-oriented modeling language for convex optimization, including the Second-Order Cone Optimization(SOCopt) required to minimize Expected Quadratic Shortfall(EQS) problem, which is not supported by other solvers in PortfolioAnalytics. Hence, CVXR is a significant extension of PortfolioAnalytics.

The purpose of this vignette is to demonstrate examples of optimization problems that can be solved in PortfolioAnalytics with CVXR and its many supported solvers. The problem types covered include not only Linear Programming(LP), Quadratic Programming(QP) but also Second-Order Cone Programming(SOCP). Multiple solvers supported by CVXR can be selected according to optimization types. For example, SCS and ECOS can completely cover the types of problems that ROI can deal with, such as mean-variance and ES problem. In order to better understand the functionality and use of PortfolioAnalytics, users are recommended to read the Vignette Introduction to PortfolioAnalytics first.

The R code demo_cvxrPortfolioAnalytics.R in the PortfolioAnalytics demo folder, reproduces the results in this Vignette.

2 Getting Started

2.1 Load Packages

Load the necessary packages.

library(PortfolioAnalytics)
library(CVXR)
library(data.table)
library(xts)
library(PCRA)

2.2 Solvers

The website https://cvxr.rbind.io/ shows that CVXR currently supports the use of 9 solvers, some of which are commercial (CBC, CPLEX, GUROBI, MOSEK)^[MOSEK ,GUROBI, and CPLEX require licenses to use, but free academic licenses are available for all three.] and the others are open source(GLPK, GLPK_MI, OSQP, SCS, ECOS).

The PortfolioAnalytics package provides the following two main wrapper functions for constrainted optimization of portfolios using a wide variety of methods:

optimize.portfolio(R =, portfolio =, optimize_method =)

optimize.portfolio.rebalancing(R =, portfolio =, optimize_method =, rebalance_on =, training_period =, rolling_window =)

Different solvers support different types of portfolio optimization problems, which should be specified by the argument optimize_method. The optimize_method=c("CVXR", {CVXRsolver}) argument of the function optimize.portfolio and optimize.portfolio.rebalancing allows the user to specify the solver to use with CVXR. If the argument is optimize_method="CVXR", the default solver for LP and QP type portfolio optimization problems, such as maximum mean return and minimum variance portfolio optimization is OSQP, and the default solver for SOCP type portfolio optimizations, such as "robust portfolio optimization" to control for alpha uncertainty, and Expected Quadratic Shortfall (EQS) portfolio optimization, is ECOS.

| Solver | LP | QP | SOCP | | :-- | :----: | :----: | :----: | |CBC|$\checkmark$| | | |GLPK|$\checkmark$| | | |GLPK_MI|$\checkmark$| | | |OSQP|$\colorbox{gray}{\checkmark}$|$\colorbox{gray}{\checkmark}$| | |SCS|$\checkmark$|$\checkmark$|$\checkmark$| |ECOS|$\checkmark$|$\checkmark$|$\colorbox{gray}{\checkmark}$| |CPLEX|$\checkmark$|$\checkmark$|$\checkmark$| |GUROBI|$\checkmark$|$\checkmark$|$\checkmark$| |MOSEK|$\checkmark$|$\checkmark$|$\checkmark$|

2.3 Data

The edhec data set from the PerformanceAnalytics package is used as example data for examples in Section 3 to Section 8. The edhec data set is an xts object that contains monthly returns for 13 hedge fund style indexes from 1997-01 to 2019-11. We use the edhec data of the last 5 years as the example data to show how to use the code.

data(edhec)
# Use edhec for a returns object
ret_edhec <- tail(edhec, 60)
colnames(ret_edhec) <- c("CA", "CTAG", "DS", "EM", "EMN", "ED", "FIA", 
                       "GM", "LSE", "MA", "RV", "SS", "FF")
print(head(ret_edhec, 5))
# Get a character vector of the asset names
fund_edhec <- colnames(ret_edhec)

tsPlotMP is a function of R package PCRA which can plot time series for the return data.

tsPlotMP(ret_edhec, layout = c(2, 7), main = "Time Series Plot of Edhec Return")

\begin{center} Fig 2.1 \end{center}

2.4 Optimization Problems

In this Vignette, all mean vectors and covariance matrices in the optimization formula will use standard sample based estimates. All optimization problems treated will use linear constraints unless stated otherwise. There will be one equality constraint, i.e., the full-investment constraint, and one or more inequality constraints such as the long-only and box constraints. More comprehensive constraint types can be found in the vignette Ross (2018) Introduction to PortfolioAnalytics.

This vignette will be organized by objective type and provide some visual examples.

3 Maximizing Mean Return

The objective to maximize mean return is a linear problem of the form: $$\max\limits_{w} \quad \boldsymbol{\mu}'\boldsymbol{w}$$ \begin{equation} \begin{aligned} s.t. \quad & A\boldsymbol{w} \ge b\ & B\boldsymbol{w} = c \end{aligned} \end{equation}

Where $\boldsymbol{\mu}$ is the estimated asset returns mean vector and $\boldsymbol{w}$ is the vector of portfolio weights.

3.1 Portfolio Object

The first step in setting up a model is to create the portfolio object, which contains the form of the constraints, and the objective specifications. In the following we create full-investment and box constraints specifications, and a maximum return objective specification.

# Create portfolio object
pspec_maxret <- portfolio.spec(assets = fund_edhec)
# Add constraints to the portfolio object
pspec_maxret <- add.constraint(pspec_maxret, type = "full_investment")
pspec_maxret <- add.constraint(portfolio = pspec_maxret, type = "box",
                              min = rep(0.02, 13),
                              max = c(rep(0.15, 8), rep(0.1, 5)))
# Add objective to the portfolio object
pspec_maxret <- add.objective(portfolio = pspec_maxret, 
                             type = "return", name = "mean")
pspec_maxret

3.2 Optimization

The next step is to run the optimization. Note that optimize_method = c("CVXR", {CVXRsolver}) should be specified in the function optimize.portfolio to use CVXR solvers for the optimization, or use the default solver by giving optimize_method = "CVXR". For maximizing mean return problem, which is a linear programming, the default solver is OSQP.

# Run the optimization with default solver
opt_maxret <- optimize.portfolio(R = ret_edhec, portfolio = pspec_maxret, 
                                 optimize_method = "CVXR", trace = TRUE)
opt_maxret
opt_maxret$solver

# Run the optimization with specific solver
opt_maxret_glpk <- optimize.portfolio(R = ret_edhec, portfolio = pspec_maxret, 
                                 optimize_method = c("CVXR", "GLPK"), trace = TRUE)
opt_maxret_glpk$solver

3.3 Backtesting

An out of sample backtest is run with optimize.portfolio.rebalancing. In this example, an initial training period of 36 months is used and the portfolio is rebalanced quarterly.

bt_maxret <- optimize.portfolio.rebalancing(R = ret_edhec, portfolio = pspec_maxret,
                                            optimize_method = "CVXR",
                                            rebalance_on = "quarters",
                                            training_period = 36)

The call to optimize.portfolio.rebalancing returns the bt_maxret object which is a list containing the optimal weights and objective measure at each rebalance period.

class(bt_maxret)
names(bt_maxret)

4 Minimizing Variance

The objective to minimize variance is a quadratic problem of the form: $$\min\limits_{w} \quad \boldsymbol{w}'\Sigma \boldsymbol{w}$$ subject to portfolio managers' desired constraints, where $\Sigma$ is the estimated covariance matrix of asset returns and $\boldsymbol{w}$ is the set of portfolio weights. It is a quadratic problem.

4.1 Global Minimum Variance Portfolio

4.1.1 Portfolio Object

In this example, the only constraint specified is the full investment constraint, therefore the optimization problem is solving for the Global Minimum Variance (GMV) portfolio.

# Create portfolio object
pspec_gmv <- portfolio.spec(assets = fund_edhec)
# Add full-investment constraint
pspec_gmv <- add.constraint(pspec_gmv, type = "full_investment")
# Add objective of minimizing variance
pspec_gmv <- add.objective(portfolio = pspec_gmv, type = "risk", name = "var")

4.1.2 Optimization

opt_gmv <- optimize.portfolio(ret_edhec, pspec_gmv, optimize_method = "CVXR")
opt_gmv

As this example illustrates, a global minimum variance portfolio with only a full-investment constraint can have short positions.

4.2 Long-Only and Group Constrained Minimum Variance Portfolio

Various linear inequality constraint, such as box constraints, group constraints and a target mean return constraint, can be used with GMV portfolio construction. Here we demonstrate the case of linearly constrained minimum variance portfolio.

# portfolio object
pspec_mv <- add.constraint(pspec_gmv, type = "long_only")
pspec_mv <- add.constraint(pspec_mv, type = "group",
                            groups = list(groupA=1,
                                          groupB=c(2:12),
                                          groupC=13),
                            group_min = c(0, 0.05, 0.05),
                            group_max = c(0.4, 0.8, 0.5))
pspec_mv <- add.constraint(pspec_mv, type = "return", return_target = 0.003)
pspec_mv

# optimization
opt_mv <- optimize.portfolio(ret_edhec, pspec_mv, optimize_method = "CVXR")
opt_mv

Compared to Section 4.1, the result shows that the more constraints, the larger the optimal value may be. But this example is closer to the real situation. The optimal weights show that the first group constraint is not binding, but the second one is binding with FIA plus MA at the upper bound of 0.8, and the third group constraint is binding at the lower bound with FF.

The use of an alternative to the CVXR default solver will get the same result to many significant digits. In this example we use optimize_method = c("CVXR", "ECOS"), since OSQP is the default solver, and get the very similar results.

opt_mv_ecos <- optimize.portfolio(ret_edhec, pspec_mv, optimize_method = c("CVXR", "ECOS"))
opt_mv_ecos

opt_mv$solver
opt_mv_ecos$solver

5 Maximizing Quadratic Utility

Next we demonstrate the classical quadratic utility form of Markowitz's mean-variance optimal portfolios, where the quadratic utility function is $\rm QU(\boldsymbol{w}) = \mu_p - \lambda \sigma_p^2 = \boldsymbol{\mu'w}-\lambda \boldsymbol{w'}\Sigma \boldsymbol{w}$: $$\max\limits_{w} \quad \boldsymbol{\mu}'\boldsymbol{w} - \lambda\boldsymbol{w}'\Sigma\boldsymbol{w}$$ $$s.t. \quad A\boldsymbol{w} \ge b$$ Here $\boldsymbol{\mu}$ is the vector of estimated mean asset returns, $0 \le \lambda < \inf$ is the risk aversion parameter, $\Sigma$ is the estimated covariance matrix of asset returns, and $\boldsymbol{w}$ is the vector of weights. Quadratic utility maximizes return while penalizing variance risk. The risk aversion parameter $\lambda$ controls how much portfolio variance is penalized, and when $\lambda=0$ it becomes a maximum mean return problem of Section 3, and as $\lambda \rightarrow \inf$, it becomes the global minimum variance problem of Section 4.

5.1 Portfolio Object

In this case the objectives of the portfolio should be both return and risk, and for this example we will use a risk aversion parameter $\lambda$ to be 20 by setting risk_aversion = 20.

pspec_mvo <- portfolio.spec(assets = fund_edhec)
pspec_mvo <- add.constraint(pspec_mvo, type = "full_investment")
pspec_mvo <- add.constraint(pspec_mvo, type = "long_only")
# Add objectives
pspec_mvo <- add.objective(portfolio = pspec_mvo, type = "return", name = "mean")
pspec_mvo <- add.objective(portfolio = pspec_mvo, type = "risk", name = "var",
                           risk_aversion = 20)

5.2 Optimization

The optimization result opt_mvo shows the call, optimal weights, and the objective measure. Objective measure contains quadratic utility, mean return and standard deviation.

opt_mvo <- optimize.portfolio(ret_edhec, pspec_mvo, optimize_method = "CVXR")
opt_mvo

6 Minimizing Expected Shortfall

Expected Shortfall(ES) is also called Conditional Value-at-Risk(CVaR) and Expected Tail Loss(ETL). The ES of a portfolio is \begin{equation} \begin{aligned} ES_\gamma(r_P) = ES_\gamma(\boldsymbol{w}) &= -E(r_P|r_P \le q_\gamma(\boldsymbol{w}))\ &=-E(\boldsymbol{w'r}|\boldsymbol{w'r} \le q_\gamma(\boldsymbol{w})) \end{aligned} \end{equation}

where $r_P$ is a random return of a portfolio $P$, and $\boldsymbol{r}$ is the loss return which is negative, and $q_\gamma$ is $\gamma$-quantile and $\gamma$ is usually a "tail" probability such as 0.01, 0.05, in which case ES is a tail risk measure. But one could also choose $\gamma=0.25$ or $\gamma=0.5$, in which case ES is just a "downside" risk measure, and if $\gamma>0.5$, the problem will take $1-\gamma$ as the tail probability.

It was shown by @rockafellar2000optimization that the optimal minimum ES portfolio is the result of the minimization: $$\min \limits_{\boldsymbol{w}} ES_\gamma(\boldsymbol{w}) = \min \limits_{\boldsymbol{w}, t} F_\gamma(\boldsymbol{w}, t)$$ where $$F_\gamma(\boldsymbol{w},t)=-t+\frac{1}{\gamma} \int[t-\boldsymbol{w'r}]^+ \cdot f(\boldsymbol{r})d\boldsymbol{r}$$ by replacing $q_\gamma$ with the free variable $t$, and with the discrete data the formula is: $$\hat{F}\gamma(\boldsymbol{w}, t) = -t + \frac{1}{n \cdot \gamma} \sum{i=1}^n[t-\boldsymbol{w'r_i}]^+$$ The positive part function, $[t-\boldsymbol{w'r_i}]^+$, can easily be converted to a collection of linear constraints, hence, the minimization of ES is equivalent to solving a linear programming problem.

The ES minimization problem is $$\min\limits_{\boldsymbol{w}, t} \quad -t + \gamma^{-1}E(t-\boldsymbol{w'r_i})^+$$ where $0<\gamma<1$ is the tail probability value, and $t$ is the value from which shortfalls are measured in the optimal solution. Many authors also use $p$ or $\alpha$ as the quantile, e.g., in @rockafellar2000optimization and other vignettes of PortfolioAnalytics, and use $\eta$ as the risk measure variable, e.g., in @krokhmal2007higher.

6.1 Portfolio Object

The default probability is $\gamma = 5\%$. Specific probability could be given by arguments.

pspec_es <- portfolio.spec(assets = fund_edhec)
pspec_es <- add.constraint(pspec_es, type = "full_investment")
pspec_es <- add.constraint(pspec_es, type = "long_only")
# Add objective of minimizing ES by using the default gamma
pspec_es <- add.objective(portfolio = pspec_es, type = "risk", name = "ES")
# Add objective of minimizing ES by using the specific gamma=0.1
pspec_es_1 <- add.objective(portfolio = pspec_es, type = "risk", name = "ES",
                          arguments = list(p=0.1))

6.2 Optimization

Notice that if optimize_method is not declared, the default solver is DEoptim. But we highly recommend using ROI solvers or CVXR solvers, because they provide accurate theoretical solutions for minES problem.

# GMES with default gamma=0.05
opt_es <- optimize.portfolio(ret_edhec, pspec_es, optimize_method = "CVXR")
opt_es
# GMES with specific gamma=0.1
opt_es_1 <- optimize.portfolio(ret_edhec, pspec_es_1, optimize_method = "CVXR")
opt_es_1

7 Minimizing Expected Quadratic Shortfall

Expected Quadratic Shortfall(EQS) is also called Second-Moment Coherent Risk Measure(SMCR) in some situations, for example, in @krokhmal2007higher. The objective to minimize EQS is in the form of: $$\min\limits_{\boldsymbol{w}, t} \quad -t + \gamma^{-1}||(t-\boldsymbol{w'r_i})^+||_2$$ where $\gamma$ is the tail probability and $0<\gamma<1$, $t$ is the value from which quadratic shortfalls are measured in the optimal solution. The default probability is $\gamma = 5\%$. Minimizing EQS could be incorporated into a convex problem as a second-order cone constraints, and PortfolioAnalytics uses ECOS in CVXR as the default solver for Second-Order Cone Optimization(SOCopt).

7.1 Portfolio Object

The default probability is $\gamma = 5\%$. Specified probability could be given by arguments.

pspec_eqs <- portfolio.spec(assets = fund_edhec)
pspec_eqs <- add.constraint(pspec_eqs, type = "full_investment")
pspec_eqs <- add.constraint(pspec_eqs, type = "long_only")
# Add objective of minimizing EQS
pspec_eqs <- add.objective(portfolio = pspec_eqs, type = "risk", name = "EQS",
                          arguments = list(p=0.05))

7.2 Optimization

opt_eqs <- optimize.portfolio(ret_edhec, pspec_eqs, optimize_method = "CVXR")
opt_eqs

8 Maximizing Mean Return Per Unit Risk

There are three basic types of risk measures: variance or standard deviation, ES and EQS. The problem of maximizing mean return per unit risk can be solved in a clever way by minimizing risk with a target return constraint, as is described below. For all three of these types of problems, both return and risk objectives should be used in PortfolioAnalytics. Then for each of these three optimization problems an appropriate argument needs to be given to the optimize.portfolio to specify the type of problem, as we describe below.

8.1 Maximum Sharpe Ratio Portfolios

The Sharpe Ratio of a random return $r_P$ of a portfolio $P$ is defined as: $$\frac{E(r_P) - r_f}{\sqrt{Var(r_P)}}_.$$ The problem of maximizing the Sharpe Ratio can be formulated as a quadratic problem with a budget normalization constraint. It is shown in @cornuejols2018optimization, that this optimization problem can be formulated as the equivalent optimization: \begin{equation} \begin{aligned} \mathop{minimize}\limits_{w} \quad w'\Sigma w\ s.t. \quad (\hat{\mu} - r_f\textbf{1})^Tw &= 1\ \textbf{1}^Tw &= \kappa\ \kappa &> 0 \end{aligned} \end{equation} which has a solution$(w^,\kappa^)$ with $k^ \ne 0$, and the maximized Sharpe ratio given by $\tilde{w}^ = w^/\kappa^$.

When creating the portfolio, the argument maxSR = TRUE should be specified in the function optimize.portfolio to distinguish from the mean-variance optimization. NOTE: The default argument is maxSR = FALSE since the default action for dealing with both mean and var/StdDev objectives is to maximize quadratic utility.

# Create portfolio object
pspec_sr <- portfolio.spec(assets = fund_edhec)
## Add constraints of maximizing Sharpe Ratio
pspec_sr <- add.constraint(pspec_sr, type = "full_investment")
pspec_sr <- add.constraint(pspec_sr, type = "long_only")
## Add objectives of maximizing Sharpe Ratio
pspec_sr <- add.objective(pspec_sr, type = "return", name = "mean")
pspec_sr <- add.objective(pspec_sr, type = "risk", name = "var")

# Optimization
optimize.portfolio(ret_edhec, pspec_sr, optimize_method = "CVXR", maxSR = TRUE)

8.2 Maximum ES ratio Portfolios

The ES ratio(ESratio), which is also called STARR in PortfolioAnalytics, is defined as: $$\frac{E(r_P) - r_f}{ES_{\gamma}(r_P)}$$ Similar to maximizing Sharpe Ratio, the problem maximizing the ES ratio can be formulated as a minimizing ES problem with a budget normalization constraint.

When creating the portfolio, both return and ES objectives should be given. The default $\gamma=0.05$, and it can be specified by arguments. When solving the problem, the default argument ESratio = TRUE in the function optimize.portfolio specifies the problem type. We note that this argument is equivalent to maxSTARR = TRUE, which is used in other vignettes. If one of these two arguments is specified as FALSE, the action will be to minimize ES ignoring the return objective.

# Create portfolio object
pspec_ESratio <- portfolio.spec(assets = fund_edhec)
## Add constraints of maximizing return per unit ES
pspec_ESratio <- add.constraint(pspec_ESratio, type = "full_investment")
pspec_ESratio <- add.constraint(pspec_ESratio, type = "long_only")
## Add objectives of maximizing return per unit ES
pspec_ESratio <- add.objective(pspec_ESratio, type = "return", name = "mean")
pspec_ESratio <- add.objective(pspec_ESratio, type = "risk", name = "ES",
                               arguments = list(p=0.05))

# Optimization
optimize.portfolio(ret_edhec, pspec_ESratio, optimize_method = "CVXR", ESratio = TRUE)

8.3 Maximum EQS ratio Portfolios

The EQS ratio of a random return $r_P$ of a portfolio $P$ is defined as: $$\frac{E(r_P) - r_f}{EQS_{\gamma}(r_P)}$$ Similar to maximizing Sharpe Ratio, the problem maximizing EQS ratio could be formulated as a minimizing EQS problem with a budget normalization constraint.

When creating the portfolio, both return and EQS objectives should be given. The argument EQSratio = is used to specify the problem type and the default value is EQSratio = TRUE. If EQSratio = FALSE, the action will be to minimize EQS ignoring the return objective. The default $\gamma=0.05$, and it can be specified by arguments.

# Create portfolio object
pspec_EQSratio <- portfolio.spec(assets = fund_edhec)
## Add constraints of maximizing return per unit EQS
pspec_EQSratio <- add.constraint(pspec_EQSratio, type = "full_investment")
pspec_EQSratio <- add.constraint(pspec_EQSratio, type = "long_only")
## Add objectives of maximizing return per unit EQS
pspec_EQSratio <- add.objective(pspec_EQSratio, type = "return", name = "mean")
pspec_EQSratio <- add.objective(pspec_EQSratio, type = "risk", name = "EQS",
                                arguments = list(p=0.05))

# Optimization
optimize.portfolio(ret_edhec, pspec_EQSratio, optimize_method = "CVXR", EQSratio = TRUE)

9 Performance of minEQS, minES and minVar Portfolios

CVXR solvers provide the Second-Order Cone Optimization (SOCopt) capability required to compute minEQS portfolios. In this Section we use this capability to compute a minEQS portfolio and compare its performance that of minES and minVar portfolios. For these computations, we use a CRSP$\textsuperscript{\textregistered}$ daily returns data set to generate these optimal portfolios and show their performance by plotting cumulative gross returns and efficient frontiers.^[CRSP$\textsuperscript{\textregistered}$ stands for the Center for Research in Security Prices, LLC.] Specifically, the data set consists of 30 smallcap stocks from 1993-01 to 2015-12, contained in the stocksCRSPdaily data object in the PCRA package available on CRAN. The 30 smallcaps are the 30 of the 106 smallcap stocks with the largest market capitalization.

The entire back-testing part is relatively slow, and takes most of the time in this Vignette. To give a sense of typical times, I provide the computing times for my MacBook Air with M2, 8-Core CPU, 8-Core GPU and 16-core Neural Engine. The backtesting in Sections 9.1 and Section 9.2 takes about 3 minutes each to run, and the average running time for generating efficient frontiers in Section 9.3 is about 30 seconds.

# Get daily returns of the 30 smallcap stocks
library(PCRA)
stocksCRSPdaily <- getPCRAData(dataset = "stocksCRSPdaily")

smallcapTS <- selectCRSPandSPGMI(
  periodicity = "daily",
  stockItems = c("Date", "TickerLast", "CapGroupLast", "Return"),
  factorItems = NULL,
  subsetType = "CapGroupLast",
  subsetValues = "SmallCap",
  outputType = "xts")

# find top 30 small cap stocks based on the market capitalization
smallcapDT <- factorsSPGMI[CapGroupLast == "SmallCap"]
scSize <- smallcapDT[, mean(LogMktCap), by = "TickerLast"]
names(scSize)[2] <- "Size"
scSize <- scSize[order(scSize$Size, decreasing = TRUE),]
sc30largest <- scSize[,TickerLast][1:30]

# daily return of top 30 stocks
retD_CRSP <- smallcapTS[ , sc30largest]
print(head(retD_CRSP, 3))

In the following part, we only show the time series of monthly returns of 10 CRSP$\textsuperscript{\textregistered}$ stocks in the last five years, but you can use this code to view the time series of any subset of stocks in any frequency and any time period.

# monthly return of top 30 stocks in last 5 years
ep <- endpoints(retD_CRSP, on= "months", k=1)
prod1 <- function(x){apply(x+1, 2, prod)}
retM_CRSP <- period.apply(retD_CRSP, INDEX = ep, FUN = prod1) - 1
retM_CRSP_5 <- tail(retM_CRSP, 60)

# time series plot of 10 stocks
tsPlotMP(retM_CRSP_5[, 1:10])

\begin{center} Fig 9.1 \end{center}

9.1 Backtesting with GMV, GMES, GMEQS portfolios

In this example, we use daily return of all the CRSP$\textsuperscript{\textregistered}$ 30 stocks to generate a comparative backtesting among Global Minimum Variance (GMV), Global Minimum ES (GMES) and Global Minimum EQS (GMEQS) portfolio. The strategy is to rebalance the portfolio at the end of each month with a rolling window of 500 days, and the performance of backtesting could be shown as a plot of cumulative returns and a plot of drawdown.

# test run time
start_time1 <- Sys.time()
# Generate GMV, GMES and GMEQS portfolios
pspec_sc <- portfolio.spec(assets = sc30largest)
pspec_sc <- add.constraint(pspec_sc, type = "full_investment")
pspec_sc <- add.constraint(pspec_sc, type = "long_only")

pspec_GMV <- add.objective(pspec_sc, type = "risk", name = "var")
pspec_GMES <- add.objective(pspec_sc, type = "risk", name = "ES")
pspec_GMEQS <- add.objective(pspec_sc, type = "risk", name = "EQS")

# Optimize Portfolio at Monthly Rebalancing and 500-Day Training
bt.GMV <- optimize.portfolio.rebalancing(retD_CRSP, pspec_GMV,
                                           optimize_method = "CVXR",
                                           rebalance_on = "months",
                                           training_period = 30,
                                           rolling_window = 500)
bt.ES <- optimize.portfolio.rebalancing(retD_CRSP, pspec_GMES,
                                           optimize_method = "CVXR",
                                           rebalance_on = "months",
                                           training_period = 30,
                                           rolling_window = 500)
bt.EQS <- optimize.portfolio.rebalancing(retD_CRSP, pspec_GMEQS,
                                           optimize_method = "CVXR",
                                           rebalance_on = "months",
                                           training_period = 30,
                                           rolling_window = 500)

# Extract time series of portfolio weights
wts.GMV <- extractWeights(bt.GMV)
wts.GMV <- wts.GMV[complete.cases(wts.GMV),]

wts.ES <- extractWeights(bt.ES)
wts.ES <- wts.ES[complete.cases(wts.ES),]

wts.EQS <- extractWeights(bt.EQS)
wts.EQS <- wts.EQS[complete.cases(wts.EQS),]

# Compute cumulative returns of three portfolios
GMV <- Return.rebalancing(retM_CRSP, wts.GMV)
ES <- Return.rebalancing(retM_CRSP, wts.ES)
EQS <- Return.rebalancing(retM_CRSP, wts.EQS)

# Combine GMV, ES and EQS portfolio cumulative returns
ret.comb <- na.omit(merge(GMV, ES, EQS, all=F))
names(ret.comb) <- c("GMV", "GMES", "GMEQS")

backtest.plot(ret.comb, colorSet = c("black", "darkblue", "darkgreen"), 
              ltySet = c(3, 2, 1))

\begin{center} Fig 9.2 \end{center}

end_time1 <- Sys.time()
runningtime1 <- end_time1 - start_time1
cat("The run time for Figure 9.2 is", format(round(runningtime1, 2)))

9.2 Backtesting with SR, ESratio, EQSratio portfolios

In this example, we use daily return of all the CRSP$\textsuperscript{\textregistered}$ 30 stocks to generate a comparative backtesting among Maximum Sharpe Ratio (SR), Maximum ES Ratio (ESratio) and Maximum EQS Ratio (EQSratio) portfolio. The strategy is to rebalance the portfolio at the end of each month with a rolling window of 500 days, and the performance of backtesting could be shown as a plot of cumulative returns and a plot of drawdown.

# test run time
start_time2 <- Sys.time()
# Generate GMV, GMES and GMEQS portfolios
pspec_sc_ratio <- add.objective(pspec_sc, type = "return", name = "mean")
pspec_Sr <- add.objective(pspec_sc_ratio, type = "risk", name = "var")
pspec_ESr <- add.objective(pspec_sc_ratio, type = "risk", name = "ES")
pspec_EQSr <- add.objective(pspec_sc_ratio, type = "risk", name = "EQS")

# Optimize Portfolio at Monthly Rebalancing and 500-Day Training
bt.Sr <- optimize.portfolio.rebalancing(retD_CRSP, pspec_Sr, maxSR = TRUE,
                                           optimize_method = "CVXR",
                                           rebalance_on = "months",
                                           training_period = 30,
                                           rolling_window = 500)
bt.ESr <- optimize.portfolio.rebalancing(retD_CRSP, pspec_ESr,
                                           optimize_method = "CVXR",
                                           rebalance_on = "months",
                                           training_period = 30,
                                           rolling_window = 500)
bt.EQSr <- optimize.portfolio.rebalancing(retD_CRSP, pspec_EQSr,
                                           optimize_method = "CVXR",
                                           rebalance_on = "months",
                                           training_period = 30,
                                           rolling_window = 500)

# Extract time series of portfolio weights
wts.Sr <- extractWeights(bt.Sr)
wts.Sr <- wts.Sr[complete.cases(wts.Sr),]

wts.ESr <- extractWeights(bt.ESr)
wts.ESr <- wts.ESr[complete.cases(wts.ESr),]

wts.EQSr <- extractWeights(bt.EQSr)
wts.EQSr <- wts.EQSr[complete.cases(wts.EQSr),]

# Compute cumulative returns of three portfolios
Sr <- Return.rebalancing(retM_CRSP, wts.Sr, rebalance_on = "months")
ESr <- Return.rebalancing(retM_CRSP, wts.ESr, rebalance_on = "months")
EQSr <- Return.rebalancing(retM_CRSP, wts.EQSr, rebalance_on = "months")

# Combine Sr, ESr and EQSr portfolio cumulative returns
ret.comb <- na.omit(merge(Sr, ESr, EQSr, all=F))
names(ret.comb) <- c("Sharpe ratio", "ES ratio", "EQS ratio")

backtest.plot(ret.comb, colorSet = c("black", "darkblue", "darkgreen"), 
              ltySet = c(3, 2, 1))

\begin{center} Fig 9.3 \end{center}

end_time2 <- Sys.time()
runningtime2 <- end_time2 - start_time2
cat("The run time for Figure 9.3 is", format(round(runningtime2, 2)))

9.3 Efficient Frontier

We generate efficient frontiers with mean-StdDev, mean-ES and mean-EQS portfolios using 30 stocks from CRSP$\textsuperscript{\textregistered}$ data set. For illustrative purposes, we use the last five years Considering that the data may show different properties over a long period of time, we only use the monthly return in the last 5 years to generate efficient frontiers, that is from 2011-01 to 2015-12 and defined in Section 2.3 as retM_CRSP_5. We can use create.EfficientFrontier to calculate the mean value and risk value for the frontier, then use chart.EfficientFrontier to draw the frontier.

9.3.1 Mean-StdDev Efficient Frontier

# mean-var efficient frontier
meanvar.ef <- create.EfficientFrontier(R = retM_CRSP_5, portfolio = pspec_sc, 
                                       type = "mean-StdDev")
meanvar.ef
chart.EfficientFrontier(meanvar.ef, match.col = "StdDev", type = "l",
                        chart.assets = FALSE, main = "Mean-StdDev Efficient Frontier",
                        RAR.text = "Sharpe ratio", pch = 1)

\begin{center} Fig 9.4 \end{center} Here we compute the maximum Sharpe ratio on the efficient frontier using the following set of efficient frontier mean and standard deviation values:

meanvar.ef$frontier[, 1:2]
sr = meanvar.ef$frontier[, 1]/meanvar.ef$frontier[, 2]
cat("maximum Sharpe ratio:", max(sr))
cat("mean of the maximum SR portfolio:", meanvar.ef$frontier[, 1][sr == max(sr)])
cat("StdDev of the maximum SR portfolio:", meanvar.ef$frontier[, 2][sr == max(sr)])

Note that we have introduced a method of finding the maximum Sharpe ratio portfolio in Section 8.1, which may have the values a little different from the estimated maximum Sharpe ratio calculated using the discrete efficient frontier values.

We now use that method to compute the maximum Sharpe ratio portfolio directly:

# Mean-StdDev Efficient Frontier
pspec_MV <- add.objective(pspec_sc, type = "risk", name = "var")
pspec_MV <- add.objective(portfolio = pspec_MV, type = "return", name = "mean")
opt_MV <- optimize.portfolio(retM_CRSP_5, pspec_MV, optimize_method = "CVXR", 
                              maxSR = TRUE, trace = TRUE)
opt_MV

The above direct computation of the maximum Sharpe ratio portolio's Sharpe ratio = 0.4161, agrees to three significant digits with the value obtained by searching the discrete set of efficient frontier mean and standard deviation values.

chart.EfficientFrontier(opt_MV, match.col = "StdDev", chart.assets = FALSE, 
                        main = "Mean-StdDev Efficient Frontier",
                        RAR.text = "Sharpe Ratio", pch = 1, xlim = c(0, 0.1))

\begin{center} Fig 9.5 \end{center} The theoretical maximum Sharpe ratio portfolio is very close to the result generated by the efficient frontier, and the Sharpe ratio value is same to 3 significant digits.

With different constraint types, we can create mean-StdDev efficient frontiers for multiple portfolios and overlay the plots.

pspec_sc_init <- portfolio.spec(assets = sc30largest)
pspec_sc_init <- add.constraint(pspec_sc_init, type = "full_investment")

# Portfolio with long-only constraints
pspec_sc_lo <- add.constraint(portfolio = pspec_sc_init, type = "long_only")

# Portfolio with long-only box constraints
pspec_sc_lobox <- add.constraint(portfolio = pspec_sc_init, type = "box", 
                                 min = 0.02, max = 0.1)

# Portfolio with long-short box constraints
pspec_sc_lsbox <- add.constraint(portfolio = pspec_sc_init, type = "box", 
                                 min = -0.1, max = 0.1)

# Combine the portfolios into a list
portf_list <- combine.portfolios(list(pspec_sc_lo, pspec_sc_lobox, pspec_sc_lsbox))

# Plot the efficient frontier overlay of the portfolios with varying constraints
legend_labels <- c("Long Only", "Long Only Box", "Long Short Box")
chart.EfficientFrontierOverlay(R = retM_CRSP_5, portfolio_list = portf_list, 
                               type = "mean-StdDev", match.col = "StdDev", 
                               legend.loc = "bottomright", chart.assets = FALSE,
                               legend.labels = legend_labels, cex.legend = 1,
                               labels.assets = FALSE, lwd = c(3,3,3),
                               col = c("black", "dark red", "dark green"),
                               main = "Overlay Mean-StdDev Efficient Frontiers",
                               xlim = c(0.03, 0.11), ylim = c(0.005, 0.035))

\begin{center} Fig 9.6 \end{center} The plot clearly shows that the long-short box constrained portfolio has the best performance, though it also requires shorting which may not be possible for many real-world portfolios.

9.3.2 Mean-ES Efficient Frontier

Generate the mean-ES efficient frontier:

# Mean-ES Efficient Frontier
meanetl.ef <- create.EfficientFrontier(R = retM_CRSP_5, portfolio = pspec_sc, 
                                       type = "mean-ES")
chart.EfficientFrontier(meanetl.ef, match.col = "ES", type = "l",
                        chart.assets = FALSE, main = "Mean-ES Efficient Frontier",
                        RAR.text = "ES ratio", pch = 1)

\begin{center} Fig 9.7 \end{center} Generate multiple mean-ES efficient frontiers and overlay the plots.

legend_labels <- c("Long Only ES (p=0.05)", 
                   "Long Only Box ES (p=0.05)", "Long Short Box ES (p=0.05)")
chart.EfficientFrontierOverlay(R = retM_CRSP_5, portfolio_list = portf_list, 
                               type = "mean-ES", match.col = "ES",
                               legend.loc = "bottomright", chart.assets = FALSE,
                               legend.labels = legend_labels, cex.legend = 1,
                               labels.assets = FALSE, lwd = c(3,3,3),
                               col = c("black", "dark red", "dark green"),
                               main = "Overlay Mean-ES Efficient Frontiers",
                               xlim = c(0.03, 0.17), ylim = c(0.005, 0.035))

\begin{center} Fig 9.8 \end{center} For each mean-ES efficient frontier, the left endpoint is the global minimum ES and its corresponding mean return, the right endpoint is the maximum mean return and and its corresponding ES, that is, the minimum ES when the maximum mean return is the target return.

With the same target return, the portfolio under the long-short box constraints can have a minimum ES value. It performs best, but shorting may not be possible for many real-world portfolios. The mean-ES efficient frontier of the long-only portfolio covers the widest range, that is because it can reach a larger return range than the long-box portfolio, but it is not as easy to reach the target return as the long-short portfolio, it needs to take more risks.

We could also notice that the ordering of the three constrained portfolios is the same of for the minVar efficient frontiers.

Instead of generating efficient frontiers with different constraint types, we can also generate mean-ES efficient frontiers with different tail probability $\gamma$.

# Create long-only ES portfolios with different tail probabilities
ES_05 <- add.objective(portfolio = pspec_sc_lo, type = "risk", name = "ES", 
                       arguments = list(p=0.05))

ES_10 <- add.objective(portfolio = pspec_sc_lo, type = "risk", name = "ES", 
                       arguments = list(p=0.1))

ES_15 <- add.objective(portfolio = pspec_sc_lo, type = "risk", name = "ES", 
                       arguments = list(p=0.15))

# Combine the portfolios into a list
portf_ES_list <- combine.portfolios(list(ES_05, ES_10, ES_15))

# Plot the efficient frontier overlay of the portfolios with varying tail probabilities
legend_ES_labels <- c("ES (p=0.05)", "ES (p=0.1)", "ES (p=0.15)")
chart.EfficientFrontierOverlay(R = retM_CRSP_5, portfolio_list = portf_ES_list, 
                               type = "mean-ES", match.col = "ES", 
                               legend.loc = "bottomright", chart.assets = FALSE,
                               legend.labels = legend_ES_labels, cex.legend = 1,
                               labels.assets = FALSE, lwd = c(3,3,3),
                               col = c("black", "dark red", "dark green"),
                               main = "Overlay Mean-ES Efficient Frontiers",
                               xlim = c(0.035, 0.165), ylim = c(0.005, 0.03))

\begin{center} Fig 9.9 \end{center} ES portfolio with a larger tail probability will have better performance, but this is only because larger tail probability gives a wider ES calculation range, so the expected ES value is smaller. Actually, these three lines cannot be compared because their x-lables are ES with different tail probabilities.

This plot mainly shows that the same portfolio will have different ES values under different tail probabilities. For example, the portfolio that reaches the maximum return, which is the right endpoint of each mean-ES efficient frontier. They are the same portfolio, but the ES values are different. It informs us the importance of tail probability for calculating ES. Different tail probabilities only focus on risks in different degrees, rather than changes in risk levels.

9.3.3 Mean-EQS Efficient Frontier

# Mean-EQS Efficient Frontier
meaneqs.ef <- create.EfficientFrontier(R = retM_CRSP_5, portfolio = pspec_sc, 
                                       type = "mean-EQS")
chart.EfficientFrontier(meaneqs.ef, match.col = "EQS", type = "l",
                        chart.assets = FALSE, main = "Mean-EQS Efficient Frontier",
                        RAR.text = "EQS ratio", pch = 1)

\begin{center} Fig 9.10 \end{center} Mean-EQS efficient frontier is more like a piecewise function rather than a smooth curve.

9.3.4 Efficient Frontiers Comparison among minStd, minES and minEQS portfolios

In order to calculate the value of StdDev, ES and EQS of a portfolio, we provide a new function extract_risk(R =, w =). Asset returns R and the portfolio weights w should be specified when using the function, ES_alpha and EQS_alpha are alternative arguments and the default values are 0.05.

# usage example: minStd Portfolio
minstd_port <- add.objective(pspec_sc, type = "risk", name = "StdDev")
minstd_w <- optimize.portfolio(retM_CRSP_5, minstd_port, optimize_method = "CVXR")$weight

# risk values with default alpha = 0.05
extract_risk(retM_CRSP_5, minstd_w)

# risk values with specific alpha
extract_risk(retM_CRSP_5, minstd_w, ES_alpha = 0.1, EQS_alpha = 0.1)

For comparing a specific risk type in minStd, minES and minEQS portfolios and generating multiple efficient frontiers, we provide a new function

chart.EfficientFrontierCompare(R =, portfolio =, risk_type =, match.col =, guideline =)

where risk_type is the risk to be compared, and match.col is the vector of the portfolios that participate in the comparison. If there are only two frontiers in the comparison, the default of the argument guideline is TRUE.

# example 1: Compare StdDev of minStd and minES portfolios with guideline
chart.EfficientFrontierCompare(R = retM_CRSP_5, portfolio = pspec_sc, risk_type = "StdDev",
                               match.col = c("StdDev", "ES"), lwd = c(2, 2))

\begin{center} Fig 9.11 \end{center}

# example 2: Compare ES of minStd, minES and minEQS portfolios without guideline
chart.EfficientFrontierCompare(R = retM_CRSP_5, portfolio = pspec_sc, risk_type = "ES",
                               match.col = c("StdDev", "ES", "EQS"), guideline = FALSE,
                               col = c(1,2,4), lty = c(1, 2, 4), lwd = c(2, 2, 2))

\begin{center} Fig 9.12 \end{center}

Acknowledgements

This vignette aims to show the extended functionality of PortfolioAnalytics. CVXR solvers and EQS objective in PortfolioAnalytics were implemented by me, Xinran Zhao, when I was a student in the University of Washington Department of Applied Mathematics Computational Finance and Risk Management M.S. degree program (MS-CFRM), and was supported by Google Summer of Code (GSoC 2022) in the summer of 2022 with mentors Doug Martin and Steve Murray.

In the following year, I continuously generated more functions and improved this vignette, not limited to the implementation of CVXR and EQS, but more focused on their applications, such as back-testing analysis and efficient frontiers analysis. This process is very long, and I am grateful to everyone who accompanies and inspires me.

Doug Martin provided a lot of guidance on the derivation, application, and subsequent research directions of EQS. He also provided feedback on the style and content of this vignette. I am very grateful for learning the habit of research in this project.

Steve Murray is an expert in optimization. I'm greatly benefited by his analyzing different CVXR solvers, mathematical derivation in second-order cone optimization and mean-EQS efficient frontier, and his assistance in R coding.

Yifu Kang, who is my classmate in MS-CFRM and was funded by another PortfolioAnalytics GSoC 2022 project, expanded PortfolioAnalytics in custom robust covariance matrix. We have collaborated extensively on improving existing functions and interacting with each other's functions.

Finally, this project could not have been completed without the help of Brian Peterson, who is an Author and the Maintainer of the PortfolioAnalytics package on CRAN.

Reference



braverock/PortfolioAnalytics documentation built on April 18, 2024, 4:09 a.m.