# # Basic ETL optimization, using fPortfolio and an analytic global minimum variance # portfolio for comparision. The analytic GMV was used because fPortfolio seems to be # broken when it comes to calculating a GMV portfolio with shorts allowed. # # Author: Ian Kaplan # rm(list=ls()) library("zoo") library("fPortfolio") library("PerformanceAnalytics") # # This function is passed a matrix of returns. The result is # the weights, mean return and standard deviation of the # global minimum variance portfolio, where the only # constraint is full investment (e.g., shorting is allowed). # analyticGMV = function(returns) { S = cov(returns) Sinv = solve(S) one = rep(1, ncol(S)) denom = as.numeric(t(one) %*% Sinv %*% one) w = (Sinv %*% one) / denom sd = 1 / sqrt(denom) rslt = list(w = round(t(w), 8), mu = mean(returns %*% w), sd = sd) return(rslt) } # # Calculate the weights, mean return and standard deviation of the # tangency portfolio. # analyticTangency = function(returns, rf) { mu = apply(returns, 2, mean) S = cov(returns) Sinv = solve(S) one = rep(1, ncol(S)) mu_e = mu - rf w = (Sinv %*% mu_e) / as.numeric((t(one) %*% Sinv %*% mu_e)) sd = sqrt(t(w) %*% S %*% w) rslt = list(w = round(t(w), 6), mu = sum(w * mu), sd = sd) return(rslt) } # # Given two Markowitz mean/variance optimized portfolio weights, calculate a set of # means and 5% CVaR values for the efficient frontier using the "two fund theorm". # Note that these are not CVaR optimized portfolio. The CVaR value is calculated # from the weighted return distribution. # meanETLPoints = function(returns, gmv_wts, tan_wts, alpha = seq(from=-0.1, to=1.1, by=0.05)) { y_mu = list() x_etl = list() gmv_wts_t = as.vector(gmv_wts) tan_wts_t = as.vector(tan_wts) C = cov(returns) ix = 1 mu = apply(returns, 2, mean) for (i in alpha) { w = i * gmv_wts_t + ((1-i) * tan_wts_t) y_mu[ix] = mu %*% w x_etl[ix] = cvarRisk(returns, w) ix = ix + 1 } rslt = list(mu=y_mu, etl=x_etl) return(rslt) } # # Calculate the efficient frontier for CVaR (ETL) optimized portfolios between a mininum mean # and a maximum mean. This is done by calculating portfolios for a given target return. # etlPoints = function(returns.ts, minMu, maxMu ) { spec = portfolioSpec() setType(spec) = "CVaR" constraints = "Short" nPts = 40 muRange = seq(from=minMu, to=maxMu, by=((maxMu - minMu)/nPts)) mu = list() etl = list() for (i in 1:length(muRange)) { targetMu = muRange[i] setTargetReturn( spec ) = targetMu port = efficientPortfolio(data=returns.ts, spec=spec, constraints=constraints) mu[i] = targetMu etl[i] = cvarRisk(returns.ts, getWeights(port)) } rslt = list(mu = mu, etl = etl) } # # Read in a small cap. stock return data set and convert the dates into the ISO format used for a # zoo date index. The last two columns a "market" data set and the weekly risk free rate # (T-bill rate) # smallCap.z = read.zoo(file="smallcap_weekly.csv", sep=",", header=T, format = "%m/%d/%Y") smallCapWin.z = window(smallCap.z, start="2000-01-01") market = smallCapWin.z[,"Weekvwretd"] rf = smallCapWin.z[,"WeekRiskFree"] stocks.z = smallCapWin.z[,1:(ncol(smallCap.z)-2)] stocks.ts = as.timeSeries(stocks.z) mu = apply(stocks.ts, 2, mean) spec = portfolioSpec() setType(spec) = "CVaR" setRiskFreeRate(spec) = mean(rf) constraints = "Short" # Calculate the CVaR (ETL) optimized tangency portfolio weights tan.etl = tangencyPortfolio(data=stocks.ts, spec = spec, constraints = constraints) # calculate the CVaR (ETL) optimized global minimum variance portfolio weights gmv.etl = minvariancePortfolio(data=stocks.ts, spec=spec, constraints=constraints) # wts.tan.etl = getWeights(tan.etl) names(wts.tan.etl) = colnames(stocks.ts) # wts.gmv.etl= getWeights(gmv.etl) names(wts.gmv.etl) = colnames(stocks.ts) # gmv.etl.mu = mu %*% wts.gmv.etl tan.etl.mu = mu %*% wts.tan.etl # # Calculate the analytic Markowitz tangency portfolio weights wts.tan.mv = analyticTangency( stocks.ts, mean(rf) )$w # Calculate the Markowitz global minimum variance portfolio weights t = analyticGMV( stocks.ts ) wts.gmv.mv = t$w # Calculate the efficient Markowitz frontier, except with abs(ETL) on the x-axis mvPlot = meanETLPoints(stocks.ts, wts.gmv.mv, wts.tan.mv ) mv.mu = as.numeric(mvPlot$mu) mv.etl = abs(as.numeric(mvPlot$etl)) gmv.mu = mv.mu[ which.min(mv.etl)] gmv.etl = mv.etl[ which.min(mv.etl)] # Calculate the efficient CVaR optimized frontier etlPlot = etlPoints(stocks.ts, gmv.etl.mu, tan.etl.mu ) etl.mu = as.numeric(etlPlot$mu) etl.etl = abs(as.numeric(etlPlot$etl)) tan.mu = tan.etl.mu tan.etl = cvarRisk(stocks.ts, wts.tan.etl) par(mfrow=c(1,1)) xmin = min(etl.etl, mv.etl) xmax = max(etl.etl, mv.etl) ymin = min(etl.mu, mv.mu) ymax = max(etl.mu, mv.mu) plot(x=etl.etl, y=etl.mu, xlim=c(xmin, xmax), ylim=c(ymin, ymax), type="l", col="black", ylab="Weekly Portfolio Return", xlab="abs(ETL)", lwd=2 ) lines(x = mv.etl, y=mv.mu, type="l", col="red", lwd=2) points(x = c(abs(gmv.etl), abs(tan.etl)), y = c(gmv.mu, tan.mu), col="magenta", cex=1.5, pch=15) grid(col="grey") legend(x="topleft", legend=c("ETL optimized", "MV optimized"), col=c("black", "red"), lwd=4) par(mfrow=c(1,1))