# How to compute modified-CVaR in the PerformanceAnalytics package?

My objective is to measure the modified-CVAR for a portfolio given its weights and matrix of security returns. Luckily the wonderful package PerformanceAnalytics has an ES() function that does just this.

The issue I am having is that modified-CVAR takes minutes to compute, when according to this paper (by the same authors) the algorithm should only take seconds: "...we provide a long but explicit formula for computing the derivative of mES. Although the resulting formulae are rather complex, they lend themselves to efficient translation into a simple algorithm that computes in less than a second mES and component mES, even for portfolios with a very large number of assets." (page 14)

I re-produce code from a vignette which works fine when using method = "gaussian" but not method = "modified". I have also reviewed the CRAN reference for the PerformanceAnalytics package although it is not as clear as the paper (linked above).

library(PerformanceAnalytics)

tickers = c( "VNO" , "VMC" , "WMT" , "WAG" , "DIS" , "WPO" , "WFC" , "WDC" ,
"WY" , "WHR" , "WMB" , "WEC" , "XEL" , "XRX" , "XLNX" ,"ZION" ,"MMM" ,
"ABT", "ADBE" , "AMD" , "AET" , "AFL" , "APD" , "ARG" ,"AA" , "AGN" ,
"ALTR" , "MO" , "AEP" , "AXP" , "AIG" , "AMGN" , "APC" ,"ADI" , "AON" ,
"AVP", "BHI" , "BLL" , "BAC" , "BK" , "BCR" , "BAX" , "BBT" , "BDX" ,
"BMS" , "BBY" , "BIG" , "HRB" , "BMC" , "BA" , "BMY" , "CA" , "COG" ,
"CPB" , "CAH" , "CCL" , "CAT" , "CELG" , "CNP" , "CTL" , "CEPH", "CERN" ,
"SCHW" , "CVX" , "CB" , "CI" ,"CINF" ,"CTAS" , "CSCO" , "C" , "CLF" ,
"CLX", "CMS" , "KO" , "CCE" , "CL" , "CMCSA" ,"CMA" , "CSC" , "CAG" ,
"COP" , "ED" , "CEG" ,"GLW" , "COST" , "CVH" , "CSX" , "CMI" , "CVS" ,
"DHR" , "DE")

library(quantmod)
getSymbols(tickers, from = "2000-12-01", to = "2010-12-31")
P <- NULL; seltickers <- NULL
for(ticker in tickers) {
tmp <- Cl(to.monthly(eval(parse(text=ticker))))
if(is.null(P)){ timeP = time(tmp) }
if( any( time(tmp)!=timeP )) next
else P<-cbind(P,as.numeric(tmp))
seltickers = c( seltickers , ticker )
}

P = xts(P,order.by=timeP)
colnames(P) <- seltickers
R <- diff(log(P))
R <- R[-1,]
dim(R)

mu <- colMeans(R)
sigma <- cov(R)

obj <- function(w) {
if (sum(w) == 0) {
w <- w + 1e-2
}
w <- w / sum(w)
CVaR <- ES(weights = w,
method = "gaussian",
portfolio_method = "component",
mu = mu,
sigma = sigma)
tmp1 <- CVaR$ES tmp2 <- max(CVaR$pct_contrib_ES - 0.05, 0)
out <- tmp1 + 1e3 * tmp2
return(out)
}

N <- ncol(R)
minw <- 0
maxw <- 1
lower <- rep(minw,N)
upper <- rep(maxw,N)

w<-rep(100/120 , 100)

# works
CVaR1 <- ES(weights = w, method = "gaussian",  portfolio_method = "component", mu = mu, sigma = sigma)

# takes too long
date()
CVaR4 <- ES(R = R , weights = w, method = "modified" , portfolio_method = "component" , clean = "boudt")
date()


Update:

It turns out that if you want to estimate m3 and m4 (skewness and kurtosis) you need to construct matrices with dimension (number of assets) raised to the 3rd and 4th powers. Therefore for large matrices such as the S&P 500, the memory demands are significant - my back of the envelope calculations are 25Gb. So this procedure suffers from the curse of dimensionality.

-

The "Component ES" section of ?ES says:

For the decomposition of Gaussian ES, the estimated mean and covariance matrix are needed. For the decomposition of modified ES, also estimates of the coskewness and cokurtosis matrices are needed.

The estimate of the coskewness and cokurtosis matrices are what take such a long time. You can calculate them beforehand and pass them to ES. ?ES says:

The matrices can be estimated through the functions ‘skewness.MM’ and ‘kurtosis.MM’.

but I do not see those functions in the version of PerformanceAnalytics installed on my system. ES itself uses the unexported functions M3.MM and M4.MM, so you could call them explicitly:

m3 <- PerformanceAnalytics:::M3.MM(R)
m4 <- PerformanceAnalytics:::M4.MM(R)
CVaR4 <- ES(R=R , weights=w, method="modified", portfolio_method="component",
clean="boudt", m3=m3, m4=m4)

-
fyi - I made a small fix to 2nd line. m4 <- PerformanceAnalytics:::M4.MM(R). This step does take several minutes to compute. However, this helps tremendously since I am optimizing with CVAR and I can store m3 and m4 in memory –  Quant Guy Sep 23 '11 at 14:15
@QuantGuy: thanks for the edit. I'm a bit sick and my mind is extra-foggy. ;-) –  Joshua Ulrich Sep 23 '11 at 14:51
no prob. The community has a massive debt to you for your contributions! –  Quant Guy Sep 23 '11 at 16:35
I want to second the advice to pre-calculate the moment matrices before doing optimization. The PortfolioAnalytics package always does this to minimize recalculating things that only need to be done once. You can also do things like cleaning the raw data, Ledoit-Wolf shrinkage on the estimates, etc before applying your optimization criteria. –  Brian G. Peterson Sep 23 '11 at 16:52