How to compute modified-CVaR in the PerformanceAnalytics package?

Question

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" ,
 "APA", "AAPL" , "AMAT" ,"ADM" , "T" , "ADSK" , "ADP" , "AZO" , "AVY" ,
 "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.

Answer 1

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)