knitr::opts_chunk$set( echo = TRUE, collapse = TRUE, comment = "#>" ) Sys.setlocale("LC_TIME", "English")
\hypersetup{linkcolor=black} \tableofcontents \newpage
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.
Load the necessary packages.
library(PortfolioAnalytics) library(CVXR) library(data.table) library(xts) library(PCRA)
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$|
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}
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.
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.
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
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
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)
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.
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")
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.
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
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.
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)
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
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.
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))
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
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).
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))
opt_eqs <- optimize.portfolio(ret_edhec, pspec_eqs, optimize_method = "CVXR") opt_eqs
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.
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)
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)
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)
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}
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)))
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)))
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.
# 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.
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.
# 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.
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}
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.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.