# script used to run the portfolio optimizations
# Examples to consider
# Example 1: Consider a portfolio of stocks. Full investment and long
# only (or box) constraints. Objective to minimize portfolio variance.
# Demonstrate a custom moments function to compare a sample covariance
# matrix estimate and a robust covariance matrix estimate. An alternative
# to a MCD estimate is ledoit-wolf shrinkage, DCC GARCH model,
# factor model, etc.
# Example 2: Consider a portfolio of stocks. Dollar neutral, beta
# neutral, box constraints, and leverage_exposure constraints. Objective
# to minimize portfolio StdDev. This will demonstrate some of the
# more advanced constraint types. Could also introduce position limit
# constraints here in this example.
# Example 3: Consider an allocation to hedge funds using the
# EDHEC-Risk Alternative Index as a proxy. This will be an extended
# example starting with an objective to minimize portfolio expected
# shortfall, then risk budget percent contribution limit, then equal
# risk contribution limit.
# Example 4: Consider an allocation to hedge funds using the
# EDHEC-Risk Alternative Index as a proxy.
# Option 1 for example 4
# Objective to maximize a risk adjusted return measure
# (e.g.Calmar Ratio, Sterling Ratio, Sortino Ratio, or Upside Potential
# Ratio)
# I prefer doing this option
# Option 2 for example 4
# Objective to maximize the
# fourth order expansion of the Constant Relative Risk Aversion (CRRA)
# expected utility function. Demonstrate a custom moment function and
# a custom objective function.
# Set the directory to save the optimization results
results.dir <- "optimization_results"
# mix of blue, green, and red hues
my_colors <- c("#a6cee3", "#1f78b4", "#b2df8a", "#33a02c", "#fb9a99", "#e31a1c")
# Load the packages
library(PortfolioAnalytics)
library(foreach)
library(ROI)
library(ROI.plugin.quadprog)
# for running via Rscript
library(methods)
# Source in the lwShrink function
source("R/lwShrink.R")
# Example 1 and Example 2 will use the crsp_weekly data
# Example 3 and Example 4 will use the edhec data
source("data_prep.R")
##### Example 1 #####
stocks <- colnames(equity.data)
# Specify an initial portfolio
portf.init <- portfolio.spec(stocks)
# Add constraints
# weights sum to 1
portf.minvar <- add.constraint(portf.init, type="full_investment")
# box constraints
portf.minvar <- add.constraint(portf.minvar, type="box", min=0.01, max=0.45)
# Add objective
# objective to minimize portfolio variance
portf.minvar <- add.objective(portf.minvar, type="risk", name="var")
# Backtesting parameters
# Set rebalancing frequency
rebal.freq <- "quarters"
# Training Period
training <- 400
# Trailing Period
trailing <- 250
# Run optimization
# Sample Covariance Matrix Estimate
# By default, momentFUN uses set.portfolio.moments which computes the sample
# moment estimates
cat("Example 1: running minimum variance with sample covariance matrix
estimate backtest\n")
if(file.exists(paste(results.dir, "opt.minVarSample.rda", sep="/"))){
cat("file already exists\n")
} else {
opt.minVarSample <- optimize.portfolio.rebalancing(equity.data, portf.minvar,
optimize_method="ROI",
rebalance_on=rebal.freq,
training_period=training,
trailing_periods=trailing)
cat("opt.minVarSample complete. Saving results to ", results.dir, "\n")
save(opt.minVarSample, file=paste(results.dir, "opt.minVarSample.rda", sep="/"))
}
# Custom moment function to use Ledoit-Wolf shinkage covariance matrix estimate
lw.sigma <- function(R, ...){
out <- list()
# estimate covariance matrix via robust covariance matrix estimate,
# ledoit-wolf shrinkage, GARCH, factor model, etc.
# set.seed(1234)
# out$sigma <- MASS::cov.rob(R, method="mcd", ...)$cov
out$sigma <- lwShrink(R)$cov
#print(index(last(R)))
return(out)
}
cat("Example 1: running minimum variance with Ledoit-Wolf shrinkage covariance
matrix estimate backtest\n")
if(file.exists(paste(results.dir, "opt.minVarLW.rda", sep="/"))){
cat("file already exists\n")
} else{
# Using Ledoit-Wolf Shrinkage Covariance Matrix Estimate
opt.minVarLW <- optimize.portfolio.rebalancing(equity.data, portf.minvar,
optimize_method="ROI",
momentFUN=lw.sigma,
rebalance_on=rebal.freq,
training_period=training,
trailing_periods=trailing)
cat("opt.minVarLW complete. Saving results to ", results.dir, "\n")
save(opt.minVarLW, file=paste(results.dir, "opt.minVarLW.rda", sep="/"))
}
##### Example 2 #####
portf.init <- portfolio.spec(stocks)
# weights sum to 0
portf.dn <- add.constraint(portf.init, type="weight_sum",
min_sum=-0.01, max_sum=0.01)
# Add box constraints such that no stock has weight less than -20% or
# greater than 20%
portf.dn <- add.constraint(portf.dn, type="box",
min=-0.2, max=0.2)
# Add position limit constraint such that the portfolio has a maximum
# of 20 non-zero positions
portf.dn <- add.constraint(portf.dn, type="position_limit", max_pos=20)
# Compute the betas of each stock
betas <- t(CAPM.beta(equity.data, market, Rf))
# Add factor exposure constraint to limit portfolio beta
portf.dn <- add.constraint(portf.dn, type="factor_exposure", B=betas,
lower=-0.25, upper=0.25)
# portf.dn <- add.constraint(portf.dn, type="leverage_exposure", leverage=2)
# generate random portfolios
if(file.exists(paste(results.dir, "rp.rda", sep="/"))){
cat("random portfolios already generated\n")
} else {
cat("generating random portfolios\n")
rp <- random_portfolios(portf.dn, 10000, eliminate=TRUE)
cat("random portfolios generated. Saving rp to ", results.dir, "\n")
save(rp, file=paste(results.dir, "rp.rda", sep="/"))
}
# Add objective to maximize return
portf.dn.StdDev <- add.objective(portf.dn, type="return", name="mean",
target=0.0015)
# Add objective to target a portfolio standard deviation of 2%
portf.dn.StdDev <- add.objective(portf.dn.StdDev, type="risk", name="StdDev",
target=0.02)
cat("Example 2: running dollar neutral optimization\n")
if(file.exists(paste(results.dir, "opt.dn.rda", sep="/"))){
cat("file already exists\n")
} else {
# Run optimization
opt.dn <- optimize.portfolio(equity.data, portf.dn.StdDev,
optimize_method="random", rp=rp,
trace=TRUE)
cat("opt.dn complete. Saving results to ", results.dir, "\n")
save(opt.dn, file=paste(results.dir, "opt.dn.rda", sep="/"))
}
##### Example 3 #####
# Example 3 will consider three portfolios
# - minES
# - minES with component contribution limit
# - minES with equal risk contribution
funds <- colnames(R)
portf.init <- portfolio.spec(funds)
portf.init <- add.constraint(portf.init, type="weight_sum",
min_sum=0.99, max_sum=1.01)
portf.init <- add.constraint(portf.init, type="box",
min=0.05, max=0.4)
# Set multiplier=0 so that it is calculated, but does not affect the optimization
portf.init <- add.objective(portf.init, type="return",
name="mean", multiplier=0)
# Add objective to minimize expected shortfall
portf.minES <- add.objective(portf.init, type="risk", name="ES")
# Add risk budget objective with upper limit on percentage contribution
portf.minES.RB <- add.objective(portf.minES, type="risk_budget",
name="ES", max_prisk=0.3)
# Relax the box constraint
portf.minES.RB$constraints[[2]]$max <- rep(1,ncol(R))
# print.default(portf.minES.RB$constraints[[2]])
# Add risk budget objective to minimize concentration of percentage component
# contribution to risk. Concentration is defined as the Herfindahl-Hirschman
# Index (HHI). $\sum_i x_i^2$
portf.minES.EqRB <- add.objective(portf.minES, type="risk_budget",
name="ES", min_concentration=TRUE)
# relax the box constraint
portf.minES.EqRB <- add.constraint(portf.minES.EqRB, type="box",
min=0.05, max=1, indexnum=2)
# portf.minES.RB$constraints[[2]]$max <- rep(1,ncol(R))
# print.default(portf.minES.EqRB$constraints[[2]])
# Add risk budget objective to minES portfolio with multiplier=0 so that it
# is calculated, but does not affect optimization
portf.minES <- add.objective(portf.minES, type="risk_budget",
name="ES", multiplier=0)
# Combine the portfolios so we can make a single call to
# optimize.portfolio
portf <- combine.portfolios(list(minES=portf.minES,
minES.RB=portf.minES.RB,
minES.EqRB=portf.minES.EqRB))
cat("Example 3: running minimum ES optimizations\n")
if(file.exists(paste(results.dir, "opt.minES.rda", sep="/"))){
cat("file already exists\n")
} else {
# Run the optimization
opt.minES <- optimize.portfolio(R, portf, optimize_method="DEoptim",
search_size=5000, trace=TRUE, traceDE=0,
message=TRUE)
cat("opt.minES complete. Saving results to ", results.dir, "\n")
save(opt.minES, file=paste(results.dir, "opt.minES.rda", sep="/"))
}
# Now we want to evaluate the optimization through time
# Rebalancing parameters
# Set rebalancing frequency
rebal.freq <- "quarters"
# Training Period
training <- 120
# Trailing Period
trailing <- 72
cat("Example 3: running minimum ES backtests\n")
if(file.exists(paste(results.dir, "bt.opt.minES.rda", sep="/"))){
cat("file already exists\n")
} else {
# Backtest
bt.opt.minES <- optimize.portfolio.rebalancing(R, portf,
optimize_method="DEoptim",
rebalance_on=rebal.freq,
training_period=training,
trailing_periods=trailing,
search_size=5000,
traceDE=0, message=TRUE)
cat("bt.opt.minES complete. Saving results to ", results.dir, "\n")
save(bt.opt.minES, file=paste(results.dir, "bt.opt.minES.rda", sep="/"))
}
##### Example 4 #####
# Simple function to compute the moments used in CRRA
crra.moments <- function(R, ...){
out <- list()
out$mu <- colMeans(R)
out$sigma <- cov(R)
out$m3 <- PerformanceAnalytics:::M3.MM(R)
out$m4 <- PerformanceAnalytics:::M4.MM(R)
out
}
# Fourth order expansion of CRRA expected utility
CRRA <- function(R, weights, lambda, sigma, m3, m4){
weights <- matrix(weights, ncol=1)
M2.w <- t(weights) %*% sigma %*% weights
M3.w <- t(weights) %*% m3 %*% (weights %x% weights)
M4.w <- t(weights) %*% m4 %*% (weights %x% weights %x% weights)
term1 <- 0.5 * lambda * M2.w
term2 <- (1 / 6) * lambda * (lambda + 1) * M3.w
term3 <- (1 / 24) * lambda * (lambda + 1) * (lambda + 2) * M4.w
out <- -term1 + term2 - term3
out
}
# test the CRRA function
portf.crra <- portfolio.spec(funds)
portf.crra <- add.constraint(portf.crra, type="weight_sum",
min_sum=0.99, max_sum=1.01)
portf.crra <- add.constraint(portf.crra, type="box",
min=0.05, max=0.4)
portf.crra <- add.objective(portf.crra, type="return",
name="CRRA", arguments=list(lambda=10))
# I just want these for plotting
# Set multiplier=0 so that it is calculated, but does not affect the optimization
portf.crra <- add.objective(portf.crra, type="return", name="mean", multiplier=0)
portf.crra <- add.objective(portf.crra, type="risk", name="ES", multiplier=0)
portf.crra <- add.objective(portf.crra, type="risk", name="StdDev", multiplier=0)
cat("Example 4: running maximum CRRA optimization\n")
if(file.exists(paste(results.dir, "opt.crra.rda", sep="/"))){
cat("file already exists\n")
} else {
# Run the optimization
opt.crra <- optimize.portfolio(R, portf.crra, optimize_method="DEoptim",
search_size=5000, trace=TRUE, traceDE=0,
momentFUN="crra.moments")
cat("opt.crra complete. Saving results to ", results.dir, "\n")
save(opt.crra, file=paste(results.dir, "opt.crra.rda", sep="/"))
}
cat("Example 4: running maximum CRRA backtest\n")
if(file.exists(paste(results.dir, "bt.opt.crra.rda", sep="/"))){
cat("file already exists\n")
} else {
# Run the optimization with rebalancing
bt.opt.crra <- optimize.portfolio.rebalancing(R, portf.crra,
optimize_method="DEoptim",
search_size=5000, trace=TRUE,
traceDE=0,
momentFUN="crra.moments",
rebalance_on=rebal.freq,
training_period=training,
trailing_periods=trailing)
cat("bt.opt.crra complete. Saving results to ", results.dir, "\n")
save(bt.opt.crra, file=paste(results.dir, "bt.opt.crra.rda", sep="/"))
}
##### RP Demo #####
cat("Random portfolio method comparison\n")
if(file.exists("figures/rp_plot.png") & file.exists("figures/rp_viz.rda")){
cat("file already exists\n")
} else {
portf.lo <- portfolio.spec(colnames(R))
portf.lo <- add.constraint(portf.lo, type="weight_sum",
min_sum=0.99, max_sum=1.01)
portf.lo <- add.constraint(portf.lo, type="long_only")
# Use the long only portfolio previously created
# Generate random portfolios using the 3 methods
rp1 <- random_portfolios(portf.lo, permutations=2000,
rp_method='sample')
rp2 <- random_portfolios(portf.lo, permutations=2000,
rp_method='simplex')
rp3 <- random_portfolios(portf.lo, permutations=2000,
rp_method='grid')
# Calculate the portfolio mean return and standard deviation
rp1_mean <- apply(rp1, 1, function(x) mean(R %*% x))
rp1_StdDev <- apply(rp1, 1, function(x) StdDev(R, weights=x))
rp2_mean <- apply(rp2, 1, function(x) mean(R %*% x))
rp2_StdDev <- apply(rp2, 1, function(x) StdDev(R, weights=x))
rp3_mean <- apply(rp3, 1, function(x) mean(R %*% x))
rp3_StdDev <- apply(rp3, 1, function(x) StdDev(R, weights=x))
x.assets <- StdDev(R)
y.assets <- colMeans(R)
###
require(rCharts)
# create an interactive plot using rCharts and nvd3 scatterChart
tmp1 <- data.frame(name="sample", mean=rp1_mean, sd=rp1_StdDev)
tmp2 <- data.frame(name="simplex", mean=rp2_mean, sd=rp2_StdDev)
tmp3 <- data.frame(name="grid", mean=rp3_mean, sd=rp3_StdDev)
tmp <- rbind(tmp1, tmp2, tmp3)
rp_viz <- nPlot(mean ~ sd, group="name", data=tmp, type="scatterChart")
rp_viz$xAxis(
axisLabel = 'Risk (std. dev.)'
,tickFormat = "#!d3.format('0.4f')!#"
)
rp_viz$yAxis(
axisLabel = 'Return'
,tickFormat = "#!d3.format('0.4f')!#"
)
rp_viz$chart(color = my_colors[c(2,4,6)])
#set left margin so y axis label will show up
rp_viz$chart( margin = list(left = 100) )
# rp_viz$chart(
# tooltipContent = "#!
# function(a,b,c,d) {
# //d has all the info you need
# return( '<h3>' + d.point.series + '</h3>Return: ' + d.point.y + '<br>Risk: ' + d.point.x)
# }
# !#")
####if you do not want fisheye/magnify
####let me know, and will show how to remove
####this will solve the tooltip problem
save(rp_viz, file="figures/rp_viz.rda")
###
x.lower <- min(x.assets) * 0.9
x.upper <- max(x.assets) * 1.1
y.lower <- min(y.assets) * 0.9
y.upper <- max(y.assets) * 1.1
png("figures/rp_plot.png", height = 500, width = 1000)
# plot feasible portfolios
plot(x=rp1_StdDev, y=rp1_mean, col=my_colors[2], main="Random Portfolio Methods",
ylab="mean", xlab="StdDev", xlim=c(x.lower, x.upper),
ylim=c(y.lower, y.upper))
points(x=rp2_StdDev, y=rp2_mean, col=my_colors[4], pch=2)
points(x=rp3_StdDev, y=rp3_mean, col=my_colors[6], pch=5)
points(x=x.assets, y=y.assets)
text(x=x.assets, y=y.assets, labels=colnames(R), pos=4, cex=0.8)
legend("bottomright", legend=c("sample", "simplex", "grid"),
col=my_colors[c(2,4,6)],
pch=c(1, 2, 5), bty="n")
dev.off()
}
cat("Random portfolio simplex method fev biasing\n")
if(file.exists("figures/fev_plot.png")){
cat("file already exists\n")
} else {
png("figures/fev_plot.png", height = 500, width = 1000)
fev <- 0:5
x.assets <- StdDev(R)
y.assets <- colMeans(R)
par(mfrow=c(2, 3))
for(i in 1:length(fev)){
rp <- rp_simplex(portfolio=portf.lo, permutations=2000, fev=fev[i])
tmp.mean <- apply(rp, 1, function(x) mean(R %*% x))
tmp.StdDev <- apply(rp, 1, function(x) StdDev(R=R, weights=x))
x.lower <- min(c(tmp.StdDev, x.assets)) * 0.85
x.upper <- max(c(tmp.StdDev, x.assets)) * 1.15
y.lower <- min(c(tmp.mean, y.assets)) * 0.85
y.upper <- max(c(tmp.mean, y.assets)) * 1.15
plot(x=tmp.StdDev, y=tmp.mean, main=paste("FEV =", fev[i]),
ylab="mean", xlab="StdDev", col=rgb(0, 0, 100, 50, maxColorValue=255),
xlim=c(x.lower, x.upper),
ylim=c(y.lower, y.upper))
points(x=x.assets, y=y.assets)
text(x=x.assets, y=y.assets, labels=colnames(R), pos=4, cex=0.8)
}
par(mfrow=c(1,1))
dev.off()
}
# # Calculate the turnover per period
# turnover.rebalancing <- function(object){
# weights <- extractWeights(object)
# n <- nrow(weights)
# out <- vector("numeric", n)
# out[1] <- NA
# for(i in 2:n){
# out[i] <- out[i] <- sum(abs(as.numeric(weights[i,]) - as.numeric(weights[i-1,])))
# }
# xts(out, index(weights))
# }
#
# # Calculate the diversification per period
# diversification.rebalancing <- function(object){
# weights <- extractWeights(object)
# n <- nrow(weights)
# out <- vector("numeric", n)
# for(i in 1:n){
# out[i] <- 1 - sum(weights[i,]^2)
# }
# xts(out, index(weights))
# }
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.