risk_budget_optimization.Rnw

\documentclass[a4paper]{article}
\usepackage[round]{natbib}
\usepackage{bm}
\usepackage{verbatim}
\usepackage[latin1]{inputenc}

\usepackage{url}

\let\proglang=\textsf
\newcommand{\pkg}[1]{{\fontseries{b}\selectfont #1}}
\newcommand{\R}[1]{{\fontseries{b}\selectfont #1}}
\newcommand{\email}[1]{\href{mailto:#1}{\normalfont\texttt{#1}}}
\newcommand{\E}{\mathsf{E}}
\newcommand{\VAR}{\mathsf{VAR}}
\newcommand{\COV}{\mathsf{COV}}
\newcommand{\Prob}{\mathsf{P}}

\renewcommand{\topfraction}{0.85}
\renewcommand{\textfraction}{0.1}
\renewcommand{\baselinestretch}{1.5}
\setlength{\textwidth}{15cm} \setlength{\textheight}{22cm} \topmargin-1cm \evensidemargin0.5cm \oddsidemargin0.5cm

\usepackage[latin1]{inputenc}
% or whatever

\usepackage{lmodern}
\usepackage[T1]{fontenc}
% Or whatever. Note that the encoding and the font should match. If T1
% does not look nice, try deleting the line with the fontenc.

% \VignetteIndexEntry{Portfolio Optimization with CVaR budgets in PortfolioAnalytics}

\begin{document}
\SweaveOpts{concordance=TRUE}

\title{Portfolio Optimization with CVaR budgets in PortfolioAnalytics}
\author{Kris Boudt, Peter Carl and Brian Peterson }
\date{June 1, 2010}

\maketitle
\tableofcontents


\bigskip

\section{General information}

Risk budgets are a central tool to estimate and manage the portfolio risk allocation. They decompose total portfolio risk into the risk contribution of each position.  \citet{BoudtCarlPeterson2010}  propose several portfolio allocation strategies that use an appropriate transformation of the portfolio Conditional Value at Risk (CVaR) budget as an objective or constraint in the portfolio optimization problem.  This document explains how risk allocation optimized portfolios can be obtained under general constraints in the \verb"PortfolioAnalytics" package of \citet{PortfolioAnalytics}. 

\verb"PortfolioAnalytics" is designed to provide numerical solutions for portfolio problems with complex constraints and objective sets comprised of any R function. It can e.g.~construct portfolios that minimize a risk objective with (possibly non-linear) per-asset constraints on returns and drawdowns \citep{CarlPetersonBoudt2010}. The generality of possible constraints and objectives is a distinctive characteristic of the package with respect to RMetrics \verb"fPortfolio" of \citet{fPortfolioBook}. For standard Markowitz optimization problems, use of \verb"fPortfolio" rather than \verb"PortfolioAnalytics" is recommended.

\verb"PortfolioAnalytics"  solves the following type of problem
\begin{equation}  \min_w g(w) \ \  s.t. \ \ 
\left\{ \begin{array}{l} h_1(w)\leq 0 \\ \vdots \\  h_q(w)\leq 0. \end{array} \right. \label{optimproblem}\end{equation} \verb"PortfolioAnalytics"   first merges the objective function and constraints into a penalty augmented objective function  
\begin{equation} L(w) = g(w) + \mbox{penalty}\sum_{i=1}^q \lambda_i \max(h_i(w),0), \label{eq:constrainedobj} \end{equation}
where $\lambda_i$ is a multiplier to tune the relative importance of the constraints.  The default values of penalty and $\lambda_i$ (called \verb"multiplier" in \verb"PortfolioAnalytics") are 10000 and 1, respectively. 

The minimum of this function is found through the \emph{Differential Evolution} (DE) algorithm of \citet{StornPrice1997} and ported to R by \citet{MullenArdiaGilWindoverCline2009}. DE is known for remarkable performance regarding continuous numerical problems \citep{PriceStornLampinen2006}. It has recently been advocated for optimizing portfolios under non-convex settings by \citet{Ardia2010} and \citet{Yollin2009}, among others. We use the R implementation of DE in the \verb"DEoptim" package of \citet{DEoptim}. 

The latest version of the \verb"PortfolioAnalytics" package can be downloaded from R-forge through the following command:
\begin{verbatim}
install.packages("PortfolioAnalytics", repos="http://R-Forge.R-project.org")
\end{verbatim}

Its principal functions are:
\begin{itemize}
\item \verb"portfolio.spec(assets)": the portfolio specification starts with creating a \verb"portfolio" object with information about the assets. The first argument \verb"assets" is either a number indicating the number of portfolio assets or a vector holding the names of the assets. The \verb"portfolio" object is a list holding the constraints and objectives.

\item \verb"add.constraint(portfolio, type)": Constraints are added to the \verb"portfolio" object by the function \verb"add.constraint". Basic constraint types include leverage constraints that specify the sum of the weights have to be between \verb"min_sum" and \verb"max_sum" and box constraints where the asset weights have to be between \verb"min" and \verb"max".

\item \verb"add.objective(portfolio, type, name)": New objectives are added to the \verb"portfolio" objected with the function \verb"add.objective". Many common risk budget objectives and constraints are prespecified and can be identified by specifying the \verb"type" and \verb"name".

\item \verb"constrained_objective(w, R, portfolio)": given the portfolio weight and return data, it evaluates the penalty augmented objective function in (\ref{eq:constrainedobj}).

\item \verb"optimize.portfolio(R, portfolio)": this function returns the portfolio weight that solves the problem in (\ref{optimproblem}). {\it R} is the multivariate return series of the portfolio components. 

\item \verb"optimize.portfolio.rebalancing(R, portfolio, rebalance_on, trailing_periods)": this function solves the multiperiod optimization problem. It returns for each rebalancing period the optimal weights and allows the estimation sample to be either from inception or a moving window.

\end{itemize}

Next we illustrate these functions on monthly return data for bond, US equity, international equity and commodity indices, which are the first 4 series
in the dataset \verb"indexes". The first step is to load the package \verb"PortfolioAnalytics"  and the dataset. An important first note is that some of the functions (especially \verb" optimize.portfolio.rebalancing") requires the dataset to be a \verb"xts" object \citep{xts}.


<<echo=FALSE>>=
options(width=80)
@

<<echo=TRUE>>=
library(PortfolioAnalytics)
library(DEoptim)
library(robustbase)
data(indexes)
class(indexes)
indexes <- indexes[,1:4]
head(indexes,2)
tail(indexes,2)
@

In what follows, we first illustrate the construction of the penalty augmented objective function. Then we present the code for solving the optimization problem. 

\section{Setting of the objective function}

\subsection{Weight constraints}
  
<<echo=TRUE, tidy=FALSE>>=
# Create the portfolio specification object
Wcons <- portfolio.spec( assets = colnames(indexes) )
# Add box constraints
Wcons <- add.constraint( portfolio=Wcons, type='box', min = 0, max=1 )
# Add the full investment constraint that specifies the weights must sum to 1.
Wcons <- add.constraint( portfolio=Wcons, type="full_investment")
@

Given the weight constraints, we can call the value of the function to be minimized. We consider the case of no violation and a case of violation. By default, \verb"normalize=TRUE" which means that if the sum of weights exceeds \verb"max_sum", the weight vector is normalized by multiplying it with \verb"sum(weights)/max_sum" such that the weights evaluated in the objective function satisfy the \verb"max_sum" constraint.    
<<echo=TRUE, tidy=FALSE>>=
constrained_objective( w = rep(1/4,4) , R = indexes, portfolio = Wcons) 
constrained_objective( w = rep(1/3,4) , R = indexes, portfolio = Wcons) 
constrained_objective( w = rep(1/3,4) , R = indexes, portfolio = Wcons, 
                       normalize=FALSE) 
@

The latter value can be recalculated as penalty times the weight violation, that is: $10000 \times 1/3.$

\subsection{Minimum CVaR objective function}

Suppose now we want to find the portfolio that minimizes the 95\% portfolio CVaR subject to the weight constraints listed above.  

<<echo=TRUE, tidy=FALSE>>=
ObjSpec = add.objective( portfolio = Wcons , type="risk",name="CVaR",
                         arguments=list(p=0.95), enabled=TRUE)   
@

The value of the objective function is: 
<<echo=TRUE, tidy=FALSE>>=
constrained_objective( w = rep(1/4,4) , R = indexes, portfolio = ObjSpec) 
@

This is the CVaR of the equal-weight portfolio as computed by the function \verb"ES" in the \verb"PerformanceAnalytics" package of \citet{ Carl2007}
<<echo=TRUE, tidy=FALSE>>=
library(PerformanceAnalytics)
out<-ES(indexes, weights = rep(1/4,4),p=0.95, 
        portfolio_method="component")
out$MES
@
All arguments in the function \verb"ES" can be passed on through \verb"arguments". E.g. to reduce the impact of extremes on the portfolio results, it is recommended to winsorize the data using the option clean="boudt".

<<echo=TRUE, tidy=FALSE>>=
out<-ES(indexes, weights = rep(1/4,4),p=0.95, clean="boudt", 
        portfolio_method="component")
out$MES
@

For the formulation of the objective function, this implies setting: 
<<echo=TRUE, tidy=FALSE>>=
ObjSpec = add.objective( portfolio = Wcons , type="risk",name="CVaR",
                         arguments=list(p=0.95,clean="boudt"), enabled=TRUE)   
constrained_objective( w = rep(1/4,4) , R = indexes[,1:4] , portfolio = ObjSpec)
@

An additional argument that is not available for the moment in \verb"ES" is to estimate the conditional covariance matrix through the constant conditional correlation model of \citet{Bollerslev90}. 

For the formulation of the objective function, this implies setting: 
<<echo=TRUE, tidy=FALSE>>=
ObjSpec = add.objective( portfolio = Wcons , type="risk",name="CVaR",
                         arguments=list(p=0.95,clean="boudt"), 
                         enabled=TRUE, garch=TRUE)   
constrained_objective( w = rep(1/4,4) , R = indexes[,1:4] , portfolio = ObjSpec)
@

\subsection{Minimum CVaR concentration objective function}

Add the minimum 95\% CVaR concentration objective to the objective function: 
<<echo=TRUE, tidy=FALSE>>=
ObjSpec = add.objective( portfolio = Wcons , type="risk_budget_objective",
                         name="CVaR", arguments=list(p=0.95, clean="boudt"), 
                         min_concentration=TRUE, enabled=TRUE)   
@

The value of the objective function is: 
<<echo=TRUE, tidy=FALSE>>=
constrained_objective( w = rep(1/4,4) , R = indexes, 
                       portfolio = ObjSpec, trace=TRUE)
@

We can verify that this is effectively the largest CVaR contribution of that portfolio as follows:
<<echo=TRUE, tidy=FALSE>>=
ES(indexes[,1:4],weights = rep(1/4,4),p=0.95,clean="boudt", 
   portfolio_method="component")
@

\subsection{Risk allocation constraints}

We see that in the equal-weight portfolio, the international equities and commodities investment cause more than 30\% of total risk. We could specify as a constraint that no asset can contribute more than 30\% to total portfolio risk with the argument \verb"max_prisk=0.3". This involves the construction of the following objective function:

<<echo=TRUE, tidy=FALSE>>=
ObjSpec = add.objective( portfolio = Wcons , type="risk_budget_objective", 
                         name="CVaR", max_prisk = 0.3, 
                         arguments=list(p=0.95,clean="boudt"), enabled=TRUE) 
constrained_objective( w = rep(1/4,4) , R = indexes, portfolio = ObjSpec) 
@

This value corresponds to the penalty parameter which has by default the value of 10000 times the exceedances: $ 10000*(0.045775103+0.054685023)\approx 1004.601.$

\section{Optimization}

The penalty augmented objective function is minimized through Differential Evolution. Two parameters are crucial in tuning the optimization: \verb"search_size" and \verb"itermax". The optimization routine 
\begin{enumerate}
\item First creates the initial generation of \verb"NP = search_size/itermax" guesses for the optimal value of the parameter vector, using the \verb"random_portfolios" function generating random weights satisfying the weight constraints. 
\item Then DE evolves over this population of candidate solutions using alteration and selection operators in order to minimize the objective function. It restarts \verb"itermax" times.
\end{enumerate}
It is important that \verb"search_size/itermax" is high enough. It is generally recommended that this ratio is at least ten times the length of the weight vector. For more details on the use of DE strategy in portfolio allocation, we refer the
reader to \citet{Ardia2010}. 

\subsection{Minimum CVaR portfolio under an upper 40\% CVaR allocation constraint} 

The portfolio object and functions needed to obtain the minimum CVaR portfolio under an upper 40\% CVaR allocation objective are the following:
<<echo=TRUE, tidy=FALSE>>=
# Create the portfolio specification object
ObjSpec <- portfolio.spec(assets=colnames(indexes[,1:4]))
# Add box constraints
ObjSpec <- add.constraint(portfolio=ObjSpec, type='box', min = 0, max=1)
# Add the full investment constraint that specifies the weights must sum to 1.
ObjSpec <- add.constraint(portfolio=ObjSpec, type="weight_sum", 
                          min_sum=0.99, max_sum=1.01)
# Add objective to minimize CVaR
ObjSpec <- add.objective(portfolio=ObjSpec, type="risk", name="CVaR", 
                         arguments=list(p=0.95, clean="boudt"))
# Add objective for an upper 40% CVaR allocation
ObjSpec <- add.objective(portfolio=ObjSpec, type="risk_budget_objective",
                         name="CVaR", max_prisk=0.4, 
                         arguments=list(p=0.95, clean="boudt"))
@

After the call to these functions it starts to explore the feasible space iteratively and is shown in the output. Iterations are given as intermediate output and by default every iteration will be printed. We set \verb"traceDE=5" to print every 5 iterations and \verb"itermax=50" for a maximum of 50 iterations.

<<echo=TRUE, tidy=FALSE>>=
set.seed(1234)
out <- optimize.portfolio(R=indexes, portfolio=ObjSpec, 
                          optimize_method="DEoptim", search_size=2000,
                          traceDE=5, itermax=50, trace=TRUE)
print(out)
@


If \verb"trace=TRUE" in \verb"optimize.portfolio", additional output from the DEoptim solver is included in the \verb"out" object created by \verb"optimize.portfolio". The additional elements in the output are \verb"DEoptim_objective_results" and \verb"DEoutput". The \verb"DEoutput" element contains output from the function \verb"DEoptim". The \verb"DEoptim_objective_results" element contains the weights, value of the objective measures, and other data at each iteration.

<<echo=TRUE>>=
names(out)
# View the DEoptim_objective_results information at the last iteration
out$DEoptim_objective_results[[length(out$DEoptim_objective_results)]]

# Extract stats from the out object into a matrix
xtract <- extractStats(out)
dim(xtract)
head(xtract)
@

It can be seen from the charts that although US Bonds has a higher weight allocation, the percentage contribution to risk is the lowest of all four indexes.

<<tidy=FALSE, fig.height=4, fig.width=6>>=
plot.new()
chart.Weights(out)
@

<<tidy=FALSE, fig.height=4, fig.width=6>>=
plot.new()
chart.RiskBudget(out, risk.type="pct_contrib", col="blue", pch=18)
@


\subsection{Minimum CVaR concentration portfolio} 

The functions needed to obtain the minimum CVaR concentration portfolio are the following:

<<echo=TRUE, tidy=FALSE>>=
# Create the portfolio specification object
ObjSpec <- portfolio.spec(assets=colnames(indexes))
# Add box constraints
ObjSpec <- add.constraint(portfolio=ObjSpec, type='box', min = 0, max=1)
# Add the full investment constraint that specifies the weights must sum to 1.
ObjSpec <- add.constraint(portfolio=ObjSpec, type="weight_sum", 
                          min_sum=0.99, max_sum=1.01)
# Add objective for min CVaR concentration
ObjSpec <- add.objective(portfolio=ObjSpec, type="risk_budget_objective",
                         name="CVaR", arguments=list(p=0.95, clean="boudt"),
                         min_concentration=TRUE)

set.seed(1234)
out <- optimize.portfolio(R=indexes, portfolio=ObjSpec, 
                          optimize_method="DEoptim", search_size=5000,
                          itermax=50, traceDE=5, trace=TRUE)
@


This portfolio has the near equal risk contribution characteristic: 
<<echo=TRUE, tidy=FALSE>>=
print(out)

# Verify results with ES function
ES(indexes[,1:4], weights=out$weights, p=0.95, clean="boudt", 
   portfolio_method="component")
@

The 95\% CVaR percent contribution to risk is near equal for all four indexes. The neighbor portfolios can be plotted to view other near optimal portfolios. Alternatively, the contribution to risk in absolute terms can plotted by setting \verb"risk.type="absolute".

<<tidy=FALSE, fig.height=4.5, fig.width=6>>=
plot.new()
chart.RiskBudget(out, neighbors=25, risk.type="pct_contrib", 
                 col="blue", pch=18)
@


\subsection{Dynamic optimization}

Dynamic rebalancing of the risk budget optimized portfolio is possible through the function \verb"optimize.portfolio.rebalancing". Additional arguments are \verb"rebalance_on" which indicates the rebalancing frequency (years, quarters, months). The estimation is either done from inception (\verb"trailing_periods=0") or through moving window estimation, where each window has \verb"trailing_periods" observations. The minimum number of observations in the estimation sample is specified by \verb"training_period". Its default value is 36, which corresponds to three years for monthly data.

As an example, consider the minimum CVaR concentration portfolio, with estimation from inception and monthly rebalancing. Since we require a minimum estimation length of total number of observations -1, we can optimize the portfolio only for the last two months.  

<<echo=TRUE, tidy=FALSE>>=
set.seed(1234)
out <- optimize.portfolio.rebalancing(R=indexes, portfolio=ObjSpec, 
                                      optimize_method="DEoptim", search_size=5000,
                                      rebalance_on="quarters", 
                                      training_period=nrow(indexes)-12,
                                      traceDE=0)
@


The output of \verb"optimize.portfolio.rebalancing" in the \verb"opt_rebalancing" slot is a list of objects created by \verb"optimize.portfolio", one for each rebalancing period. 

<<echo=TRUE>>=
names(out)
names(out$opt_rebalancing[[1]])
out
@

The \verb"summary" method provides a brief output of the optimization result along with return and risk measures.
<<echo=TRUE>>=
opt.summary <- summary(out)
names(opt.summary)
opt.summary
@

The optimal weights for each rebalancing period can be extracted fron the object with \verb"extractWeights" and are charted with \verb"chart.Weights".

<<echo=TRUE>>=
extractWeights(out)
plot.new()
chart.Weights(out, colorset=bluemono)
@


Also, the value of the objective function at each rebalancing date is extracted with \verb"extractObjectiveMeasures".
<<echo=TRUE>>=
head(extractObjectiveMeasures(out))
@


The first and last observation from the estimation sample:
<<echo=TRUE>>=
out$opt_rebalancing[[1]]$data_summary
out$opt_rebalancing[[2]]$data_summary
@

The component contribution to risk at each rebalance date can be charted with \verb"chart.RiskBudget". The component contribution to risk in absolute or percentage.
<<echo=TRUE>>=
plot.new()
chart.RiskBudget(out, match.col="CVaR", risk.type="percentage", col=bluemono)
@

<<echo=TRUE>>=
plot.new()
chart.RiskBudget(out, match.col="CVaR", risk.type="absolute", col=bluemono)
@

Of course, DE is a stochastic optimizer and typically will only find a near-optimal solution that depends on the seed. The function \verb"optimize.portfolio.parallel" in \verb"PortfolioAnalytics" allows to run an arbitrary number of portfolio sets in parallel in order to develop "confidence bands" around your solution. It is based on Revolution's \verb"foreach" package \citep{foreach}.  

\bibliographystyle{abbrvnat}
\bibliography{PA}


\end{document}